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