00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <vector>
00012 #include <assert.h>
00013 #include <functional>
00014 #include <algorithm>
00015 #include "MdcGeom/Constants.h"
00016 #include "TrkBase/TrkRecoTrk.h"
00017 #include "TrkBase/TrkHitOnTrk.h"
00018 #include "TrkBase/TrkFundHit.h"
00019 #include "TrkBase/TrkDifTraj.h"
00020 #include "TrkBase/TrkExchangePar.h"
00021 #include "TrkBase/TrkErrCode.h"
00022 #include "TrkBase/TrkExtInterface.h"
00023 #include "TrkBase/TrkFitStatus.h"
00024 #include "TrkBase/TrkRepIter.h"
00025 #include "CLHEP/Vector/ThreeVector.h"
00026 #include "MdcRecoUtil/Pdt.h"
00027 #include "MdcRecoUtil/BesPointErr.h"
00028 #include "MdcRecoUtil/BesVectorErr.h"
00029 #include "MdcRecoUtil/DifPoint.h"
00030 #include "MdcRecoUtil/DifVector.h"
00031 #include "MdcRecoUtil/DifIndepPar.h"
00032 #include "TrkBase/TrkRep.h"
00033 #include "TrkBase/TrkContext.h"
00034 #include "MdcRecoUtil/PdtPid.h"
00035 #include "boost/shared_ptr.hpp"
00036 #include "BField/BField.h"
00037 using std::ostream;
00038
00039
00040 class TrkRecoTrkImpl {
00041 public:
00042 friend class TrkRecoTrk;
00043 private:
00044 TrkRecoTrkImpl()
00045 : _reps(PdtPid::nPidType,boost::shared_ptr<TrkRep>((TrkRep *)0)),
00046 _hitInterfaces(PdtPid::nPidType,boost::shared_ptr<TrkHitList>((TrkHitList *)0))
00047 { }
00048
00049
00050
00051
00052 typedef std::vector<boost::shared_ptr<TrkRep> > repList_t;
00053 typedef repList_t::const_iterator repConstIter;
00054 typedef repList_t::iterator repIter;
00055 repList_t _reps;
00056 typedef std::vector<boost::shared_ptr<TrkHitList> > hitList_t;
00057 typedef hitList_t::iterator hitListIter;
00058 hitList_t _hitInterfaces;
00059
00060 };
00061
00062 TrkRecoTrk::TrkRecoTrk(PdtPid::PidType defaultPart, const TrkContext& ctext,
00063 double t0) :
00064 _impl(new TrkRecoTrkImpl),
00065 _id(ctext.getId()),
00066 _fitNumber(PdtPid::nPidType,(int)0),
00067 _defaultType(defaultPart),
00068 _trackT0(t0),
00069 _bField( ctext.bField() )
00070 {
00071
00072 unsigned i=0;
00073 for (TrkRecoTrkImpl::hitListIter iface = _impl->_hitInterfaces.begin();
00074 iface!= _impl->_hitInterfaces.end(); ++iface) {
00075 iface->reset( new TrkHitList(this, (PdtPid::PidType)i++) );
00076 }
00077 }
00078
00079
00080 TrkRecoTrk::TrkRecoTrk(PdtPid::PidType defaultPart,long idnum,double t0) :
00081 _impl(new TrkRecoTrkImpl),
00082 _id(idnum,0),
00083 _fitNumber(PdtPid::nPidType,(int)0),
00084 _defaultType(defaultPart),
00085 _trackT0(t0),
00086 _bField(0)
00087 {
00088
00089 unsigned i=0;
00090 for (TrkRecoTrkImpl::hitListIter iface = _impl->_hitInterfaces.begin();
00091 iface!= _impl->_hitInterfaces.end(); ++iface) {
00092 iface->reset( new TrkHitList(this, (PdtPid::PidType)i++) );
00093 }
00094 }
00095
00096
00097 TrkRecoTrk::TrkRecoTrk(const TrkRecoTrk& rhs) :
00098 _impl(new TrkRecoTrkImpl),
00099 _id(rhs._id.idManager()),
00100 _fitNumber(PdtPid::nPidType,(int)0),
00101 _storage(rhs._storage),
00102 _trackT0( rhs._trackT0),
00103 _bField(rhs._bField)
00104 {
00105 _defaultType = rhs.defaultType();
00106 unsigned i=0;
00107 for (TrkRecoTrkImpl::hitListIter iface = _impl->_hitInterfaces.begin();
00108 iface!= _impl->_hitInterfaces.end(); ++iface) {
00109 iface->reset( new TrkHitList(this, (PdtPid::PidType)i++) );
00110 }
00111 copyReps(rhs);
00112 }
00113
00114 TrkRecoTrk::~TrkRecoTrk()
00115 {
00116 delete _impl;
00117 }
00118
00119 const TrkRecoTrk&
00120 TrkRecoTrk::operator=(const TrkRecoTrk &right)
00121 {
00122 if (&right == this) return *this;
00123 _trackT0 = right._trackT0;
00124 _defaultType = right.defaultType();
00125 copyReps(right);
00126 _bField = right._bField;
00127
00128 _id.setNewValue(right._id);
00129 _storage = right._storage;
00130 return *this;
00131 }
00132
00133 const TrkId&
00134 TrkRecoTrk::id() const
00135 {
00136 return _id;
00137 }
00138
00139 double
00140 TrkRecoTrk::trackT0() const
00141 {
00142 return _trackT0;
00143 }
00144
00145 PdtPid::PidType
00146 TrkRecoTrk::whichFit(PdtPid::PidType hypo) const
00147 {
00148 if (_impl->_reps[hypo].get() == 0) {
00149 return hypo;
00150 }
00151
00152
00153 return _impl->_reps[hypo]->particleType();
00154 }
00155
00156 int
00157 TrkRecoTrk::fitNumber(PdtPid::PidType hypo) const
00158 {
00159 PdtPid::PidType used = whichFit(hypo);
00160 if (used == PdtPid::null) return -1;
00161 int index = used;
00162 return _fitNumber[index];
00163 }
00164
00165 void
00166 TrkRecoTrk::print(ostream& ostr) const
00167 {
00168 ostr << "Trk: " << id() << " def: "
00169 << Pdt::lookup(defaultType())->name()
00170 << " fitNumber:" << fitNumber(defaultType());
00171 const TrkFit* daFit = fitResult();
00172 const TrkHitList* daList = hits();
00173 const TrkFitStatus* daStatus = status();
00174 if (daList != 0) {
00175 ostr << " nhit: " << daList->nHit();
00176 }
00177 if (daFit != 0) {
00178 ostr << " nactive: " << daFit->nActive() << " chisq: " << daFit->chisq();
00179 }
00180 if (daStatus != 0) {
00181 ostr << " 3-d: " << (daStatus->is2d() == 0);
00182 }
00183 ostr << " t0: " << _trackT0 << "\n";
00184 if (daFit != 0) {
00185 TrkExchangePar par = daFit->helix(0.);
00186 ostr << "phi0: " << par.phi0()
00187 << " om: "<< par.omega()
00188 << " d0: " << par.d0() << " z0: " << par.z0()
00189 << " ct: " << par.tanDip();
00190 }
00191 ostr << std::endl;
00192 }
00193
00194 void
00195 TrkRecoTrk::printAll(ostream& ostr) const
00196 {
00197
00198 ostr << "Trk: " << id() << " def: "
00199 << Pdt::lookup(defaultType())->name()
00200 << " fitNumber:" << fitNumber(defaultType());
00201 const TrkFit* daFit = fitResult();
00202 const TrkHitList* daList = hits();
00203 const TrkFitStatus* daStatus = status();
00204 if (daList != 0) {
00205 ostr << " nhit: " << daList->nHit();
00206 }
00207 if (daFit != 0) {
00208 ostr << " nactive: " << daFit->nActive() << " chisq: " << daFit->chisq();
00209 }
00210 if (daStatus != 0) {
00211 ostr << " 3-d: " << (daStatus->is2d() == 0);
00212 }
00213 ostr << " t0: " << _trackT0 << "\n";
00214 getRep(defaultType())->printAll(ostr);
00215
00216 ostr << std::endl;
00217 }
00218
00219 TrkErrCode
00220 TrkRecoTrk::addFit(PdtPid::PidType hypo,bool fit)
00221 {
00222
00223
00224
00225 if (hits() == 0) {
00226
00227 return TrkErrCode(TrkErrCode::fail, 11,
00228 "TrkRecoTrk::addFit(): cannot add a fit to this track.");
00229 }
00230 if (whichFit(hypo) == hypo) {
00231 return TrkErrCode(TrkErrCode::succeed, 11,
00232 "TrkRecoTrk::addFit(): requested fit already exists.");
00233 }
00234 _impl->_reps[hypo].reset( _impl->_reps[defaultType()]->cloneNewHypo(hypo) );
00235 TrkErrCode fitErr(TrkErrCode::succeed, 1);
00236 if (fit && !_impl->_reps[hypo]->fitCurrent()) {
00237 fitErr = _impl->_reps[hypo]->fit();
00238 }
00239 ++_fitNumber[hypo];
00240 return fitErr;
00241 }
00242
00243 void
00244 TrkRecoTrk::resetT0(double t)
00245 {
00246 _trackT0 = t;
00247 updateReps();
00248 }
00249
00250 void
00251 TrkRecoTrk::updateReps()
00252 {
00253 std::pair<TrkRepIter,TrkRepIter> x = uniqueReps();
00254 std::for_each(x.first,x.second,std::mem_fun_ref(&TrkRep::updateHots));
00255 std::transform(_fitNumber.begin(),_fitNumber.end(),
00256 _fitNumber.begin(),
00257 std::bind2nd(std::plus<int>(),1));
00258 }
00259
00260 bool
00261 TrkRecoTrk::operator==(const TrkRecoTrk &other) const
00262 {
00263 return _id == other._id;
00264 }
00265
00266 bool
00267 TrkRecoTrk::operator<(const TrkRecoTrk &other) const
00268 {
00269 return _id < other._id;
00270 }
00271
00272
00273
00274
00275
00276 TrkRep*
00277 TrkRecoTrk::getRep(PdtPid::PidType hypo)
00278 {
00279 assert(hypo>=PdtPid::electron && hypo<= PdtPid::proton);
00280 TrkRep* theRep = _impl->_reps[hypo].get();
00281
00282 if(hypo == defaultType())assert(0 != theRep);
00283 return theRep;
00284 }
00285
00286 const TrkRep*
00287 TrkRecoTrk::getRep(PdtPid::PidType hypo) const
00288 {
00289 assert(hypo>=PdtPid::electron && hypo<= PdtPid::proton);
00290 const TrkRep* theRep = _impl->_reps[hypo].get();
00291 if(hypo == defaultType())assert(0 != theRep);
00292 return theRep;
00293 }
00294
00295 void
00296 TrkRecoTrk::changeDefault(PdtPid::PidType newHypo)
00297 {
00298 if (newHypo == defaultType()) return;
00299 assert(whichFit(newHypo) != PdtPid::null);
00300
00301 TrkHotList *oldList = getRep(defaultType())->hotList();
00302 std::for_each(oldList->begin(),
00303 oldList->end(),
00304 std::mem_fun_ref(&TrkHitOnTrk::setUnusedHit));
00305 assert(getRep(newHypo) != 0);
00306 TrkHotList *newList= getRep(newHypo)->hotList();
00307 std::for_each(newList->begin(),
00308 newList->end(),
00309 std::mem_fun_ref(&TrkHitOnTrk::setUsedHit));
00310 _defaultType = newHypo;
00311 }
00312
00313 void
00314 TrkRecoTrk::copyReps(const TrkRecoTrk& rhs)
00315 {
00316 TrkRecoTrkImpl::repIter lhs = _impl->_reps.begin();
00317 for (TrkRecoTrkImpl::repIter i=rhs._impl->_reps.begin();i!=rhs._impl->_reps.end();++i,++lhs) {
00318 TrkRecoTrkImpl::repIter j=std::find(rhs._impl->_reps.begin(),i,*i);
00319 if (j==i) {
00320 lhs->reset( (*i)->clone(this) );
00321 (*lhs)->setValid((*i)->fitValid());
00322 (*lhs)->setCurrent((*i)->fitCurrent());
00323 } else {
00324 *lhs = *(_impl->_reps.begin()+(j-rhs._impl->_reps.begin()));
00325 }
00326 }
00327 assert(_fitNumber.size()==rhs._fitNumber.size());
00328 std::copy(rhs._fitNumber.begin(),rhs._fitNumber.end(),_fitNumber.begin());
00329 }
00330
00331 void
00332 TrkRecoTrk::setRep(TrkRep* r)
00333 {
00334
00335
00336 std::fill(_impl->_reps.begin(),_impl->_reps.end(),boost::shared_ptr<TrkRep>(r));
00337 std::transform(_fitNumber.begin(),_fitNumber.end(),
00338 _fitNumber.begin(),
00339 std::bind2nd(std::plus<int>(),1));
00340 }
00341
00342 void
00343 TrkRecoTrk::repointHypo(PdtPid::PidType hypo, PdtPid::PidType fit)
00344 {
00345
00346 if (fit == hypo || getRep(fit) == getRep(hypo)) return;
00347
00348 if (hypo == defaultType()) {
00349 std::cout << "ErrMsg(error) "<<
00350 "TrkRecoTrk: can't make default hypothesis point at different fit"
00351 << std::endl;
00352 return;
00353 }
00354 _impl->_reps[hypo] = _impl->_reps[fit];
00355 }
00356
00357 bool
00358 TrkRecoTrk::attach(TrkExtInterface& interface, PdtPid::PidType hypo) const
00359 {
00360 const TrkRep* rp = getRep(hypo);
00361 return rp!=0?interface.attach(rp):0;
00362 }
00363
00364
00365 bool
00366 TrkRecoTrk::attach(TrkExtInterface& interface, PdtPid::PidType hypo)
00367 {
00368 TrkRep* rp = getRep(hypo);
00369 return rp!=0?interface.attach(rp):0;
00370 }
00371
00372 TrkHitList*
00373 TrkRecoTrk::hits(PdtPid::PidType hypo)
00374 {
00375 const TrkRep* rp = getRep(hypo);
00376 return rp==0?0:_impl->_hitInterfaces[hypo].get();
00377 }
00378
00379 const TrkHitList*
00380 TrkRecoTrk::hits(PdtPid::PidType hypo) const
00381 {
00382 const TrkRep* rp = getRep(hypo);
00383 return rp==0?0:_impl->_hitInterfaces[hypo].get();
00384 }
00385
00386 const TrkFit*
00387 TrkRecoTrk::fitResult() const
00388 {
00389 return fitResult(defaultType());
00390 }
00391
00392 const TrkFit*
00393 TrkRecoTrk::fitResult(PdtPid::PidType hypo) const
00394 {
00395 const TrkRep* rp = getRep(hypo);
00396 return rp==0?0:(rp->fitValid()?rp:0);
00397 }
00398
00399 const TrkFitStatus*
00400 TrkRecoTrk::status() const
00401 {
00402 return status(defaultType());
00403 }
00404
00405 const TrkFitStatus*
00406 TrkRecoTrk::status(PdtPid::PidType hypo) const
00407 {
00408 return getRep(hypo);
00409 }
00410
00411 TrkFitStatus*
00412 TrkRecoTrk::status()
00413 {
00414 return status(defaultType());
00415 }
00416
00417 TrkFitStatus*
00418 TrkRecoTrk::status(PdtPid::PidType hypo)
00419 {
00420 return getRep(hypo);
00421 }
00422
00423 void
00424 TrkRecoTrk::setFitNumber(PdtPid::PidType hypo, int newNumber)
00425 {
00426 _fitNumber[hypo] = newNumber;
00427 }
00428
00429 void
00430 TrkRecoTrk::addHypoTo(TrkRep* newRep, PdtPid::PidType hypo)
00431 {
00432 _impl->_reps[hypo].reset( newRep );
00433 }
00434
00435 void
00436 TrkRecoTrk::setIdManager(TrkIdManager* idMan)
00437 {
00438 _id.setIdManager(idMan);
00439 }
00440
00441 ostream&
00442 operator<<(ostream& os, const TrkRecoTrk& tk)
00443 {
00444 tk.print(os);
00445 return os;
00446 }
00447
00448 void
00449 TrkRecoTrk::setBField(const BField* field)
00450 {
00451 _bField = field;
00452 }
00453
00454 void
00455 TrkRecoTrk::markForStore(PdtPid::PidType hypo,double fltlen,
00456 const char* listname) {
00457
00458 PdtPid::PidType realhypo = whichFit(hypo);
00459
00460 if(getRep(realhypo)!= 0 && getRep(realhypo)->fitValid())
00461
00462 _storage[std::string(listname)].insert(TrkStoreHypo(realhypo,fltlen));
00463 else
00464
00465 std::cout<<"ErrMsg(error) "<< "Invalid fits cannot be marked for storage" << std::endl;
00466 }
00467
00468 const std::set<TrkStoreHypo>&
00469 TrkRecoTrk::storageRequests(const char* listname) const {
00470 static std::set<TrkStoreHypo> empty;
00471 std::map<std::string,std::set<TrkStoreHypo> >::const_iterator foundit =
00472 _storage.find(std::string(listname));
00473 if(foundit != _storage.end())
00474 return foundit->second;
00475 else
00476 return empty;
00477 }
00478
00479 void
00480 TrkRecoTrk::clearStorageRequests(const char* listname) {
00481 _storage[std::string(listname)].clear();
00482 }
00483
00484 void
00485 TrkRecoTrk::storageLists(std::set<std::string>& storage) const {
00486
00487 storage.clear();
00488
00489 std::map<std::string,std::set<TrkStoreHypo> >::const_iterator miter = _storage.begin();
00490 while(miter != _storage.end()){
00491 storage.insert(miter->first);
00492 miter++;
00493 }
00494 }
00495
00496 TrkHotList*
00497 TrkRecoTrk::hots(PdtPid::PidType hypo)
00498 {
00499 TrkRep* rp = getRep(hypo);
00500 return rp==0?0:rp->hotList();
00501 }
00502
00503 const TrkHotList*
00504 TrkRecoTrk::hots(PdtPid::PidType hypo) const
00505 {
00506 const TrkRep* rp = getRep(hypo);
00507 return rp==0?0:rp->hotList();
00508 }
00509
00510 std::pair<TrkRepIter,TrkRepIter>
00511 TrkRecoTrk::allReps() const
00512 {
00513 typedef std::vector<TrkRep*> RPL;
00514 boost::shared_ptr<RPL> x( new RPL );
00515
00516
00517 for (TrkRecoTrkImpl::repConstIter i=_impl->_reps.begin();i!=_impl->_reps.end();++i) {
00518 x->push_back(i->get());
00519
00520
00521
00522
00523 }
00524
00525 return std::make_pair(TrkRepIter(x,0),TrkRepIter(x,x->size()));
00526 }
00527
00528 std::pair<TrkRepIter,TrkRepIter>
00529 TrkRecoTrk::uniqueReps() const
00530 {
00531 typedef std::vector<TrkRep*> RPL;
00532 boost::shared_ptr<RPL> x( new RPL );
00533 for (TrkRecoTrkImpl::repConstIter i=_impl->_reps.begin();i!=_impl->_reps.end();++i) {
00534 if (std::find(x->begin(),x->end(),i->get())==x->end()) x->push_back(i->get());
00535
00536 }
00537 return std::make_pair(TrkRepIter(x,0),TrkRepIter(x,x->size()));
00538 }