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