Belle II Software  release-08-01-10
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 
85 using namespace CLHEP;
86 using namespace Belle2;
87 using 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 
94 const double
95 M_PI2 = 2. * M_PI;
96 
97 const double
98 M_PI4 = 4. * M_PI;
99 
100 const double
101 M_PI8 = 8. * M_PI;
102 
103 const double
104 Helix::ConstantAlpha = 222.376063;
105 
106 const std::string Helix::invalidhelix("Invalid Helix");
107 HepVector Helix::ms_amin(5, 0), Helix::ms_amax(5, 0);
108 bool Helix::ms_check_range(false);
109 bool Helix::ms_throw_exception(false);
110 bool Helix::ms_print_debug(false);
111 
112 bool Helix::set_exception(bool t)
113 {
114  return ms_throw_exception = t;
115 }
116 
117 bool Helix::set_print(bool t)
118 {
119  return ms_print_debug = t;
120 }
121 
122 
123 void 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 
132 Helix::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 
150 Helix::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 
166 Helix::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 
192 Helix::~Helix()
193 {
194 }
195 
196 HepPoint3D
197 Helix::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 
215 double*
216 Helix::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 
234 HepPoint3D
235 Helix::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 
258 Hep3Vector
259 Helix::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 
280 Hep3Vector
281 Helix::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 
305 HepLorentzVector
306 Helix::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 
330 HepLorentzVector
331 Helix::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 
357 HepLorentzVector
358 Helix::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 
392 const HepPoint3D&
393 Helix::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 
466 void
467 Helix::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 
480 Helix&
481 Helix::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 
507 void
508 Helix::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 
583 HepMatrix
584 Helix::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 
653 HepMatrix
654 Helix::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 
709 HepMatrix
710 Helix::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 
748 HepMatrix
749 Helix::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 
797 HepMatrix
798 Helix::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 
871 void
872 Helix::ignoreErrorMatrix(void)
873 {
874  m_matrixValid = false;
875  m_Ea *= 0.;
876 }
877 
878 void
879 Helix::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 
894 void
895 Helix::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 
917 void 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::Helix
Helix(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
Constructor with pivot, helix parameter a, and its error matrix.
Definition: Helix.cc:132
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::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::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::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