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