00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "MdcGeom/Constants.h"
00010 #include "TrkFitter/TrkHelixRep.h"
00011 #include "BField/BField.h"
00012 #include "TrkBase/TrkSimpTraj.h"
00013 #include "TrkBase/TrkExchangePar.h"
00014 #include "MdcGeom/BesAngle.h"
00015 #include "ProxyDict/IfdIntKey.h"
00016 #include "TrkBase/TrkRecoTrk.h"
00017 using std::endl;
00018 using std::ostream;
00019
00020
00021 TrkHelixRep::TrkHelixRep(const TrkExchangePar& inPar,
00022 TrkRecoTrk* trk,
00023 PdtPid::PidType hypo,
00024 const TrkHotList *hots) :
00025 TrkSimpleRep(hots, trk, hypo),
00026 _traj(inPar)
00027 {}
00028
00029
00030
00031 TrkHelixRep::TrkHelixRep(const TrkExchangePar& inPar, TrkRecoTrk* trk,
00032 PdtPid::PidType hypo, int nact, int nsv, int ndc, double chi2,
00033 double stFndRng, double endFndRng) :
00034 TrkSimpleRep(trk, hypo, nact, nsv, ndc, chi2, stFndRng, endFndRng),
00035 _traj(inPar.params(), inPar.covariance())
00036 {}
00037
00038
00039 TrkHelixRep::TrkHelixRep(const TrkHelixRep& old, TrkRecoTrk* trk, PdtPid::PidType hypo) :
00040 TrkSimpleRep(old, trk, hypo), _traj(old._traj)
00041 {}
00042
00043 TrkHelixRep::~TrkHelixRep()
00044 {}
00045
00046 TrkHelixRep*
00047 TrkHelixRep::clone(TrkRecoTrk* theTrack) const
00048 {
00049 TrkHelixRep* newRep = new TrkHelixRep(*this, theTrack, this->particleType());
00050 newRep->setCurrent(fitCurrent());
00051 return newRep;
00052 }
00053
00054 TrkHelixRep*
00055 TrkHelixRep::cloneNewHypo(PdtPid::PidType hypo)
00056 {
00057 TrkHelixRep* newRep = new TrkHelixRep(*this, parentTrack(), hypo);
00058 newRep->setValid(fitValid());
00059 newRep->setCurrent(fitCurrent());
00060 return newRep;
00061 }
00062
00063 TrkExchangePar
00064 TrkHelixRep::helix(double ) const
00065 {
00066 HepVector outPar = _traj.parameters()->parameter();
00067 outPar[HelixTraj::phi0Index] = BesAngle(outPar[HelixTraj::phi0Index]).rad();
00068 return TrkExchangePar(outPar, _traj.parameters()->covariance());
00069 }
00070
00071 void
00072 TrkHelixRep::print(ostream& ostr) const
00073 {
00074 ostr << "TrkHelixRep "
00075 << " q: "<<charge()
00076 << " phi0: " << BesAngle(_traj.phi0()).rad()
00077 << " om: "<<_traj.omega()
00078 << " pt: "<<pt()
00079 << " p: "<<momentum()
00080 << " d0: " << _traj.d0()
00081 << " z0: " << _traj.z0()
00082 << " ct: " << _traj.tanDip()
00083 << "parent track: "<<parentTrack()->id() << endl;
00084 }
00085
00086 void
00087 TrkHelixRep::printAll(ostream& ostr) const
00088 {
00089 print(ostr);
00090 }
00091
00092 TrkSimpTraj&
00093 TrkHelixRep::simpTraj()
00094 {
00095 return _traj;
00096 }
00097
00098 const TrkSimpTraj&
00099 TrkHelixRep::simpTraj() const
00100 {
00101 return _traj;
00102 }
00103
00104 TrkDifTraj&
00105 TrkHelixRep::traj()
00106 {
00107 return _traj;
00108 }
00109
00110 const TrkDifTraj&
00111 TrkHelixRep::traj() const
00112 {
00113 return _traj;
00114 }
00115
00116
00117 const IfdKey&
00118 TrkHelixRep::myKey() const
00119 {
00120 static IfdIntKey _theKey(3133);
00121 return _theKey;
00122 }
00123
00124
00125 bool
00126 TrkHelixRep::resid(const TrkHitOnTrk *h, double& residual, double& residErr,
00127 bool exclude) const
00128 {
00129 if (this != h->getParentRep() ) return false;
00130 static HepVector derivs; double deltaChi;
00131 TrkErrCode errCode = h->getFitStuff( derivs, deltaChi );
00132
00133 if (errCode.failure()) return false;
00134 if (!exclude || !h->isActive()) {
00135 HepSymMatrix c = _traj.parameters()->covariance().similarityT(derivs);
00136 double tErr2 = c(1,1);
00137 double rms = h->hitRms();
00138 residual=h->residual();
00139 double e2=rms*rms+(h->isActive()?-1:1)*tErr2;
00140 if (e2<0) return false;
00141 residErr=sqrt(e2);
00142 } else {
00143 assert(h->isActive());
00144
00145
00146 HepSymMatrix w=_traj.parameters()->covariance();
00147 int ier;
00148 w.invert(ier);
00149 if (ier!=0) return false;
00150 HepVector p0 = _traj.parameters()->parameter();
00151 HepVector p1 = w*p0;
00152 p1 -= deltaChi*derivs;
00153 for (unsigned i=0; i<p0.num_row(); ++i) for (unsigned j=0; j<i+1 ;++j) {
00154 w.fast(i+1,j+1) -= derivs[i]*derivs[j];
00155 }
00156 w.invert(ier);
00157 if (ier!=0) return false;
00158 residual=h->residual()+h->hitRms()*dot(derivs,w*p1-p0);
00159 residErr=h->hitRms()*sqrt(w(1,1)+1);
00160 }
00161 return true;
00162 }