CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx
Go to the documentation of this file.
1//
2// $Id: Helix.cxx,v 1.8 2020/05/11 02:38:28 wangll Exp $
3//
4// Class Helix
5//
6// Author Date comments
7// Y.Ohnishi 03/01/1997 original version
8// Y.Ohnishi 06/03/1997 updated
9// Y.Iwasaki 17/02/1998 BFILED removed, func. name changed, func. added
10// J.Tanaka 06/12/1998 add some utilities.
11// Y.Iwasaki 07/07/1998 cache added to speed up
12//
13#include <iostream>
14#include <math.h>
15#include <float.h>
16#include "TrackUtil/Helix.h"
17#include "CLHEP/Matrix/Matrix.h"
18#include "GaudiKernel/StatusCode.h"
19#include "GaudiKernel/IInterface.h"
20#include "GaudiKernel/Kernel.h"
21#include "GaudiKernel/Service.h"
22#include "GaudiKernel/ISvcLocator.h"
23#include "GaudiKernel/SvcFactory.h"
24#include "GaudiKernel/IDataProviderSvc.h"
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/MsgStream.h"
27#include "GaudiKernel/SmartDataPtr.h"
28#include "GaudiKernel/IMessageSvc.h"
29
30// const double
31// Helix::m_BFIELD = 10.0; // KG
32// const double
33// Helix::m_ALPHA = 222.376063; // = 10000. / 2.99792458 / BFIELD
34// now Helix::m_ALPHA = 333.564095;
35
36const double
37M_PI2 = 2. * M_PI;
38
39const double
40M_PI4 = 4. * M_PI;
41
42const double
43M_PI8 = 8. * M_PI;
44
45const double Helix::ConstantAlpha = 333.564095;
46
47Helix::Helix()
48: m_matrixValid(true)
49{
50 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
51 if(scmgn!=StatusCode::SUCCESS) {
52 std::cout<< "Unable to open Magnetic field service"<<std::endl;
53 }
54 m_bField = -10000*(m_pmgnIMF->getReferField());//zhangr 2012-09-06 for MagnField
55 //m_bField = 10000*(m_pmgnIMF->getReferField());
56 m_alpha = 10000. / 2.99792458 / m_bField;
57
58 HepPoint3D pivot(0,0,0);
59 HepVector a(5,0);
60 HepSymMatrix Ea(5,0);
61 m_pivot = pivot;
62 m_a = a;
63 m_Ea = Ea;
64 updateCache();
65}
66
67Helix::Helix(const HepPoint3D & pivot,
68 const HepVector & a,
69 const HepSymMatrix & Ea)
70 : //m_bField(-10.0),
71 //m_alpha(-333.564095),
72 m_pivot(pivot),
73 m_a(a),
74 m_matrixValid(true),
75 m_Ea(Ea) {
76 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
77 if(scmgn!=StatusCode::SUCCESS) {
78 std::cout<< "Unable to open Magnetic field service"<<std::endl;
79 }
80 m_bField = -10000*(m_pmgnIMF->getReferField());//zhangr 2012-09-06 for MagnField
81 //m_bField = 10000*(m_pmgnIMF->getReferField());
82 m_alpha = 10000. / 2.99792458 / m_bField;
83 // m_alpha = 10000. / 2.99792458 / m_bField;
84 // m_alpha = 333.564095;
85 updateCache();
86}
87
88Helix::Helix(const HepPoint3D & pivot,
89 const HepVector & a)
90 : //m_bField(-10.0),
91 //m_alpha(-333.564095),
92 m_pivot(pivot),
93 m_a(a),
94 m_matrixValid(false),
95 m_Ea(HepSymMatrix(5,0)) {
96 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
97 if(scmgn!=StatusCode::SUCCESS) {
98 // log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
99 std::cout<< "Unable to open Magnetic field service"<<std::endl;
100 }
101 m_bField = -10000*(m_pmgnIMF->getReferField());//zhangr 2012-09-06 for MagnField
102 //m_bField = 10000*(m_pmgnIMF->getReferField());
103 m_alpha = 10000. / 2.99792458 / m_bField;
104 // m_alpha = 333.564095;
105 //cout<<"MdcFastTrakAlg:: bField,alpha: "<<m_bField<<" , "<<m_alpha<<endl;
106 updateCache();
107}
108
109Helix::Helix(const HepPoint3D & position,
110 const Hep3Vector & momentum,
111 double charge)
112 : //m_bField(-10.0),
113 //m_alpha(-333.564095),
114 m_pivot(position),
115 m_a(HepVector(5,0)),
116 m_matrixValid(false),
117 m_Ea(HepSymMatrix(5,0)) {
118 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
119 if(scmgn!=StatusCode::SUCCESS) {
120 // log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
121 std::cout<< "Unable to open Magnetic field service"<<std::endl;
122 }
123 m_bField = -10000*(m_pmgnIMF->getReferField());//zhangr 2012-09-06 for MagnField
124 //m_bField = 10000*(m_pmgnIMF->getReferField());
125 m_alpha = 10000. / 2.99792458 / m_bField;
126
127 m_a[0] = 0.;
128 m_a[1] = fmod(atan2(- momentum.x(), momentum.y())
129 + M_PI4, M_PI2);
130 m_a[3] = 0.;
131 double perp(momentum.perp());
132 if (perp != 0.0) {
133 m_a[2] = charge / perp;
134 m_a[4] = momentum.z() / perp;
135 }
136 else {
137 m_a[2] = charge * (DBL_MAX);
138 if (momentum.z() >= 0) {
139 m_a[4] = (DBL_MAX);
140 } else {
141 m_a[4] = -(DBL_MAX);
142 }
143 }
144 // m_alpha = 333.564095;
145 updateCache();
146}
147
148Helix::~Helix() {
149}
150
151HepPoint3D
152Helix::x(double phi) const {
153 //
154 // Calculate position (x,y,z) along helix.
155 //
156 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
157 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
158 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
159 //
160
161 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
162 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
163 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
164
165 return HepPoint3D(x, y, z);
166}
167
168double *
169Helix::x(double phi, double p[3]) const {
170 //
171 // Calculate position (x,y,z) along helix.
172 //
173 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
174 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
175 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
176 //
177
178 p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
179 p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
180 p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
181
182 return p;
183}
184
185HepPoint3D
186Helix::x(double phi, HepSymMatrix & Ex) const {
187 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
188 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
189 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
190
191 //
192 // Calculate position error matrix.
193 // Ex(phi) = (@x/@a)(Ea)(@x/@a)^T, phi is deflection angle to specify the
194 // point to be calcualted.
195 //
196 // HepMatrix dXDA(3, 5, 0);
197 // dXDA = delXDelA(phi);
198 // Ex.assign(dXDA * m_Ea * dXDA.T());
199
200 if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
201 else Ex = m_Ea;
202
203 return HepPoint3D(x, y, z);
204}
205
206Hep3Vector
207Helix::momentum(double phi) const {
208 //
209 // Calculate momentum.
210 //
211 // Pt = | 1/kappa | (GeV/c)
212 //
213 // Px = -Pt * sin(phi0 + phi)
214 // Py = Pt * cos(phi0 + phi)
215 // Pz = Pt * tan(lambda)
216 //
217
218 double pt = fabs(m_pt);
219 double px = - pt * sin(m_ac[1] + phi);
220 double py = pt * cos(m_ac[1] + phi);
221 double pz = pt * m_ac[4];
222
223 return Hep3Vector(px, py, pz);
224}
225
226Hep3Vector
227Helix::momentum(double phi, HepSymMatrix & Em) const {
228 //
229 // Calculate momentum.
230 //
231 // Pt = | 1/kappa | (GeV/c)
232 //
233 // Px = -Pt * sin(phi0 + phi)
234 // Py = Pt * cos(phi0 + phi)
235 // Pz = Pt * tan(lambda)
236 //
237
238 double pt = fabs(m_pt);
239 double px = - pt * sin(m_ac[1] + phi);
240 double py = pt * cos(m_ac[1] + phi);
241 double pz = pt * m_ac[4];
242
243 if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
244 else Em = m_Ea;
245
246 return Hep3Vector(px, py, pz);
247}
248
249HepLorentzVector
250Helix::momentum(double phi, double mass) const {
251 //
252 // Calculate momentum.
253 //
254 // Pt = | 1/kappa | (GeV/c)
255 //
256 // Px = -Pt * sin(phi0 + phi)
257 // Py = Pt * cos(phi0 + phi)
258 // Pz = Pt * tan(lambda)
259 //
260 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
261
262 double pt = fabs(m_pt);
263 double px = - pt * sin(m_ac[1] + phi);
264 double py = pt * cos(m_ac[1] + phi);
265 double pz = pt * m_ac[4];
266 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
267
268 return HepLorentzVector(px, py, pz, E);
269}
270
271
272HepLorentzVector
273Helix::momentum(double phi, double mass, HepSymMatrix & Em) const {
274 //
275 // Calculate momentum.
276 //
277 // Pt = | 1/kappa | (GeV/c)
278 //
279 // Px = -Pt * sin(phi0 + phi)
280 // Py = Pt * cos(phi0 + phi)
281 // Pz = Pt * tan(lambda)
282 //
283 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
284
285 double pt = fabs(m_pt);
286 double px = - pt * sin(m_ac[1] + phi);
287 double py = pt * cos(m_ac[1] + phi);
288 double pz = pt * m_ac[4];
289 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
290
291 if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi,mass));
292 else Em = m_Ea;
293
294 return HepLorentzVector(px, py, pz, E);
295}
296
297HepLorentzVector
298Helix::momentum(double phi,
299 double mass,
300 HepPoint3D & x,
301 HepSymMatrix & Emx) const {
302 //
303 // Calculate momentum.
304 //
305 // Pt = | 1/kappa | (GeV/c)
306 //
307 // Px = -Pt * sin(phi0 + phi)
308 // Py = Pt * cos(phi0 + phi)
309 // Pz = Pt * tan(lambda)
310 //
311 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
312
313 double pt = fabs(m_pt);
314 double px = - pt * sin(m_ac[1] + phi);
315 double py = pt * cos(m_ac[1] + phi);
316 double pz = pt * m_ac[4];
317 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
318
319 x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
320 x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
321 x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);
322
323 if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi,mass));
324 else Emx = m_Ea;
325
326 return HepLorentzVector(px, py, pz, E);
327}
328
329
330const HepPoint3D &
331Helix::pivot(const HepPoint3D & newPivot) {
332 const double & dr = m_ac[0];
333 const double & phi0 = m_ac[1];
334 const double & kappa = m_ac[2];
335 const double & dz = m_ac[3];
336 const double & tanl = m_ac[4];
337
338 double rdr = dr + m_r;
339 double phi = fmod(phi0 + M_PI4, M_PI2);
340 double csf0 = cos(phi);
341 double snf0 = (1. - csf0) * (1. + csf0);
342 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
343 if(phi > M_PI) snf0 = - snf0;
344
345 double xc = m_pivot.x() + rdr * csf0;
346 double yc = m_pivot.y() + rdr * snf0;
347 double csf, snf;
348 if(m_r != 0.0) {
349 csf = (xc - newPivot.x()) / m_r;
350 snf = (yc - newPivot.y()) / m_r;
351 double anrm = sqrt(csf * csf + snf * snf);
352 if(anrm != 0.0) {
353 csf /= anrm;
354 snf /= anrm;
355 phi = atan2(snf, csf);
356 } else {
357 csf = 1.0;
358 snf = 0.0;
359 phi = 0.0;
360 }
361 } else {
362 csf = 1.0;
363 snf = 0.0;
364 phi = 0.0;
365 }
366 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
367 if(phid > M_PI) phid = phid - M_PI2;
368 double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
369 * csf
370 + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
371 double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();
372
373 HepVector ap(5);
374 ap[0] = drp;
375 ap[1] = fmod(phi + M_PI4, M_PI2);
376 ap[2] = kappa;
377 ap[3] = dzp;
378 ap[4] = tanl;
379
380 // if (m_matrixValid) m_Ea.assign(delApDelA(ap) * m_Ea * delApDelA(ap).T());
381 if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));
382
383 m_a = ap;
384 m_pivot = newPivot;
385
386 //...Are these needed?...iw...
387 updateCache();
388 return m_pivot;
389}
390
391void
392Helix::set(const HepPoint3D & pivot,
393 const HepVector & a,
394 const HepSymMatrix & Ea) {
395 m_pivot = pivot;
396 m_a = a;
397 m_Ea = Ea;
398 m_matrixValid = true;
399 updateCache();
400}
401
402Helix &
403Helix::operator = (const Helix & i) {
404 if (this == & i) return * this;
405
406 m_bField = i.m_bField;
407 m_alpha = i.m_alpha;
408 m_pivot = i.m_pivot;
409 m_a = i.m_a;
410 m_Ea = i.m_Ea;
411 m_matrixValid = i.m_matrixValid;
412
413 m_center = i.m_center;
414 m_cp = i.m_cp;
415 m_sp = i.m_sp;
416 m_pt = i.m_pt;
417 m_r = i.m_r;
418 m_ac[0] = i.m_ac[0];
419 m_ac[1] = i.m_ac[1];
420 m_ac[2] = i.m_ac[2];
421 m_ac[3] = i.m_ac[3];
422 m_ac[4] = i.m_ac[4];
423
424 return * this;
425}
426
427Helix::Helix(const Helix& i)
428{
429 m_bField = i.m_bField;
430 m_alpha = i.m_alpha;
431 m_pivot = i.m_pivot;
432 m_a = i.m_a;
433 m_Ea = i.m_Ea;
434 m_matrixValid = i.m_matrixValid;
435
436 m_center = i.m_center;
437 m_cp = i.m_cp;
438 m_sp = i.m_sp;
439 m_pt = i.m_pt;
440 m_r = i.m_r;
441 m_ac[0] = i.m_ac[0];
442 m_ac[1] = i.m_ac[1];
443 m_ac[2] = i.m_ac[2];
444 m_ac[3] = i.m_ac[3];
445 m_ac[4] = i.m_ac[4];
446}
447
448
449void
450Helix::updateCache(void) {
451 //
452 // Calculate Helix center( xc, yc ).
453 //
454 // xc = x0 + (dr + (alpha / kappa)) * cos(phi0) (cm)
455 // yc = y0 + (dr + (alpha / kappa)) * sin(phi0) (cm)
456 //
457
458 m_ac[0] = m_a[0];
459 m_ac[1] = m_a[1];
460 m_ac[2] = m_a[2];
461 m_ac[3] = m_a[3];
462 m_ac[4] = m_a[4];
463
464 m_cp = cos(m_ac[1]);
465 m_sp = sin(m_ac[1]);
466 if (m_ac[2] != 0.0) {
467 m_pt = 1. / m_ac[2];
468 m_r = m_alpha / m_ac[2];
469 }
470 else {
471 m_pt = (DBL_MAX);
472 m_r = (DBL_MAX);
473 }
474
475 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
476 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
477 m_center.setX(x);
478 m_center.setY(y);
479 m_center.setZ(0.);
480}
481
482HepMatrix
483Helix::delApDelA(const HepVector & ap) const {
484 //
485 // Calculate Jacobian (@ap/@a)
486 // Vector ap is new helix parameters and a is old helix parameters.
487 //
488
489 HepMatrix dApDA(5,5,0);
490
491 const double & dr = m_ac[0];
492 const double & phi0 = m_ac[1];
493 const double & cpa = m_ac[2];
494 //const double & dz = m_ac[3];
495 const double & tnl = m_ac[4];
496
497 double drp = ap[0];
498 double phi0p = ap[1];
499 //double cpap = ap[2];
500 //double dzp = ap[3];
501 //double tnlp = ap[4];
502
503 double rdr = m_r + dr;
504 double rdrpr;
505 if ((m_r + drp) != 0.0) {
506 rdrpr = 1. / (m_r + drp);
507 } else {
508 rdrpr = (DBL_MAX);
509 }
510 // double csfd = cos(phi0)*cos(phi0p) + sin(phi0)*sin(phi0p);
511 // double snfd = cos(phi0)*sin(phi0p) - sin(phi0)*cos(phi0p);
512 double csfd = cos(phi0p - phi0);
513 double snfd = sin(phi0p - phi0);
514 double phid = fmod(phi0p - phi0 + M_PI8, M_PI2);
515 if (phid > M_PI) phid = phid - M_PI2;
516
517 dApDA[0][0] = csfd;
518 dApDA[0][1] = rdr*snfd;
519 if(cpa!=0.0) {
520 dApDA[0][2] = (m_r/cpa)*( 1.0 - csfd );
521 } else {
522 dApDA[0][2] = (DBL_MAX);
523 }
524
525 dApDA[1][0] = - rdrpr*snfd;
526 dApDA[1][1] = rdr*rdrpr*csfd;
527 if(cpa!=0.0) {
528 dApDA[1][2] = (m_r/cpa)*rdrpr*snfd;
529 } else {
530 dApDA[1][2] = (DBL_MAX);
531 }
532
533 dApDA[2][2] = 1.0;
534
535 dApDA[3][0] = m_r*rdrpr*tnl*snfd;
536 dApDA[3][1] = m_r*tnl*(1.0 - rdr*rdrpr*csfd);
537 if(cpa!=0.0) {
538 dApDA[3][2] = (m_r/cpa)*tnl*(phid - m_r*rdrpr*snfd);
539 } else {
540 dApDA[3][2] = (DBL_MAX);
541 }
542 dApDA[3][3] = 1.0;
543 dApDA[3][4] = - m_r*phid;
544
545 dApDA[4][4] = 1.0;
546
547 return dApDA;
548}
549
550HepMatrix
551Helix::delXDelA(double phi) const {
552 //
553 // Calculate Jacobian (@x/@a)
554 // Vector a is helix parameters and phi is internal parameter
555 // which specifys the point to be calculated for Ex(phi).
556 //
557
558 HepMatrix dXDA(3,5,0);
559
560 const double & dr = m_ac[0];
561 const double & phi0 = m_ac[1];
562 const double & cpa = m_ac[2];
563 //const double & dz = m_ac[3];
564 const double & tnl = m_ac[4];
565
566 double cosf0phi = cos(phi0 + phi);
567 double sinf0phi = sin(phi0 + phi);
568
569 dXDA[0][0] = m_cp;
570 dXDA[0][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
571 if(cpa!=0.0) {
572 dXDA[0][2] = - (m_r / cpa) * (m_cp - cosf0phi);
573 } else {
574 dXDA[0][2] = (DBL_MAX);
575 }
576 // dXDA[0][3] = 0.0;
577 // dXDA[0][4] = 0.0;
578
579 dXDA[1][0] = m_sp;
580 dXDA[1][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
581 if(cpa!=0.0) {
582 dXDA[1][2] = - (m_r / cpa) * (m_sp - sinf0phi);
583 } else {
584 dXDA[1][2] = (DBL_MAX);
585 }
586 // dXDA[1][3] = 0.0;
587 // dXDA[1][4] = 0.0;
588
589 // dXDA[2][0] = 0.0;
590 // dXDA[2][1] = 0.0;
591 if(cpa!=0.0) {
592 dXDA[2][2] = (m_r / cpa) * tnl * phi;
593 } else {
594 dXDA[2][2] = (DBL_MAX);
595 }
596 dXDA[2][3] = 1.0;
597 dXDA[2][4] = - m_r * phi;
598
599 return dXDA;
600}
601
602
603
604HepMatrix
605Helix::delMDelA(double phi) const {
606 //
607 // Calculate Jacobian (@m/@a)
608 // Vector a is helix parameters and phi is internal parameter.
609 // Vector m is momentum.
610 //
611
612 HepMatrix dMDA(3,5,0);
613
614 const double & phi0 = m_ac[1];
615 const double & cpa = m_ac[2];
616 const double & tnl = m_ac[4];
617
618 double cosf0phi = cos(phi0+phi);
619 double sinf0phi = sin(phi0+phi);
620
621 double rho;
622 if(cpa != 0.)rho = 1./cpa;
623 else rho = (DBL_MAX);
624
625 double charge = 1.;
626 if(cpa < 0.)charge = -1.;
627
628 dMDA[0][1] = -fabs(rho)*cosf0phi;
629 dMDA[0][2] = charge*rho*rho*sinf0phi;
630
631 dMDA[1][1] = -fabs(rho)*sinf0phi;
632 dMDA[1][2] = -charge*rho*rho*cosf0phi;
633
634 dMDA[2][2] = -charge*rho*rho*tnl;
635 dMDA[2][4] = fabs(rho);
636
637 return dMDA;
638}
639
640
641HepMatrix
642Helix::del4MDelA(double phi, double mass) const {
643 //
644 // Calculate Jacobian (@4m/@a)
645 // Vector a is helix parameters and phi is internal parameter.
646 // Vector 4m is 4 momentum.
647 //
648
649 HepMatrix d4MDA(4,5,0);
650
651 double phi0 = m_ac[1];
652 double cpa = m_ac[2];
653 double tnl = m_ac[4];
654
655 double cosf0phi = cos(phi0+phi);
656 double sinf0phi = sin(phi0+phi);
657
658 double rho;
659 if(cpa != 0.)rho = 1./cpa;
660 else rho = (DBL_MAX);
661
662 double charge = 1.;
663 if(cpa < 0.)charge = -1.;
664
665 double E = sqrt(rho*rho*(1.+tnl*tnl)+mass*mass);
666
667 d4MDA[0][1] = -fabs(rho)*cosf0phi;
668 d4MDA[0][2] = charge*rho*rho*sinf0phi;
669
670 d4MDA[1][1] = -fabs(rho)*sinf0phi;
671 d4MDA[1][2] = -charge*rho*rho*cosf0phi;
672
673 d4MDA[2][2] = -charge*rho*rho*tnl;
674 d4MDA[2][4] = fabs(rho);
675
676 if (cpa != 0.0 && E != 0.0) {
677 d4MDA[3][2] = (-1.-tnl*tnl)/(cpa*cpa*cpa*E);
678 d4MDA[3][4] = tnl/(cpa*cpa*E);
679 } else {
680 d4MDA[3][2] = (DBL_MAX);
681 d4MDA[3][4] = (DBL_MAX);
682 }
683 return d4MDA;
684}
685
686
687HepMatrix
688Helix::del4MXDelA(double phi, double mass) const {
689 //
690 // Calculate Jacobian (@4mx/@a)
691 // Vector a is helix parameters and phi is internal parameter.
692 // Vector 4xm is 4 momentum and position.
693 //
694
695 HepMatrix d4MXDA(7,5,0);
696
697 const double & dr = m_ac[0];
698 const double & phi0 = m_ac[1];
699 const double & cpa = m_ac[2];
700 //const double & dz = m_ac[3];
701 const double & tnl = m_ac[4];
702
703 double cosf0phi = cos(phi0+phi);
704 double sinf0phi = sin(phi0+phi);
705
706 double rho;
707 if(cpa != 0.)rho = 1./cpa;
708 else rho = (DBL_MAX);
709
710 double charge = 1.;
711 if(cpa < 0.)charge = -1.;
712
713 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
714
715 d4MXDA[0][1] = - fabs(rho) * cosf0phi;
716 d4MXDA[0][2] = charge * rho * rho * sinf0phi;
717
718 d4MXDA[1][1] = - fabs(rho) * sinf0phi;
719 d4MXDA[1][2] = - charge * rho * rho * cosf0phi;
720
721 d4MXDA[2][2] = - charge * rho * rho * tnl;
722 d4MXDA[2][4] = fabs(rho);
723
724 if (cpa != 0.0 && E != 0.0) {
725 d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
726 d4MXDA[3][4] = tnl / (cpa * cpa * E);
727 } else {
728 d4MXDA[3][2] = (DBL_MAX);
729 d4MXDA[3][4] = (DBL_MAX);
730 }
731
732 d4MXDA[4][0] = m_cp;
733 d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
734 if (cpa != 0.0) {
735 d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
736 } else {
737 d4MXDA[4][2] = (DBL_MAX);
738 }
739
740 d4MXDA[5][0] = m_sp;
741 d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
742 if (cpa != 0.0) {
743 d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);
744
745 d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
746 } else {
747 d4MXDA[5][2] = (DBL_MAX);
748
749 d4MXDA[6][2] = (DBL_MAX);
750 }
751
752 d4MXDA[6][3] = 1.;
753 d4MXDA[6][4] = - m_r * phi;
754
755 return d4MXDA;
756}
757
758void
759Helix::ignoreErrorMatrix(void) {
760 m_matrixValid = false;
761 m_Ea *= 0.;
762}
763
764double Helix::IntersectCylinder(double r) const
765{
766 double m_rad = radius();
767 double l = center().perp();
768
769 double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
770
771 if(cosPhi < -1 || cosPhi > 1) return 0;
772
773 double dPhi = center().phi() - acos(cosPhi) - phi0();
774
775 if(dPhi < -M_PI) dPhi += 2 * M_PI;
776
777 return dPhi;
778}
779
780double Helix::flightArc(HepPoint3D& hit) const
781{
782 /*
783 double rc = center().perp();
784 //double l = center().perp(hit);
785 double l = sqrt( (center().x()-hit.x()) * (center().x()-hit.x()) + (center().y()-hit.y()) * (center().y()-hit.y()) );
786 double rHit = hit.perp();
787 double cosPhi = (rc*rc + l * l - rHit * rHit) / (2 * rc * l);
788 //double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
789 if(cosPhi < -1 || cosPhi > 1) return 0;
790 double dPhi = acos(cosPhi);// [0, pi)
791 //double phiTrkFlt = IntersectCylinder(rHit);
792 //double arcLength = rc*dPhi;
793 double xHit = hit.getHitPoint().x();
794 double yHit = hit.getHitPoint().y();
795 double xCenter = center().x();
796 double yCenter = center().y();
797 double leftOrRight = xHit*yCenter - xCenter*yHit;
798 */
799 double dphi = dPhi(hit);
800 double arcLength = fabs(m_r*dphi);
801 return arcLength;
802}
803
804double Helix::flightArc(double r) const
805{
806 double phi_turn = IntersectCylinder(r);
807 double arcLength = m_r*phi_turn;
808 arcLength = fabs(arcLength);
809 return arcLength;
810}
811
812double Helix::flightLength(HepPoint3D& hit) const
813{
814 /*
815 double rc = center().perp();
816 //double l = center().perp(hit);
817 double l = sqrt( (center().x()-hit.x()) * (center().x()-hit.x()) + (center().y()-hit.y()) * (center().y()-hit.y()) );
818 double rHit = hit.perp();
819 double cosPhi = (rc*rc + l * l - rHit * rHit) / (2 * rc * l);
820 //double cosPhi = (m_rad * m_rad + l * l - r * r) / (2 * m_rad * l);
821 if(cosPhi < -1 || cosPhi > 1) return 0;
822 double dPhi = acos(cosPhi);
823 //double dPhi = center().phi() - acos(cosPhi) - phi0();
824 //double phiTrkFlt = IntersectCylinder(rHit);
825 double arcLength = rc*dPhi;
826 */
827
828 double arcLength = flightArc(hit);
829
830 double cosLambda = 1/sqrt(1+m_ac[4]*m_ac[4]);
831 double flightLength = arcLength/cosLambda;
832 return flightLength;
833}
834
835double Helix::dPhi(HepPoint3D& hit) const
836{
837 double isClockWise = m_r/fabs(m_r);
838 double x_cw = center().x()-hit.x();
839 double y_cw = center().y()-hit.y();
840 x_cw = isClockWise*x_cw;
841 y_cw = isClockWise*y_cw;
842 double phi_cw = atan2(y_cw,x_cw);
843 double dphi = phi_cw-phi0();
844 if(isClockWise>0)
845 {
846 while(dphi>0) dphi-=2*M_PI;
847 while(dphi<-2*M_PI) dphi+=2*M_PI;
848 }
849 else {
850 while(dphi<0) dphi+=2*M_PI;
851 while(dphi>2*M_PI) dphi-=2*M_PI;
852 }
853 return dphi;
854}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
momentum
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
Definition: BesBdkRc/BesBdkRc-00-00-02/src/fortran/HepEvt.h:23
mass
double mass
Definition: CosmicGenerator.cxx:138
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
DBL_MAX
#define DBL_MAX
Definition: KalFitAlg.h:13
m_pt
NTuple::Item< double > m_pt
Definition: MdcHistItem.h:76
M_PI
#define M_PI
Definition: TConstant.h:4
Helix::delApDelA
HepMatrix delApDelA(const HepVector &ap) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:483
Helix::flightLength
double flightLength(HepPoint3D &hit) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:812
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:212
Helix::set
void set(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
sets helix pivot position, parameters, and error matrix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:392
Helix::ignoreErrorMatrix
void ignoreErrorMatrix(void)
unsets error matrix. Error calculations will be ignored after this function call until an error matri...
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:759
Helix::operator=
Helix & operator=(const Helix &)
Copy operator.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:403
Helix::del4MXDelA
HepMatrix del4MXDelA(double phi, double mass) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:688
Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:207
Helix::ConstantAlpha
static const double ConstantAlpha
Constant alpha for uniform field.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:180
Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:278
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:152
Helix::del4MDelA
HepMatrix del4MDelA(double phi, double mass) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:642
Helix::radius
double radius(void) const
returns radious of helix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:224
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:236
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:272
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:218
IMagneticFieldSvc::getReferField
virtual double getReferField()=0