TrkMomCalculator Class Reference

#include <TrkMomCalculator.h>

List of all members.

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)


Detailed Description

Definition at line 36 of file TrkMomCalculator.h.


Constructor & Destructor Documentation

TrkMomCalculator::TrkMomCalculator (  ) 

Definition at line 41 of file TrkMomCalculator.cxx.

00041                                    {
00042 //------------------------------------------------------------------------
00043 }

TrkMomCalculator::~TrkMomCalculator (  )  [virtual]

Definition at line 36 of file TrkMomCalculator.cxx.

00036                                     {
00037 //------------------------------------------------------------------------
00038 }


Member Function Documentation

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 }


Generated on Tue Nov 29 23:36:15 2016 for BOSS_7.0.2 by  doxygen 1.4.7