00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "TrkBase/TrkSimpTraj.h"
00018 #include "TrkBase/TrkMomCalculator.h"
00019 #include "TrkBase/TrkMomVisitor.h"
00020 #include "TrkBase/TrkExchangePar.h"
00021 #include "TrkBase/NeutParams.h"
00022 #include "BField/BField.h"
00023 #include "CLHEP/Vector/ThreeVector.h"
00024 #include "CLHEP/Matrix/Vector.h"
00025 #include "CLHEP/Geometry/Point3D.h"
00026 #include "MdcRecoUtil/DifPoint.h"
00027 #include "MdcRecoUtil/DifNumber.h"
00028 #include "MdcRecoUtil/DifVector.h"
00029 #include "MdcRecoUtil/DifIndepPar.h"
00030 #include "MdcRecoUtil/BesVectorErr.h"
00031 #include "TrkBase/HelixTraj.h"
00032
00033 using std::endl;
00034
00035
00036 TrkMomCalculator::~TrkMomCalculator() {
00037
00038 }
00039
00040
00041 TrkMomCalculator::TrkMomCalculator() {
00042
00043 }
00044
00045
00046 double
00047 TrkMomCalculator::ptMom(const TrkSimpTraj& theTraj, const BField&
00048 theField, double fltlen) {
00049
00050
00051 TrkMomVisitor theVisitor(theTraj);
00052
00053 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) {
00054
00055
00056
00057 return calcCurvPtMom(theTraj.direction(fltlen),
00058 theTraj.curvature(fltlen), theField);
00059
00060 } else if (theVisitor.neut() != 0) {
00061
00062
00063 double sindip = theTraj.direction(fltlen).z();
00064 double arg = 1.0-sindip*sindip;
00065 if (arg < 0.0) arg = 0.0;
00066 double cosdip = sqrt(arg);
00067 double ptot = theTraj.parameters()->parameter()[NeutParams::_p];
00068
00069 return cosdip * ptot;
00070
00071 } else {
00072
00073
00074 return 999999.99;
00075
00076 }
00077
00078 }
00079
00080
00081 Hep3Vector
00082 TrkMomCalculator::vecMom(const TrkSimpTraj& theTraj, const BField&
00083 theField, double fltlen) {
00084
00085 TrkMomVisitor theVisitor(theTraj);
00086
00087 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) {
00088
00089
00090
00091
00092 return calcCurvVecMom(theTraj.direction(fltlen),
00093 theTraj.curvature(fltlen), theField);
00094
00095 } else if (theVisitor.neut() != 0) {
00096
00097
00098 Hep3Vector theMom = theTraj.direction(fltlen);
00099 theMom.setMag(theTraj.parameters()->parameter()[NeutParams::_p]);
00100 return theMom;
00101
00102 } else {
00103
00104
00105 return Hep3Vector(999, 999, 999);
00106
00107 }
00108 }
00109
00110
00111 BesVectorErr
00112 TrkMomCalculator::errMom(const TrkSimpTraj& theTraj, const BField&
00113 theField, double fltlen) {
00114
00115
00116 TrkMomVisitor theVisitor(theTraj);
00117
00118 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) {
00119
00120
00121
00122 return calcCurvErrMom(theTraj, theField, fltlen);
00123
00124 } else if (theVisitor.neut() != 0) {
00125
00126
00127 return calcNeutErrMom(theTraj, theField, fltlen);
00128
00129 } else {
00130
00131
00132
00133 BesError theErr(3);
00134 return BesVectorErr(Hep3Vector(999, 999, 999), theErr);
00135 }
00136 }
00137
00138
00139 int
00140 TrkMomCalculator::charge(const TrkSimpTraj& theTraj, const BField&
00141 theField, double fltlen) {
00142
00143
00144 TrkMomVisitor theVisitor(theTraj);
00145
00146 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) {
00147
00148
00149
00150
00151 bool plus = false;
00152 HelixTraj::ParIndex parInd = HelixTraj::omegaIndex;
00153 int index = parInd;
00154 plus = (theField.bFieldNominal() < 0.0 &&
00155 theTraj.parameters()->parameter()[index] > 0.0)
00156 ||
00157 (theField.bFieldNominal() > 0.0 &&
00158 theTraj.parameters()->parameter()[index] <
00159 0.0);
00160 return ( plus ? 1 : -1 );
00161
00162
00163
00164
00165
00166 } else if (theVisitor.neut() != 0) {
00167
00168
00169 return 0;
00170
00171 } else {
00172
00173
00174 return 0;
00175 }
00176 }
00177
00178
00179 HepMatrix
00180 TrkMomCalculator::posmomCov(const TrkSimpTraj& theTraj,const BField& theField,
00181 double fltlen) {
00182
00183
00184 TrkMomVisitor theVisitor(theTraj);
00185
00186 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) {
00187
00188
00189
00190 return calcCurvPosmomCov(theTraj, theField, fltlen);
00191
00192 } else if (theVisitor.neut() != 0) {
00193
00194
00195 return calcNeutPosmomCov(theTraj, theField, fltlen);
00196
00197 } else {
00198
00199
00200
00201 return HepMatrix(3,3,0);
00202 }
00203 }
00204
00205
00206 void
00207 TrkMomCalculator::getAllCovs(const TrkSimpTraj& theTraj,
00208 const BField& theField,
00209 double fltlen,
00210 HepSymMatrix& xxCov,
00211 HepSymMatrix& ppCov,
00212 HepMatrix& xpCov) {
00213
00214
00215 TrkMomVisitor theVisitor(theTraj);
00216
00217 if (theVisitor.helix() != 0) {
00218
00219
00220 calcCurvAllCovs(theTraj,theField,fltlen,xxCov,ppCov,xpCov);
00221
00222 } else if (theVisitor.circle() != 0){
00223
00224
00225
00226 calcCurvAllCovsOLD(theTraj,theField,fltlen,xxCov,ppCov,xpCov);
00227
00228 } else if (theVisitor.neut() != 0) {
00229
00230
00231 calcNeutAllCovs(theTraj,theField,fltlen,xxCov,ppCov,xpCov);
00232
00233 } else {
00234
00235
00236
00237 int i,j;
00238 for(i=0;i<3;i++)
00239 for(j=i;j<3;j++)
00240 {
00241 xxCov[i][j]=0;
00242 ppCov[i][j]=0;
00243 xpCov[i][j]=0;
00244 xpCov[j][i]=0;
00245 }
00246 }
00247
00248 }
00249
00250
00251 void
00252 TrkMomCalculator::getAllWeights(const TrkSimpTraj& theTraj,
00253 const BField& theField,
00254 double fltlen,
00255 HepVector& pos,
00256 HepVector& mom,
00257 HepSymMatrix& xxWeight,
00258 HepSymMatrix& ppWeight,
00259 HepMatrix& xpWeight) {
00260
00261
00262 TrkMomVisitor theVisitor(theTraj);
00263
00264 if (theVisitor.helix() != 0){
00265
00266
00267
00268 calcCurvAllWeights(theTraj,theField,fltlen,
00269 pos,mom,xxWeight,ppWeight,xpWeight);
00270
00271 } else if (theVisitor.circle() != 0) {
00272
00273 calcCurvAllWeightsOLD(theTraj,theField,fltlen,
00274 pos,mom,xxWeight,ppWeight,xpWeight);
00275
00276 } else if (theVisitor.neut() != 0) {
00277
00278
00279 calcNeutAllWeights(theTraj,theField,fltlen,
00280 pos,mom,xxWeight,ppWeight,xpWeight);
00281
00282 } else {
00283
00284
00285
00286 int i,j;
00287 for(i=0;i<3;i++)
00288 for(j=i;j<3;j++) {
00289 xxWeight[i][j]=0;
00290 ppWeight[i][j]=0;
00291 xpWeight[i][j]=0;
00292 xpWeight[j][i]=0;
00293 }
00294
00295 for(i=0;i<3;i++) {
00296 pos[i]=0;
00297 mom[i]=0;
00298 }
00299 }
00300 }
00301
00302
00303 double
00304 TrkMomCalculator::calcCurvPtMom(const Hep3Vector& direction,
00305 double curvature,
00306 const BField& theField) {
00307
00308 Hep3Vector theMomVec = calcCurvVecMom(direction, curvature,
00309 theField);
00310 return theMomVec.perp();
00311
00312 }
00313
00314
00315 Hep3Vector
00316 TrkMomCalculator::calcCurvVecMom(const Hep3Vector& direction,
00317 double curvature,
00318 const BField& theField) {
00319
00320 double sindip = direction.z();
00321 double arg = 1.0-sindip*sindip;
00322 if (arg < 0.0) arg = 0.0;
00323 double cosdip = sqrt(arg);
00324 double momMag =
00325 fabs( BField::cmTeslaToGeVc*theField.bFieldNominal()*
00326 cosdip/curvature );
00327 Hep3Vector momVec = direction;
00328 momVec.setMag(momMag);
00329
00330 return momVec;
00331
00332 }
00333
00334
00335 BesVectorErr
00336 TrkMomCalculator::calcCurvErrMom(const TrkSimpTraj& theTraj,
00337 const BField& theField,
00338 double fltlen) {
00339
00340
00341 DifPoint PosDif;
00342 DifVector DirDif;
00343 DifVector delDirDif;
00344 DifVector MomDif(0., 0., 0.);
00345
00346 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
00347 if (delDirDif.length() == 0.) {
00348 }
00349 else {
00350 DifNumber sindip=DirDif.z;
00351 DifNumber arg = 1.0-sindip*sindip;
00352 if (arg.number() < 0.0) {arg.setNumber(0.0);}
00353 DifNumber cosdip = sqrt(arg);
00354
00355 DifNumber momMag =
00356 BField::cmTeslaToGeVc*theField.bFieldNominal()*cosdip /
00357 delDirDif.length();
00358 momMag.absolute();
00359
00360 MomDif = DirDif * momMag;
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405 }
00406 BesError symErr(MomDif.errorMatrix(
00407 MomDif.x.indepPar()->covariance()));
00408 return BesVectorErr(Hep3Vector(MomDif.x.number(), MomDif.y.number(),
00409 MomDif.z.number()), symErr);
00410 }
00411
00412
00413 BesVectorErr
00414 TrkMomCalculator::calcNeutErrMom(const TrkSimpTraj& theTraj,
00415 const BField& theField,
00416 double fltlen) {
00417
00418
00419 DifPoint posDif;
00420 DifVector dirDif;
00421 DifVector delDirDif;
00422
00423 theTraj.getDFInfo(fltlen, posDif, dirDif, delDirDif);
00424
00425
00426 DifVector momDif = dirDif;
00427
00428
00429
00430 momDif *= theTraj.parameters()->difPar(NeutParams::_p + 1);
00431
00432
00433 BesError symErr(momDif.errorMatrix(
00434 momDif.x.indepPar()->covariance()));
00435
00436
00437 return BesVectorErr(Hep3Vector(momDif.x.number(), momDif.y.number(),
00438 momDif.z.number()), symErr);
00439 }
00440
00441
00442
00443
00444
00445
00446 int
00447 TrkMomCalculator::calcCurvCharge(const Hep3Vector& direction,
00448 double curvature,
00449 const BField& theField) {
00450
00451
00452 Hep3Vector momVec =
00453 calcCurvVecMom(direction, curvature, theField);
00454
00455 if (theField.bFieldNominal() > 0.) {
00456 return -nearestInt(momVec.mag() * curvature / theField.bFieldNominal());
00457 } else {
00458 return nearestInt(momVec.mag() * curvature / theField.bFieldNominal());
00459 }
00460 }
00461
00462
00463 int
00464 TrkMomCalculator::nearestInt(double floatPt) {
00465
00466
00467 if (floatPt > 0.) {
00468 return (int)(floatPt+0.5);
00469 } else {
00470 return (int)(floatPt-0.5);
00471 }
00472 }
00473
00474
00475
00476
00477
00478
00479
00480 HepMatrix
00481 TrkMomCalculator::calcCurvPosmomCov(const TrkSimpTraj& theTraj,
00482 const BField& theField,
00483 double fltlen) {
00484
00485
00486 DifPoint PosDif;
00487 DifVector DirDif;
00488 DifVector delDirDif;
00489 DifVector MomDif(0., 0., 0.);
00490
00491 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
00492 if (delDirDif.length() == 0.) {
00493 }
00494 else {
00495 DifNumber sindip=DirDif.z;
00496 DifNumber arg = 1.0-sindip*sindip;
00497 if (arg.number() < 0.0) {arg.setNumber(0.0);}
00498 DifNumber cosdip = sqrt(arg);
00499
00500
00501 DifNumber momMag =
00502 BField::cmTeslaToGeVc*theField.bFieldNominal()*cosdip /
00503 delDirDif.length();
00504 momMag.absolute();
00505
00506 MomDif = DirDif * momMag;
00507
00508 }
00509
00510
00511 HepMatrix xpCov(3,3);
00512 const HepSymMatrix& nnCov=MomDif.x.indepPar()->covariance();
00513 xpCov(1,1)=correlation(PosDif.x,MomDif.x,nnCov);
00514 xpCov(1,2)=correlation(PosDif.x,MomDif.y,nnCov);
00515 xpCov(1,3)=correlation(PosDif.x,MomDif.z,nnCov);
00516 xpCov(2,1)=correlation(PosDif.y,MomDif.x,nnCov);
00517 xpCov(2,2)=correlation(PosDif.y,MomDif.y,nnCov);
00518 xpCov(2,3)=correlation(PosDif.y,MomDif.z,nnCov);
00519 xpCov(3,1)=correlation(PosDif.z,MomDif.x,nnCov);
00520 xpCov(3,2)=correlation(PosDif.z,MomDif.y,nnCov);
00521 xpCov(3,3)=correlation(PosDif.z,MomDif.z,nnCov);
00522
00523 return xpCov;
00524 }
00525
00526
00527 HepMatrix
00528 TrkMomCalculator::calcNeutPosmomCov(const TrkSimpTraj& theTraj,
00529 const BField& theField,
00530 double fltlen) {
00531
00532
00533 DifPoint PosDif;
00534 DifVector DirDif;
00535 DifVector delDirDif;
00536
00537 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
00538
00539
00540 DifVector MomDif = DirDif;
00541
00542
00543
00544 MomDif *= theTraj.parameters()->difPar(NeutParams::_p + 1);
00545
00546
00547 HepMatrix xpCov(3,3);
00548 const HepSymMatrix& nnCov=MomDif.x.indepPar()->covariance();
00549 xpCov(1,1)=correlation(PosDif.x,MomDif.x,nnCov);
00550 xpCov(1,2)=correlation(PosDif.x,MomDif.y,nnCov);
00551 xpCov(1,3)=correlation(PosDif.x,MomDif.z,nnCov);
00552 xpCov(2,1)=correlation(PosDif.y,MomDif.x,nnCov);
00553 xpCov(2,2)=correlation(PosDif.y,MomDif.y,nnCov);
00554 xpCov(2,3)=correlation(PosDif.y,MomDif.z,nnCov);
00555 xpCov(3,1)=correlation(PosDif.z,MomDif.x,nnCov);
00556 xpCov(3,2)=correlation(PosDif.z,MomDif.y,nnCov);
00557 xpCov(3,3)=correlation(PosDif.z,MomDif.z,nnCov);
00558
00559 return HepMatrix(3,3,0);
00560 }
00561
00562 bool TrkMomCalculator::weightToCov(const HepSymMatrix& inXX,const HepSymMatrix& inPP,const HepMatrix& inXP,
00563 HepSymMatrix& outXX,HepSymMatrix& outPP,HepMatrix& outXP){
00564 assert(inXX.num_row()==outXX.num_row());
00565 assert(inPP.num_row()==outPP.num_row());
00566 assert(inXP.num_row()==outXP.num_row());
00567 assert(inXP.num_col()==outXP.num_col());
00568 assert(inXX.num_row()==inXP.num_row());
00569 assert(inPP.num_row()==inXP.num_col());
00570 int status;
00571 HepSymMatrix aInv = inXX.inverse(status);
00572 if(status)return false;
00573 HepSymMatrix beta = inPP-aInv.similarityT(inXP);
00574 outPP = beta.inverse(status);
00575 if(status)return false;
00576 outXP = -aInv*inXP*outPP;
00577 HepMatrix alpha(aInv-aInv*inXP*outXP.T());
00578 outXX.assign(alpha);
00579 return true;
00580 }
00581
00582
00583 void
00584 TrkMomCalculator::calcCurvAllCovs(const TrkSimpTraj& theTraj,
00585 const BField& theField,
00586 double fltlen,
00587 HepSymMatrix& xxCov,
00588 HepSymMatrix& ppCov,
00589 HepMatrix& xpCov) {
00590
00591
00592 const HepVector& v = theTraj.parameters()->parameter();
00593 const HepSymMatrix& m = theTraj.parameters()->covariance();
00594
00595
00596
00597
00598 double s = v[TrkExchangePar::ex_tanDip];
00599 double omega = v[TrkExchangePar::ex_omega];
00600 double d0 = v[TrkExchangePar::ex_d0];
00601
00602 double phi0 = v[TrkExchangePar::ex_phi0];
00603 double cosDip = 1/sqrt(1.0+s*s);
00604 double l = fltlen*cosDip ;
00605
00606
00607 double dlds = -fltlen*s*(cosDip*cosDip*cosDip) ;
00608 double phi = phi0 + omega*l;
00609 double sinphi0 = sin(phi0);
00610 double cosphi0 = cos(phi0);
00611 double sinphi = sin(phi);
00612 double cosphi = cos(phi);
00613 double pt = fabs( BField::cmTeslaToGeVc * theField.bFieldNominal() / omega );
00614 double r = 1.0/omega;
00615
00616
00617 double d_x_d0 = -sinphi0;
00618 double d_x_phi0 = r*cosphi - (r+d0)*cosphi0;
00619 double d_x_omega = -r*r*sinphi + r*r*sinphi0 + l*r*cosphi;
00620 double d_x_tanDip = cosphi*dlds;
00621 double d_y_d0 = cosphi0;
00622 double d_y_phi0 = r*sinphi - (r+d0)*sinphi0;
00623 double d_y_omega = r*r*cosphi - r*r*cosphi0 + l*r*sinphi;
00624 double d_y_tanDip = sinphi*dlds;
00625 double d_z_z0 = 1.0;
00626 double d_z_tanDip = l + dlds*s;
00627 double d_px_phi0 = -pt*sinphi;
00628 double d_px_omega = -pt*cosphi/omega - pt*l*sinphi;
00629 double d_px_tanDip = -pt*omega*sinphi*dlds;
00630 double d_py_phi0 = pt*cosphi;
00631 double d_py_omega = -pt*sinphi/omega + pt*l*cosphi;
00632 double d_py_tanDip = pt*omega*cosphi*dlds;
00633 double d_pz_omega = -pt*s/omega;
00634 double d_pz_tanDip = pt;
00635
00636
00637 double m_d0_d0 = m[TrkExchangePar::ex_d0][TrkExchangePar::ex_d0];
00638 double m_phi0_d0 = m[TrkExchangePar::ex_phi0][TrkExchangePar::ex_d0];
00639 double m_phi0_phi0 = m[TrkExchangePar::ex_phi0][TrkExchangePar::ex_phi0];
00640 double m_omega_d0 = m[TrkExchangePar::ex_omega][TrkExchangePar::ex_d0];
00641 double m_omega_phi0 = m[TrkExchangePar::ex_omega][TrkExchangePar::ex_phi0];
00642 double m_omega_omega = m[TrkExchangePar::ex_omega][TrkExchangePar::ex_omega];
00643 double m_z0_d0 = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_d0];
00644 double m_z0_phi0 = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_phi0];
00645 double m_z0_omega = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_omega];
00646 double m_z0_z0 = m[TrkExchangePar::ex_z0][TrkExchangePar::ex_z0];
00647 double m_tanDip_d0 = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_d0];
00648 double m_tanDip_phi0 = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_phi0];
00649 double m_tanDip_omega = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_omega];
00650 double m_tanDip_z0 = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_z0];
00651 double m_tanDip_tanDip = m[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_tanDip];
00652
00653
00654 xxCov(1,1) =
00655 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 ) +
00656 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 ) +
00657 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 ) +
00658 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 ) ;
00659 xxCov(2,1) =
00660 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 ) +
00661 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 ) +
00662 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 ) +
00663 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 ) ;
00664 xxCov(2,2) =
00665 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 ) +
00666 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 ) +
00667 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 ) +
00668 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 ) ;
00669 xxCov(3,1) =
00670 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 ) +
00671 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 ) ;
00672 xxCov(3,2) =
00673 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 ) +
00674 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 ) ;
00675 xxCov(3,3) =
00676 d_z_z0* ( d_z_z0*m_z0_z0 + d_z_tanDip*m_tanDip_z0 ) +
00677 d_z_tanDip* ( d_z_z0*m_tanDip_z0 + d_z_tanDip*m_tanDip_tanDip ) ;
00678
00679 xpCov(1,1) =
00680 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 ) +
00681 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 ) +
00682 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 ) ;
00683 xpCov(2,1) =
00684 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 ) +
00685 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 ) +
00686 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 ) ;
00687 xpCov(3,1) =
00688 d_px_phi0* ( d_z_z0*m_z0_phi0 + d_z_tanDip*m_tanDip_phi0 ) +
00689 d_px_omega* ( d_z_z0*m_z0_omega + d_z_tanDip*m_tanDip_omega ) +
00690 d_px_tanDip* ( d_z_z0*m_tanDip_z0 + d_z_tanDip*m_tanDip_tanDip ) ;
00691 xpCov(1,2) =
00692 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 ) +
00693 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 ) +
00694 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 ) ;
00695 xpCov(2,2) =
00696 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 ) +
00697 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 ) +
00698 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 ) ;
00699 xpCov(3,2) =
00700 d_py_phi0* ( d_z_z0*m_z0_phi0 + d_z_tanDip*m_tanDip_phi0 ) +
00701 d_py_omega* ( d_z_z0*m_z0_omega + d_z_tanDip*m_tanDip_omega ) +
00702 d_py_tanDip* ( d_z_z0*m_tanDip_z0 + d_z_tanDip*m_tanDip_tanDip ) ;
00703 xpCov(1,3) =
00704 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 ) +
00705 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 ) ;
00706 xpCov(2,3) =
00707 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 ) +
00708 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 ) ;
00709 xpCov(3,3) =
00710 d_pz_omega* ( d_z_z0*m_z0_omega + d_z_tanDip*m_tanDip_omega ) +
00711 d_pz_tanDip* ( d_z_z0*m_tanDip_z0 + d_z_tanDip*m_tanDip_tanDip ) ;
00712
00713 ppCov(1,1) =
00714 d_px_phi0* ( d_px_phi0*m_phi0_phi0 + d_px_omega*m_omega_phi0 + d_px_tanDip*m_tanDip_phi0 ) +
00715 d_px_omega* ( d_px_phi0*m_omega_phi0 + d_px_omega*m_omega_omega + d_px_tanDip*m_tanDip_omega ) +
00716 d_px_tanDip* ( d_px_phi0*m_tanDip_phi0 + d_px_omega*m_tanDip_omega + d_px_tanDip*m_tanDip_tanDip ) ;
00717 ppCov(2,1) =
00718 d_py_phi0* ( d_px_phi0*m_phi0_phi0 + d_px_omega*m_omega_phi0 + d_px_tanDip*m_tanDip_phi0 ) +
00719 d_py_omega* ( d_px_phi0*m_omega_phi0 + d_px_omega*m_omega_omega + d_px_tanDip*m_tanDip_omega ) +
00720 d_py_tanDip* ( d_px_phi0*m_tanDip_phi0 + d_px_omega*m_tanDip_omega + d_px_tanDip*m_tanDip_tanDip ) ;
00721 ppCov(2,2) =
00722 d_py_phi0* ( d_py_phi0*m_phi0_phi0 + d_py_omega*m_omega_phi0 + d_py_tanDip*m_tanDip_phi0 ) +
00723 d_py_omega* ( d_py_phi0*m_omega_phi0 + d_py_omega*m_omega_omega + d_py_tanDip*m_tanDip_omega ) +
00724 d_py_tanDip* ( d_py_phi0*m_tanDip_phi0 + d_py_omega*m_tanDip_omega + d_py_tanDip*m_tanDip_tanDip ) ;
00725 ppCov(3,1) =
00726 d_pz_omega* ( d_px_phi0*m_omega_phi0 + d_px_omega*m_omega_omega + d_px_tanDip*m_tanDip_omega ) +
00727 d_pz_tanDip* ( d_px_phi0*m_tanDip_phi0 + d_px_omega*m_tanDip_omega + d_px_tanDip*m_tanDip_tanDip ) ;
00728 ppCov(3,2) =
00729 d_pz_omega* ( d_py_phi0*m_omega_phi0 + d_py_omega*m_omega_omega + d_py_tanDip*m_tanDip_omega ) +
00730 d_pz_tanDip* ( d_py_phi0*m_tanDip_phi0 + d_py_omega*m_tanDip_omega + d_py_tanDip*m_tanDip_tanDip ) ;
00731 ppCov(3,3) =
00732 d_pz_omega* ( d_pz_omega*m_omega_omega + d_pz_tanDip*m_tanDip_omega ) +
00733 d_pz_tanDip* ( d_pz_omega*m_tanDip_omega + d_pz_tanDip*m_tanDip_tanDip ) ;
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766 }
00767
00768
00769 void
00770 TrkMomCalculator::calcCurvAllCovsOLD(const TrkSimpTraj& theTraj,
00771 const BField& theField,
00772 double fltlen,
00773 HepSymMatrix& xxCov,
00774 HepSymMatrix& ppCov,
00775 HepMatrix& xpCov) {
00776
00777
00778 DifPoint PosDif;
00779 DifVector DirDif;
00780 DifVector delDirDif;
00781 DifVector MomDif(0., 0., 0.);
00782
00783 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
00784 if (delDirDif.length() != 0) {
00785 DifNumber sindip=DirDif.z;
00786 DifNumber arg = 1.0-sindip*sindip;
00787 if (arg.number() < 0.0) {arg.setNumber(0.0);}
00788 DifNumber cosdip = sqrt(arg);
00789
00790 DifNumber momMag =
00791 BField::cmTeslaToGeVc*theField.bFieldNominal()*cosdip /
00792 delDirDif.length();
00793 momMag.absolute();
00794
00795 MomDif = DirDif * momMag;
00796
00797 }
00798
00799
00800
00801 xxCov.assign(PosDif.errorMatrix(PosDif.x.indepPar()->covariance()));
00802
00803
00804 ppCov.assign(MomDif.errorMatrix(MomDif.x.indepPar()->covariance()));
00805
00806
00807 const HepSymMatrix& nnCov=MomDif.x.indepPar()->covariance();
00808 xpCov(1,1)=correlation(PosDif.x,MomDif.x,nnCov);
00809 xpCov(1,2)=correlation(PosDif.x,MomDif.y,nnCov);
00810 xpCov(1,3)=correlation(PosDif.x,MomDif.z,nnCov);
00811 xpCov(2,1)=correlation(PosDif.y,MomDif.x,nnCov);
00812 xpCov(2,2)=correlation(PosDif.y,MomDif.y,nnCov);
00813 xpCov(2,3)=correlation(PosDif.y,MomDif.z,nnCov);
00814 xpCov(3,1)=correlation(PosDif.z,MomDif.x,nnCov);
00815 xpCov(3,2)=correlation(PosDif.z,MomDif.y,nnCov);
00816 xpCov(3,3)=correlation(PosDif.z,MomDif.z,nnCov);
00817
00818 }
00819
00820
00821 void
00822 TrkMomCalculator::calcNeutAllCovs(const TrkSimpTraj& theTraj,
00823 const BField& theField,
00824 double fltlen,
00825 HepSymMatrix& xxCov,
00826 HepSymMatrix& ppCov,
00827 HepMatrix& xpCov) {
00828
00829 HepVector p0(3),x0(3);
00830 calcNeutAllCovs(theTraj,theField,fltlen,x0,p0,xxCov,ppCov,xpCov);
00831
00832 }
00833
00834
00835 void
00836 TrkMomCalculator::calcNeutAllCovs(const TrkSimpTraj& theTraj,
00837 const BField& theField,
00838 double fltlen,
00839 HepVector& x0,HepVector& p0,
00840 HepSymMatrix& xxCov,
00841 HepSymMatrix& ppCov,
00842 HepMatrix& xpCov) {
00843
00844
00845 DifPoint PosDif;
00846 DifVector DirDif;
00847 DifVector delDirDif;
00848
00849 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
00850
00851
00852 DifVector MomDif = DirDif;
00853
00854
00855
00856 MomDif *= theTraj.parameters()->difPar(NeutParams::_p + 1);
00857
00858
00859 p0[0] = MomDif.x.number();
00860 p0[1] = MomDif.y.number();
00861 p0[2] = MomDif.z.number();
00862 x0[0] = PosDif.x.number();
00863 x0[1] = PosDif.y.number();
00864 x0[2] = PosDif.z.number();
00865
00866
00867 xxCov.assign(PosDif.errorMatrix(PosDif.x.indepPar()->covariance()));
00868
00869
00870 ppCov.assign(MomDif.errorMatrix(MomDif.x.indepPar()->covariance()));
00871
00872
00873 const HepSymMatrix& nnCov=MomDif.x.indepPar()->covariance();
00874 xpCov(1,1)=correlation(PosDif.x,MomDif.x,nnCov);
00875 xpCov(1,2)=correlation(PosDif.x,MomDif.y,nnCov);
00876 xpCov(1,3)=correlation(PosDif.x,MomDif.z,nnCov);
00877 xpCov(2,1)=correlation(PosDif.y,MomDif.x,nnCov);
00878 xpCov(2,2)=correlation(PosDif.y,MomDif.y,nnCov);
00879 xpCov(2,3)=correlation(PosDif.y,MomDif.z,nnCov);
00880 xpCov(3,1)=correlation(PosDif.z,MomDif.x,nnCov);
00881 xpCov(3,2)=correlation(PosDif.z,MomDif.y,nnCov);
00882 xpCov(3,3)=correlation(PosDif.z,MomDif.z,nnCov);
00883
00884 }
00885
00886
00887 void
00888 TrkMomCalculator::calcCurvAllWeights(const TrkSimpTraj& theTraj,
00889 const BField& theField,
00890 double fltlen,
00891 HepVector& pos,
00892 HepVector& mom,
00893 HepSymMatrix& xxWeight,
00894 HepSymMatrix& ppWeight,
00895 HepMatrix& xpWeight) {
00896
00897 const HepVector& v = theTraj.parameters()->parameter();
00898 const HepSymMatrix& w = theTraj.parameters()->weightMatrix();
00899
00900
00901
00902
00903 double s = v[TrkExchangePar::ex_tanDip];
00904 double omega = v[TrkExchangePar::ex_omega];
00905 double d0 = v[TrkExchangePar::ex_d0];
00906 double z0 = v[TrkExchangePar::ex_z0];
00907 double phi0 = v[TrkExchangePar::ex_phi0];
00908 double l = fltlen / sqrt(1.0+s*s) ;
00909
00910 double phi = phi0 + omega*l;
00911 double sinphi0 = sin(phi0);
00912 double cosphi0 = cos(phi0);
00913 double sinphi = sin(phi);
00914 double cosphi = cos(phi);
00915 double C = BField::cmTeslaToGeVc * theField.bFieldNominal();
00916 double q(1.0);
00917 -omega>0 ? q=1.0: q=-1.0;
00918 double qC = q*C;
00919 double pt = fabs( -qC / omega );
00920 double r = 1.0/omega;
00921
00922
00923 pos(1) = r*sinphi - (r + d0)*sinphi0;
00924 pos(2) = -r*cosphi + (r + d0)*cosphi0;
00925 pos(3) = z0 + l*s;
00926 mom(1) = pt*cosphi;
00927 mom(2) = pt*sinphi;
00928 mom(3) = pt*s;
00929
00930
00931
00932
00933 if ( (1+d0*omega)==0.0 ){
00934 calcCurvAllWeightsOLD(theTraj,theField,fltlen,
00935 pos,mom,xxWeight,ppWeight,xpWeight);
00936 return;
00937 }
00938
00939 double dinv_x_d0 = -sinphi0;
00940 double dinv_x_phi0 = -omega * cosphi0 / (1 + d0 * omega);
00941 double dinv_x_z0 = -s * cosphi0 / (1 + d0 * omega);
00942
00943 double dinv_y_d0 = cosphi0;
00944 double dinv_y_phi0 = -omega * sinphi0 / (1 + d0 * omega);
00945 double dinv_y_z0 = -s * sinphi0 / (1 + d0 * omega);
00946
00947 double dinv_z_z0 = 1;
00948
00949 double dinv_px_d0 = (cosphi - cosphi0) / qC;
00950 double dinv_px_phi0 = omega * sinphi0 / (1 + d0 * omega) / qC;
00951 double dinv_px_omega = omega * omega * cosphi / qC;
00952 double dinv_px_z0 = -s * ( sinphi*(1 + d0 * omega) - sinphi0 ) / (qC * (1 + d0 * omega));
00953 double dinv_px_tanDip = omega * cosphi * s / qC;
00954
00955 double dinv_py_d0 = (sinphi - sinphi0) / qC;
00956 double dinv_py_phi0 = -omega * cosphi0 / (1 + d0 * omega) / qC;
00957 double dinv_py_omega = omega * omega * sinphi / qC;
00958 double dinv_py_z0 = s * ( cosphi*(1 + d0 * omega) - cosphi0) / (qC * (1 + d0 * omega));
00959 double dinv_py_tanDip = omega * sinphi * s / qC;
00960
00961 double dinv_pz_z0 = l * omega / qC;
00962 double dinv_pz_tanDip = -omega / qC;
00963
00964
00965 double w_d0_d0 = w[TrkExchangePar::ex_d0][TrkExchangePar::ex_d0];
00966 double w_phi0_d0 = w[TrkExchangePar::ex_phi0][TrkExchangePar::ex_d0];
00967 double w_phi0_phi0 = w[TrkExchangePar::ex_phi0][TrkExchangePar::ex_phi0];
00968 double w_omega_d0 = w[TrkExchangePar::ex_omega][TrkExchangePar::ex_d0];
00969 double w_omega_phi0 = w[TrkExchangePar::ex_omega][TrkExchangePar::ex_phi0];
00970 double w_omega_omega = w[TrkExchangePar::ex_omega][TrkExchangePar::ex_omega];
00971 double w_z0_d0 = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_d0];
00972 double w_z0_phi0 = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_phi0];
00973 double w_z0_omega = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_omega];
00974 double w_z0_z0 = w[TrkExchangePar::ex_z0][TrkExchangePar::ex_z0];
00975 double w_tanDip_d0 = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_d0];
00976 double w_tanDip_phi0 = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_phi0];
00977 double w_tanDip_omega = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_omega];
00978 double w_tanDip_z0 = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_z0];
00979 double w_tanDip_tanDip = w[TrkExchangePar::ex_tanDip][TrkExchangePar::ex_tanDip];
00980
00981
00982 xxWeight(1,1) =
00983 dinv_x_d0* ( dinv_x_d0*w_d0_d0 +dinv_x_phi0*w_phi0_d0 +dinv_x_z0*w_z0_d0 ) +
00984 dinv_x_phi0* ( dinv_x_d0*w_phi0_d0 +dinv_x_phi0*w_phi0_phi0 +dinv_x_z0*w_z0_phi0 ) +
00985 dinv_x_z0* ( dinv_x_d0*w_z0_d0 +dinv_x_phi0*w_z0_phi0 +dinv_x_z0*w_z0_z0 ) ;
00986 xxWeight(2,1) =
00987 dinv_y_d0* ( dinv_x_d0*w_d0_d0 +dinv_x_phi0*w_phi0_d0 +dinv_x_z0*w_z0_d0 ) +
00988 dinv_y_phi0* ( dinv_x_d0*w_phi0_d0 +dinv_x_phi0*w_phi0_phi0 +dinv_x_z0*w_z0_phi0 ) +
00989 dinv_y_z0* ( dinv_x_d0*w_z0_d0 +dinv_x_phi0*w_z0_phi0 +dinv_x_z0*w_z0_z0 ) ;
00990 xxWeight(2,2) =
00991 dinv_y_d0* ( dinv_y_d0*w_d0_d0 +dinv_y_phi0*w_phi0_d0 +dinv_y_z0*w_z0_d0 ) +
00992 dinv_y_phi0* ( dinv_y_d0*w_phi0_d0 +dinv_y_phi0*w_phi0_phi0 +dinv_y_z0*w_z0_phi0 ) +
00993 dinv_y_z0* ( dinv_y_d0*w_z0_d0 +dinv_y_phi0*w_z0_phi0 +dinv_y_z0*w_z0_z0 ) ;
00994 xxWeight(3,1) =
00995 dinv_z_z0* ( dinv_x_d0*w_z0_d0 +dinv_x_phi0*w_z0_phi0 +dinv_x_z0*w_z0_z0 ) ;
00996 xxWeight(3,2) =
00997 dinv_z_z0* ( dinv_y_d0*w_z0_d0 +dinv_y_phi0*w_z0_phi0 +dinv_y_z0*w_z0_z0 ) ;
00998 xxWeight(3,3) =
00999 dinv_z_z0* ( dinv_z_z0*w_z0_z0 ) ;
01000
01001 xpWeight(1,1) =
01002 dinv_px_d0* ( dinv_x_d0*w_d0_d0 +dinv_x_phi0*w_phi0_d0 +dinv_x_z0*w_z0_d0 ) +
01003 dinv_px_phi0* ( dinv_x_d0*w_phi0_d0 +dinv_x_phi0*w_phi0_phi0 +dinv_x_z0*w_z0_phi0 ) +
01004 dinv_px_omega* ( dinv_x_d0*w_omega_d0 +dinv_x_phi0*w_omega_phi0 +dinv_x_z0*w_z0_omega ) +
01005 dinv_px_z0* ( dinv_x_d0*w_z0_d0 +dinv_x_phi0*w_z0_phi0 +dinv_x_z0*w_z0_z0 ) +
01006 dinv_px_tanDip* ( dinv_x_d0*w_tanDip_d0 +dinv_x_phi0*w_tanDip_phi0 +dinv_x_z0*w_tanDip_z0 ) ;
01007 xpWeight(2,1) =
01008 dinv_px_d0* ( dinv_y_d0*w_d0_d0 +dinv_y_phi0*w_phi0_d0 +dinv_y_z0*w_z0_d0 ) +
01009 dinv_px_phi0* ( dinv_y_d0*w_phi0_d0 +dinv_y_phi0*w_phi0_phi0 +dinv_y_z0*w_z0_phi0 ) +
01010 dinv_px_omega* ( dinv_y_d0*w_omega_d0 +dinv_y_phi0*w_omega_phi0 +dinv_y_z0*w_z0_omega ) +
01011 dinv_px_z0* ( dinv_y_d0*w_z0_d0 +dinv_y_phi0*w_z0_phi0 +dinv_y_z0*w_z0_z0 ) +
01012 dinv_px_tanDip* ( dinv_y_d0*w_tanDip_d0 +dinv_y_phi0*w_tanDip_phi0 +dinv_y_z0*w_tanDip_z0 ) ;
01013 xpWeight(3,1) =
01014 dinv_px_d0* ( dinv_z_z0*w_z0_d0 ) +
01015 dinv_px_phi0* ( dinv_z_z0*w_z0_phi0 ) +
01016 dinv_px_omega* ( dinv_z_z0*w_z0_omega ) +
01017 dinv_px_z0* ( dinv_z_z0*w_z0_z0 ) +
01018 dinv_px_tanDip* ( dinv_z_z0*w_tanDip_z0 ) ;
01019 xpWeight(1,2) =
01020 dinv_py_d0* ( dinv_x_d0*w_d0_d0 +dinv_x_phi0*w_phi0_d0 +dinv_x_z0*w_z0_d0 ) +
01021 dinv_py_phi0* ( dinv_x_d0*w_phi0_d0 +dinv_x_phi0*w_phi0_phi0 +dinv_x_z0*w_z0_phi0 ) +
01022 dinv_py_omega* ( dinv_x_d0*w_omega_d0 +dinv_x_phi0*w_omega_phi0 +dinv_x_z0*w_z0_omega ) +
01023 dinv_py_z0* ( dinv_x_d0*w_z0_d0 +dinv_x_phi0*w_z0_phi0 +dinv_x_z0*w_z0_z0 ) +
01024 dinv_py_tanDip* ( dinv_x_d0*w_tanDip_d0 +dinv_x_phi0*w_tanDip_phi0 +dinv_x_z0*w_tanDip_z0 ) ;
01025 xpWeight(2,2) =
01026 dinv_py_d0* ( dinv_y_d0*w_d0_d0 +dinv_y_phi0*w_phi0_d0 +dinv_y_z0*w_z0_d0 ) +
01027 dinv_py_phi0* ( dinv_y_d0*w_phi0_d0 +dinv_y_phi0*w_phi0_phi0 +dinv_y_z0*w_z0_phi0 ) +
01028 dinv_py_omega* ( dinv_y_d0*w_omega_d0 +dinv_y_phi0*w_omega_phi0 +dinv_y_z0*w_z0_omega ) +
01029 dinv_py_z0* ( dinv_y_d0*w_z0_d0 +dinv_y_phi0*w_z0_phi0 +dinv_y_z0*w_z0_z0 ) +
01030 dinv_py_tanDip* ( dinv_y_d0*w_tanDip_d0 +dinv_y_phi0*w_tanDip_phi0 +dinv_y_z0*w_tanDip_z0 ) ;
01031 xpWeight(3,2) =
01032 dinv_py_d0* ( dinv_z_z0*w_z0_d0 ) +
01033 dinv_py_phi0* ( dinv_z_z0*w_z0_phi0 ) +
01034 dinv_py_omega* ( dinv_z_z0*w_z0_omega ) +
01035 dinv_py_z0* ( dinv_z_z0*w_z0_z0 ) +
01036 dinv_py_tanDip* ( dinv_z_z0*w_tanDip_z0 ) ;
01037 xpWeight(1,3) =
01038 dinv_pz_z0* ( dinv_x_d0*w_z0_d0 +dinv_x_phi0*w_z0_phi0 +dinv_x_z0*w_z0_z0 ) +
01039 dinv_pz_tanDip* ( dinv_x_d0*w_tanDip_d0 +dinv_x_phi0*w_tanDip_phi0 +dinv_x_z0*w_tanDip_z0 ) ;
01040 xpWeight(2,3) =
01041 dinv_pz_z0* ( dinv_y_d0*w_z0_d0 +dinv_y_phi0*w_z0_phi0 +dinv_y_z0*w_z0_z0 ) +
01042 dinv_pz_tanDip* ( dinv_y_d0*w_tanDip_d0 +dinv_y_phi0*w_tanDip_phi0 +dinv_y_z0*w_tanDip_z0 ) ;
01043 xpWeight(3,3) =
01044 dinv_pz_z0* ( dinv_z_z0*w_z0_z0 ) +
01045 dinv_pz_tanDip* ( dinv_z_z0*w_tanDip_z0 ) ;
01046
01047
01048 ppWeight(1,1) =
01049 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 ) +
01050 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 ) +
01051 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 ) +
01052 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 ) +
01053 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 ) ;
01054 ppWeight(2,1) =
01055 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 ) +
01056 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 ) +
01057 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 ) +
01058 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 ) +
01059 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 ) ;
01060 ppWeight(2,2) =
01061 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 ) +
01062 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 ) +
01063 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 ) +
01064 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 ) +
01065 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 ) ;
01066 ppWeight(3,1) =
01067 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 ) +
01068 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 ) ;
01069 ppWeight(3,2) =
01070 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 ) +
01071 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 ) ;
01072 ppWeight(3,3) =
01073 dinv_pz_z0* ( dinv_pz_z0*w_z0_z0 +dinv_pz_tanDip*w_tanDip_z0 ) +
01074 dinv_pz_tanDip* ( dinv_pz_z0*w_tanDip_z0 +dinv_pz_tanDip*w_tanDip_tanDip ) ;
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106 }
01107
01108
01109 void
01110 TrkMomCalculator::calcCurvAllWeightsOLD(const TrkSimpTraj& theTraj,
01111 const BField& theField,
01112 double fltlen,
01113 HepVector& pos,
01114 HepVector& mom,
01115 HepSymMatrix& xxWeight,
01116 HepSymMatrix& ppWeight,
01117 HepMatrix& xpWeight) {
01118
01119
01120 DifPoint PosDif;
01121 DifVector DirDif;
01122 DifVector delDirDif;
01123 DifNumber momMag;
01124 DifVector MomDif;
01125
01126 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
01127 if (delDirDif.length() != 0) {
01128 DifNumber sindip=DirDif.z;
01129 DifNumber arg = 1.0-sindip*sindip;
01130 if (arg.number() < 0.0) {arg.setNumber(0.0);}
01131 DifNumber cosdip = sqrt(arg);
01132
01133 momMag =
01134 BField::cmTeslaToGeVc*theField.bFieldNominal()*cosdip /
01135 delDirDif.length();
01136 momMag.absolute();
01137
01138 MomDif = DirDif * momMag;
01139
01140 }
01141
01142
01143
01144
01145
01146
01147
01148
01149 HepMatrix Jx_n(PosDif.jacobian());
01150 HepMatrix Jp_n(MomDif.jacobian());
01151
01152 int i,j;
01153 HepMatrix Jxp_ns(6,6);
01154
01155 for(i=0;i<3;i++)
01156 for(j=0;j<5;j++)
01157 {
01158 Jxp_ns[i ][j]=Jx_n[i][j];
01159 Jxp_ns[i+3][j]=Jp_n[i][j];
01160 }
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171 Jxp_ns[0][5]=DirDif.x.number();
01172 Jxp_ns[1][5]=DirDif.y.number();
01173 Jxp_ns[2][5]=DirDif.z.number();
01174
01175 Jxp_ns[3][5]=momMag.number()*delDirDif.x.number();
01176 Jxp_ns[4][5]=momMag.number()*delDirDif.y.number();
01177 Jxp_ns[5][5]=momMag.number()*delDirDif.z.number();
01178
01179
01180
01181
01182
01183
01184
01185
01186 int invStatus;
01187
01188 Jxp_ns.invert(invStatus);
01189
01190 HepMatrix Jn_x(5,3);
01191 HepMatrix Jn_p(5,3);
01192
01193 for(i=0;i<5;i++)
01194 for(j=0;j<3;j++)
01195 {
01196 Jn_x[i][j]=Jxp_ns[i][j ];
01197 Jn_p[i][j]=Jxp_ns[i][j+3];
01198 }
01199
01200 HepSymMatrix Wnn(PosDif.x.indepPar()->covariance());
01201
01202 Wnn.invert(invStatus);
01203
01204
01205
01206
01207 xxWeight = Wnn.similarityT(Jn_x);
01208 ppWeight = Wnn.similarityT(Jn_p);
01209 xpWeight = Jn_x.T()*Wnn*Jn_p;
01210
01211 pos[0]=PosDif.x.number();
01212 pos[1]=PosDif.y.number();
01213 pos[2]=PosDif.z.number();
01214
01215 mom[0]=MomDif.x.number();
01216 mom[1]=MomDif.y.number();
01217 mom[2]=MomDif.z.number();
01218
01219 }
01220
01221
01222 void
01223 TrkMomCalculator::calcNeutAllWeights(const TrkSimpTraj& theTraj,
01224 const BField& theField,
01225 double fltlen,
01226 HepVector& pos,
01227 HepVector& mom,
01228 HepSymMatrix& xxWeight,
01229 HepSymMatrix& ppWeight,
01230 HepMatrix& xpWeight) {
01231
01232 DifPoint PosDif;
01233 DifVector DirDif;
01234 DifVector delDirDif;
01235 DifNumber momMag;
01236
01237 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif);
01238
01239
01240 DifVector MomDif = DirDif;
01241
01242 MomDif *= theTraj.parameters()->difPar(NeutParams::_p + 1);
01243
01244 HepMatrix Jx_n(PosDif.jacobian());
01245 HepMatrix Jp_n(MomDif.jacobian());
01246
01247 int i,j;
01248 HepMatrix Jxp_ns(6,6);
01249
01250 for(i=0;i<3;i++)
01251 for(j=0;j<6;j++)
01252 {
01253 Jxp_ns[i ][j]=Jx_n[i][j];
01254 Jxp_ns[i+3][j]=Jp_n[i][j];
01255 }
01256 int invStatus;
01257
01258 Jxp_ns.invert(invStatus);
01259
01260 HepMatrix Jn_x(5,3);
01261 HepMatrix Jn_p(5,3);
01262
01263 for(i=0;i<5;i++)
01264 for(j=0;j<3;j++)
01265 {
01266 Jn_x[i][j]=Jxp_ns[i][j ];
01267 Jn_p[i][j]=Jxp_ns[i][j+3];
01268 }
01269
01270 HepSymMatrix Wnn(PosDif.x.indepPar()->covariance().sub(1,5));
01271
01272 Wnn.invert(invStatus);
01273 xxWeight = Wnn.similarityT(Jn_x);
01274 ppWeight = Wnn.similarityT(Jn_p);
01275 xpWeight = Jn_x.T()*Wnn*Jn_p;
01276
01277 pos[0]=PosDif.x.number();
01278 pos[1]=PosDif.y.number();
01279 pos[2]=PosDif.z.number();
01280
01281 mom[0]=MomDif.x.number();
01282 mom[1]=MomDif.y.number();
01283 mom[2]=MomDif.z.number();
01284 }