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