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