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