00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "MdcGeom/Constants.h"
00019 #include <assert.h>
00020 #include <algorithm>
00021 #include <iostream>
00022 #include "TrkBase/TrkRep.h"
00023 #include "MdcRecoUtil/Pdt.h"
00024 #include "MdcRecoUtil/PdtEntry.h"
00025 #include "TrkBase/TrkDifTraj.h"
00026 #include "TrkBase/TrkHotListFull.h"
00027 #include "TrkBase/TrkHotListEmpty.h"
00028 #include "TrkBase/TrkHitOnTrk.h"
00029 #include "TrkBase/TrkRecoTrk.h"
00030 #include "TrkBase/TrkFunctors.h"
00031 #include "TrkBase/TrkErrCode.h"
00032 #include "MdcRecoUtil/DifPoint.h"
00033 #include "MdcRecoUtil/DifVector.h"
00034 #include "MdcRecoUtil/DifIndepPar.h"
00035 #include "ProbTools/ChisqConsistency.h"
00036 #include "ProxyDict/IfdIntKey.h"
00037 #include "TrkBase/TrkExchangePar.h"
00038 using std::cout;
00039 using std::endl;
00040
00041 TrkRep::TrkRep(TrkRecoTrk* trk, PdtPid::PidType hypo,bool createHotList)
00042 : _hotList( createHotList?new TrkHotListFull:0 )
00043 {
00044 init(trk, hypo);
00045 }
00046
00047 TrkRep::TrkRep(const TrkHotList& hotlist, TrkRecoTrk* trk,
00048 PdtPid::PidType hypo)
00049 :_hotList( hotlist.clone(TrkBase::Functors::cloneHot(this)) )
00050 {
00051 init(trk, hypo);
00052 }
00053
00054 TrkRep::TrkRep(TrkHotList& hotlist, TrkRecoTrk* trk,
00055 PdtPid::PidType hypo, bool stealHots)
00056 {
00057 init(trk, hypo);
00058 if (!stealHots) {
00059 _hotList.reset( hotlist.clone(TrkBase::Functors::cloneHot(this)) );
00060 } else {
00061 _hotList.reset( new TrkHotListFull(hotlist,setParent(this)) );
00062 }
00063 }
00064
00065 TrkRep::TrkRep(const TrkHotList* hotlist, TrkRecoTrk* trk,
00066 PdtPid::PidType hypo)
00067 {
00068
00069
00070 init(trk, hypo);
00071 _hotList.reset( hotlist!=0?
00072 hotlist->clone(TrkBase::Functors::cloneHot(this)):
00073 new TrkHotListFull );
00074 }
00075
00076 TrkRep::TrkRep(TrkHotList* hotlist, TrkRecoTrk* trk,
00077 PdtPid::PidType hypo,bool takeownership)
00078 {
00079 init(trk, hypo);
00080 if (!takeownership) {
00081 _hotList.reset( hotlist!=0?
00082 hotlist->clone(TrkBase::Functors::cloneHot(this)):
00083 new TrkHotListFull );
00084 } else {
00085 assert(hotlist!=0);
00086 _hotList.reset( hotlist->resetParent(setParent(this)) );
00087 }
00088 }
00089
00090
00091 TrkRep::TrkRep(TrkRecoTrk* trk, PdtPid::PidType hypo, int nact, int nsv,
00092 int ndc, double stFndRng, double endFndRng)
00093 :_hotList(new TrkHotListEmpty(nact, nsv, ndc, stFndRng, endFndRng))
00094 {
00095
00096
00097 init(trk, hypo);
00098 }
00099
00100
00101 TrkRep::TrkRep(const TrkRep& oldRep, TrkRecoTrk* trk, PdtPid::PidType hypo) :
00102 TrkFitStatus(oldRep)
00103 {
00104
00105
00106 init(trk, hypo);
00107
00108 }
00109
00110 TrkRep&
00111 TrkRep::operator= (const TrkRep& right)
00112 {
00113 if(&right != this){
00114 init(right._parentTrack, right._partHypo);
00115 _hotList.reset( right._hotList->clone(this) );
00116 TrkFitStatus::operator=(right);
00117 }
00118 return *this;
00119 }
00120
00121 void
00122 TrkRep::init(TrkRecoTrk* trk, PdtPid::PidType hypo)
00123 {
00124 _parentTrack = trk;
00125 _partHypo = hypo;
00126 _betainv = -999999.;
00127 }
00128
00129 TrkRep::~TrkRep()
00130 {
00131 }
00132
00133 bool
00134 TrkRep::operator== (const TrkRep& rhs)
00135 {
00136 return (&rhs == this);
00137 }
00138
00139 void
00140 TrkRep::addHot(TrkHitOnTrk *newHot)
00141 {
00142 if (newHot->isActive()) setCurrent(false);
00143 hotList()->append(newHot);
00144 }
00145
00146 void
00147 TrkRep::removeHot(TrkHitOnTrk *theHot)
00148 {
00149 if(theHot->isActive()) setCurrent(false);
00150 hotList()->remove(theHot);
00151 }
00152
00153 void
00154 TrkRep::activateHot(TrkHitOnTrk* hot)
00155 {
00156 if(!hot->isActive()){
00157
00158 if(this == hot->getParentRep()){
00159 setCurrent(false);
00160
00161 hot->setActive(true);
00162 }
00163 }
00164 }
00165
00166 void
00167 TrkRep::deactivateHot(TrkHitOnTrk* hot)
00168 {
00169 if(hot->isActive()){
00170
00171 if(this == hot->getParentRep()){
00172 setCurrent(false);
00173
00174 hot->setActive(false);
00175 }
00176 }
00177 }
00178
00179 HepPoint3D
00180 TrkRep::position(double fltL) const
00181 {
00182 return traj().position(fltL);
00183 }
00184
00185 Hep3Vector
00186 TrkRep::direction(double fltL) const
00187 {
00188 return traj().direction(fltL);
00189 }
00190
00191 double
00192 TrkRep::arrivalTime(double fltL) const
00193 {
00194 static double cinv = 1./Constants::c;
00195
00196 if (_betainv < 0.0) {
00197 double mass2 = Pdt::lookup(particleType())->mass();
00198 mass2 = mass2 * mass2;
00199 double ptot2 = momentum(0.).mag2();
00200 assert(ptot2 != 0.0);
00201 _betainv = sqrt( (ptot2 + mass2)/ ptot2);
00202 }
00203 double tof = fltL * _betainv * cinv;
00204 return trackT0() + tof;
00205 }
00206
00207 double
00208 TrkRep::trackT0() const
00209 {
00210 return parentTrack()->trackT0();
00211 }
00212
00213 BesPointErr
00214 TrkRep::positionErr(double fltL) const
00215 {
00216 static DifPoint posD;
00217 static DifVector dirD;
00218 traj().getDFInfo2(fltL, posD, dirD);
00219 HepMatrix err = posD.errorMatrix( posD.x.indepPar()->covariance() );
00220 HepPoint3D point(posD.x.number(), posD.y.number(), posD.z.number());
00221 BesError symErr(3);
00222 symErr.assign(err);
00223
00224 if (false) {
00225 #ifdef MDCPATREC_ROUTINE
00226 cout<< "ErrMsg(routine) " << "Pos "
00227 << err.num_row() << " " << err.num_col() << endl
00228 << "output:" << endl
00229
00230
00231
00232 << "x deriv: " << endl
00233 << posD.x.derivatives() << endl
00234 << "y deriv: " << endl
00235 << posD.y.derivatives() << endl
00236 << endl;
00237 #endif
00238
00239
00240 Hep3Vector pointDir(point.x(), point.y());
00241 double dirMag = pointDir.mag();
00242 double dist = 5.e-3;
00243 double delx = dist * point.x() / dirMag;
00244 double dely = dist * point.y() / dirMag;
00245 int ierr = 0;
00246 HepMatrix weight = err.inverse(ierr);
00247 double chisq = weight(1,1) * delx * delx +
00248 2 * weight(2,1) * delx * dely +
00249 weight(2,2) * dely * dely;
00250 #ifdef MDCPATREC_DEBUG
00251 cout << point << endl;
00252 cout << symErr << endl;
00253 cout << "delta: " << delx << " " << dely << endl;
00254 cout << "chisq: " << chisq << endl;
00255 #endif
00256 double phi0 = helix(fltL).phi0();
00257 delx = dist * cos(phi0);
00258 dely = dist * sin(phi0);
00259 chisq = weight(1,1) * delx * delx +
00260 2 * weight(2,1) * delx * dely +
00261 weight(2,2) * dely * dely;
00262 #ifdef MDCPATREC_DEBUG
00263 cout << "delta: " << delx << " " << dely << endl;
00264 cout << "chisq: " << chisq << endl;
00265 cout << endl << endl;
00266 #endif
00267 }
00268 return BesPointErr(point, symErr);
00269 }
00270
00271 BesVectorErr
00272 TrkRep::directionErr(double fltL) const
00273 {
00274 static DifPoint posD;
00275 static DifVector dirD;
00276 traj().getDFInfo2(fltL, posD, dirD);
00277 BesError symErr(3);
00278 symErr.assign( dirD.errorMatrix( dirD.x.indepPar()->covariance() ));
00279 Hep3Vector dir(dirD.x.number(), dirD.y.number(), dirD.z.number());
00280 return BesVectorErr(dir, symErr);
00281 }
00282
00283 double
00284 TrkRep::startValidRange() const
00285 {
00286 return traj().lowRange();
00287 }
00288
00289 double
00290 TrkRep::endValidRange() const
00291 {
00292 return traj().hiRange();
00293 }
00294
00295 double
00296 TrkRep::startFoundRange() const
00297 {
00298 return hotList()->startFoundRange();
00299 }
00300
00301 double
00302 TrkRep::endFoundRange() const
00303 {
00304 return hotList()->endFoundRange();
00305 }
00306
00307 PdtPid::PidType
00308 TrkRep::particleType() const
00309 {
00310 return _partHypo;
00311 }
00312
00313 const IfdKey&
00314 TrkRep::myKey() const
00315 {
00316
00317
00318 static IfdIntKey _theKey(0);
00319 return _theKey;
00320 }
00321
00322 void
00323 TrkRep::updateHots()
00324 {
00325 setCurrent(false);
00326 hotList()->updateHots();
00327 }
00328
00329 int
00330 TrkRep::nActive() const
00331 {
00332 return hotList()->nActive();
00333 }
00334
00335 int
00336 TrkRep::nSvt() const
00337 {
00338 return hotList()->nSvt();
00339 }
00340
00341 int
00342 TrkRep::nMdc() const
00343 {
00344 return hotList()->nMdc();
00345 }
00346
00347 bool
00348 TrkRep::resid(const TrkHitOnTrk *h,
00349 double& residual, double& residErr,
00350 bool exclude) const
00351 {
00352 assert (h != 0);
00353 if (h->parentRep() != this) return false;
00354 if (!h->hasResidual()) return false;
00355 if (exclude) return false;
00356
00357 residual=h->residual();
00358 residErr=h->hitRms();
00359 return true;
00360 }
00361
00362 ChisqConsistency
00363 TrkRep::chisqConsistency() const {
00364 if(fitValid())
00365 return ChisqConsistency(chisq(),nDof());
00366 else
00367 return ChisqConsistency();
00368 }