TrkRecoTrk.cxx

Go to the documentation of this file.

// File and Version Information:
//      $Id: TrkRecoTrk.cxx,v 1.9 2008/09/23 00:58:03 zhangy Exp $
//
// Description:
//      implementation for TrkRecoTrk
//
//
// Author List: Stephen Schaffner
//
//------------------------------------------------------------------------
#include <vector>
#include <assert.h>
#include <functional>
#include <algorithm>
#include "MdcGeom/Constants.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TrkBase/TrkHitOnTrk.h"
#include "TrkBase/TrkFundHit.h"
#include "TrkBase/TrkDifTraj.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkBase/TrkErrCode.h"
#include "TrkBase/TrkExtInterface.h"
#include "TrkBase/TrkFitStatus.h"
#include "TrkBase/TrkRepIter.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "MdcRecoUtil/Pdt.h"
#include "MdcRecoUtil/BesPointErr.h"
#include "MdcRecoUtil/BesVectorErr.h"
#include "MdcRecoUtil/DifPoint.h"
#include "MdcRecoUtil/DifVector.h"
#include "MdcRecoUtil/DifIndepPar.h"
#include "TrkBase/TrkRep.h"
#include "TrkBase/TrkContext.h"
#include "MdcRecoUtil/PdtPid.h"
#include "boost/shared_ptr.hpp"
#include "BField/BField.h"
using std::ostream;
 
 
class TrkRecoTrkImpl {
public:
  friend class TrkRecoTrk;
private:
  TrkRecoTrkImpl()
     : _reps(PdtPid::nPidType,boost::shared_ptr<TrkRep>((TrkRep *)0)),
       _hitInterfaces(PdtPid::nPidType,boost::shared_ptr<TrkHitList>((TrkHitList *)0))
  { }
 
  // the reps live here; owned by trk;
  // Pointers to the Reps, indexed by particle type; always has nPart entries,
  // which may point to as few as one single distinct Rep...
  typedef std::vector<boost::shared_ptr<TrkRep> > repList_t;
  typedef repList_t::const_iterator repConstIter;
  typedef repList_t::iterator repIter;
  repList_t  _reps;
  typedef std::vector<boost::shared_ptr<TrkHitList> > hitList_t;
  typedef hitList_t::iterator hitListIter;
  hitList_t _hitInterfaces;
 
};
 
TrkRecoTrk::TrkRecoTrk(PdtPid::PidType defaultPart, const TrkContext& ctext,
                       double t0) :
  _impl(new TrkRecoTrkImpl), 
  _id(ctext.getId()),
  _fitNumber(PdtPid::nPidType,(int)0),
  _defaultType(defaultPart),
  _trackT0(t0),
  _bField( ctext.bField() )
{
  // No TrkRep is defined here; must be created in appropriate FitMaker.
  unsigned i=0;
  for (TrkRecoTrkImpl::hitListIter iface = _impl->_hitInterfaces.begin();
       iface!= _impl->_hitInterfaces.end(); ++iface) {
    iface->reset( new TrkHitList(this, (PdtPid::PidType)i++) );  //cast
  }
}
 
// persistence reconstitution.  This sets a nul value for BField and IdManager
TrkRecoTrk::TrkRecoTrk(PdtPid::PidType defaultPart,long idnum,double t0) :
  _impl(new TrkRecoTrkImpl),
  _id(idnum,0),
  _fitNumber(PdtPid::nPidType,(int)0),
  _defaultType(defaultPart),
  _trackT0(t0),
  _bField(0)
{
  // No TrkRep is defined here; must be created in appropriate FitMaker.
  unsigned i=0;
  for (TrkRecoTrkImpl::hitListIter iface = _impl->_hitInterfaces.begin();
       iface!= _impl->_hitInterfaces.end(); ++iface) {
    iface->reset( new TrkHitList(this, (PdtPid::PidType)i++) );  //cast
  }
}
 
//-- Copy constructor
TrkRecoTrk::TrkRecoTrk(const TrkRecoTrk& rhs) :
  _impl(new TrkRecoTrkImpl),
  _id(rhs._id.idManager()),
  _fitNumber(PdtPid::nPidType,(int)0),
  _storage(rhs._storage),
  _trackT0( rhs._trackT0),
  _bField(rhs._bField)
{
  _defaultType = rhs.defaultType();
  unsigned i=0;
  for (TrkRecoTrkImpl::hitListIter iface = _impl->_hitInterfaces.begin();
       iface!= _impl->_hitInterfaces.end(); ++iface) {
    iface->reset( new TrkHitList(this, (PdtPid::PidType)i++) );  //cast
  }
  copyReps(rhs);
}
 
TrkRecoTrk::~TrkRecoTrk()
{
  delete _impl;
}
 
const TrkRecoTrk&
TrkRecoTrk::operator=(const TrkRecoTrk &right)
{
  if (&right == this) return *this;
  _trackT0 = right._trackT0;
  _defaultType = right.defaultType();
  copyReps(right);
  _bField = right._bField;
//  AbsEvtObj::operator=(right);
  _id.setNewValue(right._id);
  _storage = right._storage;
  return *this;
}
 
const TrkId&
TrkRecoTrk::id() const
{
  return _id;
}
 
double
TrkRecoTrk::trackT0() const
{
  return _trackT0;
}
 
PdtPid::PidType
TrkRecoTrk::whichFit(PdtPid::PidType hypo) const
{
  if (_impl->_reps[hypo].get() == 0) {
    return hypo;
  }
//  return _repPtr[hypo]->fitValid() ? _repPtr[hypo]->particleType()
//                                   : PdtPid::null;
  return _impl->_reps[hypo]->particleType();
}
 
int
TrkRecoTrk::fitNumber(PdtPid::PidType hypo) const
{
  PdtPid::PidType used = whichFit(hypo);
  if (used == PdtPid::null) return -1;
  int index = used;
  return _fitNumber[index];
}
 
void
TrkRecoTrk::print(ostream& ostr) const
{
  ostr << "Trk: " << id() << " def: "
       << Pdt::lookup(defaultType())->name()
       << " fitNumber:" << fitNumber(defaultType());
  const TrkFit* daFit = fitResult();
  const TrkHitList* daList = hits();
  const TrkFitStatus* daStatus = status();
  if (daList != 0) {
    ostr << " nhit: " << daList->nHit();
  }
  if (daFit != 0) {
    ostr << " nactive: " << daFit->nActive() << " chisq: " << daFit->chisq();
  }
  if (daStatus != 0) {
    ostr << " 3-d: " << (daStatus->is2d() == 0);
  }
  ostr << " t0: " << _trackT0 << "\n";
  if (daFit != 0) {
    TrkExchangePar par = daFit->helix(0.);
    ostr << "phi0: " << par.phi0()
       << " om: "<< par.omega()
       << " d0: " << par.d0() << " z0: " << par.z0()
       << " ct: " << par.tanDip();
  }
  ostr << std::endl;
}
 
void
TrkRecoTrk::printAll(ostream& ostr) const
{
  //  This should be expanded to print other hypotheses as well
  ostr << "Trk: " << id() << " def: "
       << Pdt::lookup(defaultType())->name()
       << " fitNumber:" << fitNumber(defaultType());
  const TrkFit* daFit = fitResult();
  const TrkHitList* daList = hits();
  const TrkFitStatus* daStatus = status();
  if (daList != 0) {
    ostr << " nhit: " << daList->nHit();
  }
  if (daFit != 0) {
    ostr << " nactive: " << daFit->nActive() << " chisq: " << daFit->chisq();
  }
  if (daStatus != 0) {
    ostr << " 3-d: " << (daStatus->is2d() == 0);
  }
  ostr << " t0: " << _trackT0 << "\n";
  getRep(defaultType())->printAll(ostr);
  //  allReps();//yzhang debug
  ostr << std::endl;
}
 
TrkErrCode
TrkRecoTrk::addFit(PdtPid::PidType hypo,bool fit)
{
  // If there is no fit, create one.  If hypo points to a fit for a different
  //   particle type, create a fit of type hypo, and point at that.  Carry
  //   out the fit if needed.
  if (hits() == 0) {
    // Unfittable rep
    return TrkErrCode(TrkErrCode::fail, 11,
                      "TrkRecoTrk::addFit(): cannot add a fit to this track.");
  }
  if (whichFit(hypo) == hypo) {
    return TrkErrCode(TrkErrCode::succeed, 11,
                      "TrkRecoTrk::addFit(): requested fit already exists.");
  }
  _impl->_reps[hypo].reset( _impl->_reps[defaultType()]->cloneNewHypo(hypo) );
  TrkErrCode fitErr(TrkErrCode::succeed, 1);
  if (fit && !_impl->_reps[hypo]->fitCurrent()) {
    fitErr = _impl->_reps[hypo]->fit();
  }
  ++_fitNumber[hypo];
  return fitErr;
}
 
void
TrkRecoTrk::resetT0(double t)
{
  _trackT0 = t;
  updateReps();
}
 
void
TrkRecoTrk::updateReps()
{
  std::pair<TrkRepIter,TrkRepIter> x = uniqueReps();
  std::for_each(x.first,x.second,std::mem_fun_ref(&TrkRep::updateHots));
  std::transform(_fitNumber.begin(),_fitNumber.end(),
                 _fitNumber.begin(),
                 std::bind2nd(std::plus<int>(),1));
}
 
bool
TrkRecoTrk::operator==(const TrkRecoTrk &other) const
{
  return _id == other._id;
}
 
bool
TrkRecoTrk::operator<(const TrkRecoTrk &other) const
{
  return _id < other._id;
}
 
//***********************
// Protected functions
//***********************
 
TrkRep*
TrkRecoTrk::getRep(PdtPid::PidType hypo)
{
  assert(hypo>=PdtPid::electron && hypo<= PdtPid::proton);
  TrkRep* theRep = _impl->_reps[hypo].get();
// insist the default rep exist
  if(hypo == defaultType())assert(0 != theRep);
  return theRep;
}
 
const TrkRep*
TrkRecoTrk::getRep(PdtPid::PidType hypo) const
{
  assert(hypo>=PdtPid::electron && hypo<= PdtPid::proton);
  const TrkRep* theRep = _impl->_reps[hypo].get();
  if(hypo == defaultType())assert(0 != theRep);
  return theRep;
}
 
void
TrkRecoTrk::changeDefault(PdtPid::PidType newHypo)
{
  if (newHypo == defaultType()) return;
  assert(whichFit(newHypo) != PdtPid::null);
 
  TrkHotList *oldList = getRep(defaultType())->hotList();
  std::for_each(oldList->begin(),
                oldList->end(),
                std::mem_fun_ref(&TrkHitOnTrk::setUnusedHit));
  assert(getRep(newHypo) != 0);
  TrkHotList *newList= getRep(newHypo)->hotList();
  std::for_each(newList->begin(),
                newList->end(),
                std::mem_fun_ref(&TrkHitOnTrk::setUsedHit));
  _defaultType = newHypo;
}
 
void
TrkRecoTrk::copyReps(const TrkRecoTrk& rhs)
{
  TrkRecoTrkImpl::repIter lhs = _impl->_reps.begin();
  for (TrkRecoTrkImpl::repIter i=rhs._impl->_reps.begin();i!=rhs._impl->_reps.end();++i,++lhs) {
      TrkRecoTrkImpl::repIter j=std::find(rhs._impl->_reps.begin(),i,*i);
      if (j==i) { // first time this one is seen
          lhs->reset( (*i)->clone(this) );
          (*lhs)->setValid((*i)->fitValid());
          (*lhs)->setCurrent((*i)->fitCurrent());
      } else {
          *lhs = *(_impl->_reps.begin()+(j-rhs._impl->_reps.begin()));
      }
  }
  assert(_fitNumber.size()==rhs._fitNumber.size());
  std::copy(rhs._fitNumber.begin(),rhs._fitNumber.end(),_fitNumber.begin());
}
 
void
TrkRecoTrk::setRep(TrkRep* r)
{
  // Sets the default rep to be r, clears out other reps., and sets
  //  non-default rep ptrs to point at default.  Increments all fit numbers.
  std::fill(_impl->_reps.begin(),_impl->_reps.end(),boost::shared_ptr<TrkRep>(r));
  std::transform(_fitNumber.begin(),_fitNumber.end(),
                 _fitNumber.begin(),
                 std::bind2nd(std::plus<int>(),1));
}
 
void
TrkRecoTrk::repointHypo(PdtPid::PidType hypo, PdtPid::PidType fit)
{
  // Do we have to do anything?
  if (fit == hypo || getRep(fit) == getRep(hypo)) return;
 
  if (hypo == defaultType()) {
    std::cout << "ErrMsg(error) "<<
      "TrkRecoTrk: can't make default hypothesis point at different fit"
                  << std::endl;
    return;
  }
  _impl->_reps[hypo] = _impl->_reps[fit];
}
 
bool
TrkRecoTrk::attach(TrkExtInterface& interface, PdtPid::PidType hypo) const
{
  const TrkRep* rp = getRep(hypo);
  return rp!=0?interface.attach(rp):0;
}
 
 
bool
TrkRecoTrk::attach(TrkExtInterface& interface, PdtPid::PidType hypo)
{
  TrkRep* rp = getRep(hypo);
  return rp!=0?interface.attach(rp):0;
}
 
TrkHitList*
TrkRecoTrk::hits(PdtPid::PidType hypo)
{
  const TrkRep* rp = getRep(hypo);
  return rp==0?0:_impl->_hitInterfaces[hypo].get();
}
 
const TrkHitList*
TrkRecoTrk::hits(PdtPid::PidType hypo) const
{
  const TrkRep* rp = getRep(hypo);
  return rp==0?0:_impl->_hitInterfaces[hypo].get();
}
 
const TrkFit*
TrkRecoTrk::fitResult() const
{
  return fitResult(defaultType());
}
 
const TrkFit*
TrkRecoTrk::fitResult(PdtPid::PidType hypo) const
{
  const TrkRep* rp = getRep(hypo);
  return rp==0?0:(rp->fitValid()?rp:0);
}
 
const TrkFitStatus*
TrkRecoTrk::status() const
{
  return status(defaultType());
}
 
const TrkFitStatus*
TrkRecoTrk::status(PdtPid::PidType hypo) const
{
  return getRep(hypo);
}
 
TrkFitStatus*
TrkRecoTrk::status()
{
  return status(defaultType());
}
 
TrkFitStatus*
TrkRecoTrk::status(PdtPid::PidType hypo)
{
  return getRep(hypo);
}
 
void
TrkRecoTrk::setFitNumber(PdtPid::PidType hypo, int newNumber)
{
  _fitNumber[hypo] = newNumber;
}
 
void
TrkRecoTrk::addHypoTo(TrkRep* newRep, PdtPid::PidType hypo)
{
  _impl->_reps[hypo].reset( newRep );
}
 
void
TrkRecoTrk::setIdManager(TrkIdManager* idMan)
{
  _id.setIdManager(idMan);
}
 
ostream&
operator<<(ostream& os, const TrkRecoTrk& tk)
{
  tk.print(os);
  return os;
}
 
void
TrkRecoTrk::setBField(const BField* field)
{
  _bField = field;
}
 
void
TrkRecoTrk::markForStore(PdtPid::PidType hypo,double fltlen,
             const char* listname) {
// first, translate to the real hypo
  PdtPid::PidType realhypo = whichFit(hypo);
// Then, make sure this hypo has a valid fit
  if(getRep(realhypo)!= 0 && getRep(realhypo)->fitValid())
// add an entry in the toStore list (if it's unique)
    _storage[std::string(listname)].insert(TrkStoreHypo(realhypo,fltlen));
  else
// It's an error to try to store invalid fits
    std::cout<<"ErrMsg(error) "<< "Invalid fits cannot be marked for storage" << std::endl;
}
 
const std::set<TrkStoreHypo>&
TrkRecoTrk::storageRequests(const char* listname) const {
  static std::set<TrkStoreHypo> empty; // empty set to return if list doesn't exist
  std::map<std::string,std::set<TrkStoreHypo> >::const_iterator foundit =
    _storage.find(std::string(listname));
  if(foundit != _storage.end())
    return foundit->second;
  else
    return empty;
}
 
void
TrkRecoTrk::clearStorageRequests(const char* listname) {
  _storage[std::string(listname)].clear();
}
 
void
TrkRecoTrk::storageLists(std::set<std::string>& storage) const {
  // clear the output set
  storage.clear();
  // iterate over all the storage requests
  std::map<std::string,std::set<TrkStoreHypo> >::const_iterator miter = _storage.begin();
  while(miter != _storage.end()){
    storage.insert(miter->first);
    miter++;
  }
}
 
TrkHotList*
TrkRecoTrk::hots(PdtPid::PidType hypo)
{
  TrkRep* rp = getRep(hypo);
  return rp==0?0:rp->hotList();
}
 
const TrkHotList*
TrkRecoTrk::hots(PdtPid::PidType hypo) const
{
  const TrkRep* rp = getRep(hypo);
  return rp==0?0:rp->hotList();
}
 
std::pair<TrkRepIter,TrkRepIter>
TrkRecoTrk::allReps() const
{
  typedef std::vector<TrkRep*> RPL;
  boost::shared_ptr<RPL> x( new RPL );
  //std::cout << " TrkRecoTrk::allReps" << std::endl;//yzhang debug
        
  for (TrkRecoTrkImpl::repConstIter i=_impl->_reps.begin();i!=_impl->_reps.end();++i) {
          x->push_back(i->get());
      /*
      i->get()->printType(std::cout);//yzhang debug
      i->get()->printAll(std::cout);//yzhang debug
      */
  }
  //std::cout << "------ " << std::endl;//yzhang debug
  return std::make_pair(TrkRepIter(x,0),TrkRepIter(x,x->size()));
}
 
std::pair<TrkRepIter,TrkRepIter>
TrkRecoTrk::uniqueReps() const
{
  typedef std::vector<TrkRep*> RPL;
  boost::shared_ptr<RPL> x( new RPL );
  for (TrkRecoTrkImpl::repConstIter i=_impl->_reps.begin();i!=_impl->_reps.end();++i) {
     if (std::find(x->begin(),x->end(),i->get())==x->end()) x->push_back(i->get());
      
  }
  return std::make_pair(TrkRepIter(x,0),TrkRepIter(x,x->size()));
}