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