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