BOSS 6.6.4.p03
BESIII Offline Software System
Loading...
Searching...
No Matches
Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx
Go to the documentation of this file.
1//
2// $Id: Helix.cxx,v 1.4 2011/05/12 10:25: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(const HepPoint3D & pivot,
48 const HepVector & a,
49 const HepSymMatrix & Ea)
50 : //m_bField(-10.0),
51 //m_alpha(-333.564095),
52 m_pivot(pivot),
53 m_a(a),
54 m_matrixValid(true),
55 m_Ea(Ea) {
56 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
57 if(scmgn!=StatusCode::SUCCESS) {
58 std::cout<< "Unable to open Magnetic field service"<<std::endl;
59 }
60 m_bField = -10000*(m_pmgnIMF->getReferField());
61 m_alpha = 10000. / 2.99792458 / m_bField;
62 // m_alpha = 10000. / 2.99792458 / m_bField;
63 // m_alpha = 333.564095;
64 updateCache();
65}
66
67Helix::Helix(const HepPoint3D & pivot,
68 const HepVector & a)
69 : //m_bField(-10.0),
70 //m_alpha(-333.564095),
71 m_pivot(pivot),
72 m_a(a),
73 m_matrixValid(false),
74 m_Ea(HepSymMatrix(5,0)) {
75 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
76 if(scmgn!=StatusCode::SUCCESS) {
77 // log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
78 std::cout<< "Unable to open Magnetic field service"<<std::endl;
79 }
80 m_bField = -10000*(m_pmgnIMF->getReferField());
81 m_alpha = 10000. / 2.99792458 / m_bField;
82 // m_alpha = 333.564095;
83 //cout<<"MdcFastTrakAlg:: bField,alpha: "<<m_bField<<" , "<<m_alpha<<endl;
84 updateCache();
85}
86
87Helix::Helix(const HepPoint3D & position,
88 const Hep3Vector & momentum,
89 double charge)
90 : //m_bField(-10.0),
91 //m_alpha(-333.564095),
92 m_pivot(position),
93 m_a(HepVector(5,0)),
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());
102 m_alpha = 10000. / 2.99792458 / m_bField;
103
104 m_a[0] = 0.;
105 m_a[1] = fmod(atan2(- momentum.x(), momentum.y())
106 + M_PI4, M_PI2);
107 m_a[3] = 0.;
108 double perp(momentum.perp());
109 if (perp != 0.0) {
110 m_a[2] = charge / perp;
111 m_a[4] = momentum.z() / perp;
112 }
113 else {
114 m_a[2] = charge * (DBL_MAX);
115 if (momentum.z() >= 0) {
116 m_a[4] = (DBL_MAX);
117 } else {
118 m_a[4] = -(DBL_MAX);
119 }
120 }
121 // m_alpha = 333.564095;
122 updateCache();
123}
124
125Helix::~Helix() {
126}
127
128HepPoint3D
129Helix::x(double phi) const {
130 //
131 // Calculate position (x,y,z) along helix.
132 //
133 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
134 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
135 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
136 //
137
138 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
139 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
140 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
141
142 return HepPoint3D(x, y, z);
143}
144
145double *
146Helix::x(double phi, double p[3]) const {
147 //
148 // Calculate position (x,y,z) along helix.
149 //
150 // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
151 // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
152 // z = z0 + dz - (alpha / kappa) * tan(lambda) * phi
153 //
154
155 p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
156 p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
157 p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
158
159 return p;
160}
161
162HepPoint3D
163Helix::x(double phi, HepSymMatrix & Ex) const {
164 double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
165 double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
166 double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;
167
168 //
169 // Calculate position error matrix.
170 // Ex(phi) = (@x/@a)(Ea)(@x/@a)^T, phi is deflection angle to specify the
171 // point to be calcualted.
172 //
173 // HepMatrix dXDA(3, 5, 0);
174 // dXDA = delXDelA(phi);
175 // Ex.assign(dXDA * m_Ea * dXDA.T());
176
177 if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
178 else Ex = m_Ea;
179
180 return HepPoint3D(x, y, z);
181}
182
183Hep3Vector
184Helix::momentum(double phi) const {
185 //
186 // Calculate momentum.
187 //
188 // Pt = | 1/kappa | (GeV/c)
189 //
190 // Px = -Pt * sin(phi0 + phi)
191 // Py = Pt * cos(phi0 + phi)
192 // Pz = Pt * tan(lambda)
193 //
194
195 double pt = fabs(m_pt);
196 double px = - pt * sin(m_ac[1] + phi);
197 double py = pt * cos(m_ac[1] + phi);
198 double pz = pt * m_ac[4];
199
200 return Hep3Vector(px, py, pz);
201}
202
203Hep3Vector
204Helix::momentum(double phi, HepSymMatrix & Em) const {
205 //
206 // Calculate momentum.
207 //
208 // Pt = | 1/kappa | (GeV/c)
209 //
210 // Px = -Pt * sin(phi0 + phi)
211 // Py = Pt * cos(phi0 + phi)
212 // Pz = Pt * tan(lambda)
213 //
214
215 double pt = fabs(m_pt);
216 double px = - pt * sin(m_ac[1] + phi);
217 double py = pt * cos(m_ac[1] + phi);
218 double pz = pt * m_ac[4];
219
220 if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
221 else Em = m_Ea;
222
223 return Hep3Vector(px, py, pz);
224}
225
226HepLorentzVector
227Helix::momentum(double phi, double mass) 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 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
238
239 double pt = fabs(m_pt);
240 double px = - pt * sin(m_ac[1] + phi);
241 double py = pt * cos(m_ac[1] + phi);
242 double pz = pt * m_ac[4];
243 double E = sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);
244
245 return HepLorentzVector(px, py, pz, E);
246}
247
248
249HepLorentzVector
250Helix::momentum(double phi, double mass, HepSymMatrix & Em) 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 if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi,mass));
269 else Em = m_Ea;
270
271 return HepLorentzVector(px, py, pz, E);
272}
273
274HepLorentzVector
275Helix::momentum(double phi,
276 double mass,
277 HepPoint3D & x,
278 HepSymMatrix & Emx) const {
279 //
280 // Calculate momentum.
281 //
282 // Pt = | 1/kappa | (GeV/c)
283 //
284 // Px = -Pt * sin(phi0 + phi)
285 // Py = Pt * cos(phi0 + phi)
286 // Pz = Pt * tan(lambda)
287 //
288 // E = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
289
290 double pt = fabs(m_pt);
291 double px = - pt * sin(m_ac[1] + phi);
292 double py = pt * cos(m_ac[1] + phi);
293 double pz = pt * m_ac[4];
294 double E = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);
295
296 x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
297 x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
298 x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);
299
300 if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi,mass));
301 else Emx = m_Ea;
302
303 return HepLorentzVector(px, py, pz, E);
304}
305
306
307const HepPoint3D &
308Helix::pivot(const HepPoint3D & newPivot) {
309 const double & dr = m_ac[0];
310 const double & phi0 = m_ac[1];
311 const double & kappa = m_ac[2];
312 const double & dz = m_ac[3];
313 const double & tanl = m_ac[4];
314
315 double rdr = dr + m_r;
316 double phi = fmod(phi0 + M_PI4, M_PI2);
317 double csf0 = cos(phi);
318 double snf0 = (1. - csf0) * (1. + csf0);
319 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
320 if(phi > M_PI) snf0 = - snf0;
321
322 double xc = m_pivot.x() + rdr * csf0;
323 double yc = m_pivot.y() + rdr * snf0;
324 double csf, snf;
325 if(m_r != 0.0) {
326 csf = (xc - newPivot.x()) / m_r;
327 snf = (yc - newPivot.y()) / m_r;
328 double anrm = sqrt(csf * csf + snf * snf);
329 if(anrm != 0.0) {
330 csf /= anrm;
331 snf /= anrm;
332 phi = atan2(snf, csf);
333 } else {
334 csf = 1.0;
335 snf = 0.0;
336 phi = 0.0;
337 }
338 } else {
339 csf = 1.0;
340 snf = 0.0;
341 phi = 0.0;
342 }
343 double phid = fmod(phi - phi0 + M_PI8, M_PI2);
344 if(phid > M_PI) phid = phid - M_PI2;
345 double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
346 * csf
347 + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
348 double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();
349
350 HepVector ap(5);
351 ap[0] = drp;
352 ap[1] = fmod(phi + M_PI4, M_PI2);
353 ap[2] = kappa;
354 ap[3] = dzp;
355 ap[4] = tanl;
356
357 // if (m_matrixValid) m_Ea.assign(delApDelA(ap) * m_Ea * delApDelA(ap).T());
358 if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));
359
360 m_a = ap;
361 m_pivot = newPivot;
362
363 //...Are these needed?...iw...
364 updateCache();
365 return m_pivot;
366}
367
368void
369Helix::set(const HepPoint3D & pivot,
370 const HepVector & a,
371 const HepSymMatrix & Ea) {
372 m_pivot = pivot;
373 m_a = a;
374 m_Ea = Ea;
375 m_matrixValid = true;
376 updateCache();
377}
378
379Helix &
380Helix::operator = (const Helix & i) {
381 if (this == & i) return * this;
382
383 m_bField = i.m_bField;
384 m_alpha = i.m_alpha;
385 m_pivot = i.m_pivot;
386 m_a = i.m_a;
387 m_Ea = i.m_Ea;
388 m_matrixValid = i.m_matrixValid;
389
390 m_center = i.m_center;
391 m_cp = i.m_cp;
392 m_sp = i.m_sp;
393 m_pt = i.m_pt;
394 m_r = i.m_r;
395 m_ac[0] = i.m_ac[0];
396 m_ac[1] = i.m_ac[1];
397 m_ac[2] = i.m_ac[2];
398 m_ac[3] = i.m_ac[3];
399 m_ac[4] = i.m_ac[4];
400
401 return * this;
402}
403
404void
405Helix::updateCache(void) {
406 //
407 // Calculate Helix center( xc, yc ).
408 //
409 // xc = x0 + (dr + (alpha / kappa)) * cos(phi0) (cm)
410 // yc = y0 + (dr + (alpha / kappa)) * sin(phi0) (cm)
411 //
412
413 m_ac[0] = m_a[0];
414 m_ac[1] = m_a[1];
415 m_ac[2] = m_a[2];
416 m_ac[3] = m_a[3];
417 m_ac[4] = m_a[4];
418
419 m_cp = cos(m_ac[1]);
420 m_sp = sin(m_ac[1]);
421 if (m_ac[2] != 0.0) {
422 m_pt = 1. / m_ac[2];
423 m_r = m_alpha / m_ac[2];
424 }
425 else {
426 m_pt = (DBL_MAX);
427 m_r = (DBL_MAX);
428 }
429
430 double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
431 double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
432 m_center.setX(x);
433 m_center.setY(y);
434 m_center.setZ(0.);
435}
436
437HepMatrix
438Helix::delApDelA(const HepVector & ap) const {
439 //
440 // Calculate Jacobian (@ap/@a)
441 // Vector ap is new helix parameters and a is old helix parameters.
442 //
443
444 HepMatrix dApDA(5,5,0);
445
446 const double & dr = m_ac[0];
447 const double & phi0 = m_ac[1];
448 const double & cpa = m_ac[2];
449 const double & dz = m_ac[3];
450 const double & tnl = m_ac[4];
451
452 double drp = ap[0];
453 double phi0p = ap[1];
454 double cpap = ap[2];
455 double dzp = ap[3];
456 double tnlp = ap[4];
457
458 double rdr = m_r + dr;
459 double rdrpr;
460 if ((m_r + drp) != 0.0) {
461 rdrpr = 1. / (m_r + drp);
462 } else {
463 rdrpr = (DBL_MAX);
464 }
465 // double csfd = cos(phi0)*cos(phi0p) + sin(phi0)*sin(phi0p);
466 // double snfd = cos(phi0)*sin(phi0p) - sin(phi0)*cos(phi0p);
467 double csfd = cos(phi0p - phi0);
468 double snfd = sin(phi0p - phi0);
469 double phid = fmod(phi0p - phi0 + M_PI8, M_PI2);
470 if (phid > M_PI) phid = phid - M_PI2;
471
472 dApDA[0][0] = csfd;
473 dApDA[0][1] = rdr*snfd;
474 if(cpa!=0.0) {
475 dApDA[0][2] = (m_r/cpa)*( 1.0 - csfd );
476 } else {
477 dApDA[0][2] = (DBL_MAX);
478 }
479
480 dApDA[1][0] = - rdrpr*snfd;
481 dApDA[1][1] = rdr*rdrpr*csfd;
482 if(cpa!=0.0) {
483 dApDA[1][2] = (m_r/cpa)*rdrpr*snfd;
484 } else {
485 dApDA[1][2] = (DBL_MAX);
486 }
487
488 dApDA[2][2] = 1.0;
489
490 dApDA[3][0] = m_r*rdrpr*tnl*snfd;
491 dApDA[3][1] = m_r*tnl*(1.0 - rdr*rdrpr*csfd);
492 if(cpa!=0.0) {
493 dApDA[3][2] = (m_r/cpa)*tnl*(phid - m_r*rdrpr*snfd);
494 } else {
495 dApDA[3][2] = (DBL_MAX);
496 }
497 dApDA[3][3] = 1.0;
498 dApDA[3][4] = - m_r*phid;
499
500 dApDA[4][4] = 1.0;
501
502 return dApDA;
503}
504
505HepMatrix
506Helix::delXDelA(double phi) const {
507 //
508 // Calculate Jacobian (@x/@a)
509 // Vector a is helix parameters and phi is internal parameter
510 // which specifys the point to be calculated for Ex(phi).
511 //
512
513 HepMatrix dXDA(3,5,0);
514
515 const double & dr = m_ac[0];
516 const double & phi0 = m_ac[1];
517 const double & cpa = m_ac[2];
518 const double & dz = m_ac[3];
519 const double & tnl = m_ac[4];
520
521 double cosf0phi = cos(phi0 + phi);
522 double sinf0phi = sin(phi0 + phi);
523
524 dXDA[0][0] = m_cp;
525 dXDA[0][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
526 if(cpa!=0.0) {
527 dXDA[0][2] = - (m_r / cpa) * (m_cp - cosf0phi);
528 } else {
529 dXDA[0][2] = (DBL_MAX);
530 }
531 // dXDA[0][3] = 0.0;
532 // dXDA[0][4] = 0.0;
533
534 dXDA[1][0] = m_sp;
535 dXDA[1][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
536 if(cpa!=0.0) {
537 dXDA[1][2] = - (m_r / cpa) * (m_sp - sinf0phi);
538 } else {
539 dXDA[1][2] = (DBL_MAX);
540 }
541 // dXDA[1][3] = 0.0;
542 // dXDA[1][4] = 0.0;
543
544 // dXDA[2][0] = 0.0;
545 // dXDA[2][1] = 0.0;
546 if(cpa!=0.0) {
547 dXDA[2][2] = (m_r / cpa) * tnl * phi;
548 } else {
549 dXDA[2][2] = (DBL_MAX);
550 }
551 dXDA[2][3] = 1.0;
552 dXDA[2][4] = - m_r * phi;
553
554 return dXDA;
555}
556
557
558
559HepMatrix
560Helix::delMDelA(double phi) const {
561 //
562 // Calculate Jacobian (@m/@a)
563 // Vector a is helix parameters and phi is internal parameter.
564 // Vector m is momentum.
565 //
566
567 HepMatrix dMDA(3,5,0);
568
569 const double & phi0 = m_ac[1];
570 const double & cpa = m_ac[2];
571 const double & tnl = m_ac[4];
572
573 double cosf0phi = cos(phi0+phi);
574 double sinf0phi = sin(phi0+phi);
575
576 double rho;
577 if(cpa != 0.)rho = 1./cpa;
578 else rho = (DBL_MAX);
579
580 double charge = 1.;
581 if(cpa < 0.)charge = -1.;
582
583 dMDA[0][1] = -fabs(rho)*cosf0phi;
584 dMDA[0][2] = charge*rho*rho*sinf0phi;
585
586 dMDA[1][1] = -fabs(rho)*sinf0phi;
587 dMDA[1][2] = -charge*rho*rho*cosf0phi;
588
589 dMDA[2][2] = -charge*rho*rho*tnl;
590 dMDA[2][4] = fabs(rho);
591
592 return dMDA;
593}
594
595
596HepMatrix
597Helix::del4MDelA(double phi, double mass) const {
598 //
599 // Calculate Jacobian (@4m/@a)
600 // Vector a is helix parameters and phi is internal parameter.
601 // Vector 4m is 4 momentum.
602 //
603
604 HepMatrix d4MDA(4,5,0);
605
606 double phi0 = m_ac[1];
607 double cpa = m_ac[2];
608 double tnl = m_ac[4];
609
610 double cosf0phi = cos(phi0+phi);
611 double sinf0phi = sin(phi0+phi);
612
613 double rho;
614 if(cpa != 0.)rho = 1./cpa;
615 else rho = (DBL_MAX);
616
617 double charge = 1.;
618 if(cpa < 0.)charge = -1.;
619
620 double E = sqrt(rho*rho*(1.+tnl*tnl)+mass*mass);
621
622 d4MDA[0][1] = -fabs(rho)*cosf0phi;
623 d4MDA[0][2] = charge*rho*rho*sinf0phi;
624
625 d4MDA[1][1] = -fabs(rho)*sinf0phi;
626 d4MDA[1][2] = -charge*rho*rho*cosf0phi;
627
628 d4MDA[2][2] = -charge*rho*rho*tnl;
629 d4MDA[2][4] = fabs(rho);
630
631 if (cpa != 0.0 && E != 0.0) {
632 d4MDA[3][2] = (-1.-tnl*tnl)/(cpa*cpa*cpa*E);
633 d4MDA[3][4] = tnl/(cpa*cpa*E);
634 } else {
635 d4MDA[3][2] = (DBL_MAX);
636 d4MDA[3][4] = (DBL_MAX);
637 }
638 return d4MDA;
639}
640
641
642HepMatrix
643Helix::del4MXDelA(double phi, double mass) const {
644 //
645 // Calculate Jacobian (@4mx/@a)
646 // Vector a is helix parameters and phi is internal parameter.
647 // Vector 4xm is 4 momentum and position.
648 //
649
650 HepMatrix d4MXDA(7,5,0);
651
652 const double & dr = m_ac[0];
653 const double & phi0 = m_ac[1];
654 const double & cpa = m_ac[2];
655 const double & dz = m_ac[3];
656 const double & tnl = m_ac[4];
657
658 double cosf0phi = cos(phi0+phi);
659 double sinf0phi = sin(phi0+phi);
660
661 double rho;
662 if(cpa != 0.)rho = 1./cpa;
663 else rho = (DBL_MAX);
664
665 double charge = 1.;
666 if(cpa < 0.)charge = -1.;
667
668 double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);
669
670 d4MXDA[0][1] = - fabs(rho) * cosf0phi;
671 d4MXDA[0][2] = charge * rho * rho * sinf0phi;
672
673 d4MXDA[1][1] = - fabs(rho) * sinf0phi;
674 d4MXDA[1][2] = - charge * rho * rho * cosf0phi;
675
676 d4MXDA[2][2] = - charge * rho * rho * tnl;
677 d4MXDA[2][4] = fabs(rho);
678
679 if (cpa != 0.0 && E != 0.0) {
680 d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
681 d4MXDA[3][4] = tnl / (cpa * cpa * E);
682 } else {
683 d4MXDA[3][2] = (DBL_MAX);
684 d4MXDA[3][4] = (DBL_MAX);
685 }
686
687 d4MXDA[4][0] = m_cp;
688 d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
689 if (cpa != 0.0) {
690 d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
691 } else {
692 d4MXDA[4][2] = (DBL_MAX);
693 }
694
695 d4MXDA[5][0] = m_sp;
696 d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
697 if (cpa != 0.0) {
698 d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);
699
700 d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
701 } else {
702 d4MXDA[5][2] = (DBL_MAX);
703
704 d4MXDA[6][2] = (DBL_MAX);
705 }
706
707 d4MXDA[6][3] = 1.;
708 d4MXDA[6][4] = - m_r * phi;
709
710 return d4MXDA;
711}
712
713void
714Helix::ignoreErrorMatrix(void) {
715 m_matrixValid = false;
716 m_Ea *= 0.;
717}
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-01/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-08/src/Helix.cxx:438
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-08/src/Helix.cxx:369
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-08/src/Helix.cxx:714
Helix::operator=
Helix & operator=(const Helix &)
Copy operator.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:380
Helix::del4MXDelA
HepMatrix del4MXDelA(double phi, double mass) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:643
Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:184
Helix::ConstantAlpha
static const double ConstantAlpha
Constant alpha for uniform field.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:171
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:129
Helix::del4MDelA
HepMatrix del4MDelA(double phi, double mass) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:597
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:227
Helix::Helix
Helix(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
Constructor with pivot, helix parameter a, and its error matrix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:47
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:209
IMagneticFieldSvc::getReferField
virtual double getReferField()=0