Belle II Software development
Helix.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9//-----------------------------------------------------------------------------
10// Description : A class to represent a track helix parameter in Belle style
11//-----------------------------------------------------------------------------
12
13#include <math.h>
14#include <float.h>
15#include "trg/cdc/Helix.h"
16#include "CLHEP/Matrix/Matrix.h"
17
18namespace Belle2 {
24// const double
25// TRGCDCHelix::m_BFIELD = 15.0; // KG
26// const double
27// TRGCDCHelix::m_ALPHA = 222.376063; // = 10000. / 2.99792458 / BFIELD
28
30 const double
31 M_PI2 = 2. * M_PI;
32
34 const double
35 M_PI4 = 4. * M_PI;
36
38 const double
39 M_PI8 = 8. * M_PI;
40
41 const double
42 TRGCDCHelix::ConstantAlpha = 222.376063;
43
44 const std::string TRGCDCHelix::invalidhelix("Invalid TRGCDCHelix");
45 CLHEP::HepVector TRGCDCHelix::ms_amin(5, 0), TRGCDCHelix::ms_amax(5, 0);
46 bool TRGCDCHelix::ms_check_range(false);
47 bool TRGCDCHelix::ms_throw_exception(false);
48 bool TRGCDCHelix::ms_print_debug(false);
49
50 bool TRGCDCHelix::set_exception(bool t)
51 {
52 return ms_throw_exception = t;
53 }
54
55 bool TRGCDCHelix::set_print(bool t)
56 {
57 return ms_print_debug = t;
58 }
59
60
61 void TRGCDCHelix::set_limits(const CLHEP::HepVector& a_min, const CLHEP::HepVector& a_max)
62 {
63 if (a_min.num_row() != 5 || a_max.num_row() != 5) return;
64 ms_amin = a_min;
65 ms_amax = a_max;
66 ms_check_range = true;
67 }
68
69
70 TRGCDCHelix::TRGCDCHelix(const HepGeom::Point3D<double>& pivot,
71 const CLHEP::HepVector& a,
72 const CLHEP::HepSymMatrix& Ea)
73 : m_matrixValid(true),
74 m_helixValid(false),
75 m_bField(15.0),
76 m_alpha(222.376063),
77 m_pivot(pivot),
78 m_a(a),
79 m_Ea(Ea)
80 {
81 // m_alpha = 10000. / 2.99792458 / m_bField;
82 // m_alpha = 222.376063;
83 if (m_a.num_row() == 5 && m_Ea.num_row() == 5) {
84 updateCache();
85 }
86 }
87
88 TRGCDCHelix::TRGCDCHelix(const HepGeom::Point3D<double>& pivot,
89 const CLHEP::HepVector& a)
90 : m_matrixValid(false),
91 m_helixValid(false),
92 m_bField(15.0),
93 m_alpha(222.376063),
94 m_pivot(pivot),
95 m_a(a),
96 m_Ea(CLHEP::HepSymMatrix(5, 0))
97 {
98 // m_alpha = 222.376063;
99 if (m_a.num_row() == 5) {
100 updateCache();
101 }
102 }
103
104 TRGCDCHelix::TRGCDCHelix(const HepGeom::Point3D<double>& position,
105 const CLHEP::Hep3Vector& momentum,
106 double charge)
107 : m_matrixValid(false),
108 m_helixValid(false),
109 m_bField(15.0),
110 m_alpha(222.376063),
111 m_pivot(position),
112 m_a(CLHEP::HepVector(5, 0)),
113 m_Ea(CLHEP::HepSymMatrix(5, 0))
114 {
115 m_a[0] = 0.;
116 m_a[3] = 0.;
117 double perp(momentum.perp());
118 if (perp != 0.0) {
119 m_a[1] = fmod(atan2(- momentum.x(), momentum.y())
120 + M_PI4, M_PI2);
121 m_a[2] = charge / perp;
122 m_a[4] = momentum.z() / perp;
123 } else {
124 m_a[2] = charge * (DBL_MAX);
125 }
126 // m_alpha = 222.376063;
127 updateCache();
128 }
129
130 TRGCDCHelix::~TRGCDCHelix()
131 {
132 }
133
134 HepGeom::Point3D<double>
135 TRGCDCHelix::x(double phi) const
136 {
137 //
138 // Calculate position (x,y,z) along helix.
139 //
140 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
141 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
142 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
143 //
144
145 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
146 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
147 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
148
149 return HepGeom::Point3D<double> (x, y, z);
150 }
151
152 double*
153 TRGCDCHelix::x(double phi, double p[3]) const
154 {
155 //
156 // Calculate position (x,y,z) along helix.
157 //
158 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
159 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
160 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
161 //
162
163 p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
164 p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
165 p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
166
167 return p;
168 }
169
170 HepGeom::Point3D<double>
171 TRGCDCHelix::x(double phi, CLHEP::HepSymMatrix& Ex) const
172 {
173
174 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
175 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
176 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
177
178 //
179 // Calculate position error matrix.
180 // Ex(phi) = (@x/@a)(Ea)(@x/@a)^T, phi is deflection angle to specify the
181 // point to be calcualted.
182 //
183 // CLHEP::HepMatrix dXDA(3, 5, 0);
184 // dXDA = delXDelA(phi);
185 // Ex.assign(dXDA * m_Ea * dXDA.T());
186
187 if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
188 else Ex = m_Ea;
189
190 return HepGeom::Point3D<double> (x, y, z);
191 }
192
193 CLHEP::Hep3Vector
194 TRGCDCHelix::momentum(double phi) const
195 {
196 //
197 // Calculate momentum.
198 //
199 // Pt = | 1/kappa | (GeV/c)
200 //
201 // Px = -Pt * sin(phi0 + phi)
202 // Py = Pt * cos(phi0 + phi)
203 // Pz = Pt * tan(lambda)
204 //
205
206 double pt = fabs(m_pt);
207 double px = - pt * sin(m_ac[1] + phi);
208 double py = pt * cos(m_ac[1] + phi);
209 double pz = pt * m_ac[4];
210
211 return CLHEP::Hep3Vector(px, py, pz);
212 }
213
214 CLHEP::Hep3Vector
215 TRGCDCHelix::momentum(double phi, CLHEP::HepSymMatrix& Em) const
216 {
217 //
218 // Calculate momentum.
219 //
220 // Pt = | 1/kappa | (GeV/c)
221 //
222 // Px = -Pt * sin(phi0 + phi)
223 // Py = Pt * cos(phi0 + phi)
224 // Pz = Pt * tan(lambda)
225 //
226
227 double pt = fabs(m_pt);
228 double px = - pt * sin(m_ac[1] + phi);
229 double py = pt * cos(m_ac[1] + phi);
230 double pz = pt * m_ac[4];
231
232 if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
233 else Em = m_Ea;
234
235 return CLHEP::Hep3Vector(px, py, pz);
236 }
237
238 CLHEP::HepLorentzVector
239 TRGCDCHelix::momentum(double phi, double mass) const
240 {
241 //
242 // Calculate momentum.
243 //
244 // Pt = | 1/kappa | (GeV/c)
245 //
246 // Px = -Pt * sin(phi0 + phi)
247 // Py = Pt * cos(phi0 + phi)
248 // Pz = Pt * tan(lambda)
249 //
250 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
251
252 double pt = fabs(m_pt);
253 double px = - pt * sin(m_ac[1] + phi);
254 double py = pt * cos(m_ac[1] + phi);
255 double pz = pt * m_ac[4];
256 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
257
258 return CLHEP::HepLorentzVector(px, py, pz, E);
259 }
260
261
262 CLHEP::HepLorentzVector
263 TRGCDCHelix::momentum(double phi, double mass, CLHEP::HepSymMatrix& Em) const
264 {
265 //
266 // Calculate momentum.
267 //
268 // Pt = | 1/kappa | (GeV/c)
269 //
270 // Px = -Pt * sin(phi0 + phi)
271 // Py = Pt * cos(phi0 + phi)
272 // Pz = Pt * tan(lambda)
273 //
274 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
275
276 double pt = fabs(m_pt);
277 double px = - pt * sin(m_ac[1] + phi);
278 double py = pt * cos(m_ac[1] + phi);
279 double pz = pt * m_ac[4];
280 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
281
282 if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi, mass));
283 else Em = m_Ea;
284
285 return CLHEP::HepLorentzVector(px, py, pz, E);
286 }
287
288 CLHEP::HepLorentzVector
289 TRGCDCHelix::momentum(double phi,
290 double mass,
291 HepGeom::Point3D<double>& x,
292 CLHEP::HepSymMatrix& Emx) const
293 {
294
295 //
296 // Calculate momentum.
297 //
298 // Pt = | 1/kappa | (GeV/c)
299 //
300 // Px = -Pt * sin(phi0 + phi)
301 // Py = Pt * cos(phi0 + phi)
302 // Pz = Pt * tan(lambda)
303 //
304 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
305
306 double pt = fabs(m_pt);
307 double px = - pt * sin(m_ac[1] + phi);
308 double py = pt * cos(m_ac[1] + phi);
309 double pz = pt * m_ac[4];
310 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
311
312 x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
313 x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
314 x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);
315
316 if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi, mass));
317 else Emx = m_Ea;
318
319 return CLHEP::HepLorentzVector(px, py, pz, E);
320 }
321
322
323 const HepGeom::Point3D<double>&
324 TRGCDCHelix::pivot(const HepGeom::Point3D<double>& newPivot)
325 {
326#if defined(BELLE_DEBUG)
327 try {
328#endif
329 const double& dr = m_ac[0];
330 const double& phi0 = m_ac[1];
331 const double& kappa = m_ac[2];
332 const double& dz = m_ac[3];
333 const double& tanl = m_ac[4];
334
335 double rdr = dr + m_r;
336 double phi = fmod(phi0 + M_PI4, M_PI2);
337 double csf0 = cos(phi);
338 double snf0 = (1. - csf0) * (1. + csf0);
339 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
340 if (phi > M_PI) snf0 = - snf0;
341
342 double xc = m_pivot.x() + rdr * csf0;
343 double yc = m_pivot.y() + rdr * snf0;
344 double csf, snf;
345 if (m_r != 0.0) {
346 csf = (xc - newPivot.x()) / m_r;
347 snf = (yc - newPivot.y()) / m_r;
348 double anrm = sqrt(csf * csf + snf * snf);
349 if (anrm != 0.0) {
350 csf /= anrm;
351 snf /= anrm;
352 phi = atan2(snf, csf);
353 } else {
354 csf = 1.0;
355 snf = 0.0;
356 phi = 0.0;
357 }
358 } else {
359 csf = 1.0;
360 snf = 0.0;
361 phi = 0.0;
362 }
363 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
364 if (phid > M_PI) phid = phid - M_PI2;
365 double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
366 * csf
367 + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
368 double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();
369
370 CLHEP::HepVector ap(5);
371 ap[0] = drp;
372 ap[1] = fmod(phi + M_PI4, M_PI2);
373 ap[2] = kappa;
374 ap[3] = dzp;
375 ap[4] = tanl;
376
377 // if (m_matrixValid) m_Ea.assign(delApDelA(ap) * m_Ea * delApDelA(ap).T());
378 if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));
379
380 m_a = ap;
381 m_pivot = newPivot;
382
383 //...Are these needed?...iw...
384 updateCache();
385 return m_pivot;
386#if defined(BELLE_DEBUG)
387 } catch (...) {
388 m_helixValid = false;
389 if (ms_throw_exception) throw invalidhelix;
390 }
391#endif
392 return m_pivot;
393 }
394
395 void
396 TRGCDCHelix::set(const HepGeom::Point3D<double>& pivot,
397 const CLHEP::HepVector& a,
398 const CLHEP::HepSymMatrix& Ea)
399 {
400 m_pivot = pivot;
401 m_a = a;
402 m_Ea = Ea;
403 m_matrixValid = true;
404 m_helixValid = false;
405 updateCache();
406 }
407
408 TRGCDCHelix&
409 TRGCDCHelix::operator = (const TRGCDCHelix& i)
410 {
411 if (this == & i) return * this;
412
413 m_bField = i.m_bField;
414 m_alpha = i.m_alpha;
415 m_pivot = i.m_pivot;
416 m_a = i.m_a;
417 m_Ea = i.m_Ea;
418 m_matrixValid = i.m_matrixValid;
419 m_helixValid = i.m_helixValid;
420
421 m_center = i.m_center;
422 m_cp = i.m_cp;
423 m_sp = i.m_sp;
424 m_pt = i.m_pt;
425 m_r = i.m_r;
426 m_ac[0] = i.m_ac[0];
427 m_ac[1] = i.m_ac[1];
428 m_ac[2] = i.m_ac[2];
429 m_ac[3] = i.m_ac[3];
430 m_ac[4] = i.m_ac[4];
431
432 return * this;
433 }
434
435 void
436 TRGCDCHelix::updateCache(void)
437 {
438#if defined(BELLE_DEBUG)
439 checkValid();
440 if (m_helixValid) {
441#endif
442 //
443 // Calculate TRGCDCHelix center( xc, yc ).
444 //
445 // xc = x0 + (dr + (alpha / kappa)) * cos(phi0) (cm)
446 // yc = y0 + (dr + (alpha / kappa)) * sin(phi0) (cm)
447 //
448 m_ac[0] = m_a[0];
449 m_ac[1] = m_a[1];
450 m_ac[2] = m_a[2];
451 m_ac[3] = m_a[3];
452 m_ac[4] = m_a[4];
453
454 m_cp = cos(m_ac[1]);
455 m_sp = sin(m_ac[1]);
456 if (m_ac[2] != 0.0) {
457 if (m_ac[2] == DBL_MAX || m_ac[2] == (-DBL_MAX)) {
458 m_pt = m_r = 0;
459 return;
460 } else {
461 m_pt = 1. / m_ac[2];
462 m_r = m_alpha / m_ac[2];
463 }
464 } else {
465 m_pt = (DBL_MAX);
466 m_r = (DBL_MAX);
467 return;
468 }
469
470 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
471 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
472 m_center.setX(x);
473 m_center.setY(y);
474 m_center.setZ(0.);
475#if defined(BELLE_DEBUG)
476 } else {
477 m_ac[0] = m_a[0];
478 m_ac[1] = m_a[1];
479 m_ac[2] = m_a[2];
480 m_ac[3] = m_a[3];
481 m_ac[4] = m_a[4];
482
483 m_cp = cos(m_ac[1]);
484 m_sp = sin(m_ac[1]);
485 if (m_ac[2] != 0.0) {
486 if (m_ac[2] == DBL_MAX || m_ac[2] == (-DBL_MAX)) {
487 m_pt = m_r = 0;
488 return;
489 } else {
490 m_pt = 1. / m_ac[2];
491 m_r = m_alpha / m_ac[2];
492 }
493 } else {
494 m_pt = (DBL_MAX);
495 m_r = (DBL_MAX);
496 return;
497 }
498
499 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
500 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
501 m_center.setX(x);
502 m_center.setY(y);
503 m_center.setZ(0.);
504 }
505#endif
506 }
507
508 CLHEP::HepMatrix
509 TRGCDCHelix::delApDelA(const CLHEP::HepVector& ap) const
510 {
511 //
512 // Calculate Jacobian (@ap/@a)
513 // CLHEP::HepVector ap is new helix parameters and a is old helix parameters.
514 //
515
516 CLHEP::HepMatrix dApDA(5, 5, 0);
517
518 const double& dr = m_ac[0];
519 const double& phi0 = m_ac[1];
520 const double& cpa = m_ac[2];
521//cnv const double & dz = m_ac[3];
522 const double& tnl = m_ac[4];
523
524 double drp = ap[0];
525 double phi0p = ap[1];
526//cnv double cpap = ap[2];
527//cnv double dzp = ap[3];
528//cnv double tnlp = ap[4];
529
530 double rdr = m_r + dr;
531 double rdrpr;
532 if ((m_r + drp) != 0.0) {
533 rdrpr = 1. / (m_r + drp);
534 } else {
535 rdrpr = (DBL_MAX);
536 }
537 // double csfd = cos(phi0)*cos(phi0p) + sin(phi0)*sin(phi0p);
538 // double snfd = cos(phi0)*sin(phi0p) - sin(phi0)*cos(phi0p);
539 double csfd = cos(phi0p - phi0);
540 double snfd = sin(phi0p - phi0);
541 double phid = fmod(phi0p - phi0 + M_PI8, M_PI2);
542 if (phid > M_PI) phid = phid - M_PI2;
543
544 dApDA[0][0] = csfd;
545 dApDA[0][1] = rdr * snfd;
546 if (cpa != 0.0) {
547 dApDA[0][2] = (m_r / cpa) * (1.0 - csfd);
548 } else {
549 dApDA[0][2] = (DBL_MAX);
550 }
551
552 dApDA[1][0] = - rdrpr * snfd;
553 dApDA[1][1] = rdr * rdrpr * csfd;
554 if (cpa != 0.0) {
555 dApDA[1][2] = (m_r / cpa) * rdrpr * snfd;
556 } else {
557 dApDA[1][2] = (DBL_MAX);
558 }
559
560 dApDA[2][2] = 1.0;
561
562 dApDA[3][0] = m_r * rdrpr * tnl * snfd;
563 dApDA[3][1] = m_r * tnl * (1.0 - rdr * rdrpr * csfd);
564 if (cpa != 0.0) {
565 dApDA[3][2] = (m_r / cpa) * tnl * (phid - m_r * rdrpr * snfd);
566 } else {
567 dApDA[3][2] = (DBL_MAX);
568 }
569 dApDA[3][3] = 1.0;
570 dApDA[3][4] = - m_r * phid;
571
572 dApDA[4][4] = 1.0;
573
574 return dApDA;
575 }
576
577 CLHEP::HepMatrix
578 TRGCDCHelix::delXDelA(double phi) const
579 {
580 //
581 // Calculate Jacobian (@x/@a)
582 // CLHEP::HepVector a is helix parameters and phi is internal parameter
583 // which specifys the point to be calculated for Ex(phi).
584 //
585
586 CLHEP::HepMatrix dXDA(3, 5, 0);
587
588 const double& dr = m_ac[0];
589 const double& phi0 = m_ac[1];
590 const double& cpa = m_ac[2];
591//cnv const double & dz = m_ac[3];
592 const double& tnl = m_ac[4];
593
594 double cosf0phi = cos(phi0 + phi);
595 double sinf0phi = sin(phi0 + phi);
596
597 dXDA[0][0] = m_cp;
598 dXDA[0][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
599 if (cpa != 0.0) {
600 dXDA[0][2] = - (m_r / cpa) * (m_cp - cosf0phi);
601 } else {
602 dXDA[0][2] = (DBL_MAX);
603 }
604 // dXDA[0][3] = 0.0;
605 // dXDA[0][4] = 0.0;
606
607 dXDA[1][0] = m_sp;
608 dXDA[1][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
609 if (cpa != 0.0) {
610 dXDA[1][2] = - (m_r / cpa) * (m_sp - sinf0phi);
611 } else {
612 dXDA[1][2] = (DBL_MAX);
613 }
614 // dXDA[1][3] = 0.0;
615 // dXDA[1][4] = 0.0;
616
617 // dXDA[2][0] = 0.0;
618 // dXDA[2][1] = 0.0;
619 if (cpa != 0.0) {
620 dXDA[2][2] = (m_r / cpa) * tnl * phi;
621 } else {
622 dXDA[2][2] = (DBL_MAX);
623 }
624 dXDA[2][3] = 1.0;
625 dXDA[2][4] = - m_r * phi;
626
627 return dXDA;
628 }
629
630
631
632 CLHEP::HepMatrix
633 TRGCDCHelix::delMDelA(double phi) const
634 {
635 //
636 // Calculate Jacobian (@m/@a)
637 // CLHEP::HepVector a is helix parameters and phi is internal parameter.
638 // CLHEP::HepVector m is momentum.
639 //
640
641 CLHEP::HepMatrix dMDA(3, 5, 0);
642
643 const double& phi0 = m_ac[1];
644 const double& cpa = m_ac[2];
645 const double& tnl = m_ac[4];
646
647 double cosf0phi = cos(phi0 + phi);
648 double sinf0phi = sin(phi0 + phi);
649
650 double rho;
651 if (cpa != 0.)rho = 1. / cpa;
652 else rho = (DBL_MAX);
653
654 double charge = 1.;
655 if (cpa < 0.)charge = -1.;
656
657 dMDA[0][1] = -fabs(rho) * cosf0phi;
658 dMDA[0][2] = charge * rho * rho * sinf0phi;
659
660 dMDA[1][1] = -fabs(rho) * sinf0phi;
661 dMDA[1][2] = -charge * rho * rho * cosf0phi;
662
663 dMDA[2][2] = -charge * rho * rho * tnl;
664 dMDA[2][4] = fabs(rho);
665
666 return dMDA;
667 }
668
669
670 CLHEP::HepMatrix
671 TRGCDCHelix::del4MDelA(double phi, double mass) const
672 {
673
674 //
675 // Calculate Jacobian (@4m/@a)
676 // CLHEP::HepVector a is helix parameters and phi is internal parameter.
677 // CLHEP::HepVector 4m is 4 momentum.
678 //
679
680 CLHEP::HepMatrix d4MDA(4, 5, 0);
681
682 double phi0 = m_ac[1];
683 double cpa = m_ac[2];
684 double tnl = m_ac[4];
685
686 double cosf0phi = cos(phi0 + phi);
687 double sinf0phi = sin(phi0 + phi);
688
689 double rho;
690 if (cpa != 0.)rho = 1. / cpa;
691 else rho = (DBL_MAX);
692
693 double charge = 1.;
694 if (cpa < 0.)charge = -1.;
695
696 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
697
698 d4MDA[0][1] = -fabs(rho) * cosf0phi;
699 d4MDA[0][2] = charge * rho * rho * sinf0phi;
700
701 d4MDA[1][1] = -fabs(rho) * sinf0phi;
702 d4MDA[1][2] = -charge * rho * rho * cosf0phi;
703
704 d4MDA[2][2] = -charge * rho * rho * tnl;
705 d4MDA[2][4] = fabs(rho);
706
707 if (cpa != 0.0 && E != 0.0) {
708 d4MDA[3][2] = (-1. - tnl * tnl) / (cpa * cpa * cpa * E);
709 d4MDA[3][4] = tnl / (cpa * cpa * E);
710 } else {
711 d4MDA[3][2] = (DBL_MAX);
712 d4MDA[3][4] = (DBL_MAX);
713 }
714 return d4MDA;
715 }
716
717
718 CLHEP::HepMatrix
719 TRGCDCHelix::del4MXDelA(double phi, double mass) const
720 {
721
722 //
723 // Calculate Jacobian (@4mx/@a)
724 // CLHEP::HepVector a is helix parameters and phi is internal parameter.
725 // CLHEP::HepVector 4xm is 4 momentum and position.
726 //
727
728 CLHEP::HepMatrix d4MXDA(7, 5, 0);
729
730 const double& dr = m_ac[0];
731 const double& phi0 = m_ac[1];
732 const double& cpa = m_ac[2];
733//cnv const double & dz = m_ac[3];
734 const double& tnl = m_ac[4];
735
736 double cosf0phi = cos(phi0 + phi);
737 double sinf0phi = sin(phi0 + phi);
738
739 double rho;
740 if (cpa != 0.)rho = 1. / cpa;
741 else rho = (DBL_MAX);
742
743 double charge = 1.;
744 if (cpa < 0.)charge = -1.;
745
746 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
747
748 d4MXDA[0][1] = - fabs(rho) * cosf0phi;
749 d4MXDA[0][2] = charge * rho * rho * sinf0phi;
750
751 d4MXDA[1][1] = - fabs(rho) * sinf0phi;
752 d4MXDA[1][2] = - charge * rho * rho * cosf0phi;
753
754 d4MXDA[2][2] = - charge * rho * rho * tnl;
755 d4MXDA[2][4] = fabs(rho);
756
757 if (cpa != 0.0 && E != 0.0) {
758 d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
759 d4MXDA[3][4] = tnl / (cpa * cpa * E);
760 } else {
761 d4MXDA[3][2] = (DBL_MAX);
762 d4MXDA[3][4] = (DBL_MAX);
763 }
764
765 d4MXDA[4][0] = m_cp;
766 d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
767 if (cpa != 0.0) {
768 d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
769 } else {
770 d4MXDA[4][2] = (DBL_MAX);
771 }
772
773 d4MXDA[5][0] = m_sp;
774 d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
775 if (cpa != 0.0) {
776 d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);
777
778 d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
779 } else {
780 d4MXDA[5][2] = (DBL_MAX);
781
782 d4MXDA[6][2] = (DBL_MAX);
783 }
784
785 d4MXDA[6][3] = 1.;
786 d4MXDA[6][4] = - m_r * phi;
787
788 return d4MXDA;
789 }
790
791 void
792 TRGCDCHelix::ignoreErrorMatrix(void)
793 {
794 m_matrixValid = false;
795 m_Ea *= 0.;
796 }
797
798 void
799 TRGCDCHelix::debugPrint(void) const
800 {
801
802 const double dr = m_a[0];
803 const double phi0 = m_a[1];
804 const double cpa = m_a[2];
805 const double dz = m_a[3];
806 const double tnl = m_a[4];
807 if (ms_print_debug) {
808 std::cout << "TRGCDCHelix::dr = " << dr << " phi0 = " << phi0 << " cpa = " << cpa
809 << " dz = " << dz << " tnl = " << tnl << std::endl;
810 std::cout << " pivot = " << m_pivot << std::endl;
811 }
812 }
813
814 void
815 TRGCDCHelix::checkValid(void)
816 {
817 if (!ms_check_range) return;
818//cnv const double adr = abs(m_a[0]);
819 const double adr = fabs(m_a[0]);
820//cnv const double acpa = abs(m_a[2]);
821 const double acpa = fabs(m_a[2]);
822 if (!(adr >= ms_amin[0] && adr <= ms_amax[0])) {
823 m_helixValid = false;
824 } else if (!(acpa >= ms_amin[2] && acpa <= ms_amax[2])) {
825 m_helixValid = false;
826 } else {
827 m_helixValid = true;
828 }
829 if (!m_helixValid) {
830 if (m_a[0] != 0.0 || m_a[1] != 0.0 || m_a[2] != 0.0 ||
831 m_a[3] != 0.0 || m_a[4] != 0.0)
832 std::cout << "something wrong" << std::endl;
833 }
834
835 }
836
838} // namespace Belle2
839
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36
Belle2::TRGCDCHelix
TRGCDCHelix parameter class.
Definition: Helix.h:34
Belle2::TRGCDCHelix::m_center
HepGeom::Point3D< double > m_center
caches
Definition: Helix.h:222
Belle2::TRGCDCHelix::m_a
CLHEP::HepVector m_a
a HepVector parameter
Definition: Helix.h:216
Belle2::TRGCDCHelix::m_Ea
CLHEP::HepSymMatrix m_Ea
Ea HepSymMatrix parameter.
Definition: Helix.h:218
Belle2::TRGCDCHelix::m_pivot
HepGeom::Point3D< double > m_pivot
pivot
Definition: Helix.h:214
Belle2::TRGCDCHelix::m_matrixValid
bool m_matrixValid
matrix validity
Definition: Helix.h:206
Belle2::TRGCDCHelix::m_ac
double m_ac[5]
caches
Definition: Helix.h:232
Belle2::TRGCDCHelix::m_helixValid
bool m_helixValid
helix validity
Definition: Helix.h:208
Belle2::TRGCDCHelix::m_cp
double m_cp
caches
Definition: Helix.h:224
Belle2::TRGCDCHelix::m_alpha
double m_alpha
alpha parameter
Definition: Helix.h:212
Belle2::TRGCDCHelix::m_sp
double m_sp
caches
Definition: Helix.h:226
Belle2::TRGCDCHelix::m_bField
double m_bField
magnetic field
Definition: Helix.h:210
Belle2::TRGCDCHelix::m_r
double m_r
caches
Definition: Helix.h:230
Belle2::TRGCDCHelix::m_pt
double m_pt
caches
Definition: Helix.h:228
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::TRGCDCHelix::set
void set(const HepGeom::Point3D< double > &pivot, const CLHEP::HepVector &a, const CLHEP::HepSymMatrix &Ea)
sets helix pivot position, parameters, and error matrix.
Definition: Helix.cc:396
Belle2::TRGCDCHelix::ms_amax
static CLHEP::HepVector ms_amax
limits for helix parameters
Definition: Helix.h:165
Belle2::TRGCDCHelix::ms_print_debug
static bool ms_print_debug
print debug info or not
Definition: Helix.h:169
Belle2::TRGCDCHelix::operator=
TRGCDCHelix & operator=(const TRGCDCHelix &)
Copy operator.
Definition: Helix.cc:409
Belle2::TRGCDCHelix::checkValid
void checkValid(void)
check validity
Definition: Helix.cc:815
Belle2::TRGCDCHelix::x
HepGeom::Point3D< double > x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Helix.cc:135
Belle2::TRGCDCHelix::set_limits
static void set_limits(const CLHEP::HepVector &a_min, const CLHEP::HepVector &a_max)
set limits for helix parameters
Definition: Helix.cc:61
Belle2::TRGCDCHelix::ignoreErrorMatrix
void ignoreErrorMatrix(void)
unsets error matrix. Error calculations will be ignored after this function call until an error matri...
Definition: Helix.cc:792
Belle2::TRGCDCHelix::Ea
const CLHEP::HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Helix.h:320
Belle2::M_PI8
const double M_PI8
8*PI
Definition: Helix.cc:39
Belle2::TRGCDCHelix::set_exception
static bool set_exception(bool)
set to throw exception or not
Definition: Helix.cc:50
Belle2::M_PI2
const double M_PI2
2*PI
Definition: Helix.cc:31
Belle2::TRGCDCHelix::updateCache
void updateCache(void)
update Caches
Definition: Helix.cc:436
Belle2::TRGCDCHelix::tanl
double tanl(void) const
returns tanl.
Definition: Helix.h:299
Belle2::TRGCDCHelix::phi0
double phi0(void) const
returns phi0.
Definition: Helix.h:278
Belle2::TRGCDCHelix::ms_throw_exception
static bool ms_throw_exception
throw exception or not
Definition: Helix.h:171
Belle2::TRGCDCHelix::a
const CLHEP::HepVector & a(void) const
returns helix parameters.
Definition: Helix.h:313
Belle2::TRGCDCHelix::ms_check_range
static bool ms_check_range
range in checked or not
Definition: Helix.h:167
Belle2::TRGCDCHelix::del4MXDelA
CLHEP::HepMatrix del4MXDelA(double phi, double mass) const
Mathmatical functions.
Definition: Helix.cc:719
Belle2::TRGCDCHelix::dr
double dr(void) const
returns dr.
Definition: Helix.h:271
Belle2::TRGCDCHelix::debugPrint
void debugPrint(void) const
print debug info
Definition: Helix.cc:799
Belle2::TRGCDCHelix::delApDelA
CLHEP::HepMatrix delApDelA(const CLHEP::HepVector &ap) const
Mathmatical functions.
Definition: Helix.cc:509
Belle2::TRGCDCHelix::dz
double dz(void) const
returns dz.
Definition: Helix.h:292
Belle2::TRGCDCHelix::ConstantAlpha
static const double ConstantAlpha
Constant alpha for uniform field.
Definition: Helix.h:202
Belle2::TRGCDCHelix::invalidhelix
static const std::string invalidhelix
string of invalid helix
Definition: Helix.h:235
Belle2::TRGCDCHelix::~TRGCDCHelix
virtual ~TRGCDCHelix()
Destructor.
Definition: Helix.cc:130
Belle2::TRGCDCHelix::ms_amin
static CLHEP::HepVector ms_amin
limits for helix parameters
Definition: Helix.h:163
Belle2::TRGCDCHelix::delXDelA
CLHEP::HepMatrix delXDelA(double phi) const
Mathmatical functions.
Definition: Helix.cc:578
Belle2::TRGCDCHelix::del4MDelA
CLHEP::HepMatrix del4MDelA(double phi, double mass) const
Mathmatical functions.
Definition: Helix.cc:671
Belle2::TRGCDCHelix::delMDelA
CLHEP::HepMatrix delMDelA(double phi) const
Mathmatical functions.
Definition: Helix.cc:633
Belle2::TRGCDCHelix::kappa
double kappa(void) const
returns kappa.
Definition: Helix.h:285
Belle2::M_PI4
const double M_PI4
4*PI
Definition: Helix.cc:35
Belle2::TRGCDCHelix::set_print
static bool set_print(bool)
set to print debug info or not
Definition: Helix.cc:55
Belle2::TRGCDCHelix::momentum
CLHEP::Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Helix.cc:194
Belle2::TRGCDCHelix::pivot
const HepGeom::Point3D< double > & pivot(void) const
returns pivot position.
Definition: Helix.h:250
Belle2::TRGCDCHelix::TRGCDCHelix
TRGCDCHelix(const HepGeom::Point3D< double > &pivot, const CLHEP::HepVector &a, const CLHEP::HepSymMatrix &Ea)
Constructor with pivot, helix parameter a, and its error matrix.
Definition: Helix.cc:70
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17