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