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