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