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// $Id: Helix.cc 10002 2007-02-26 06:56:17Z katayama $
11//
12// The local helix-parameter shorthands below deliberately carry the same names
13// as the accessors returning them.
14// cppcheck-suppress-file shadowFunction
15//
16// $Log$
17// Revision 1.19 2002/01/03 11:05:06 katayama
18// Point3D and other header files are cleaned
19//
20// Revision 1.18 2001/05/07 21:06:37 yiwasaki
21// <float.h> included for linux
22//
23// Revision 1.17 2001/04/25 02:55:39 yiwasaki
24// cache m_ac[5] added
25//
26// Revision 1.16 2001/04/18 12:09:19 katayama
27// optimized
28//
29// Revision 1.15 2000/01/26 10:22:53 yiwasaki
30// copy operator bug fix(reported by M.Yokoyama)
31//
32// Revision 1.14 1999/11/23 10:29:22 yiwasaki
33// static cosnt double Helix::ConstantAlpha added
34//
35// Revision 1.13 1999/06/15 01:49:41 yiwasaki
36// constructor bug fix
37//
38// Revision 1.12 1999/06/15 01:44:53 yiwasaki
39// minor change
40//
41// Revision 1.11 1999/06/06 07:34:27 yiwasaki
42// i/o functions to set/get mag. field added
43//
44// Revision 1.10 1999/05/11 23:26:59 yiwasaki
45// option to ignore error calculations added
46//
47// Revision 1.9 1998/11/06 08:51:19 katayama
48// protect 0div
49//
50// Revision 1.8 1998/07/08 04:35:50 jtanaka
51// add some members to the constructor by Iwasaki-san and Ozaki-san
52//
53// Revision 1.7 1998/06/18 10:27:20 katayama
54// Added several new functions from Tanaka san
55//
56// Revision 1.6 1998/02/24 01:07:59 yiwasaki
57// minor fix
58//
59// Revision 1.5 1998/02/24 00:31:52 yiwasaki
60// bug fix
61//
62// Revision 1.4 1998/02/20 02:22:23 yiwasaki
63// for new helix class
64//
65// Revision 1.3 1998/01/22 08:38:07 hitoshi
66// Fixed a bug in fid cal. (found at Osaka).
67//
68// Revision 1.2 1997/09/19 00:35:52 katayama
69// Added Id and Log
70//
71//
72//
73// Class Helix
74//
75// Author Date comments
76// Y.Ohnishi 03/01/1997 original version
77// Y.Ohnishi 06/03/1997 updated
78// Y.Iwasaki 17/02/1998 BFILED removed, func. name changed, func. added
79// J.Tanaka 06/12/1998 add some utilities.
80// Y.Iwasaki 07/07/1998 cache added to speed up
81//
82
83#include <math.h>
84#include <float.h>
85
86#include <cdc/simulation/Helix.h>
87#include <CLHEP/Matrix/Matrix.h>
88
89using namespace CLHEP;
90using namespace Belle2;
91using namespace CDC;
92
93// const double
94// Helix::m_BFIELD = 15.0; // KG
95// const double
96// Helix::m_ALPHA = 222.376063; // = 10000. / 2.99792458 / BFIELD
97
98const double
99M_PI2 = 2. * M_PI;
100
101const double
102M_PI4 = 4. * M_PI;
103
104const double
105M_PI8 = 8. * M_PI;
106
107const double
108Helix::ConstantAlpha = 222.376063;
109
110const std::string Helix::invalidhelix("Invalid Helix");
111HepVector Helix::ms_amin(5, 0), Helix::ms_amax(5, 0);
112bool Helix::ms_check_range(false);
113bool Helix::ms_throw_exception(false);
114bool Helix::ms_print_debug(false);
115
116bool Helix::set_exception(bool t)
117{
118 return ms_throw_exception = t;
119}
120
121bool Helix::set_print(bool t)
122{
123 return ms_print_debug = t;
124}
125
126
127void Helix::set_limits(const HepVector& a_min, const HepVector& a_max)
128{
129 if (a_min.num_row() != 5 || a_max.num_row() != 5) return;
130 ms_amin = a_min;
131 ms_amax = a_max;
132 ms_check_range = true;
133}
134
135
136Helix::Helix(const HepPoint3D& pivot,
137 const HepVector& a,
138 const HepSymMatrix& Ea)
139 : m_matrixValid(true),
140 m_helixValid(false),
141 m_bField(15.0),
142 m_alpha(222.376063),
143 m_pivot(pivot),
144 m_a(a),
145 m_Ea(Ea)
146{
147 // m_alpha = 10000. / 2.99792458 / m_bField;
148 // m_alpha = 222.376063;
149 if (m_a.num_row() == 5 && m_Ea.num_row() == 5) {
150 updateCache();
151 }
152}
153
154Helix::Helix(const HepPoint3D& pivot,
155 const HepVector& a)
156 : m_matrixValid(false),
157 m_helixValid(false),
158 m_bField(15.0),
159 m_alpha(222.376063),
160 m_pivot(pivot),
161 m_a(a),
162 m_Ea(HepSymMatrix(5, 0))
163{
164 // m_alpha = 222.376063;
165 if (m_a.num_row() == 5) {
166 updateCache();
167 }
168}
169
170Helix::Helix(const HepPoint3D& position,
171 const Hep3Vector& momentum,
172 double charge)
173 : m_matrixValid(false),
174 m_helixValid(false),
175 m_bField(15.0),
176 m_alpha(222.376063),
177 m_pivot(position),
178 m_a(HepVector(5, 0)),
179 m_Ea(HepSymMatrix(5, 0))
180{
181 m_a[0] = 0.;
182 m_a[3] = 0.;
183 double perp(momentum.perp());
184 if (perp != 0.0) {
185 m_a[1] = fmod(atan2(- momentum.x(), momentum.y())
186 + M_PI4, M_PI2);
187 m_a[2] = charge / perp;
188 m_a[4] = momentum.z() / perp;
189 } else {
190 m_a[2] = charge * (DBL_MAX);
191 }
192 // m_alpha = 222.376063;
193 updateCache();
194}
195
196Helix::~Helix()
197{
198}
199
200HepPoint3D
201Helix::x(double phi) const
202{
203 DEBUG_HELIX;
204 //
205 // Calculate position (x,y,z) along helix.
206 //
207 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
208 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
209 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
210 //
211
212 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
213 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
214 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
215
216 return HepPoint3D(x, y, z);
217}
218
219double*
220Helix::x(double phi, double p[3]) const
221{
222 DEBUG_HELIX;
223 //
224 // Calculate position (x,y,z) along helix.
225 //
226 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
227 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
228 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
229 //
230
231 p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
232 p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
233 p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
234
235 return p;
236}
237
238HepPoint3D
239Helix::x(double phi, HepSymMatrix& Ex) const
240{
241 DEBUG_HELIX;
242
243 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
244 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
245 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
246
247 //
248 // Calculate position error matrix.
249 // Ex(phi) = (@x/@a)(Ea)(@x/@a)^T, phi is deflection angle to specify the
250 // point to be calculated.
251 //
252 // HepMatrix dXDA(3, 5, 0);
253 // dXDA = delXDelA(phi);
254 // Ex.assign(dXDA * m_Ea * dXDA.T());
255
256 if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
257 else Ex = m_Ea;
258
259 return HepPoint3D(x, y, z);
260}
261
262Hep3Vector
263Helix::momentum(double phi) const
264{
265 DEBUG_HELIX;
266 //
267 // Calculate momentum.
268 //
269 // Pt = | 1/kappa | (GeV/c)
270 //
271 // Px = -Pt * sin(phi0 + phi)
272 // Py = Pt * cos(phi0 + phi)
273 // Pz = Pt * tan(lambda)
274 //
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
281 return Hep3Vector(px, py, pz);
282}
283
284Hep3Vector
285Helix::momentum(double phi, HepSymMatrix& Em) const
286{
287 DEBUG_HELIX;
288 //
289 // Calculate momentum.
290 //
291 // Pt = | 1/kappa | (GeV/c)
292 //
293 // Px = -Pt * sin(phi0 + phi)
294 // Py = Pt * cos(phi0 + phi)
295 // Pz = Pt * tan(lambda)
296 //
297
298 double pt = fabs(m_pt);
299 double px = - pt * sin(m_ac[1] + phi);
300 double py = pt * cos(m_ac[1] + phi);
301 double pz = pt * m_ac[4];
302
303 if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
304 else Em = m_Ea;
305
306 return Hep3Vector(px, py, pz);
307}
308
309HepLorentzVector
310Helix::momentum(double phi, double mass) const
311{
312 DEBUG_HELIX;
313 //
314 // Calculate momentum.
315 //
316 // Pt = | 1/kappa | (GeV/c)
317 //
318 // Px = -Pt * sin(phi0 + phi)
319 // Py = Pt * cos(phi0 + phi)
320 // Pz = Pt * tan(lambda)
321 //
322 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
323
324 double pt = fabs(m_pt);
325 double px = - pt * sin(m_ac[1] + phi);
326 double py = pt * cos(m_ac[1] + phi);
327 double pz = pt * m_ac[4];
328 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
329
330 return HepLorentzVector(px, py, pz, E);
331}
332
333
334HepLorentzVector
335Helix::momentum(double phi, double mass, HepSymMatrix& Em) const
336{
337 DEBUG_HELIX;
338 //
339 // Calculate momentum.
340 //
341 // Pt = | 1/kappa | (GeV/c)
342 //
343 // Px = -Pt * sin(phi0 + phi)
344 // Py = Pt * cos(phi0 + phi)
345 // Pz = Pt * tan(lambda)
346 //
347 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
348
349 double pt = fabs(m_pt);
350 double px = - pt * sin(m_ac[1] + phi);
351 double py = pt * cos(m_ac[1] + phi);
352 double pz = pt * m_ac[4];
353 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
354
355 if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi, mass));
356 else Em = m_Ea;
357
358 return HepLorentzVector(px, py, pz, E);
359}
360
361HepLorentzVector
362Helix::momentum(double phi,
363 double mass,
364 HepPoint3D& x,
365 HepSymMatrix& Emx) const
366{
367 DEBUG_HELIX;
368 //
369 // Calculate momentum.
370 //
371 // Pt = | 1/kappa | (GeV/c)
372 //
373 // Px = -Pt * sin(phi0 + phi)
374 // Py = Pt * cos(phi0 + phi)
375 // Pz = Pt * tan(lambda)
376 //
377 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
378
379 double pt = fabs(m_pt);
380 double px = - pt * sin(m_ac[1] + phi);
381 double py = pt * cos(m_ac[1] + phi);
382 double pz = pt * m_ac[4];
383 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
384
385 x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
386 x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
387 x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);
388
389 if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi, mass));
390 else Emx = m_Ea;
391
392 return HepLorentzVector(px, py, pz, E);
393}
394
395
396const HepPoint3D&
397Helix::pivot(const HepPoint3D& newPivot)
398{
399 DEBUG_HELIX;
400#if defined(BELLE_DEBUG)
401 try {
402#endif
403 const double& dr = m_ac[0];
404 const double& phi0 = m_ac[1];
405 const double& kappa = m_ac[2];
406 const double& dz = m_ac[3];
407 const double& tanl = m_ac[4];
408
409 double rdr = dr + m_r;
410 double phi = fmod(phi0 + M_PI4, M_PI2);
411 double csf0 = cos(phi);
412 double snf0 = (1. - csf0) * (1. + csf0);
413 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
414 if (phi > M_PI) snf0 = - snf0;
415
416 double xc = m_pivot.x() + rdr * csf0;
417 double yc = m_pivot.y() + rdr * snf0;
418 double csf, snf;
419 if (m_r != 0.0) {
420 csf = (xc - newPivot.x()) / m_r;
421 snf = (yc - newPivot.y()) / m_r;
422 double anrm = sqrt(csf * csf + snf * snf);
423 if (anrm != 0.0) {
424 csf /= anrm;
425 snf /= anrm;
426 phi = atan2(snf, csf);
427 } else {
428 csf = 1.0;
429 snf = 0.0;
430 phi = 0.0;
431 }
432 } else {
433 csf = 1.0;
434 snf = 0.0;
435 phi = 0.0;
436 }
437 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
438 if (phid > M_PI) phid = phid - M_PI2;
439 double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
440 * csf
441 + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
442 double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();
443
444 HepVector ap(5);
445 ap[0] = drp;
446 ap[1] = fmod(phi + M_PI4, M_PI2);
447 ap[2] = kappa;
448 ap[3] = dzp;
449 ap[4] = tanl;
450
451 // if (m_matrixValid) m_Ea.assign(delApDelA(ap) * m_Ea * delApDelA(ap).T());
452 if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));
453
454 m_a = ap;
455 m_pivot = newPivot;
456
457 //...Are these needed?...iw...
458 updateCache();
459 return m_pivot;
460#if defined(BELLE_DEBUG)
461 } catch (...) {
462 m_helixValid = false;
463 DEBUG_PRINT;
464 if (ms_throw_exception) throw invalidhelix;
465 }
466#endif
467 return m_pivot;
468}
469
470void
471Helix::set(const HepPoint3D& pivot,
472 const HepVector& a,
473 const HepSymMatrix& Ea)
474{
475 m_pivot = pivot;
476 m_a = a;
477 m_Ea = Ea;
478 m_matrixValid = true;
479 m_helixValid = false;
480 updateCache();
481 DEBUG_HELIX;
482}
483
484Helix&
485// cppcheck-suppress operatorEqVarError ; m_ac is copied element by element below
486Helix::operator = (const Helix& i)
487{
488 if (this == & i) return * this;
489 DEBUG_HELIX;
490
491 m_bField = i.m_bField;
492 m_alpha = i.m_alpha;
493 m_pivot = i.m_pivot;
494 m_a = i.m_a;
495 m_Ea = i.m_Ea;
496 m_matrixValid = i.m_matrixValid;
497
498 m_center = i.m_center;
499 m_cp = i.m_cp;
500 m_sp = i.m_sp;
501 m_pt = i.m_pt;
502 m_r = i.m_r;
503 m_ac[0] = i.m_ac[0];
504 m_ac[1] = i.m_ac[1];
505 m_ac[2] = i.m_ac[2];
506 m_ac[3] = i.m_ac[3];
507 m_ac[4] = i.m_ac[4];
508
509 return * this;
510}
511
512void
513Helix::updateCache(void)
514{
515
516#if defined(BELLE_DEBUG)
517 checkValid();
518 if (m_helixValid) {
519#endif
520
521 //
522 // Calculate Helix center( xc, yc ).
523 //
524 // xc = x0 + (dr + (alpha / kappa)) * cos(phi0) (cm)
525 // yc = y0 + (dr + (alpha / kappa)) * sin(phi0) (cm)
526 //
527
528 m_ac[0] = m_a[0];
529 m_ac[1] = m_a[1];
530 m_ac[2] = m_a[2];
531 m_ac[3] = m_a[3];
532 m_ac[4] = m_a[4];
533
534 m_cp = cos(m_ac[1]);
535 m_sp = sin(m_ac[1]);
536 if (m_ac[2] != 0.0) {
537 if (m_ac[2] == DBL_MAX || m_ac[2] == (-DBL_MAX)) {
538 m_pt = m_r = 0;
539 return;
540 } else {
541 m_pt = 1. / m_ac[2];
542 m_r = m_alpha / m_ac[2];
543 }
544 } else {
545 m_pt = (DBL_MAX);
546 m_r = (DBL_MAX);
547 return;
548 }
549
550 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
551 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
552 m_center.setX(x);
553 m_center.setY(y);
554 m_center.setZ(0.);
555#if defined(BELLE_DEBUG)
556 } else {
557 m_ac[0] = m_a[0];
558 m_ac[1] = m_a[1];
559 m_ac[2] = m_a[2];
560 m_ac[3] = m_a[3];
561 m_ac[4] = m_a[4];
562
563 m_cp = cos(m_ac[1]);
564 m_sp = sin(m_ac[1]);
565 if (m_ac[2] != 0.0) {
566 if (m_ac[2] == DBL_MAX || m_ac[2] == (-DBL_MAX)) {
567 m_pt = m_r = 0;
568 return;
569 } else {
570 m_pt = 1. / m_ac[2];
571 m_r = m_alpha / m_ac[2];
572 }
573 } else {
574 m_pt = (DBL_MAX);
575 m_r = (DBL_MAX);
576 return;
577 }
578
579 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
580 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
581 m_center.setX(x);
582 m_center.setY(y);
583 m_center.setZ(0.);
584 }
585#endif
586}
587
588HepMatrix
589Helix::delApDelA(const HepVector& ap) const
590{
591 DEBUG_HELIX;
592 //
593 // Calculate Jacobian (@ap/@a)
594 // Vector ap is new helix parameters and a is old helix parameters.
595 //
596
597 HepMatrix dApDA(5, 5, 0);
598
599 const double& dr = m_ac[0];
600 const double& phi0 = m_ac[1];
601 const double& cpa = m_ac[2];
602 //const double & dz = m_ac[3];
603 const double& tnl = m_ac[4];
604
605 double drp = ap[0];
606 double phi0p = ap[1];
607 //double cpap = ap[2];
608 //double dzp = ap[3];
609 //double tnlp = ap[4];
610
611 double rdr = m_r + dr;
612 double rdrpr;
613 if ((m_r + drp) != 0.0) {
614 rdrpr = 1. / (m_r + drp);
615 } else {
616 rdrpr = (DBL_MAX);
617 }
618 // double csfd = cos(phi0)*cos(phi0p) + sin(phi0)*sin(phi0p);
619 // double snfd = cos(phi0)*sin(phi0p) - sin(phi0)*cos(phi0p);
620 double csfd = cos(phi0p - phi0);
621 double snfd = sin(phi0p - phi0);
622 double phid = fmod(phi0p - phi0 + M_PI8, M_PI2);
623 if (phid > M_PI) phid = phid - M_PI2;
624
625 dApDA[0][0] = csfd;
626 dApDA[0][1] = rdr * snfd;
627 if (cpa != 0.0) {
628 dApDA[0][2] = (m_r / cpa) * (1.0 - csfd);
629 } else {
630 dApDA[0][2] = (DBL_MAX);
631 }
632
633 dApDA[1][0] = - rdrpr * snfd;
634 dApDA[1][1] = rdr * rdrpr * csfd;
635 if (cpa != 0.0) {
636 dApDA[1][2] = (m_r / cpa) * rdrpr * snfd;
637 } else {
638 dApDA[1][2] = (DBL_MAX);
639 }
640
641 dApDA[2][2] = 1.0;
642
643 dApDA[3][0] = m_r * rdrpr * tnl * snfd;
644 dApDA[3][1] = m_r * tnl * (1.0 - rdr * rdrpr * csfd);
645 if (cpa != 0.0) {
646 dApDA[3][2] = (m_r / cpa) * tnl * (phid - m_r * rdrpr * snfd);
647 } else {
648 dApDA[3][2] = (DBL_MAX);
649 }
650 dApDA[3][3] = 1.0;
651 dApDA[3][4] = - m_r * phid;
652
653 dApDA[4][4] = 1.0;
654
655 return dApDA;
656}
657
658HepMatrix
659Helix::delXDelA(double phi) const
660{
661 DEBUG_HELIX;
662 //
663 // Calculate Jacobian (@x/@a)
664 // Vector a is helix parameters and phi is internal parameter
665 // which specifys the point to be calculated for Ex(phi).
666 //
667
668 HepMatrix dXDA(3, 5, 0);
669
670 const double& dr = m_ac[0];
671 const double& phi0 = m_ac[1];
672 const double& cpa = m_ac[2];
673 //const double & dz = m_ac[3];
674 const double& tnl = m_ac[4];
675
676 double cosf0phi = cos(phi0 + phi);
677 double sinf0phi = sin(phi0 + phi);
678
679 dXDA[0][0] = m_cp;
680 dXDA[0][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
681 if (cpa != 0.0) {
682 dXDA[0][2] = - (m_r / cpa) * (m_cp - cosf0phi);
683 } else {
684 dXDA[0][2] = (DBL_MAX);
685 }
686 // dXDA[0][3] = 0.0;
687 // dXDA[0][4] = 0.0;
688
689 dXDA[1][0] = m_sp;
690 dXDA[1][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
691 if (cpa != 0.0) {
692 dXDA[1][2] = - (m_r / cpa) * (m_sp - sinf0phi);
693 } else {
694 dXDA[1][2] = (DBL_MAX);
695 }
696 // dXDA[1][3] = 0.0;
697 // dXDA[1][4] = 0.0;
698
699 // dXDA[2][0] = 0.0;
700 // dXDA[2][1] = 0.0;
701 if (cpa != 0.0) {
702 dXDA[2][2] = (m_r / cpa) * tnl * phi;
703 } else {
704 dXDA[2][2] = (DBL_MAX);
705 }
706 dXDA[2][3] = 1.0;
707 dXDA[2][4] = - m_r * phi;
708
709 return dXDA;
710}
711
712
713
714HepMatrix
715Helix::delMDelA(double phi) const
716{
717 DEBUG_HELIX;
718 //
719 // Calculate Jacobian (@m/@a)
720 // Vector a is helix parameters and phi is internal parameter.
721 // Vector m is momentum.
722 //
723
724 HepMatrix dMDA(3, 5, 0);
725
726 const double& phi0 = m_ac[1];
727 const double& cpa = m_ac[2];
728 const double& tnl = m_ac[4];
729
730 double cosf0phi = cos(phi0 + phi);
731 double sinf0phi = sin(phi0 + phi);
732
733 double rho;
734 if (cpa != 0.)rho = 1. / cpa;
735 else rho = (DBL_MAX);
736
737 double charge = 1.;
738 if (cpa < 0.)charge = -1.;
739
740 dMDA[0][1] = -fabs(rho) * cosf0phi;
741 dMDA[0][2] = charge * rho * rho * sinf0phi;
742
743 dMDA[1][1] = -fabs(rho) * sinf0phi;
744 dMDA[1][2] = -charge * rho * rho * cosf0phi;
745
746 dMDA[2][2] = -charge * rho * rho * tnl;
747 dMDA[2][4] = fabs(rho);
748
749 return dMDA;
750}
751
752
753HepMatrix
754Helix::del4MDelA(double phi, double mass) const
755{
756 DEBUG_HELIX;
757
758 //
759 // Calculate Jacobian (@4m/@a)
760 // Vector a is helix parameters and phi is internal parameter.
761 // Vector 4m is 4 momentum.
762 //
763
764 HepMatrix d4MDA(4, 5, 0);
765
766 double phi0 = m_ac[1];
767 double cpa = m_ac[2];
768 double tnl = m_ac[4];
769
770 double cosf0phi = cos(phi0 + phi);
771 double sinf0phi = sin(phi0 + phi);
772
773 double rho;
774 if (cpa != 0.)rho = 1. / cpa;
775 else rho = (DBL_MAX);
776
777 double charge = 1.;
778 if (cpa < 0.)charge = -1.;
779
780 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
781
782 d4MDA[0][1] = -fabs(rho) * cosf0phi;
783 d4MDA[0][2] = charge * rho * rho * sinf0phi;
784
785 d4MDA[1][1] = -fabs(rho) * sinf0phi;
786 d4MDA[1][2] = -charge * rho * rho * cosf0phi;
787
788 d4MDA[2][2] = -charge * rho * rho * tnl;
789 d4MDA[2][4] = fabs(rho);
790
791 if (cpa != 0.0 && E != 0.0) {
792 d4MDA[3][2] = (-1. - tnl * tnl) / (cpa * cpa * cpa * E);
793 d4MDA[3][4] = tnl / (cpa * cpa * E);
794 } else {
795 d4MDA[3][2] = (DBL_MAX);
796 d4MDA[3][4] = (DBL_MAX);
797 }
798 return d4MDA;
799}
800
801
802HepMatrix
803Helix::del4MXDelA(double phi, double mass) const
804{
805 DEBUG_HELIX;
806
807 //
808 // Calculate Jacobian (@4mx/@a)
809 // Vector a is helix parameters and phi is internal parameter.
810 // Vector 4xm is 4 momentum and position.
811 //
812
813 HepMatrix d4MXDA(7, 5, 0);
814
815 const double& dr = m_ac[0];
816 const double& phi0 = m_ac[1];
817 const double& cpa = m_ac[2];
818 //const double & dz = m_ac[3];
819 const double& tnl = m_ac[4];
820
821 double cosf0phi = cos(phi0 + phi);
822 double sinf0phi = sin(phi0 + phi);
823
824 double rho;
825 if (cpa != 0.)rho = 1. / cpa;
826 else rho = (DBL_MAX);
827
828 double charge = 1.;
829 if (cpa < 0.)charge = -1.;
830
831 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
832
833 d4MXDA[0][1] = - fabs(rho) * cosf0phi;
834 d4MXDA[0][2] = charge * rho * rho * sinf0phi;
835
836 d4MXDA[1][1] = - fabs(rho) * sinf0phi;
837 d4MXDA[1][2] = - charge * rho * rho * cosf0phi;
838
839 d4MXDA[2][2] = - charge * rho * rho * tnl;
840 d4MXDA[2][4] = fabs(rho);
841
842 if (cpa != 0.0 && E != 0.0) {
843 d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
844 d4MXDA[3][4] = tnl / (cpa * cpa * E);
845 } else {
846 d4MXDA[3][2] = (DBL_MAX);
847 d4MXDA[3][4] = (DBL_MAX);
848 }
849
850 d4MXDA[4][0] = m_cp;
851 d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
852 if (cpa != 0.0) {
853 d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
854 } else {
855 d4MXDA[4][2] = (DBL_MAX);
856 }
857
858 d4MXDA[5][0] = m_sp;
859 d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
860 if (cpa != 0.0) {
861 d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);
862
863 d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
864 } else {
865 d4MXDA[5][2] = (DBL_MAX);
866
867 d4MXDA[6][2] = (DBL_MAX);
868 }
869
870 d4MXDA[6][3] = 1.;
871 d4MXDA[6][4] = - m_r * phi;
872
873 return d4MXDA;
874}
875
876void
877Helix::ignoreErrorMatrix(void)
878{
879 m_matrixValid = false;
880 m_Ea *= 0.;
881}
882
883void
884// cppcheck-suppress unusedPrivateFunction ; only called through DEBUG_HELIX/DEBUG_PRINT
885Helix::debugPrint(void) const
886{
887
888 const double dr = m_a[0];
889 const double phi0 = m_a[1];
890 const double cpa = m_a[2];
891 const double dz = m_a[3];
892 const double tnl = m_a[4];
893 if (ms_print_debug) {
894 std::cout << "Helix::dr = " << dr << " phi0 = " << phi0 << " cpa = " << cpa
895 << " dz = " << dz << " tnl = " << tnl << std::endl;
896 std::cout << " pivot = " << m_pivot << std::endl;
897 }
898}
899
900void
901// cppcheck-suppress unusedPrivateFunction ; only called through DEBUG_HELIX/DEBUG_PRINT
902Helix::checkValid(void)
903{
904
905 if (!ms_check_range) return;
906 const double adr = fabs(m_a[0]);
907 const double acpa = fabs(m_a[2]);
908 if (!(adr >= ms_amin[0] && adr <= ms_amax[0])) {
909 m_helixValid = false;
910 } else if (!(acpa >= ms_amin[2] && acpa <= ms_amax[2])) {
911 m_helixValid = false;
912 } else {
913 m_helixValid = true;
914 }
915 if (!m_helixValid) {
916 if (m_a[0] != 0.0 || m_a[1] != 0.0 || m_a[2] != 0.0 ||
917 m_a[3] != 0.0 || m_a[4] != 0.0) {
918 DEBUG_PRINT;
919 }
920 }
921
922}
923
924// cppcheck-suppress unusedPrivateFunction ; only called through DEBUG_HELIX/DEBUG_PRINT
925void Helix::debugHelix(void) const
926{
927 if (!m_helixValid) {
928 if (ms_check_range) {DEBUG_PRINT;}
929 if (ms_throw_exception) { throw invalidhelix;}
930 }
931}
932
E
R E
internal precision of FFTW codelets
Definition DiscreteCosineTransform_31points.cc:36
Belle2::CDC::Helix
Helix parameter class.
Definition Helix.h:48
Belle2::CDC::Helix::ms_print_debug
static bool ms_print_debug
Debug option flag.
Definition Helix.h:237
Belle2::CDC::Helix::checkValid
void checkValid(void)
Check whether helix parameters is valid or not.
Definition Helix.cc:902
Belle2::CDC::Helix::delMDelA
HepMatrix delMDelA(double phi) const
DM/DA.
Definition Helix.cc:715
Belle2::CDC::Helix::del4MDelA
HepMatrix del4MDelA(double phi, double mass) const
DM4/DA.
Definition Helix.cc:754
Belle2::CDC::Helix::Helix
Helix(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
Constructor with pivot, helix parameter a, and its error matrix.
Definition Helix.cc:136
Belle2::CDC::Helix::ignoreErrorMatrix
void ignoreErrorMatrix(void)
Unsets error matrix.
Definition Helix.cc:877
Belle2::CDC::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Helix.h:356
Belle2::CDC::Helix::delXDelA
HepMatrix delXDelA(double phi) const
DX/DA.
Definition Helix.cc:659
Belle2::CDC::Helix::m_center
HepPoint3D m_center
Cache of the center position of Helix.
Definition Helix.h:305
Belle2::CDC::Helix::Ea
const HepSymMatrix & Ea(void) const
Returns error matrix.
Definition Helix.h:436
Belle2::CDC::Helix::~Helix
virtual ~Helix()
Destructor.
Definition Helix.cc:196
Belle2::CDC::Helix::set_exception
static bool set_exception(bool)
set exception
Definition Helix.cc:116
Belle2::CDC::Helix::del4MXDelA
HepMatrix del4MXDelA(double phi, double mass) const
DMX4/DA.
Definition Helix.cc:803
Belle2::CDC::Helix::updateCache
void updateCache(void)
updateCache
Definition Helix.cc:513
Belle2::CDC::Helix::m_matrixValid
bool m_matrixValid
True: matrix valid, False: matrix not valid.
Definition Helix.h:288
Belle2::CDC::Helix::m_ac
double m_ac[5]
Cache of the helix parameter.
Definition Helix.h:315
Belle2::CDC::Helix::tanl
double tanl(void) const
Return helix parameter tangent lambda.
Definition Helix.h:412
Belle2::CDC::Helix::m_helixValid
bool m_helixValid
True: helix valid, False: helix not valid.
Definition Helix.h:290
Belle2::CDC::Helix::phi0
double phi0(void) const
Return helix parameter phi0.
Definition Helix.h:388
Belle2::CDC::Helix::ms_throw_exception
static bool ms_throw_exception
Throw exception flag.
Definition Helix.h:239
Belle2::CDC::Helix::ms_check_range
static bool ms_check_range
Check the helix parameter's range.
Definition Helix.h:235
Belle2::CDC::Helix::set
void set(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
Sets helix pivot position, parameters, and error matrix.
Definition Helix.cc:471
Belle2::CDC::Helix::a
const HepVector & a(void) const
Returns helix parameters.
Definition Helix.h:428
Belle2::CDC::Helix::dr
double dr(void) const
Return helix parameter dr.
Definition Helix.h:380
Belle2::CDC::Helix::debugHelix
void debugHelix(void) const
Debug Helix.
Definition Helix.cc:925
Belle2::CDC::Helix::debugPrint
void debugPrint(void) const
Print the helix parameters to stdout.
Definition Helix.cc:885
Belle2::CDC::Helix::ms_amin
static HepVector ms_amin
minimum limit of Helix parameter a
Definition Helix.h:231
Belle2::CDC::Helix::delApDelA
HepMatrix delApDelA(const HepVector &ap) const
DAp/DA.
Definition Helix.cc:589
Belle2::CDC::Helix::m_a
HepVector m_a
Helix parameter.
Definition Helix.h:298
Belle2::CDC::Helix::dz
double dz(void) const
Return helix parameter dz.
Definition Helix.h:404
Belle2::CDC::Helix::ConstantAlpha
static const double ConstantAlpha
Constant alpha for uniform field.
Definition Helix.h:283
Belle2::CDC::Helix::ms_amax
static HepVector ms_amax
maxiimum limit of Helix parameter a
Definition Helix.h:233
Belle2::CDC::Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition Helix.cc:263
Belle2::CDC::Helix::invalidhelix
static const std::string invalidhelix
String "Invalid Helix".
Definition Helix.h:318
Belle2::CDC::Helix::operator=
Helix & operator=(const Helix &)
Copy operator.
Definition Helix.cc:486
Belle2::CDC::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Helix.cc:201
Belle2::CDC::Helix::m_cp
double m_cp
Cache of the cos phi0.
Definition Helix.h:307
Belle2::CDC::Helix::m_alpha
double m_alpha
10000.0/(speed of light)/B.
Definition Helix.h:294
Belle2::CDC::Helix::m_Ea
HepSymMatrix m_Ea
Error of the helix parameter.
Definition Helix.h:300
Belle2::CDC::Helix::m_sp
double m_sp
Cache of the sin phi0.
Definition Helix.h:309
Belle2::CDC::Helix::m_pivot
HepPoint3D m_pivot
Pivot.
Definition Helix.h:296
Belle2::CDC::Helix::m_bField
double m_bField
Magnetic field, assuming uniform Bz in the unit of kG.
Definition Helix.h:292
Belle2::CDC::Helix::m_r
double m_r
Cache of the r.
Definition Helix.h:313
Belle2::CDC::Helix::kappa
double kappa(void) const
Return helix parameter kappa.
Definition Helix.h:396
Belle2::CDC::Helix::m_pt
double m_pt
Cache of the pt.
Definition Helix.h:311
Belle2::CDC::Helix::set_limits
static void set_limits(const HepVector &a_min, const HepVector &a_max)
set limit for parameter "a"
Definition Helix.cc:127
Belle2::CDC::Helix::set_print
static bool set_print(bool)
Set print option for debugging.
Definition Helix.cc:121
Belle2::Helix::Helix
Helix()
Constructor initializing all perigee parameters to zero.
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2::M_PI8
const double M_PI8
8*PI
Definition Helix.cc:39
Belle2::M_PI2
const double M_PI2
2*PI
Definition Helix.cc:31
Belle2::M_PI4
const double M_PI4
4*PI
Definition Helix.cc:35
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17