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