BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TrkMomCalculator.cxx
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// File and Version Information:
3// $Id: TrkMomCalculator.cxx,v 1.2 2007/12/07 07:04:14 zhangy Exp $
4// Description: See the .h file
5//
6//
7// Environment:
8// Software developed for the BaBar Detector at the SLAC B-Factory.
9//
10// Author(s): Justin Albert, Steve Schaffner
11// Add functions needed for Vertexing: Mario Bondioli, Riccardo Faccini, Eugenio Paoloni.
12// Add fast Helix->PX methods: Aaron Roodman
13//
14//------------------------------------------------------------------------
15
16// #include "BaBar/BaBar.h"
17#include "TrkBase/TrkMomCalculator.h"
18#include "CLHEP/Geometry/Point3D.h"
19#include "CLHEP/Matrix/Vector.h"
20#include "CLHEP/Vector/ThreeVector.h"
21#include "MdcRecoUtil/BesVectorErr.h"
22#include "MdcRecoUtil/DifIndepPar.h"
23#include "MdcRecoUtil/DifNumber.h"
24#include "MdcRecoUtil/DifPoint.h"
25#include "MdcRecoUtil/DifVector.h"
26#include "PatBField/BField.h"
27#include "TrkBase/HelixTraj.h"
28#include "TrkBase/NeutParams.h"
29#include "TrkBase/TrkExchangePar.h"
30#include "TrkBase/TrkMomVisitor.h"
31#include "TrkBase/TrkSimpTraj.h"
32// #include "ErrLogger/ErrLog.h"
33using std::endl;
34
35using MdcPatRec::BField;
36
37//------------------------------------------------------------------------
38TrkMomCalculator::~TrkMomCalculator() {
39 //------------------------------------------------------------------------
40}
41
42//------------------------------------------------------------------------
43TrkMomCalculator::TrkMomCalculator() {
44 //------------------------------------------------------------------------
45}
46
47//------------------------------------------------------------------------
48double TrkMomCalculator::ptMom( const TrkSimpTraj& theTraj, const BField& theField,
49 double fltlen ) {
50 //------------------------------------------------------------------------
51
52 TrkMomVisitor theVisitor( theTraj );
53
54 if ( theVisitor.helix() != 0 || theVisitor.circle() != 0 )
55 {
56
57 // treat as curve, calculate Pt accordingly (see calcCurvMom function
58 // below)
59 return calcCurvPtMom( theTraj.direction( fltlen ), theTraj.curvature( fltlen ), theField );
60 }
61 else if ( theVisitor.neut() != 0 )
62 {
63
64 // treat as neutral particle (with the ptot as a parameter)
65 double sindip = theTraj.direction( fltlen ).z();
66 double arg = 1.0 - sindip * sindip;
67 if ( arg < 0.0 ) arg = 0.0;
68 double cosdip = sqrt( arg );
69 double ptot = theTraj.parameters()->parameter()[NeutParams::_p];
70
71 return cosdip * ptot;
72 }
73 else
74 {
75
76 // particle must be a plain line--no way to calculate momentum
77 return 999999.99;
78 }
79}
80
81//------------------------------------------------------------------------
82Hep3Vector TrkMomCalculator::vecMom( const TrkSimpTraj& theTraj, const BField& theField,
83 double fltlen ) {
84 //------------------------------------------------------------------------
85 TrkMomVisitor theVisitor( theTraj );
86
87 if ( theVisitor.helix() != 0 || theVisitor.circle() != 0 )
88 {
89
90 // treat as curve, calculate VecMom accordingly (see calcCurvMom function
91 // below)
92
93 return calcCurvVecMom( theTraj.direction( fltlen ), theTraj.curvature( fltlen ),
94 theField );
95 }
96 else if ( theVisitor.neut() != 0 )
97 {
98
99 // treat as neutral particle (with the pt as a parameter)
100 Hep3Vector theMom = theTraj.direction( fltlen );
101 theMom.setMag( theTraj.parameters()->parameter()[NeutParams::_p] );
102 return theMom;
103 }
104 else
105 {
106
107 // particle must be a plain line--no way to calculate momentum
108 return Hep3Vector( 999, 999, 999 );
109 }
110}
111
112//------------------------------------------------------------------------
113BesVectorErr TrkMomCalculator::errMom( const TrkSimpTraj& theTraj, const BField& theField,
114 double fltlen ) {
115 //------------------------------------------------------------------------
116
117 TrkMomVisitor theVisitor( theTraj );
118
119 if ( theVisitor.helix() != 0 || theVisitor.circle() != 0 )
120 {
121
122 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function
123 // below)
124 return calcCurvErrMom( theTraj, theField, fltlen );
125 }
126 else if ( theVisitor.neut() != 0 )
127 {
128
129 // treat as neutral particle, same as curve in this case
130 return calcNeutErrMom( theTraj, theField, fltlen );
131 }
132 else
133 {
134
135 // particle must be a plain line--no way to calculate momentum or err
136 // The matrix is initialized to zero (see BesError constructor)
137 BesError theErr( 3 );
138 return BesVectorErr( Hep3Vector( 999, 999, 999 ), theErr );
139 }
140}
141
142//------------------------------------------------------------------------
143int TrkMomCalculator::charge( const TrkSimpTraj& theTraj, const BField& theField,
144 double fltlen ) {
145 //------------------------------------------------------------------------
146
147 TrkMomVisitor theVisitor( theTraj );
148
149 if ( theVisitor.helix() != 0 || theVisitor.circle() != 0 )
150 {
151
152 // treat as curve, calculate Pt accordingly (see calcCurvMom function
153 // below)
154
155 bool plus = false;
156 HelixTraj::ParIndex parInd = HelixTraj::omegaIndex;
157 int index = parInd;
158 plus =
159 ( theField.bFieldNominal() < 0.0 && theTraj.parameters()->parameter()[index] > 0.0 ) ||
160 ( theField.bFieldNominal() > 0.0 && theTraj.parameters()->parameter()[index] < 0.0 );
161 return ( plus ? 1 : -1 );
162
163 // return calcCurvCharge(
164 // theTraj.direction(fltlen),
165 // theTraj.curvature(fltlen), theField);
166 }
167 else if ( theVisitor.neut() != 0 )
168 {
169
170 // treat as neutral particle, so charge is zero
171 return 0;
172 }
173 else
174 {
175
176 // particle must be a plain line--take charge as zero
177 return 0;
178 }
179}
180
181//------------------------------------------------------------------------
182HepMatrix TrkMomCalculator::posmomCov( const TrkSimpTraj& theTraj, const BField& theField,
183 double fltlen ) {
184 //------------------------------------------------------------------------
185
186 TrkMomVisitor theVisitor( theTraj );
187
188 if ( theVisitor.helix() != 0 || theVisitor.circle() != 0 )
189 {
190
191 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function
192 // below)
193 return calcCurvPosmomCov( theTraj, theField, fltlen );
194 }
195 else if ( theVisitor.neut() != 0 )
196 {
197
198 // treat as neutral particle, same as curve in this case
199 return calcNeutPosmomCov( theTraj, theField, fltlen );
200 }
201 else
202 {
203
204 // particle must be a plain line--no way to calculate momentum or err
205 // The matrix is initialized to zero (see BesError constructor)
206 return HepMatrix( 3, 3, 0 );
207 }
208}
209
210//------------------------------------------------------------------------
211void TrkMomCalculator::getAllCovs( const TrkSimpTraj& theTraj, const BField& theField,
212 double fltlen, HepSymMatrix& xxCov, HepSymMatrix& ppCov,
213 HepMatrix& xpCov ) {
214 //------------------------------------------------------------------------
215
216 TrkMomVisitor theVisitor( theTraj );
217
218 if ( theVisitor.helix() != 0 )
219 {
220
221 // fast inline calculation...
222 calcCurvAllCovs( theTraj, theField, fltlen, xxCov, ppCov, xpCov );
223 }
224 else if ( theVisitor.circle() != 0 )
225 {
226
227 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function
228 // below)
229 calcCurvAllCovsOLD( theTraj, theField, fltlen, xxCov, ppCov, xpCov );
230 }
231 else if ( theVisitor.neut() != 0 )
232 {
233
234 // treat as neutral particle, same as curve in this case
235 calcNeutAllCovs( theTraj, theField, fltlen, xxCov, ppCov, xpCov );
236 }
237 else
238 {
239
240 // particle must be a plain line--no way to calculate momentum or err
241 // The matrix is initialized to zero (see BesError constructor)
242 int i, j;
243 for ( i = 0; i < 3; i++ )
244 for ( j = i; j < 3; j++ )
245 {
246 xxCov[i][j] = 0;
247 ppCov[i][j] = 0;
248 xpCov[i][j] = 0;
249 xpCov[j][i] = 0;
250 }
251 }
252}
253
254//------------------------------------------------------------------------
255void TrkMomCalculator::getAllWeights( const TrkSimpTraj& theTraj, const BField& theField,
256 double fltlen, HepVector& pos, HepVector& mom,
257 HepSymMatrix& xxWeight, HepSymMatrix& ppWeight,
258 HepMatrix& xpWeight ) {
259 //------------------------------------------------------------------------
260
261 TrkMomVisitor theVisitor( theTraj );
262
263 if ( theVisitor.helix() != 0 )
264 {
265
266 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function
267 // below)
268 calcCurvAllWeights( theTraj, theField, fltlen, pos, mom, xxWeight, ppWeight, xpWeight );
269 }
270 else if ( theVisitor.circle() != 0 )
271 {
272
273 calcCurvAllWeightsOLD( theTraj, theField, fltlen, pos, mom, xxWeight, ppWeight, xpWeight );
274 }
275 else if ( theVisitor.neut() != 0 )
276 {
277
278 // treat as neutral particle, same as curve in this case
279 calcNeutAllWeights( theTraj, theField, fltlen, pos, mom, xxWeight, ppWeight, xpWeight );
280 }
281 else
282 {
283
284 // particle must be a plain line--no way to calculate momentum or err
285 // temporary: initialize everything to 0
286 int i, j;
287 for ( i = 0; i < 3; i++ )
288 for ( j = i; j < 3; j++ )
289 {
290 xxWeight[i][j] = 0;
291 ppWeight[i][j] = 0;
292 xpWeight[i][j] = 0;
293 xpWeight[j][i] = 0;
294 }
295
296 for ( i = 0; i < 3; i++ )
297 {
298 pos[i] = 0;
299 mom[i] = 0;
300 }
301 }
302}
303
304//------------------------------------------------------------------------
305double TrkMomCalculator::calcCurvPtMom( const Hep3Vector& direction, double curvature,
306 const BField& theField ) {
307 //------------------------------------------------------------------------
308 Hep3Vector theMomVec = calcCurvVecMom( direction, curvature, theField );
309 return theMomVec.perp();
310}
311
312//------------------------------------------------------------------------
313Hep3Vector TrkMomCalculator::calcCurvVecMom( const Hep3Vector& direction, double curvature,
314 const BField& theField ) {
315 //------------------------------------------------------------------------
316 double sindip = direction.z();
317 double arg = 1.0 - sindip * sindip;
318 if ( arg < 0.0 ) arg = 0.0;
319 double cosdip = sqrt( arg );
320 double momMag =
321 fabs( BField::cmTeslaToGeVc * theField.bFieldNominal() * cosdip / curvature );
322 Hep3Vector momVec = direction;
323 momVec.setMag( momMag );
324
325 return momVec;
326}
327
328//------------------------------------------------------------------------
329BesVectorErr TrkMomCalculator::calcCurvErrMom( const TrkSimpTraj& theTraj,
330 const BField& theField, double fltlen ) {
331 //------------------------------------------------------------------------
332
333 DifPoint PosDif;
334 DifVector DirDif;
335 DifVector delDirDif;
336 DifVector MomDif( 0., 0., 0. );
337
338 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
339 if ( delDirDif.length() == 0. ) {}
340 else
341 {
342 DifNumber sindip = DirDif.z;
343 DifNumber arg = 1.0 - sindip * sindip;
344 if ( arg.number() < 0.0 ) { arg.setNumber( 0.0 ); }
345 DifNumber cosdip = sqrt( arg );
346
347 DifNumber momMag =
348 BField::cmTeslaToGeVc * theField.bFieldNominal() * cosdip / delDirDif.length();
349 momMag.absolute();
350
351 MomDif = DirDif * momMag;
352 /*
353 if (ErrLogging(debugging) && 0) {
354 HepMatrix e1 = DirDif.errorMatrix(MomDif.x.indepPar()->covariance());
355 double e2 = momMag.error(MomDif.x.indepPar()->covariance());
356 HepMatrix e3 = MomDif.errorMatrix(MomDif.x.indepPar()->covariance());
357 const HepMatrix& c1 = MomDif.x.indepPar()->covariance();
358
359 std::cout<<"ErrMsg(debugging) "<< endl << endl << "Start" << endl
360 << "param cov: " << endl
361 << c1(1,1) << endl
362 << c1(2,1) << " " << c1(2,2) << endl
363 << c1(3,1) << " " << c1(3,2) << " " << c1(3,3) << endl
364 << c1(4,1) << " " << c1(4,2) << " " << c1(4,3)
365 << " " << c1(4,4) << endl
366 << c1(5,1) << " " << c1(5,2) << " " << c1(5,3)
367 << " " << c1(5,4) << " " << c1(5,5) << endl
368
369 << "Dir " << e1.num_row() << " " << e1.num_col() << endl
370 << "output:" << endl
371 << e1(1,1) << endl
372 << e1(2,1) << " " << e1(2,2) << endl
373 << e1(3,1) << " " << e1(3,2) << " " << e1(3,3) << endl
374 << "deriv: " << endl
375 << DirDif.x.derivatives() << endl
376 << endl
377
378 << "Mag " << endl
379 << "output:" << endl
380 << e2 << endl
381 << "deriv: " << endl
382 << momMag.derivatives() << endl
383 << endl
384
385 << "Momdif " << e3.num_row() << " " << e3.num_col() << endl
386 << "output:" << endl
387 << e3(1,1) << endl
388 << e3(2,1) << " " << e3(2,2) << endl
389 << e3(3,1) << " " << e3(3,2) << " " << e3(3,3) << endl
390 << "deriv: " << endl
391 << MomDif.x.derivatives() << endl
392 << endl
393
394 << "End" << endl << std::endl;
395 } */
396 }
397 BesError symErr( MomDif.errorMatrix( MomDif.x.indepPar()->covariance() ) );
398 return BesVectorErr( Hep3Vector( MomDif.x.number(), MomDif.y.number(), MomDif.z.number() ),
399 symErr );
400}
401
402//------------------------------------------------------------------------
403BesVectorErr TrkMomCalculator::calcNeutErrMom( const TrkSimpTraj& theTraj,
404 const BField& theField, double fltlen ) {
405 //------------------------------------------------------------------------
406
407 DifPoint posDif;
408 DifVector dirDif;
409 DifVector delDirDif;
410
411 theTraj.getDFInfo( fltlen, posDif, dirDif, delDirDif );
412
413 // set the momentum's direction, and then its magnitude
414 DifVector momDif = dirDif;
415
416 // The + 1 is necessary because DifIndepPar starts counting at 1, whereas
417 // the enum for NeutParams starts at 0. Stupid, ain't it?
418 momDif *= theTraj.parameters()->difPar( NeutParams::_p + 1 );
419
420 // now create an error matrix
421 BesError symErr( momDif.errorMatrix( momDif.x.indepPar()->covariance() ) );
422
423 // get the result in the correct form
424 return BesVectorErr( Hep3Vector( momDif.x.number(), momDif.y.number(), momDif.z.number() ),
425 symErr );
426}
427
428// The following functions may be used at a later date (if and when we
429// ever want to drop the assumption that all recon'ed charges are 0, +1,
430// or -1):
431
432//------------------------------------------------------------------------
433int TrkMomCalculator::calcCurvCharge( const Hep3Vector& direction, double curvature,
434 const BField& theField ) {
435 //------------------------------------------------------------------------
436
437 Hep3Vector momVec = calcCurvVecMom( direction, curvature, theField );
438
439 if ( theField.bFieldNominal() > 0. )
440 { return -nearestInt( momVec.mag() * curvature / theField.bFieldNominal() ); }
441 else { return nearestInt( momVec.mag() * curvature / theField.bFieldNominal() ); }
442}
443
444//------------------------------------------------------------------------
445int TrkMomCalculator::nearestInt( double floatPt ) {
446 //------------------------------------------------------------------------
447
448 if ( floatPt > 0. ) { return (int)( floatPt + 0.5 ); }
449 else { return (int)( floatPt - 0.5 ); }
450}
451
452/*
453 * These function were added on 7/18/98 to accomplish
454 * vertexing interface (M.Bondioli)
455 */
456
457//------------------------------------------------------------------------
458HepMatrix TrkMomCalculator::calcCurvPosmomCov( const TrkSimpTraj& theTraj,
459 const BField& theField, double fltlen ) {
460 //------------------------------------------------------------------------
461
462 DifPoint PosDif;
463 DifVector DirDif;
464 DifVector delDirDif;
465 DifVector MomDif( 0., 0., 0. );
466
467 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
468 if ( delDirDif.length() == 0. ) {}
469 else
470 {
471 DifNumber sindip = DirDif.z;
472 DifNumber arg = 1.0 - sindip * sindip;
473 if ( arg.number() < 0.0 ) { arg.setNumber( 0.0 ); }
474 DifNumber cosdip = sqrt( arg );
475
476 DifNumber momMag =
477 BField::cmTeslaToGeVc * theField.bFieldNominal() * cosdip / delDirDif.length();
478 momMag.absolute();
479
480 MomDif = DirDif * momMag;
481 }
482
483 // computes the correlation among position and momentum
484 HepMatrix xpCov( 3, 3 );
485 const HepSymMatrix& nnCov = MomDif.x.indepPar()->covariance();
486 xpCov( 1, 1 ) = correlation( PosDif.x, MomDif.x, nnCov );
487 xpCov( 1, 2 ) = correlation( PosDif.x, MomDif.y, nnCov );
488 xpCov( 1, 3 ) = correlation( PosDif.x, MomDif.z, nnCov );
489 xpCov( 2, 1 ) = correlation( PosDif.y, MomDif.x, nnCov );
490 xpCov( 2, 2 ) = correlation( PosDif.y, MomDif.y, nnCov );
491 xpCov( 2, 3 ) = correlation( PosDif.y, MomDif.z, nnCov );
492 xpCov( 3, 1 ) = correlation( PosDif.z, MomDif.x, nnCov );
493 xpCov( 3, 2 ) = correlation( PosDif.z, MomDif.y, nnCov );
494 xpCov( 3, 3 ) = correlation( PosDif.z, MomDif.z, nnCov );
495
496 return xpCov;
497}
498
499//------------------------------------------------------------------------
500HepMatrix TrkMomCalculator::calcNeutPosmomCov( const TrkSimpTraj& theTraj,
501 const BField& theField, double fltlen ) {
502 //------------------------------------------------------------------------
503
504 DifPoint PosDif;
505 DifVector DirDif;
506 DifVector delDirDif;
507
508 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
509
510 // set the momentum's direction, and then its magnitude
511 DifVector MomDif = DirDif;
512
513 // The + 1 is necessary because DifIndepPar starts counting at 1, whereas
514 // the enum for NeutParams starts at 0. Stupid, ain't it?
515 MomDif *= theTraj.parameters()->difPar( NeutParams::_p + 1 );
516
517 // computes the correlation among position and momentum
518 HepMatrix xpCov( 3, 3 );
519 const HepSymMatrix& nnCov = MomDif.x.indepPar()->covariance();
520 xpCov( 1, 1 ) = correlation( PosDif.x, MomDif.x, nnCov );
521 xpCov( 1, 2 ) = correlation( PosDif.x, MomDif.y, nnCov );
522 xpCov( 1, 3 ) = correlation( PosDif.x, MomDif.z, nnCov );
523 xpCov( 2, 1 ) = correlation( PosDif.y, MomDif.x, nnCov );
524 xpCov( 2, 2 ) = correlation( PosDif.y, MomDif.y, nnCov );
525 xpCov( 2, 3 ) = correlation( PosDif.y, MomDif.z, nnCov );
526 xpCov( 3, 1 ) = correlation( PosDif.z, MomDif.x, nnCov );
527 xpCov( 3, 2 ) = correlation( PosDif.z, MomDif.y, nnCov );
528 xpCov( 3, 3 ) = correlation( PosDif.z, MomDif.z, nnCov );
529
530 return HepMatrix( 3, 3, 0 );
531}
532
533bool TrkMomCalculator::weightToCov( const HepSymMatrix& inXX, const HepSymMatrix& inPP,
534 const HepMatrix& inXP, HepSymMatrix& outXX,
535 HepSymMatrix& outPP, HepMatrix& outXP ) {
536 assert( inXX.num_row() == outXX.num_row() );
537 assert( inPP.num_row() == outPP.num_row() );
538 assert( inXP.num_row() == outXP.num_row() );
539 assert( inXP.num_col() == outXP.num_col() );
540 assert( inXX.num_row() == inXP.num_row() );
541 assert( inPP.num_row() == inXP.num_col() );
542 int status;
543 HepSymMatrix aInv = inXX.inverse( status );
544 if ( status ) return false;
545 HepSymMatrix beta = inPP - aInv.similarityT( inXP );
546 outPP = beta.inverse( status );
547 if ( status ) return false;
548 outXP = -aInv * inXP * outPP;
549 HepMatrix alpha( aInv - aInv * inXP * outXP.T() );
550 outXX.assign( alpha );
551 return true;
552}
553
554//------------------------------------------------------------------------
555void TrkMomCalculator::calcCurvAllCovs( const TrkSimpTraj& theTraj, const BField& theField,
556 double fltlen, HepSymMatrix& xxCov,
557 HepSymMatrix& ppCov, HepMatrix& xpCov ) {
558 //------------------------------------------------------------------------
559
560 const HepVector& v = theTraj.parameters()->parameter();
561 const HepSymMatrix& m = theTraj.parameters()->covariance();
562
563 // fast inlined conversion from Helix to PX representation
564
565 // track parameters
566 double s = v[TrkExchangePar::ex_tanDip];
567 double omega = v[TrkExchangePar::ex_omega];
568 double d0 = v[TrkExchangePar::ex_d0];
569 // double z0 = v[TrkExchangePar::ex_z0];
570 double phi0 = v[TrkExchangePar::ex_phi0];
571 double cosDip = 1 / sqrt( 1.0 + s * s );
572 double l = fltlen * cosDip;
573
574 // calculate some sin,cos etc..
575 double dlds = -fltlen * s * ( cosDip * cosDip * cosDip );
576 double phi = phi0 + omega * l;
577 double sinphi0 = sin( phi0 );
578 double cosphi0 = cos( phi0 );
579 double sinphi = sin( phi );
580 double cosphi = cos( phi );
581 double pt = fabs( BField::cmTeslaToGeVc * theField.bFieldNominal() / omega );
582 double r = 1.0 / omega;
583
584 // Calculate derivatives for Jacobian matrix
585 double d_x_d0 = -sinphi0;
586 double d_x_phi0 = r * cosphi - ( r + d0 ) * cosphi0;
587 double d_x_omega = -r * r * sinphi + r * r * sinphi0 + l * r * cosphi;
588 double d_x_tanDip = cosphi * dlds;
589 double d_y_d0 = cosphi0;
590 double d_y_phi0 = r * sinphi - ( r + d0 ) * sinphi0;
591 double d_y_omega = r * r * cosphi - r * r * cosphi0 + l * r * sinphi;
592 double d_y_tanDip = sinphi * dlds;
593 double d_z_z0 = 1.0;
594 double d_z_tanDip = l + dlds * s;
595 double d_px_phi0 = -pt * sinphi;
596 double d_px_omega = -pt * cosphi / omega - pt * l * sinphi;
597 double d_px_tanDip = -pt * omega * sinphi * dlds;
598 double d_py_phi0 = pt * cosphi;
599 double d_py_omega = -pt * sinphi / omega + pt * l * cosphi;
600 double d_py_tanDip = pt * omega * cosphi * dlds;
601 double d_pz_omega = -pt * s / omega;
602 double d_pz_tanDip = pt;
603
604 // Fill temporary variables for m
605 double m_d0_d0 = m[TrkExchangePar::ex_d0][TrkExchangePar::ex_d0];
606 double m_phi0_d0 = m[TrkExchangePar::ex_phi0][TrkExchangePar::ex_d0];
607 double m_phi0_phi0 = m[TrkExchangePar::ex_phi0][TrkExchangePar::ex_phi0];
608 double m_omega_d0 = m[TrkExchangePar::ex_omega][TrkExchangePar::ex_d0];
609 double m_omega_phi0 = m[TrkExchangePar::ex_omega][TrkExchangePar::ex_phi0];
610 double m_omega_omega = m[TrkExchangePar::ex_omega][TrkExchangePar::ex_omega];
611 double m_z0_d0 = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_d0];
612 double m_z0_phi0 = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_phi0];
613 double m_z0_omega = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_omega];
614 double m_z0_z0 = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_z0];
615 double m_tanDip_d0 = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_d0];
616 double m_tanDip_phi0 = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_phi0];
617 double m_tanDip_omega = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_omega];
618 double m_tanDip_z0 = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_z0];
619 double m_tanDip_tanDip = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_tanDip];
620
621 // Fill xxCov,xpCov,ppCov - nb. this code generated by script writecov.py
622 xxCov( 1, 1 ) = d_x_d0 * ( d_x_d0 * m_d0_d0 + d_x_phi0 * m_phi0_d0 + d_x_omega * m_omega_d0 +
623 d_x_tanDip * m_tanDip_d0 ) +
624 d_x_phi0 * ( d_x_d0 * m_phi0_d0 + d_x_phi0 * m_phi0_phi0 +
625 d_x_omega * m_omega_phi0 + d_x_tanDip * m_tanDip_phi0 ) +
626 d_x_omega * ( d_x_d0 * m_omega_d0 + d_x_phi0 * m_omega_phi0 +
627 d_x_omega * m_omega_omega + d_x_tanDip * m_tanDip_omega ) +
628 d_x_tanDip * ( d_x_d0 * m_tanDip_d0 + d_x_phi0 * m_tanDip_phi0 +
629 d_x_omega * m_tanDip_omega + d_x_tanDip * m_tanDip_tanDip );
630 xxCov( 2, 1 ) = d_y_d0 * ( d_x_d0 * m_d0_d0 + d_x_phi0 * m_phi0_d0 + d_x_omega * m_omega_d0 +
631 d_x_tanDip * m_tanDip_d0 ) +
632 d_y_phi0 * ( d_x_d0 * m_phi0_d0 + d_x_phi0 * m_phi0_phi0 +
633 d_x_omega * m_omega_phi0 + d_x_tanDip * m_tanDip_phi0 ) +
634 d_y_omega * ( d_x_d0 * m_omega_d0 + d_x_phi0 * m_omega_phi0 +
635 d_x_omega * m_omega_omega + d_x_tanDip * m_tanDip_omega ) +
636 d_y_tanDip * ( d_x_d0 * m_tanDip_d0 + d_x_phi0 * m_tanDip_phi0 +
637 d_x_omega * m_tanDip_omega + d_x_tanDip * m_tanDip_tanDip );
638 xxCov( 2, 2 ) = d_y_d0 * ( d_y_d0 * m_d0_d0 + d_y_phi0 * m_phi0_d0 + d_y_omega * m_omega_d0 +
639 d_y_tanDip * m_tanDip_d0 ) +
640 d_y_phi0 * ( d_y_d0 * m_phi0_d0 + d_y_phi0 * m_phi0_phi0 +
641 d_y_omega * m_omega_phi0 + d_y_tanDip * m_tanDip_phi0 ) +
642 d_y_omega * ( d_y_d0 * m_omega_d0 + d_y_phi0 * m_omega_phi0 +
643 d_y_omega * m_omega_omega + d_y_tanDip * m_tanDip_omega ) +
644 d_y_tanDip * ( d_y_d0 * m_tanDip_d0 + d_y_phi0 * m_tanDip_phi0 +
645 d_y_omega * m_tanDip_omega + d_y_tanDip * m_tanDip_tanDip );
646 xxCov( 3, 1 ) = d_z_z0 * ( d_x_d0 * m_z0_d0 + d_x_phi0 * m_z0_phi0 + d_x_omega * m_z0_omega +
647 d_x_tanDip * m_tanDip_z0 ) +
648 d_z_tanDip * ( d_x_d0 * m_tanDip_d0 + d_x_phi0 * m_tanDip_phi0 +
649 d_x_omega * m_tanDip_omega + d_x_tanDip * m_tanDip_tanDip );
650 xxCov( 3, 2 ) = d_z_z0 * ( d_y_d0 * m_z0_d0 + d_y_phi0 * m_z0_phi0 + d_y_omega * m_z0_omega +
651 d_y_tanDip * m_tanDip_z0 ) +
652 d_z_tanDip * ( d_y_d0 * m_tanDip_d0 + d_y_phi0 * m_tanDip_phi0 +
653 d_y_omega * m_tanDip_omega + d_y_tanDip * m_tanDip_tanDip );
654 xxCov( 3, 3 ) = d_z_z0 * ( d_z_z0 * m_z0_z0 + d_z_tanDip * m_tanDip_z0 ) +
655 d_z_tanDip * ( d_z_z0 * m_tanDip_z0 + d_z_tanDip * m_tanDip_tanDip );
656
657 xpCov( 1, 1 ) = d_px_phi0 * ( d_x_d0 * m_phi0_d0 + d_x_phi0 * m_phi0_phi0 +
658 d_x_omega * m_omega_phi0 + d_x_tanDip * m_tanDip_phi0 ) +
659 d_px_omega * ( d_x_d0 * m_omega_d0 + d_x_phi0 * m_omega_phi0 +
660 d_x_omega * m_omega_omega + d_x_tanDip * m_tanDip_omega ) +
661 d_px_tanDip * ( d_x_d0 * m_tanDip_d0 + d_x_phi0 * m_tanDip_phi0 +
662 d_x_omega * m_tanDip_omega + d_x_tanDip * m_tanDip_tanDip );
663 xpCov( 2, 1 ) = d_px_phi0 * ( d_y_d0 * m_phi0_d0 + d_y_phi0 * m_phi0_phi0 +
664 d_y_omega * m_omega_phi0 + d_y_tanDip * m_tanDip_phi0 ) +
665 d_px_omega * ( d_y_d0 * m_omega_d0 + d_y_phi0 * m_omega_phi0 +
666 d_y_omega * m_omega_omega + d_y_tanDip * m_tanDip_omega ) +
667 d_px_tanDip * ( d_y_d0 * m_tanDip_d0 + d_y_phi0 * m_tanDip_phi0 +
668 d_y_omega * m_tanDip_omega + d_y_tanDip * m_tanDip_tanDip );
669 xpCov( 3, 1 ) = d_px_phi0 * ( d_z_z0 * m_z0_phi0 + d_z_tanDip * m_tanDip_phi0 ) +
670 d_px_omega * ( d_z_z0 * m_z0_omega + d_z_tanDip * m_tanDip_omega ) +
671 d_px_tanDip * ( d_z_z0 * m_tanDip_z0 + d_z_tanDip * m_tanDip_tanDip );
672 xpCov( 1, 2 ) = d_py_phi0 * ( d_x_d0 * m_phi0_d0 + d_x_phi0 * m_phi0_phi0 +
673 d_x_omega * m_omega_phi0 + d_x_tanDip * m_tanDip_phi0 ) +
674 d_py_omega * ( d_x_d0 * m_omega_d0 + d_x_phi0 * m_omega_phi0 +
675 d_x_omega * m_omega_omega + d_x_tanDip * m_tanDip_omega ) +
676 d_py_tanDip * ( d_x_d0 * m_tanDip_d0 + d_x_phi0 * m_tanDip_phi0 +
677 d_x_omega * m_tanDip_omega + d_x_tanDip * m_tanDip_tanDip );
678 xpCov( 2, 2 ) = d_py_phi0 * ( d_y_d0 * m_phi0_d0 + d_y_phi0 * m_phi0_phi0 +
679 d_y_omega * m_omega_phi0 + d_y_tanDip * m_tanDip_phi0 ) +
680 d_py_omega * ( d_y_d0 * m_omega_d0 + d_y_phi0 * m_omega_phi0 +
681 d_y_omega * m_omega_omega + d_y_tanDip * m_tanDip_omega ) +
682 d_py_tanDip * ( d_y_d0 * m_tanDip_d0 + d_y_phi0 * m_tanDip_phi0 +
683 d_y_omega * m_tanDip_omega + d_y_tanDip * m_tanDip_tanDip );
684 xpCov( 3, 2 ) = d_py_phi0 * ( d_z_z0 * m_z0_phi0 + d_z_tanDip * m_tanDip_phi0 ) +
685 d_py_omega * ( d_z_z0 * m_z0_omega + d_z_tanDip * m_tanDip_omega ) +
686 d_py_tanDip * ( d_z_z0 * m_tanDip_z0 + d_z_tanDip * m_tanDip_tanDip );
687 xpCov( 1, 3 ) = d_pz_omega * ( d_x_d0 * m_omega_d0 + d_x_phi0 * m_omega_phi0 +
688 d_x_omega * m_omega_omega + d_x_tanDip * m_tanDip_omega ) +
689 d_pz_tanDip * ( d_x_d0 * m_tanDip_d0 + d_x_phi0 * m_tanDip_phi0 +
690 d_x_omega * m_tanDip_omega + d_x_tanDip * m_tanDip_tanDip );
691 xpCov( 2, 3 ) = d_pz_omega * ( d_y_d0 * m_omega_d0 + d_y_phi0 * m_omega_phi0 +
692 d_y_omega * m_omega_omega + d_y_tanDip * m_tanDip_omega ) +
693 d_pz_tanDip * ( d_y_d0 * m_tanDip_d0 + d_y_phi0 * m_tanDip_phi0 +
694 d_y_omega * m_tanDip_omega + d_y_tanDip * m_tanDip_tanDip );
695 xpCov( 3, 3 ) = d_pz_omega * ( d_z_z0 * m_z0_omega + d_z_tanDip * m_tanDip_omega ) +
696 d_pz_tanDip * ( d_z_z0 * m_tanDip_z0 + d_z_tanDip * m_tanDip_tanDip );
697
698 ppCov( 1, 1 ) = d_px_phi0 * ( d_px_phi0 * m_phi0_phi0 + d_px_omega * m_omega_phi0 +
699 d_px_tanDip * m_tanDip_phi0 ) +
700 d_px_omega * ( d_px_phi0 * m_omega_phi0 + d_px_omega * m_omega_omega +
701 d_px_tanDip * m_tanDip_omega ) +
702 d_px_tanDip * ( d_px_phi0 * m_tanDip_phi0 + d_px_omega * m_tanDip_omega +
703 d_px_tanDip * m_tanDip_tanDip );
704 ppCov( 2, 1 ) = d_py_phi0 * ( d_px_phi0 * m_phi0_phi0 + d_px_omega * m_omega_phi0 +
705 d_px_tanDip * m_tanDip_phi0 ) +
706 d_py_omega * ( d_px_phi0 * m_omega_phi0 + d_px_omega * m_omega_omega +
707 d_px_tanDip * m_tanDip_omega ) +
708 d_py_tanDip * ( d_px_phi0 * m_tanDip_phi0 + d_px_omega * m_tanDip_omega +
709 d_px_tanDip * m_tanDip_tanDip );
710 ppCov( 2, 2 ) = d_py_phi0 * ( d_py_phi0 * m_phi0_phi0 + d_py_omega * m_omega_phi0 +
711 d_py_tanDip * m_tanDip_phi0 ) +
712 d_py_omega * ( d_py_phi0 * m_omega_phi0 + d_py_omega * m_omega_omega +
713 d_py_tanDip * m_tanDip_omega ) +
714 d_py_tanDip * ( d_py_phi0 * m_tanDip_phi0 + d_py_omega * m_tanDip_omega +
715 d_py_tanDip * m_tanDip_tanDip );
716 ppCov( 3, 1 ) = d_pz_omega * ( d_px_phi0 * m_omega_phi0 + d_px_omega * m_omega_omega +
717 d_px_tanDip * m_tanDip_omega ) +
718 d_pz_tanDip * ( d_px_phi0 * m_tanDip_phi0 + d_px_omega * m_tanDip_omega +
719 d_px_tanDip * m_tanDip_tanDip );
720 ppCov( 3, 2 ) = d_pz_omega * ( d_py_phi0 * m_omega_phi0 + d_py_omega * m_omega_omega +
721 d_py_tanDip * m_tanDip_omega ) +
722 d_pz_tanDip * ( d_py_phi0 * m_tanDip_phi0 + d_py_omega * m_tanDip_omega +
723 d_py_tanDip * m_tanDip_tanDip );
724 ppCov( 3, 3 ) =
725 d_pz_omega * ( d_pz_omega * m_omega_omega + d_pz_tanDip * m_tanDip_omega ) +
726 d_pz_tanDip * ( d_pz_omega * m_tanDip_omega + d_pz_tanDip * m_tanDip_tanDip );
727
728 // code added for testing...
729 // HepSymMatrix xxCovOld(3),ppCovOld(3);
730 // HepMatrix xpCovOld(3,3);
731
732 // calcCurvAllCovsOLD(theTraj,theField,fltlen,
733 // xxCovOld,ppCovOld,xpCovOld);
734
735 // if (_HdiffCov==0){
736 // _manager = gblEnv->getGen()->ntupleManager();
737 // _HdiffCov = _manager->histogram("log10(|Diff|) old - new Cov",100,-20.,0.);
738 // _HfdiffCov = _manager->histogram("log10(|fDiff|) (old - new)/old Cov",100,-20.,0.);
739 // }
740
741 // for (int i=1;i<=3;i++){
742 // for (int j=i;j<=3;j++){
743 // _HdiffCov->accumulate(log10(fabs(xxCovOld(i,j)-xxCov(i,j))));
744 // _HdiffCov->accumulate(log10(fabs(ppCovOld(i,j)-ppCov(i,j))));
745 // _HfdiffCov->accumulate(log10 ( fabs( (xxCovOld(i,j)-xxCov(i,j))/xxCovOld(i,j)) ) );
746 // _HfdiffCov->accumulate(log10 ( fabs( (ppCovOld(i,j)-ppCov(i,j))/ppCovOld(i,j)) ) );
747 // }
748 // }
749 // for (int i=1;i<=3;i++){
750 // for (int j=1;j<=3;j++){
751 // _HdiffCov->accumulate(log10(fabs(xpCovOld(i,j)-xpCov(i,j))));
752 // _HfdiffCov->accumulate(log10 ( fabs( (xpCovOld(i,j)-xpCov(i,j))/xpCovOld(i,j)) ) );
753 // }
754 // }
755}
756
757//------------------------------------------------------------------------
758void TrkMomCalculator::calcCurvAllCovsOLD( const TrkSimpTraj& theTraj, const BField& theField,
759 double fltlen, HepSymMatrix& xxCov,
760 HepSymMatrix& ppCov, HepMatrix& xpCov ) {
761 //------------------------------------------------------------------------
762
763 DifPoint PosDif;
764 DifVector DirDif;
765 DifVector delDirDif;
766 DifVector MomDif( 0., 0., 0. );
767
768 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
769 if ( delDirDif.length() != 0 )
770 {
771 DifNumber sindip = DirDif.z;
772 DifNumber arg = 1.0 - sindip * sindip;
773 if ( arg.number() < 0.0 ) { arg.setNumber( 0.0 ); }
774 DifNumber cosdip = sqrt( arg );
775
776 DifNumber momMag =
777 BField::cmTeslaToGeVc * theField.bFieldNominal() * cosdip / delDirDif.length();
778 momMag.absolute();
779
780 MomDif = DirDif * momMag;
781 }
782
783 // computes position covariances
784
785 xxCov.assign( PosDif.errorMatrix( PosDif.x.indepPar()->covariance() ) );
786
787 // computes momentum covariances
788 ppCov.assign( MomDif.errorMatrix( MomDif.x.indepPar()->covariance() ) );
789
790 // computes correlations
791 const HepSymMatrix& nnCov = MomDif.x.indepPar()->covariance();
792 xpCov( 1, 1 ) = correlation( PosDif.x, MomDif.x, nnCov );
793 xpCov( 1, 2 ) = correlation( PosDif.x, MomDif.y, nnCov );
794 xpCov( 1, 3 ) = correlation( PosDif.x, MomDif.z, nnCov );
795 xpCov( 2, 1 ) = correlation( PosDif.y, MomDif.x, nnCov );
796 xpCov( 2, 2 ) = correlation( PosDif.y, MomDif.y, nnCov );
797 xpCov( 2, 3 ) = correlation( PosDif.y, MomDif.z, nnCov );
798 xpCov( 3, 1 ) = correlation( PosDif.z, MomDif.x, nnCov );
799 xpCov( 3, 2 ) = correlation( PosDif.z, MomDif.y, nnCov );
800 xpCov( 3, 3 ) = correlation( PosDif.z, MomDif.z, nnCov );
801}
802
803//------------------------------------------------------------------------
804void TrkMomCalculator::calcNeutAllCovs( const TrkSimpTraj& theTraj, const BField& theField,
805 double fltlen, HepSymMatrix& xxCov,
806 HepSymMatrix& ppCov, HepMatrix& xpCov ) {
807 //------------------------------------------------------------------------
808 HepVector p0( 3 ), x0( 3 );
809 calcNeutAllCovs( theTraj, theField, fltlen, x0, p0, xxCov, ppCov, xpCov );
810}
811
812//------------------------------------------------------------------------
813void TrkMomCalculator::calcNeutAllCovs( const TrkSimpTraj& theTraj, const BField& theField,
814 double fltlen, HepVector& x0, HepVector& p0,
815 HepSymMatrix& xxCov, HepSymMatrix& ppCov,
816 HepMatrix& xpCov ) {
817 //------------------------------------------------------------------------
818
819 DifPoint PosDif;
820 DifVector DirDif;
821 DifVector delDirDif;
822
823 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
824
825 // set the momentum's direction, and then its magnitude
826 DifVector MomDif = DirDif;
827
828 // The + 1 is necessary because DifIndepPar starts counting at 1, whereas
829 // the enum for NeutParams starts at 0. Stupid, ain't it?
830 MomDif *= theTraj.parameters()->difPar( NeutParams::_p + 1 );
831
832 // fill momenta and positions
833 p0[0] = MomDif.x.number();
834 p0[1] = MomDif.y.number();
835 p0[2] = MomDif.z.number();
836 x0[0] = PosDif.x.number();
837 x0[1] = PosDif.y.number();
838 x0[2] = PosDif.z.number();
839
840 // computes momentum covariances
841 xxCov.assign( PosDif.errorMatrix( PosDif.x.indepPar()->covariance() ) );
842
843 // computes momentum covariances
844 ppCov.assign( MomDif.errorMatrix( MomDif.x.indepPar()->covariance() ) );
845
846 // computes correlations
847 const HepSymMatrix& nnCov = MomDif.x.indepPar()->covariance();
848 xpCov( 1, 1 ) = correlation( PosDif.x, MomDif.x, nnCov );
849 xpCov( 1, 2 ) = correlation( PosDif.x, MomDif.y, nnCov );
850 xpCov( 1, 3 ) = correlation( PosDif.x, MomDif.z, nnCov );
851 xpCov( 2, 1 ) = correlation( PosDif.y, MomDif.x, nnCov );
852 xpCov( 2, 2 ) = correlation( PosDif.y, MomDif.y, nnCov );
853 xpCov( 2, 3 ) = correlation( PosDif.y, MomDif.z, nnCov );
854 xpCov( 3, 1 ) = correlation( PosDif.z, MomDif.x, nnCov );
855 xpCov( 3, 2 ) = correlation( PosDif.z, MomDif.y, nnCov );
856 xpCov( 3, 3 ) = correlation( PosDif.z, MomDif.z, nnCov );
857}
858
859//------------------------------------------------------------------------
860void TrkMomCalculator::calcCurvAllWeights( const TrkSimpTraj& theTraj, const BField& theField,
861 double fltlen, HepVector& pos, HepVector& mom,
862 HepSymMatrix& xxWeight, HepSymMatrix& ppWeight,
863 HepMatrix& xpWeight ) {
864 //------------------------------------------------------------------------
865 const HepVector& v = theTraj.parameters()->parameter();
866 const HepSymMatrix& w = theTraj.parameters()->weightMatrix();
867
868 // fast inlined conversion from Helix to PX representation
869
870 // track parameters
871 double s = v[TrkExchangePar::ex_tanDip];
872 double omega = v[TrkExchangePar::ex_omega];
873 double d0 = v[TrkExchangePar::ex_d0];
874 double z0 = v[TrkExchangePar::ex_z0];
875 double phi0 = v[TrkExchangePar::ex_phi0];
876 double l = fltlen / sqrt( 1.0 + s * s );
877
878 double phi = phi0 + omega * l;
879 double sinphi0 = sin( phi0 );
880 double cosphi0 = cos( phi0 );
881 double sinphi = sin( phi );
882 double cosphi = cos( phi );
883 double C = BField::cmTeslaToGeVc * theField.bFieldNominal();
884 double q( 1.0 );
885 -omega > 0 ? q = 1.0 : q = -1.0;
886 double qC = q * C;
887 double pt = fabs( -qC / omega );
888 double r = 1.0 / omega;
889
890 // calculate position and momentum
891 pos( 1 ) = r * sinphi - ( r + d0 ) * sinphi0;
892 pos( 2 ) = -r * cosphi + ( r + d0 ) * cosphi0;
893 pos( 3 ) = z0 + l * s;
894 mom( 1 ) = pt * cosphi;
895 mom( 2 ) = pt * sinphi;
896 mom( 3 ) = pt * s;
897
898 // inverse of the jacobian matrix - see helix.mws Maple worksheet
899
900 // protect against divide by 0.
901 if ( ( 1 + d0 * omega ) == 0.0 )
902 {
903 calcCurvAllWeightsOLD( theTraj, theField, fltlen, pos, mom, xxWeight, ppWeight, xpWeight );
904 return;
905 }
906
907 double dinv_x_d0 = -sinphi0;
908 double dinv_x_phi0 = -omega * cosphi0 / ( 1 + d0 * omega );
909 double dinv_x_z0 = -s * cosphi0 / ( 1 + d0 * omega );
910
911 double dinv_y_d0 = cosphi0;
912 double dinv_y_phi0 = -omega * sinphi0 / ( 1 + d0 * omega );
913 double dinv_y_z0 = -s * sinphi0 / ( 1 + d0 * omega );
914
915 double dinv_z_z0 = 1;
916
917 double dinv_px_d0 = ( cosphi - cosphi0 ) / qC;
918 double dinv_px_phi0 = omega * sinphi0 / ( 1 + d0 * omega ) / qC;
919 double dinv_px_omega = omega * omega * cosphi / qC;
920 double dinv_px_z0 =
921 -s * ( sinphi * ( 1 + d0 * omega ) - sinphi0 ) / ( qC * ( 1 + d0 * omega ) );
922 double dinv_px_tanDip = omega * cosphi * s / qC;
923
924 double dinv_py_d0 = ( sinphi - sinphi0 ) / qC;
925 double dinv_py_phi0 = -omega * cosphi0 / ( 1 + d0 * omega ) / qC;
926 double dinv_py_omega = omega * omega * sinphi / qC;
927 double dinv_py_z0 =
928 s * ( cosphi * ( 1 + d0 * omega ) - cosphi0 ) / ( qC * ( 1 + d0 * omega ) );
929 double dinv_py_tanDip = omega * sinphi * s / qC;
930
931 double dinv_pz_z0 = l * omega / qC;
932 double dinv_pz_tanDip = -omega / qC;
933
934 // local variables for the weight matrix
935 double w_d0_d0 = w[TrkExchangePar::ex_d0][TrkExchangePar::ex_d0];
936 double w_phi0_d0 = w[TrkExchangePar::ex_phi0][TrkExchangePar::ex_d0];
937 double w_phi0_phi0 = w[TrkExchangePar::ex_phi0][TrkExchangePar::ex_phi0];
938 double w_omega_d0 = w[TrkExchangePar::ex_omega][TrkExchangePar::ex_d0];
939 double w_omega_phi0 = w[TrkExchangePar::ex_omega][TrkExchangePar::ex_phi0];
940 double w_omega_omega = w[TrkExchangePar::ex_omega][TrkExchangePar::ex_omega];
941 double w_z0_d0 = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_d0];
942 double w_z0_phi0 = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_phi0];
943 double w_z0_omega = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_omega];
944 double w_z0_z0 = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_z0];
945 double w_tanDip_d0 = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_d0];
946 double w_tanDip_phi0 = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_phi0];
947 double w_tanDip_omega = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_omega];
948 double w_tanDip_z0 = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_z0];
949 double w_tanDip_tanDip = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_tanDip];
950
951 // calculate the Weight matrix from dinv (the indices are xpWeight(ip,ix) ) (used writewgt.py
952 // script)
953 xxWeight( 1, 1 ) =
954 dinv_x_d0 * ( dinv_x_d0 * w_d0_d0 + dinv_x_phi0 * w_phi0_d0 + dinv_x_z0 * w_z0_d0 ) +
955 dinv_x_phi0 *
956 ( dinv_x_d0 * w_phi0_d0 + dinv_x_phi0 * w_phi0_phi0 + dinv_x_z0 * w_z0_phi0 ) +
957 dinv_x_z0 * ( dinv_x_d0 * w_z0_d0 + dinv_x_phi0 * w_z0_phi0 + dinv_x_z0 * w_z0_z0 );
958 xxWeight( 2, 1 ) =
959 dinv_y_d0 * ( dinv_x_d0 * w_d0_d0 + dinv_x_phi0 * w_phi0_d0 + dinv_x_z0 * w_z0_d0 ) +
960 dinv_y_phi0 *
961 ( dinv_x_d0 * w_phi0_d0 + dinv_x_phi0 * w_phi0_phi0 + dinv_x_z0 * w_z0_phi0 ) +
962 dinv_y_z0 * ( dinv_x_d0 * w_z0_d0 + dinv_x_phi0 * w_z0_phi0 + dinv_x_z0 * w_z0_z0 );
963 xxWeight( 2, 2 ) =
964 dinv_y_d0 * ( dinv_y_d0 * w_d0_d0 + dinv_y_phi0 * w_phi0_d0 + dinv_y_z0 * w_z0_d0 ) +
965 dinv_y_phi0 *
966 ( dinv_y_d0 * w_phi0_d0 + dinv_y_phi0 * w_phi0_phi0 + dinv_y_z0 * w_z0_phi0 ) +
967 dinv_y_z0 * ( dinv_y_d0 * w_z0_d0 + dinv_y_phi0 * w_z0_phi0 + dinv_y_z0 * w_z0_z0 );
968 xxWeight( 3, 1 ) =
969 dinv_z_z0 * ( dinv_x_d0 * w_z0_d0 + dinv_x_phi0 * w_z0_phi0 + dinv_x_z0 * w_z0_z0 );
970 xxWeight( 3, 2 ) =
971 dinv_z_z0 * ( dinv_y_d0 * w_z0_d0 + dinv_y_phi0 * w_z0_phi0 + dinv_y_z0 * w_z0_z0 );
972 xxWeight( 3, 3 ) = dinv_z_z0 * ( dinv_z_z0 * w_z0_z0 );
973
974 xpWeight( 1, 1 ) =
975 dinv_px_d0 * ( dinv_x_d0 * w_d0_d0 + dinv_x_phi0 * w_phi0_d0 + dinv_x_z0 * w_z0_d0 ) +
976 dinv_px_phi0 *
977 ( dinv_x_d0 * w_phi0_d0 + dinv_x_phi0 * w_phi0_phi0 + dinv_x_z0 * w_z0_phi0 ) +
978 dinv_px_omega *
979 ( dinv_x_d0 * w_omega_d0 + dinv_x_phi0 * w_omega_phi0 + dinv_x_z0 * w_z0_omega ) +
980 dinv_px_z0 * ( dinv_x_d0 * w_z0_d0 + dinv_x_phi0 * w_z0_phi0 + dinv_x_z0 * w_z0_z0 ) +
981 dinv_px_tanDip *
982 ( dinv_x_d0 * w_tanDip_d0 + dinv_x_phi0 * w_tanDip_phi0 + dinv_x_z0 * w_tanDip_z0 );
983 xpWeight( 2, 1 ) =
984 dinv_px_d0 * ( dinv_y_d0 * w_d0_d0 + dinv_y_phi0 * w_phi0_d0 + dinv_y_z0 * w_z0_d0 ) +
985 dinv_px_phi0 *
986 ( dinv_y_d0 * w_phi0_d0 + dinv_y_phi0 * w_phi0_phi0 + dinv_y_z0 * w_z0_phi0 ) +
987 dinv_px_omega *
988 ( dinv_y_d0 * w_omega_d0 + dinv_y_phi0 * w_omega_phi0 + dinv_y_z0 * w_z0_omega ) +
989 dinv_px_z0 * ( dinv_y_d0 * w_z0_d0 + dinv_y_phi0 * w_z0_phi0 + dinv_y_z0 * w_z0_z0 ) +
990 dinv_px_tanDip *
991 ( dinv_y_d0 * w_tanDip_d0 + dinv_y_phi0 * w_tanDip_phi0 + dinv_y_z0 * w_tanDip_z0 );
992 xpWeight( 3, 1 ) =
993 dinv_px_d0 * ( dinv_z_z0 * w_z0_d0 ) + dinv_px_phi0 * ( dinv_z_z0 * w_z0_phi0 ) +
994 dinv_px_omega * ( dinv_z_z0 * w_z0_omega ) + dinv_px_z0 * ( dinv_z_z0 * w_z0_z0 ) +
995 dinv_px_tanDip * ( dinv_z_z0 * w_tanDip_z0 );
996 xpWeight( 1, 2 ) =
997 dinv_py_d0 * ( dinv_x_d0 * w_d0_d0 + dinv_x_phi0 * w_phi0_d0 + dinv_x_z0 * w_z0_d0 ) +
998 dinv_py_phi0 *
999 ( dinv_x_d0 * w_phi0_d0 + dinv_x_phi0 * w_phi0_phi0 + dinv_x_z0 * w_z0_phi0 ) +
1000 dinv_py_omega *
1001 ( dinv_x_d0 * w_omega_d0 + dinv_x_phi0 * w_omega_phi0 + dinv_x_z0 * w_z0_omega ) +
1002 dinv_py_z0 * ( dinv_x_d0 * w_z0_d0 + dinv_x_phi0 * w_z0_phi0 + dinv_x_z0 * w_z0_z0 ) +
1003 dinv_py_tanDip *
1004 ( dinv_x_d0 * w_tanDip_d0 + dinv_x_phi0 * w_tanDip_phi0 + dinv_x_z0 * w_tanDip_z0 );
1005 xpWeight( 2, 2 ) =
1006 dinv_py_d0 * ( dinv_y_d0 * w_d0_d0 + dinv_y_phi0 * w_phi0_d0 + dinv_y_z0 * w_z0_d0 ) +
1007 dinv_py_phi0 *
1008 ( dinv_y_d0 * w_phi0_d0 + dinv_y_phi0 * w_phi0_phi0 + dinv_y_z0 * w_z0_phi0 ) +
1009 dinv_py_omega *
1010 ( dinv_y_d0 * w_omega_d0 + dinv_y_phi0 * w_omega_phi0 + dinv_y_z0 * w_z0_omega ) +
1011 dinv_py_z0 * ( dinv_y_d0 * w_z0_d0 + dinv_y_phi0 * w_z0_phi0 + dinv_y_z0 * w_z0_z0 ) +
1012 dinv_py_tanDip *
1013 ( dinv_y_d0 * w_tanDip_d0 + dinv_y_phi0 * w_tanDip_phi0 + dinv_y_z0 * w_tanDip_z0 );
1014 xpWeight( 3, 2 ) =
1015 dinv_py_d0 * ( dinv_z_z0 * w_z0_d0 ) + dinv_py_phi0 * ( dinv_z_z0 * w_z0_phi0 ) +
1016 dinv_py_omega * ( dinv_z_z0 * w_z0_omega ) + dinv_py_z0 * ( dinv_z_z0 * w_z0_z0 ) +
1017 dinv_py_tanDip * ( dinv_z_z0 * w_tanDip_z0 );
1018 xpWeight( 1, 3 ) =
1019 dinv_pz_z0 * ( dinv_x_d0 * w_z0_d0 + dinv_x_phi0 * w_z0_phi0 + dinv_x_z0 * w_z0_z0 ) +
1020 dinv_pz_tanDip *
1021 ( dinv_x_d0 * w_tanDip_d0 + dinv_x_phi0 * w_tanDip_phi0 + dinv_x_z0 * w_tanDip_z0 );
1022 xpWeight( 2, 3 ) =
1023 dinv_pz_z0 * ( dinv_y_d0 * w_z0_d0 + dinv_y_phi0 * w_z0_phi0 + dinv_y_z0 * w_z0_z0 ) +
1024 dinv_pz_tanDip *
1025 ( dinv_y_d0 * w_tanDip_d0 + dinv_y_phi0 * w_tanDip_phi0 + dinv_y_z0 * w_tanDip_z0 );
1026 xpWeight( 3, 3 ) =
1027 dinv_pz_z0 * ( dinv_z_z0 * w_z0_z0 ) + dinv_pz_tanDip * ( dinv_z_z0 * w_tanDip_z0 );
1028
1029 ppWeight( 1, 1 ) =
1030 dinv_px_d0 *
1031 ( dinv_px_d0 * w_d0_d0 + dinv_px_phi0 * w_phi0_d0 + dinv_px_omega * w_omega_d0 +
1032 dinv_px_z0 * w_z0_d0 + dinv_px_tanDip * w_tanDip_d0 ) +
1033 dinv_px_phi0 * ( dinv_px_d0 * w_phi0_d0 + dinv_px_phi0 * w_phi0_phi0 +
1034 dinv_px_omega * w_omega_phi0 + dinv_px_z0 * w_z0_phi0 +
1035 dinv_px_tanDip * w_tanDip_phi0 ) +
1036 dinv_px_omega * ( dinv_px_d0 * w_omega_d0 + dinv_px_phi0 * w_omega_phi0 +
1037 dinv_px_omega * w_omega_omega + dinv_px_z0 * w_z0_omega +
1038 dinv_px_tanDip * w_tanDip_omega ) +
1039 dinv_px_z0 *
1040 ( dinv_px_d0 * w_z0_d0 + dinv_px_phi0 * w_z0_phi0 + dinv_px_omega * w_z0_omega +
1041 dinv_px_z0 * w_z0_z0 + dinv_px_tanDip * w_tanDip_z0 ) +
1042 dinv_px_tanDip * ( dinv_px_d0 * w_tanDip_d0 + dinv_px_phi0 * w_tanDip_phi0 +
1043 dinv_px_omega * w_tanDip_omega + dinv_px_z0 * w_tanDip_z0 +
1044 dinv_px_tanDip * w_tanDip_tanDip );
1045 ppWeight( 2, 1 ) =
1046 dinv_py_d0 *
1047 ( dinv_px_d0 * w_d0_d0 + dinv_px_phi0 * w_phi0_d0 + dinv_px_omega * w_omega_d0 +
1048 dinv_px_z0 * w_z0_d0 + dinv_px_tanDip * w_tanDip_d0 ) +
1049 dinv_py_phi0 * ( dinv_px_d0 * w_phi0_d0 + dinv_px_phi0 * w_phi0_phi0 +
1050 dinv_px_omega * w_omega_phi0 + dinv_px_z0 * w_z0_phi0 +
1051 dinv_px_tanDip * w_tanDip_phi0 ) +
1052 dinv_py_omega * ( dinv_px_d0 * w_omega_d0 + dinv_px_phi0 * w_omega_phi0 +
1053 dinv_px_omega * w_omega_omega + dinv_px_z0 * w_z0_omega +
1054 dinv_px_tanDip * w_tanDip_omega ) +
1055 dinv_py_z0 *
1056 ( dinv_px_d0 * w_z0_d0 + dinv_px_phi0 * w_z0_phi0 + dinv_px_omega * w_z0_omega +
1057 dinv_px_z0 * w_z0_z0 + dinv_px_tanDip * w_tanDip_z0 ) +
1058 dinv_py_tanDip * ( dinv_px_d0 * w_tanDip_d0 + dinv_px_phi0 * w_tanDip_phi0 +
1059 dinv_px_omega * w_tanDip_omega + dinv_px_z0 * w_tanDip_z0 +
1060 dinv_px_tanDip * w_tanDip_tanDip );
1061 ppWeight( 2, 2 ) =
1062 dinv_py_d0 *
1063 ( dinv_py_d0 * w_d0_d0 + dinv_py_phi0 * w_phi0_d0 + dinv_py_omega * w_omega_d0 +
1064 dinv_py_z0 * w_z0_d0 + dinv_py_tanDip * w_tanDip_d0 ) +
1065 dinv_py_phi0 * ( dinv_py_d0 * w_phi0_d0 + dinv_py_phi0 * w_phi0_phi0 +
1066 dinv_py_omega * w_omega_phi0 + dinv_py_z0 * w_z0_phi0 +
1067 dinv_py_tanDip * w_tanDip_phi0 ) +
1068 dinv_py_omega * ( dinv_py_d0 * w_omega_d0 + dinv_py_phi0 * w_omega_phi0 +
1069 dinv_py_omega * w_omega_omega + dinv_py_z0 * w_z0_omega +
1070 dinv_py_tanDip * w_tanDip_omega ) +
1071 dinv_py_z0 *
1072 ( dinv_py_d0 * w_z0_d0 + dinv_py_phi0 * w_z0_phi0 + dinv_py_omega * w_z0_omega +
1073 dinv_py_z0 * w_z0_z0 + dinv_py_tanDip * w_tanDip_z0 ) +
1074 dinv_py_tanDip * ( dinv_py_d0 * w_tanDip_d0 + dinv_py_phi0 * w_tanDip_phi0 +
1075 dinv_py_omega * w_tanDip_omega + dinv_py_z0 * w_tanDip_z0 +
1076 dinv_py_tanDip * w_tanDip_tanDip );
1077 ppWeight( 3, 1 ) =
1078 dinv_pz_z0 *
1079 ( dinv_px_d0 * w_z0_d0 + dinv_px_phi0 * w_z0_phi0 + dinv_px_omega * w_z0_omega +
1080 dinv_px_z0 * w_z0_z0 + dinv_px_tanDip * w_tanDip_z0 ) +
1081 dinv_pz_tanDip * ( dinv_px_d0 * w_tanDip_d0 + dinv_px_phi0 * w_tanDip_phi0 +
1082 dinv_px_omega * w_tanDip_omega + dinv_px_z0 * w_tanDip_z0 +
1083 dinv_px_tanDip * w_tanDip_tanDip );
1084 ppWeight( 3, 2 ) =
1085 dinv_pz_z0 *
1086 ( dinv_py_d0 * w_z0_d0 + dinv_py_phi0 * w_z0_phi0 + dinv_py_omega * w_z0_omega +
1087 dinv_py_z0 * w_z0_z0 + dinv_py_tanDip * w_tanDip_z0 ) +
1088 dinv_pz_tanDip * ( dinv_py_d0 * w_tanDip_d0 + dinv_py_phi0 * w_tanDip_phi0 +
1089 dinv_py_omega * w_tanDip_omega + dinv_py_z0 * w_tanDip_z0 +
1090 dinv_py_tanDip * w_tanDip_tanDip );
1091 ppWeight( 3, 3 ) =
1092 dinv_pz_z0 * ( dinv_pz_z0 * w_z0_z0 + dinv_pz_tanDip * w_tanDip_z0 ) +
1093 dinv_pz_tanDip * ( dinv_pz_z0 * w_tanDip_z0 + dinv_pz_tanDip * w_tanDip_tanDip );
1094
1095 // code added for testing...
1096 // HepVector posOld(3),momOld(3);
1097 // HepSymMatrix xxWeightOld(3),ppWeightOld(3);
1098 // HepMatrix xpWeightOld(3,3);
1099
1100 // calcCurvAllWeightsOLD(theTraj,theField,fltlen,
1101 // posOld,momOld,xxWeightOld,ppWeightOld,xpWeightOld);
1102
1103 // if (_Hdiff==0){
1104 // _manager = gblEnv->getGen()->ntupleManager();
1105 // _Hdiff = _manager->histogram("log10(|Diff|) old - new",100,-20.,0.);
1106 // _Hfdiff = _manager->histogram("log10(|fDiff|) (old - new)/old",100,-20.,0.);
1107 // }
1108
1109 // for (int i=1;i<=3;i++){
1110 // for (int j=i;j<=3;j++){
1111 // _Hdiff->accumulate(log10(fabs(xxWeightOld(i,j)-xxWeight(i,j))));
1112 // _Hdiff->accumulate(log10(fabs(ppWeightOld(i,j)-ppWeight(i,j))));
1113 // _Hfdiff->accumulate(log10 ( fabs( (xxWeightOld(i,j)-xxWeight(i,j))/xxWeightOld(i,j))
1114 // ) ); _Hfdiff->accumulate(log10 ( fabs(
1115 // (ppWeightOld(i,j)-ppWeight(i,j))/ppWeightOld(i,j)) ) );
1116 // }
1117 // }
1118 // for (int i=1;i<=3;i++){
1119 // for (int j=1;j<=3;j++){
1120 // _Hdiff->accumulate(log10(fabs(xpWeightOld(i,j)-xpWeight(i,j))));
1121 // _Hfdiff->accumulate(log10 ( fabs( (xpWeightOld(i,j)-xpWeight(i,j))/xpWeightOld(i,j))
1122 // ) );
1123 // }
1124 // }
1125}
1126
1127//------------------------------------------------------------------------
1128void TrkMomCalculator::calcCurvAllWeightsOLD( const TrkSimpTraj& theTraj,
1129 const BField& theField, double fltlen,
1130 HepVector& pos, HepVector& mom,
1131 HepSymMatrix& xxWeight, HepSymMatrix& ppWeight,
1132 HepMatrix& xpWeight ) {
1133 //------------------------------------------------------------------------
1134
1135 DifPoint PosDif;
1136 DifVector DirDif;
1137 DifVector delDirDif;
1138 DifNumber momMag;
1139 DifVector MomDif;
1140
1141 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
1142 if ( delDirDif.length() != 0 )
1143 {
1144 DifNumber sindip = DirDif.z;
1145 DifNumber arg = 1.0 - sindip * sindip;
1146 if ( arg.number() < 0.0 ) { arg.setNumber( 0.0 ); }
1147 DifNumber cosdip = sqrt( arg );
1148
1149 momMag = BField::cmTeslaToGeVc * theField.bFieldNominal() * cosdip / delDirDif.length();
1150 momMag.absolute();
1151
1152 MomDif = DirDif * momMag;
1153 }
1154
1155 /*
1156 * start computing the inverse of the Jacobian needed:
1157 *
1158 * D(x,p)
1159 * ------
1160 * D(n)
1161 */
1162 HepMatrix Jx_n( PosDif.jacobian() );
1163 HepMatrix Jp_n( MomDif.jacobian() );
1164
1165 int i, j;
1166 HepMatrix Jxp_ns( 6, 6 );
1167
1168 for ( i = 0; i < 3; i++ )
1169 for ( j = 0; j < 5; j++ )
1170 {
1171 Jxp_ns[i][j] = Jx_n[i][j];
1172 Jxp_ns[i + 3][j] = Jp_n[i][j];
1173 }
1174
1175 /*
1176 * now we need:
1177 *
1178 * D(x,p)
1179 * ------
1180 * D(s)
1181 *
1182 */
1183
1184 Jxp_ns[0][5] = DirDif.x.number();
1185 Jxp_ns[1][5] = DirDif.y.number();
1186 Jxp_ns[2][5] = DirDif.z.number();
1187
1188 Jxp_ns[3][5] = momMag.number() * delDirDif.x.number();
1189 Jxp_ns[4][5] = momMag.number() * delDirDif.y.number();
1190 Jxp_ns[5][5] = momMag.number() * delDirDif.z.number();
1191
1192 /*
1193 * with an inversion we obtain
1194 * D(n)
1195 * --------
1196 * D(x,p)
1197 *
1198 */
1199 int invStatus;
1200
1201 Jxp_ns.invert( invStatus );
1202
1203 HepMatrix Jn_x( 5, 3 );
1204 HepMatrix Jn_p( 5, 3 );
1205
1206 for ( i = 0; i < 5; i++ )
1207 for ( j = 0; j < 3; j++ )
1208 {
1209 Jn_x[i][j] = Jxp_ns[i][j];
1210 Jn_p[i][j] = Jxp_ns[i][j + 3];
1211 }
1212 // this is the weight matrix for the helix parameters
1213 HepSymMatrix Wnn( PosDif.x.indepPar()->covariance() );
1214
1215 Wnn.invert( invStatus );
1216 /*
1217 * now we have the weight matrices
1218 *
1219 */
1220 xxWeight = Wnn.similarityT( Jn_x );
1221 ppWeight = Wnn.similarityT( Jn_p );
1222 xpWeight = Jn_x.T() * Wnn * Jn_p;
1223
1224 pos[0] = PosDif.x.number();
1225 pos[1] = PosDif.y.number();
1226 pos[2] = PosDif.z.number();
1227
1228 mom[0] = MomDif.x.number();
1229 mom[1] = MomDif.y.number();
1230 mom[2] = MomDif.z.number();
1231}
1232
1233//------------------------------------------------------------------------
1234void TrkMomCalculator::calcNeutAllWeights( const TrkSimpTraj& theTraj, const BField& theField,
1235 double fltlen, HepVector& pos, HepVector& mom,
1236 HepSymMatrix& xxWeight, HepSymMatrix& ppWeight,
1237 HepMatrix& xpWeight ) {
1238 //------------------------------------------------------------------------
1239 DifPoint PosDif;
1240 DifVector DirDif;
1241 DifVector delDirDif;
1242 DifNumber momMag;
1243
1244 theTraj.getDFInfo( fltlen, PosDif, DirDif, delDirDif );
1245
1246 // set the momentum's direction, and then its magnitude
1247 DifVector MomDif = DirDif;
1248
1249 MomDif *= theTraj.parameters()->difPar( NeutParams::_p + 1 );
1250
1251 HepMatrix Jx_n( PosDif.jacobian() );
1252 HepMatrix Jp_n( MomDif.jacobian() );
1253
1254 int i, j;
1255 HepMatrix Jxp_ns( 6, 6 );
1256
1257 for ( i = 0; i < 3; i++ )
1258 for ( j = 0; j < 6; j++ )
1259 {
1260 Jxp_ns[i][j] = Jx_n[i][j];
1261 Jxp_ns[i + 3][j] = Jp_n[i][j];
1262 }
1263 int invStatus;
1264
1265 Jxp_ns.invert( invStatus );
1266
1267 HepMatrix Jn_x( 5, 3 );
1268 HepMatrix Jn_p( 5, 3 );
1269
1270 for ( i = 0; i < 5; i++ )
1271 for ( j = 0; j < 3; j++ )
1272 {
1273 Jn_x[i][j] = Jxp_ns[i][j];
1274 Jn_p[i][j] = Jxp_ns[i][j + 3];
1275 }
1276 // this is the weight matrix for the helix parameters
1277 HepSymMatrix Wnn( PosDif.x.indepPar()->covariance().sub( 1, 5 ) );
1278
1279 Wnn.invert( invStatus );
1280 xxWeight = Wnn.similarityT( Jn_x );
1281 ppWeight = Wnn.similarityT( Jn_p );
1282 xpWeight = Jn_x.T() * Wnn * Jn_p;
1283
1284 pos[0] = PosDif.x.number();
1285 pos[1] = PosDif.y.number();
1286 pos[2] = PosDif.z.number();
1287
1288 mom[0] = MomDif.x.number();
1289 mom[1] = MomDif.y.number();
1290 mom[2] = MomDif.z.number();
1291}
correlation
double correlation(const DifNumber &a, const DifNumber &b)
Definition DifNumber.cxx:69
arg
double arg(const EvtComplex &c)
Definition EvtComplex.hh:195
alpha
double alpha
Definition EvtPhokharaDef.hh:17
w
double w
Definition EvtPhokharaDef.hh:28
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
DifIndepPar::difPar
DifNumber difPar(int i) const
Definition DifIndepPar.cxx:52
DifVector::errorMatrix
HepSymMatrix errorMatrix(const HepSymMatrix &e) const
Definition DifVector.cxx:43
DifVector::length
DifNumber length() const
DifVector::jacobian
HepMatrix jacobian() const
Definition DifVector.cxx:54
Trajectory::direction
virtual Hep3Vector direction(double) const =0
Trajectory::curvature
virtual double curvature(double) const =0
TrkDifTraj::getDFInfo
virtual void getDFInfo(double fltLen, DifPoint &pos, DifVector &direction, DifVector &delDirect) const =0
TrkMomCalculator::vecMom
static Hep3Vector vecMom(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen)
Definition TrkMomCalculator.cxx:82
TrkMomCalculator::ptMom
static double ptMom(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen)
Definition TrkMomCalculator.cxx:48
TrkMomCalculator::weightToCov
static bool weightToCov(const HepSymMatrix &inXX, const HepSymMatrix &inPP, const HepMatrix &inXP, HepSymMatrix &outXX, HepSymMatrix &outPP, HepMatrix &outXP)
Definition TrkMomCalculator.cxx:533
TrkMomCalculator::getAllWeights
static void getAllWeights(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen, HepVector &pos, HepVector &mom, HepSymMatrix &xxWeight, HepSymMatrix &ppWeight, HepMatrix &xpWeight)
Definition TrkMomCalculator.cxx:255
TrkMomCalculator::getAllCovs
static void getAllCovs(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen, HepSymMatrix &xxCov, HepSymMatrix &ppCov, HepMatrix &xpCov)
Definition TrkMomCalculator.cxx:211
TrkMomCalculator::charge
static int charge(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen)
Definition TrkMomCalculator.cxx:143
TrkMomCalculator::posmomCov
static HepMatrix posmomCov(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen)
Definition TrkMomCalculator.cxx:182
TrkMomCalculator::errMom
static BesVectorErr errMom(const TrkSimpTraj &, const MdcPatRec::BField &, double fltlen)
Definition TrkMomCalculator.cxx:113
TrkParams::weightMatrix
const HepSymMatrix & weightMatrix() const
Definition TrkParams.cxx:73