#include <TrkMomCalculator.h>
Public Member Functions | |
TrkMomCalculator () | |
virtual | ~TrkMomCalculator () |
Static Public Member Functions | |
static double | ptMom (const TrkSimpTraj &, const BField &, double fltlen) |
static Hep3Vector | vecMom (const TrkSimpTraj &, const BField &, double fltlen) |
static BesVectorErr | errMom (const TrkSimpTraj &, const BField &, double fltlen) |
static int | charge (const TrkSimpTraj &, const BField &, double fltlen) |
static HepMatrix | posmomCov (const TrkSimpTraj &, const BField &, double fltlen) |
static void | getAllCovs (const TrkSimpTraj &, const BField &, double fltlen, HepSymMatrix &xxCov, HepSymMatrix &ppCov, HepMatrix &xpCov) |
static bool | weightToCov (const HepSymMatrix &inXX, const HepSymMatrix &inPP, const HepMatrix &inXP, HepSymMatrix &outXX, HepSymMatrix &outPP, HepMatrix &outXP) |
static void | getAllWeights (const TrkSimpTraj &, const BField &, double fltlen, HepVector &pos, HepVector &mom, HepSymMatrix &xxWeight, HepSymMatrix &ppWeight, HepMatrix &xpWeight) |
Static Private Member Functions | |
static double | calcCurvPtMom (const Hep3Vector &, double curvature, const BField &) |
static Hep3Vector | calcCurvVecMom (const Hep3Vector &, double curvature, const BField &) |
static BesVectorErr | calcCurvErrMom (const TrkSimpTraj &, const BField &, double flt) |
static BesVectorErr | calcNeutErrMom (const TrkSimpTraj &, const BField &, double flt) |
static int | calcCurvCharge (const Hep3Vector &, double curvature, const BField &) |
static HepMatrix | calcCurvPosmomCov (const TrkSimpTraj &, const BField &, double fltlen) |
static HepMatrix | calcNeutPosmomCov (const TrkSimpTraj &, const BField &, double fltlen) |
static void | calcCurvAllCovs (const TrkSimpTraj &, const BField &, double fltlen, HepSymMatrix &xxCov, HepSymMatrix &ppCov, HepMatrix &xpCov) |
static void | calcCurvAllCovsOLD (const TrkSimpTraj &, const BField &, double fltlen, HepSymMatrix &xxCov, HepSymMatrix &ppCov, HepMatrix &xpCov) |
static void | calcNeutAllCovs (const TrkSimpTraj &, const BField &, double fltlen, HepSymMatrix &xxCov, HepSymMatrix &ppCov, HepMatrix &xpCov) |
static void | calcNeutAllCovs (const TrkSimpTraj &, const BField &, double fltlen, HepVector &x0, HepVector &p0, HepSymMatrix &xxCov, HepSymMatrix &ppCov, HepMatrix &xpCov) |
static void | calcCurvAllWeights (const TrkSimpTraj &, const BField &, double fltlen, HepVector &pos, HepVector &mom, HepSymMatrix &xxWeight, HepSymMatrix &ppWeight, HepMatrix &xpWeight) |
static void | calcCurvAllWeightsOLD (const TrkSimpTraj &, const BField &, double fltlen, HepVector &pos, HepVector &mom, HepSymMatrix &xxWeight, HepSymMatrix &ppWeight, HepMatrix &xpWeight) |
static void | calcNeutAllWeights (const TrkSimpTraj &, const BField &, double fltlen, HepVector &pos, HepVector &mom, HepSymMatrix &xxWeight, HepSymMatrix &ppWeight, HepMatrix &xpWeight) |
static int | nearestInt (double) |
Definition at line 36 of file TrkMomCalculator.h.
TrkMomCalculator::TrkMomCalculator | ( | ) |
TrkMomCalculator::~TrkMomCalculator | ( | ) | [virtual] |
void TrkMomCalculator::calcCurvAllCovs | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepSymMatrix & | xxCov, | |||
HepSymMatrix & | ppCov, | |||
HepMatrix & | xpCov | |||
) | [static, private] |
Definition at line 584 of file TrkMomCalculator.cxx.
References BField::bFieldNominal(), BField::cmTeslaToGeVc, cos(), TrkParams::covariance(), TrkExchangePar::ex_d0, TrkExchangePar::ex_omega, TrkExchangePar::ex_phi0, TrkExchangePar::ex_tanDip, TrkExchangePar::ex_z0, DifIndepPar::parameter(), TrkSimpTraj::parameters(), phi0, s, sin(), and v.
Referenced by getAllCovs().
00589 { 00590 //------------------------------------------------------------------------ 00591 00592 const HepVector& v = theTraj.parameters()->parameter(); 00593 const HepSymMatrix& m = theTraj.parameters()->covariance(); 00594 00595 // fast inlined conversion from Helix to PX representation 00596 00597 // track parameters 00598 double s = v[TrkExchangePar::ex_tanDip]; 00599 double omega = v[TrkExchangePar::ex_omega]; 00600 double d0 = v[TrkExchangePar::ex_d0]; 00601 //double z0 = v[TrkExchangePar::ex_z0]; 00602 double phi0 = v[TrkExchangePar::ex_phi0]; 00603 double cosDip = 1/sqrt(1.0+s*s); 00604 double l = fltlen*cosDip ; 00605 00606 // calculate some sin,cos etc.. 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 // Calculate derivatives for Jacobian matrix 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 // Fill temporary variables for m 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 // Fill xxCov,xpCov,ppCov - nb. this code generated by script writecov.py 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 // code added for testing... 00737 // HepSymMatrix xxCovOld(3),ppCovOld(3); 00738 // HepMatrix xpCovOld(3,3); 00739 00740 // calcCurvAllCovsOLD(theTraj,theField,fltlen, 00741 // xxCovOld,ppCovOld,xpCovOld); 00742 00743 00744 // if (_HdiffCov==0){ 00745 // _manager = gblEnv->getGen()->ntupleManager(); 00746 // _HdiffCov = _manager->histogram("log10(|Diff|) old - new Cov",100,-20.,0.); 00747 // _HfdiffCov = _manager->histogram("log10(|fDiff|) (old - new)/old Cov",100,-20.,0.); 00748 // } 00749 00750 // for (int i=1;i<=3;i++){ 00751 // for (int j=i;j<=3;j++){ 00752 // _HdiffCov->accumulate(log10(fabs(xxCovOld(i,j)-xxCov(i,j)))); 00753 // _HdiffCov->accumulate(log10(fabs(ppCovOld(i,j)-ppCov(i,j)))); 00754 // _HfdiffCov->accumulate(log10 ( fabs( (xxCovOld(i,j)-xxCov(i,j))/xxCovOld(i,j)) ) ); 00755 // _HfdiffCov->accumulate(log10 ( fabs( (ppCovOld(i,j)-ppCov(i,j))/ppCovOld(i,j)) ) ); 00756 // } 00757 // } 00758 // for (int i=1;i<=3;i++){ 00759 // for (int j=1;j<=3;j++){ 00760 // _HdiffCov->accumulate(log10(fabs(xpCovOld(i,j)-xpCov(i,j)))); 00761 // _HfdiffCov->accumulate(log10 ( fabs( (xpCovOld(i,j)-xpCov(i,j))/xpCovOld(i,j)) ) ); 00762 // } 00763 // } 00764 00765 00766 }
void TrkMomCalculator::calcCurvAllCovsOLD | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepSymMatrix & | xxCov, | |||
HepSymMatrix & | ppCov, | |||
HepMatrix & | xpCov | |||
) | [static, private] |
Definition at line 770 of file TrkMomCalculator.cxx.
References DifNumber::absolute(), arg(), BField::bFieldNominal(), BField::cmTeslaToGeVc, correlation(), DifIndepPar::covariance(), DifVector::errorMatrix(), TrkDifTraj::getDFInfo(), if(), DifNumber::indepPar(), DifVector::length(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by getAllCovs().
00775 { 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 // computes position covariances 00800 00801 xxCov.assign(PosDif.errorMatrix(PosDif.x.indepPar()->covariance())); 00802 00803 // computes momentum covariances 00804 ppCov.assign(MomDif.errorMatrix(MomDif.x.indepPar()->covariance())); 00805 00806 // computes correlations 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 }
void TrkMomCalculator::calcCurvAllWeights | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepVector & | pos, | |||
HepVector & | mom, | |||
HepSymMatrix & | xxWeight, | |||
HepSymMatrix & | ppWeight, | |||
HepMatrix & | xpWeight | |||
) | [static, private] |
Definition at line 888 of file TrkMomCalculator.cxx.
References BField::bFieldNominal(), EvtCyclic3::C, calcCurvAllWeightsOLD(), BField::cmTeslaToGeVc, cos(), TrkExchangePar::ex_d0, TrkExchangePar::ex_omega, TrkExchangePar::ex_phi0, TrkExchangePar::ex_tanDip, TrkExchangePar::ex_z0, DifIndepPar::parameter(), TrkSimpTraj::parameters(), phi0, boss::pos, q, s, sin(), v, w, and TrkParams::weightMatrix().
Referenced by getAllWeights().
00895 { 00896 //------------------------------------------------------------------------ 00897 const HepVector& v = theTraj.parameters()->parameter(); 00898 const HepSymMatrix& w = theTraj.parameters()->weightMatrix(); 00899 00900 // fast inlined conversion from Helix to PX representation 00901 00902 // track parameters 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 // calculate position and momentum 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 // inverse of the jacobian matrix - see helix.mws Maple worksheet 00931 00932 // protect against divide by 0. 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 // local variables for the weight matrix 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 // calculate the Weight matrix from dinv (the indices are xpWeight(ip,ix) ) (used writewgt.py script) 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 // code added for testing... 01077 // HepVector posOld(3),momOld(3); 01078 // HepSymMatrix xxWeightOld(3),ppWeightOld(3); 01079 // HepMatrix xpWeightOld(3,3); 01080 01081 // calcCurvAllWeightsOLD(theTraj,theField,fltlen, 01082 // posOld,momOld,xxWeightOld,ppWeightOld,xpWeightOld); 01083 01084 // if (_Hdiff==0){ 01085 // _manager = gblEnv->getGen()->ntupleManager(); 01086 // _Hdiff = _manager->histogram("log10(|Diff|) old - new",100,-20.,0.); 01087 // _Hfdiff = _manager->histogram("log10(|fDiff|) (old - new)/old",100,-20.,0.); 01088 // } 01089 01090 // for (int i=1;i<=3;i++){ 01091 // for (int j=i;j<=3;j++){ 01092 // _Hdiff->accumulate(log10(fabs(xxWeightOld(i,j)-xxWeight(i,j)))); 01093 // _Hdiff->accumulate(log10(fabs(ppWeightOld(i,j)-ppWeight(i,j)))); 01094 // _Hfdiff->accumulate(log10 ( fabs( (xxWeightOld(i,j)-xxWeight(i,j))/xxWeightOld(i,j)) ) ); 01095 // _Hfdiff->accumulate(log10 ( fabs( (ppWeightOld(i,j)-ppWeight(i,j))/ppWeightOld(i,j)) ) ); 01096 // } 01097 // } 01098 // for (int i=1;i<=3;i++){ 01099 // for (int j=1;j<=3;j++){ 01100 // _Hdiff->accumulate(log10(fabs(xpWeightOld(i,j)-xpWeight(i,j)))); 01101 // _Hfdiff->accumulate(log10 ( fabs( (xpWeightOld(i,j)-xpWeight(i,j))/xpWeightOld(i,j)) ) ); 01102 // } 01103 // } 01104 01105 01106 }
void TrkMomCalculator::calcCurvAllWeightsOLD | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepVector & | pos, | |||
HepVector & | mom, | |||
HepSymMatrix & | xxWeight, | |||
HepSymMatrix & | ppWeight, | |||
HepMatrix & | xpWeight | |||
) | [static, private] |
Definition at line 1110 of file TrkMomCalculator.cxx.
References DifNumber::absolute(), arg(), BField::bFieldNominal(), BField::cmTeslaToGeVc, DifIndepPar::covariance(), TrkDifTraj::getDFInfo(), genRecEmupikp::i, if(), DifNumber::indepPar(), ganga-rec::j, DifVector::jacobian(), DifVector::length(), DifNumber::number(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by calcCurvAllWeights(), and getAllWeights().
01117 { 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 * start computing the inverse of the Jacobian needed: 01144 * 01145 * D(x,p) 01146 * ------ 01147 * D(n) 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 * now we need: 01164 * 01165 * D(x,p) 01166 * ------ 01167 * D(s) 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 * with an inversion we obtain 01181 * D(n) 01182 * -------- 01183 * D(x,p) 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 // this is the weight matrix for the helix parameters 01200 HepSymMatrix Wnn(PosDif.x.indepPar()->covariance()); 01201 01202 Wnn.invert(invStatus); 01203 /* 01204 * now we have the weight matrices 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 }
int TrkMomCalculator::calcCurvCharge | ( | const Hep3Vector & | , | |
double | curvature, | |||
const BField & | ||||
) | [static, private] |
Definition at line 447 of file TrkMomCalculator.cxx.
References BField::bFieldNominal(), calcCurvVecMom(), and nearestInt().
00449 { 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 }
BesVectorErr TrkMomCalculator::calcCurvErrMom | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | flt | |||
) | [static, private] |
Definition at line 336 of file TrkMomCalculator.cxx.
References DifNumber::absolute(), arg(), BField::bFieldNominal(), BField::cmTeslaToGeVc, DifIndepPar::covariance(), DifVector::errorMatrix(), TrkDifTraj::getDFInfo(), DifNumber::indepPar(), DifVector::length(), DifNumber::number(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by errMom().
00338 { 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 if (ErrLogging(debugging) && 0) { 00363 HepMatrix e1 = DirDif.errorMatrix(MomDif.x.indepPar()->covariance()); 00364 double e2 = momMag.error(MomDif.x.indepPar()->covariance()); 00365 HepMatrix e3 = MomDif.errorMatrix(MomDif.x.indepPar()->covariance()); 00366 const HepMatrix& c1 = MomDif.x.indepPar()->covariance(); 00367 00368 std::cout<<"ErrMsg(debugging) "<< endl << endl << "Start" << endl 00369 << "param cov: " << endl 00370 << c1(1,1) << endl 00371 << c1(2,1) << " " << c1(2,2) << endl 00372 << c1(3,1) << " " << c1(3,2) << " " << c1(3,3) << endl 00373 << c1(4,1) << " " << c1(4,2) << " " << c1(4,3) 00374 << " " << c1(4,4) << endl 00375 << c1(5,1) << " " << c1(5,2) << " " << c1(5,3) 00376 << " " << c1(5,4) << " " << c1(5,5) << endl 00377 00378 << "Dir " << e1.num_row() << " " << e1.num_col() << endl 00379 << "output:" << endl 00380 << e1(1,1) << endl 00381 << e1(2,1) << " " << e1(2,2) << endl 00382 << e1(3,1) << " " << e1(3,2) << " " << e1(3,3) << endl 00383 << "deriv: " << endl 00384 << DirDif.x.derivatives() << endl 00385 << endl 00386 00387 << "Mag " << endl 00388 << "output:" << endl 00389 << e2 << endl 00390 << "deriv: " << endl 00391 << momMag.derivatives() << endl 00392 << endl 00393 00394 << "Momdif " << e3.num_row() << " " << e3.num_col() << endl 00395 << "output:" << endl 00396 << e3(1,1) << endl 00397 << e3(2,1) << " " << e3(2,2) << endl 00398 << e3(3,1) << " " << e3(3,2) << " " << e3(3,3) << endl 00399 << "deriv: " << endl 00400 << MomDif.x.derivatives() << endl 00401 << endl 00402 00403 << "End" << endl << std::endl; 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 }
HepMatrix TrkMomCalculator::calcCurvPosmomCov | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static, private] |
Definition at line 481 of file TrkMomCalculator.cxx.
References DifNumber::absolute(), arg(), BField::bFieldNominal(), BField::cmTeslaToGeVc, correlation(), DifIndepPar::covariance(), TrkDifTraj::getDFInfo(), DifNumber::indepPar(), DifVector::length(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by posmomCov().
00483 { 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 // computes the correlation among position and momentum 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 }
double TrkMomCalculator::calcCurvPtMom | ( | const Hep3Vector & | , | |
double | curvature, | |||
const BField & | ||||
) | [static, private] |
Definition at line 304 of file TrkMomCalculator.cxx.
References calcCurvVecMom().
Referenced by ptMom().
00306 { 00307 //------------------------------------------------------------------------ 00308 Hep3Vector theMomVec = calcCurvVecMom(direction, curvature, 00309 theField); 00310 return theMomVec.perp(); 00311 00312 }
Hep3Vector TrkMomCalculator::calcCurvVecMom | ( | const Hep3Vector & | , | |
double | curvature, | |||
const BField & | ||||
) | [static, private] |
Definition at line 316 of file TrkMomCalculator.cxx.
References arg(), BField::bFieldNominal(), and BField::cmTeslaToGeVc.
Referenced by calcCurvCharge(), calcCurvPtMom(), and vecMom().
00318 { 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 }
void TrkMomCalculator::calcNeutAllCovs | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepVector & | x0, | |||
HepVector & | p0, | |||
HepSymMatrix & | xxCov, | |||
HepSymMatrix & | ppCov, | |||
HepMatrix & | xpCov | |||
) | [static, private] |
Definition at line 836 of file TrkMomCalculator.cxx.
References NeutParams::_p, correlation(), DifIndepPar::covariance(), DifIndepPar::difPar(), DifVector::errorMatrix(), TrkDifTraj::getDFInfo(), DifNumber::indepPar(), DifNumber::number(), TrkSimpTraj::parameters(), DifVector::x, DifVector::y, and DifVector::z.
00842 { 00843 //------------------------------------------------------------------------ 00844 00845 DifPoint PosDif; 00846 DifVector DirDif; 00847 DifVector delDirDif; 00848 00849 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif); 00850 00851 // set the momentum's direction, and then its magnitude 00852 DifVector MomDif = DirDif; 00853 00854 // The + 1 is necessary because DifIndepPar starts counting at 1, whereas 00855 // the enum for NeutParams starts at 0. Stupid, ain't it? 00856 MomDif *= theTraj.parameters()->difPar(NeutParams::_p + 1); 00857 00858 // fill momenta and positions 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 // computes momentum covariances 00867 xxCov.assign(PosDif.errorMatrix(PosDif.x.indepPar()->covariance())); 00868 00869 // computes momentum covariances 00870 ppCov.assign(MomDif.errorMatrix(MomDif.x.indepPar()->covariance())); 00871 00872 // computes correlations 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 }
void TrkMomCalculator::calcNeutAllCovs | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepSymMatrix & | xxCov, | |||
HepSymMatrix & | ppCov, | |||
HepMatrix & | xpCov | |||
) | [static, private] |
Definition at line 822 of file TrkMomCalculator.cxx.
Referenced by getAllCovs().
00827 { 00828 //------------------------------------------------------------------------ 00829 HepVector p0(3),x0(3); 00830 calcNeutAllCovs(theTraj,theField,fltlen,x0,p0,xxCov,ppCov,xpCov); 00831 00832 }
void TrkMomCalculator::calcNeutAllWeights | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepVector & | pos, | |||
HepVector & | mom, | |||
HepSymMatrix & | xxWeight, | |||
HepSymMatrix & | ppWeight, | |||
HepMatrix & | xpWeight | |||
) | [static, private] |
Definition at line 1223 of file TrkMomCalculator.cxx.
References NeutParams::_p, DifIndepPar::covariance(), DifIndepPar::difPar(), TrkDifTraj::getDFInfo(), genRecEmupikp::i, DifNumber::indepPar(), ganga-rec::j, DifVector::jacobian(), DifNumber::number(), TrkSimpTraj::parameters(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by getAllWeights().
01230 { 01231 //------------------------------------------------------------------------ 01232 DifPoint PosDif; 01233 DifVector DirDif; 01234 DifVector delDirDif; 01235 DifNumber momMag; 01236 01237 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif); 01238 01239 // set the momentum's direction, and then its magnitude 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 // this is the weight matrix for the helix parameters 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 }
BesVectorErr TrkMomCalculator::calcNeutErrMom | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | flt | |||
) | [static, private] |
Definition at line 414 of file TrkMomCalculator.cxx.
References NeutParams::_p, DifIndepPar::covariance(), DifIndepPar::difPar(), DifVector::errorMatrix(), TrkDifTraj::getDFInfo(), DifNumber::indepPar(), DifNumber::number(), TrkSimpTraj::parameters(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by errMom().
00416 { 00417 //------------------------------------------------------------------------ 00418 00419 DifPoint posDif; 00420 DifVector dirDif; 00421 DifVector delDirDif; 00422 00423 theTraj.getDFInfo(fltlen, posDif, dirDif, delDirDif); 00424 00425 // set the momentum's direction, and then its magnitude 00426 DifVector momDif = dirDif; 00427 00428 // The + 1 is necessary because DifIndepPar starts counting at 1, whereas 00429 // the enum for NeutParams starts at 0. Stupid, ain't it? 00430 momDif *= theTraj.parameters()->difPar(NeutParams::_p + 1); 00431 00432 // now create an error matrix 00433 BesError symErr(momDif.errorMatrix( 00434 momDif.x.indepPar()->covariance())); 00435 00436 // get the result in the correct form 00437 return BesVectorErr(Hep3Vector(momDif.x.number(), momDif.y.number(), 00438 momDif.z.number()), symErr); 00439 }
HepMatrix TrkMomCalculator::calcNeutPosmomCov | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static, private] |
Definition at line 528 of file TrkMomCalculator.cxx.
References NeutParams::_p, correlation(), DifIndepPar::covariance(), DifIndepPar::difPar(), TrkDifTraj::getDFInfo(), DifNumber::indepPar(), TrkSimpTraj::parameters(), DifVector::x, DifVector::y, and DifVector::z.
Referenced by posmomCov().
00530 { 00531 //------------------------------------------------------------------------ 00532 00533 DifPoint PosDif; 00534 DifVector DirDif; 00535 DifVector delDirDif; 00536 00537 theTraj.getDFInfo(fltlen, PosDif, DirDif, delDirDif); 00538 00539 // set the momentum's direction, and then its magnitude 00540 DifVector MomDif = DirDif; 00541 00542 // The + 1 is necessary because DifIndepPar starts counting at 1, whereas 00543 // the enum for NeutParams starts at 0. Stupid, ain't it? 00544 MomDif *= theTraj.parameters()->difPar(NeutParams::_p + 1); 00545 00546 // computes the correlation among position and momentum 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 }
int TrkMomCalculator::charge | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static] |
Definition at line 140 of file TrkMomCalculator.cxx.
References BField::bFieldNominal(), TrkMomVisitor::circle(), TrkMomVisitor::helix(), TrkMomVisitor::neut(), HelixTraj::omegaIndex, DifIndepPar::parameter(), and TrkSimpTraj::parameters().
Referenced by TrkSimpleRep::charge().
00141 { 00142 //------------------------------------------------------------------------ 00143 00144 TrkMomVisitor theVisitor(theTraj); 00145 00146 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) { 00147 00148 // treat as curve, calculate Pt accordingly (see calcCurvMom function 00149 // below) 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 // return calcCurvCharge( 00163 // theTraj.direction(fltlen), 00164 // theTraj.curvature(fltlen), theField); 00165 00166 } else if (theVisitor.neut() != 0) { 00167 00168 // treat as neutral particle, so charge is zero 00169 return 0; 00170 00171 } else { 00172 00173 // particle must be a plain line--take charge as zero 00174 return 0; 00175 } 00176 }
BesVectorErr TrkMomCalculator::errMom | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static] |
Definition at line 112 of file TrkMomCalculator.cxx.
References calcCurvErrMom(), calcNeutErrMom(), TrkMomVisitor::circle(), TrkMomVisitor::helix(), and TrkMomVisitor::neut().
Referenced by TrkSimpleRep::momentumErr(), and TrkCompTrk::momentumErr().
00113 { 00114 //------------------------------------------------------------------------ 00115 00116 TrkMomVisitor theVisitor(theTraj); 00117 00118 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) { 00119 00120 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function 00121 // below) 00122 return calcCurvErrMom(theTraj, theField, fltlen); 00123 00124 } else if (theVisitor.neut() != 0) { 00125 00126 // treat as neutral particle, same as curve in this case 00127 return calcNeutErrMom(theTraj, theField, fltlen); 00128 00129 } else { 00130 00131 // particle must be a plain line--no way to calculate momentum or err 00132 // The matrix is initialized to zero (see BesError constructor) 00133 BesError theErr(3); 00134 return BesVectorErr(Hep3Vector(999, 999, 999), theErr); 00135 } 00136 }
void TrkMomCalculator::getAllCovs | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepSymMatrix & | xxCov, | |||
HepSymMatrix & | ppCov, | |||
HepMatrix & | xpCov | |||
) | [static] |
Definition at line 207 of file TrkMomCalculator.cxx.
References calcCurvAllCovs(), calcCurvAllCovsOLD(), calcNeutAllCovs(), TrkMomVisitor::circle(), TrkMomVisitor::helix(), genRecEmupikp::i, ganga-rec::j, and TrkMomVisitor::neut().
Referenced by TrkSimpleRep::getAllCovs(), and TrkCompTrk::getAllCovs().
00212 { 00213 //------------------------------------------------------------------------ 00214 00215 TrkMomVisitor theVisitor(theTraj); 00216 00217 if (theVisitor.helix() != 0) { 00218 00219 // fast inline calculation... 00220 calcCurvAllCovs(theTraj,theField,fltlen,xxCov,ppCov,xpCov); 00221 00222 } else if (theVisitor.circle() != 0){ 00223 00224 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function 00225 // below) 00226 calcCurvAllCovsOLD(theTraj,theField,fltlen,xxCov,ppCov,xpCov); 00227 00228 } else if (theVisitor.neut() != 0) { 00229 00230 // treat as neutral particle, same as curve in this case 00231 calcNeutAllCovs(theTraj,theField,fltlen,xxCov,ppCov,xpCov); 00232 00233 } else { 00234 00235 // particle must be a plain line--no way to calculate momentum or err 00236 // The matrix is initialized to zero (see BesError constructor) 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 }
void TrkMomCalculator::getAllWeights | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen, | |||
HepVector & | pos, | |||
HepVector & | mom, | |||
HepSymMatrix & | xxWeight, | |||
HepSymMatrix & | ppWeight, | |||
HepMatrix & | xpWeight | |||
) | [static] |
Definition at line 252 of file TrkMomCalculator.cxx.
References calcCurvAllWeights(), calcCurvAllWeightsOLD(), calcNeutAllWeights(), TrkMomVisitor::circle(), TrkMomVisitor::helix(), genRecEmupikp::i, ganga-rec::j, and TrkMomVisitor::neut().
Referenced by TrkSimpleRep::getAllWeights(), and TrkCompTrk::getAllWeights().
00259 { 00260 //------------------------------------------------------------------------ 00261 00262 TrkMomVisitor theVisitor(theTraj); 00263 00264 if (theVisitor.helix() != 0){ 00265 00266 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function 00267 // below) 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 // treat as neutral particle, same as curve in this case 00279 calcNeutAllWeights(theTraj,theField,fltlen, 00280 pos,mom,xxWeight,ppWeight,xpWeight); 00281 00282 } else { 00283 00284 // particle must be a plain line--no way to calculate momentum or err 00285 // temporary: initialize everything to 0 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 }
int TrkMomCalculator::nearestInt | ( | double | ) | [static, private] |
Definition at line 464 of file TrkMomCalculator.cxx.
Referenced by calcCurvCharge().
00464 { 00465 //------------------------------------------------------------------------ 00466 00467 if (floatPt > 0.) { 00468 return (int)(floatPt+0.5); 00469 } else { 00470 return (int)(floatPt-0.5); 00471 } 00472 }
HepMatrix TrkMomCalculator::posmomCov | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static] |
Definition at line 180 of file TrkMomCalculator.cxx.
References calcCurvPosmomCov(), calcNeutPosmomCov(), TrkMomVisitor::circle(), TrkMomVisitor::helix(), and TrkMomVisitor::neut().
Referenced by TrkSimpleRep::posmomCov(), and TrkCompTrk::posmomCov().
00181 { 00182 //------------------------------------------------------------------------ 00183 00184 TrkMomVisitor theVisitor(theTraj); 00185 00186 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) { 00187 00188 // treat as curve, calculate ErrMom accordingly (see calcCurvMom function 00189 // below) 00190 return calcCurvPosmomCov(theTraj, theField, fltlen); 00191 00192 } else if (theVisitor.neut() != 0) { 00193 00194 // treat as neutral particle, same as curve in this case 00195 return calcNeutPosmomCov(theTraj, theField, fltlen); 00196 00197 } else { 00198 00199 // particle must be a plain line--no way to calculate momentum or err 00200 // The matrix is initialized to zero (see BesError constructor) 00201 return HepMatrix(3,3,0); 00202 } 00203 }
double TrkMomCalculator::ptMom | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static] |
Definition at line 47 of file TrkMomCalculator.cxx.
References NeutParams::_p, arg(), calcCurvPtMom(), TrkMomVisitor::circle(), Trajectory::curvature(), Trajectory::direction(), TrkMomVisitor::helix(), TrkMomVisitor::neut(), DifIndepPar::parameter(), and TrkSimpTraj::parameters().
Referenced by TrkSimpleRep::pt(), and TrkCompTrk::pt().
00048 { 00049 //------------------------------------------------------------------------ 00050 00051 TrkMomVisitor theVisitor(theTraj); 00052 00053 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) { 00054 00055 // treat as curve, calculate Pt accordingly (see calcCurvMom function 00056 // below) 00057 return calcCurvPtMom(theTraj.direction(fltlen), 00058 theTraj.curvature(fltlen), theField); 00059 00060 } else if (theVisitor.neut() != 0) { 00061 00062 // treat as neutral particle (with the ptot as a parameter) 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 // particle must be a plain line--no way to calculate momentum 00074 return 999999.99; 00075 00076 } 00077 00078 }
Hep3Vector TrkMomCalculator::vecMom | ( | const TrkSimpTraj & | , | |
const BField & | , | |||
double | fltlen | |||
) | [static] |
Definition at line 82 of file TrkMomCalculator.cxx.
References NeutParams::_p, calcCurvVecMom(), TrkMomVisitor::circle(), Trajectory::curvature(), Trajectory::direction(), TrkMomVisitor::helix(), TrkMomVisitor::neut(), DifIndepPar::parameter(), and TrkSimpTraj::parameters().
Referenced by TrkSimpleRep::momentum(), and TrkCompTrk::momentum().
00083 { 00084 //------------------------------------------------------------------------ 00085 TrkMomVisitor theVisitor(theTraj); 00086 00087 if (theVisitor.helix() != 0 || theVisitor.circle() != 0) { 00088 00089 // treat as curve, calculate VecMom accordingly (see calcCurvMom function 00090 // below) 00091 00092 return calcCurvVecMom(theTraj.direction(fltlen), 00093 theTraj.curvature(fltlen), theField); 00094 00095 } else if (theVisitor.neut() != 0) { 00096 00097 // treat as neutral particle (with the pt as a parameter) 00098 Hep3Vector theMom = theTraj.direction(fltlen); 00099 theMom.setMag(theTraj.parameters()->parameter()[NeutParams::_p]); 00100 return theMom; 00101 00102 } else { 00103 00104 // particle must be a plain line--no way to calculate momentum 00105 return Hep3Vector(999, 999, 999); 00106 00107 } 00108 }
bool TrkMomCalculator::weightToCov | ( | const HepSymMatrix & | inXX, | |
const HepSymMatrix & | inPP, | |||
const HepMatrix & | inXP, | |||
HepSymMatrix & | outXX, | |||
HepSymMatrix & | outPP, | |||
HepMatrix & | outXP | |||
) | [static] |
Definition at line 562 of file TrkMomCalculator.cxx.
References alpha.
00563 { 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 }