00001
00002 #include <assert.h>
00003 #include <iostream>
00004 #include "TrkBase/TrkCompTrk.h"
00005 #include "TrkBase/TrkSimpTraj.h"
00006 #include "TrkBase/HelixTraj.h"
00007 #include "TrkBase/NeutTraj.h"
00008 #include "TrkBase/NeutParams.h"
00009 #include "TrkBase/TrkMomCalculator.h"
00010
00011
00012 #include "ProxyDict/IfdStrKey.h"
00013 #include "BField/BField.h"
00014 #include "TrkBase/TrkExchangePar.h"
00015 #include "TrkBase/TrkPoca.h"
00016 #include "MdcRecoUtil/BesPointErr.h"
00017 #include "MdcRecoUtil/BesVectorErr.h"
00018 #include "MdcRecoUtil/DifPoint.h"
00019 #include "MdcRecoUtil/DifVector.h"
00020 #include "MdcRecoUtil/DifIndepPar.h"
00021 #include "TrkBase/TrkHelixUtils.h"
00022 #include "MdcGeom/BesAngle.h"
00023 using std::endl;
00024 using std::ostream;
00025
00026
00027
00028 TrkCompTrk::TrkCompTrk(const BesPointErr& pos,
00029 const BesVectorErr& mom,
00030 const HepMatrix& xpCov,
00031 int charge, double chisq, int nDof, const BField* bf) :
00032 _chisq(chisq),_nDof(nDof) {
00033
00034 _bf = bf;
00035 _charge=charge;
00036 if(_charge!=0) {
00037 TrkExchangePar par1 = TrkHelixUtils::helixFromMomErr(pos,mom,xpCov,charge,bField());
00038 _traj.reset( new HelixTraj(par1.params(), par1.covariance()) );
00039 }
00040 else {
00041 _traj.reset( new NeutTraj(TrkHelixUtils::lineFromMomErr(pos,mom,xpCov,1,bField())) );
00042 }
00043 }
00044
00045
00046
00047 TrkCompTrk::TrkCompTrk(const TrkCompTrk& rhs) : _bf(rhs._bf),_chisq(rhs._chisq),_nDof(rhs._nDof) {
00048
00049 _traj.reset((TrkSimpTraj*)rhs.traj().clone());
00050 _charge=rhs.charge();
00051 }
00052
00053
00054 TrkCompTrk::~TrkCompTrk() {
00055
00056 }
00057
00058
00059 const TrkCompTrk&
00060 TrkCompTrk::operator=(const TrkCompTrk &right) {
00061
00062 if (&right == this) return *this;
00063 _bf = right._bf;
00064 _traj.reset((TrkSimpTraj*)right.traj().clone());
00065 _chisq=right.chisq();
00066 _nDof=right.nDof();
00067 _charge=right.charge();
00068 return *this;
00069 }
00070
00071
00072 int
00073 TrkCompTrk::nDof() const {
00074
00075 return _nDof;
00076 }
00077
00078
00079 double
00080 TrkCompTrk::chisq() const {
00081
00082 return _chisq;
00083 }
00084
00085
00086
00087
00088 int
00089 TrkCompTrk::charge() const {
00090
00091 return _charge;
00092 }
00093
00094
00095
00096 const TrkDifTraj&
00097 TrkCompTrk::traj() const {
00098
00099 return *_traj;
00100 }
00101
00102
00103 HepPoint3D
00104 TrkCompTrk::position( double fltL) const {
00105
00106 return traj().position(fltL);
00107 }
00108
00109
00110 Hep3Vector
00111 TrkCompTrk::direction(double fltL) const {
00112
00113 return traj().direction(fltL);
00114 }
00115
00116 Hep3Vector
00117 TrkCompTrk::momentum(double fltL) const
00118 {
00119 return TrkMomCalculator::vecMom(*_traj, bField(), fltL);
00120
00121 }
00122
00123 BesPointErr
00124 TrkCompTrk::positionErr(double fltL) const
00125 {
00126
00127 DifPoint posD;
00128 DifVector dirD;
00129 traj().getDFInfo2(fltL, posD, dirD);
00130 HepMatrix err = posD.errorMatrix( posD.x.indepPar()->covariance() );
00131 BesError symErr(3);
00132 symErr.assign(err);
00133 HepPoint3D point(posD.x.number(), posD.y.number(), posD.z.number());
00134 return BesPointErr(point, symErr);
00135 }
00136
00137 BesVectorErr
00138 TrkCompTrk::directionErr( double fltL) const
00139 {
00140
00141 DifPoint posD;
00142 DifVector dirD;
00143 traj().getDFInfo2(fltL, posD, dirD);
00144 HepMatrix err = dirD.errorMatrix( dirD.x.indepPar()->covariance() );
00145 BesError symErr(3);
00146 symErr.assign(err);
00147 Hep3Vector dir(dirD.x.number(), dirD.y.number(), dirD.z.number());
00148 return BesVectorErr(dir, symErr);
00149 }
00150
00151 BesVectorErr
00152 TrkCompTrk::momentumErr(double fltL) const
00153 {
00154 return TrkMomCalculator::errMom(*_traj, bField(), fltL);
00155 }
00156
00157 double
00158 TrkCompTrk::pt(double fltL) const
00159 {
00160 return TrkMomCalculator::ptMom(*_traj, bField(), fltL);
00161
00162
00163
00164 }
00165
00166
00167 double
00168 TrkCompTrk::startValidRange() const
00169 {
00170 return traj().lowRange();
00171 }
00172
00173 double
00174 TrkCompTrk::endValidRange() const
00175 {
00176 return traj().hiRange();
00177 }
00178
00179 void
00180 TrkCompTrk::print(ostream& ostr) const
00181 {
00182 ostr << "Traj: ";
00183 if(_charge==-1 || _charge==1) {
00184 ostr << "Charged Particle -> Helix Trajectory" << endl;
00185 HelixTraj& theTraj=(HelixTraj&)traj();
00186 ostr << " phi0=" << BesAngle(theTraj.phi0()).rad() << endl;
00187 ostr << " d0=" << theTraj.d0() << endl;
00188 ostr << " z0=" << theTraj.z0() << endl;
00189 ostr << " omega=" << theTraj.omega() << endl;
00190 ostr << " tanDip=" << theTraj.tanDip() << endl;
00191 }
00192 else {
00193 ostr << "Neutral Particle -> NeutTraj" << endl;
00194 NeutTraj& theTraj=(NeutTraj&)traj();
00195 ostr << " phi0=" << BesAngle(theTraj.params().phi0()).rad() << endl;
00196 ostr << " d0=" << theTraj.params().d0() << endl;
00197 ostr << " z0=" << theTraj.params().z0() << endl;
00198 ostr << " p=" << theTraj.params().p() << endl;
00199 ostr << " tanDip=" << theTraj.params().tanDip() << endl;
00200 }
00201 }
00202
00203 void
00204 TrkCompTrk::printAll(ostream& ostr) const
00205 {
00206 print(ostr);
00207 }
00208
00209 HepMatrix TrkCompTrk::posmomCov(double fltL) const
00210 {
00211 const BField& theField = bField();
00212 return TrkMomCalculator::posmomCov(*_traj, theField, fltL);
00213 }
00214
00215 void TrkCompTrk::getAllCovs(double fltL,
00216 HepSymMatrix& xxCov,
00217 HepSymMatrix& ppCov,
00218 HepMatrix& xpCov) const
00219 {
00220 const BField& theField = bField();
00221 TrkMomCalculator::getAllCovs(*_traj, theField, fltL,
00222 xxCov,ppCov,xpCov);
00223 }
00224
00225 void TrkCompTrk::getAllWeights(double fltL,
00226 HepVector& pos,
00227 HepVector& mom,
00228 HepSymMatrix& xxWeight,
00229 HepSymMatrix& ppWeight,
00230 HepMatrix& xpWeight) const
00231 {
00232 const BField& theField = bField();
00233 TrkMomCalculator::getAllWeights(*_traj, theField, fltL,
00234 pos,mom,xxWeight,ppWeight,xpWeight);
00235
00236 }
00237
00238
00239 void TrkCompTrk::getAllWeights(const HepPoint3D& pt,
00240 HepVector& pos,
00241 HepVector& mom,
00242 HepSymMatrix& xxWeight,
00243 HepSymMatrix& ppWeight,
00244 HepMatrix& xpWeight)const
00245 {
00246 double fltL=0;
00247 TrkPoca poca(traj(),fltL, pt);
00248 fltL = poca.flt1();
00249 getAllWeights(fltL,pos,mom,xxWeight,ppWeight,xpWeight);
00250 }