CMS 3D CMS Logo

Classes | Public Member Functions | Private Types | Private Member Functions | Private Attributes

MuonRoadTrajectoryBuilder Class Reference

#include <MuonRoadTrajectoryBuilder.h>

Inheritance diagram for MuonRoadTrajectoryBuilder:
TrajectoryBuilder

List of all members.

Classes

class  flippingPair
class  trajectory
 for the trajectory collection More...

Public Member Functions

 MuonRoadTrajectoryBuilder (const edm::ParameterSet &par, const MeasurementTracker *mt, const MagneticField *f, const Propagator *p)
 constructor from PSet and things from record
void setEvent (const edm::Event &) const
std::vector< Trajectorytrajectories (const TrajectorySeed &seed) const
void trajectories (const TrajectorySeed &seed, TrajectoryContainer &ret) const
 process the seed, in a faster manner.
 ~MuonRoadTrajectoryBuilder ()

Private Types

typedef std::list< trajectoryTrajectoryCollection
typedef flippingPair
< TrajectoryCollection
TrajectoryCollectionFPair

Private Member Functions

void checkDuplicate (TrajectoryCollection &collection) const
bool checkStep (TrajectoryCollection &collection) const
void cleanTrajectory (Trajectory &traj) const
int GatherHits (const TrajectoryStateOnSurface &step, const DetLayer *thislayer, TrajectoryCollectionFPair &Trajectories) const
void makeTrajectories (const TrajectorySeed &seed, std::vector< Trajectory > &result, int version=0) const
void makeTrajectories_0 (const TrajectorySeed &seed, std::vector< Trajectory > &result) const
void makeTrajectories_1 (const TrajectorySeed &seed, std::vector< Trajectory > &result) const
Trajectory smooth (Trajectory &) const

Private Attributes

bool theBranchonfirstlayer
bool theCarriedIPatfirstlayer
bool theCarriedIPatfirstlayerModule
std::string theCategory
 Info/Debug category "Muon|RecoMuon|MuonRoadTrajectoryBuilder".
bool theDynamicMaxNumberOfHitPerModule
const MagneticFieldtheField
bool theFirstlayer
Chi2MeasurementEstimatortheHitEstimator
unsigned int theMaxTrajectories
std::vector< unsigned int > theMaxTrajectoriesThreshold
const MeasurementTrackertheMeasurementTracker
unsigned int theMinNumberOfHitOnCandidate
unsigned int theNumberOfHitPerModule
unsigned int theNumberOfHitPerModuleDefault
std::vector< unsigned int > theNumberOfHitPerModuleThreshold
bool theOutputAllTraj
const PropagatorthePropagator
Chi2MeasurementEstimatortheRoadEstimator
KFTrajectorySmoothertheSmoother
TrajectoryStateTransform theTransformer
TrajectoryStateUpdatortheUpdator

Detailed Description

Definition at line 55 of file MuonRoadTrajectoryBuilder.h.


Member Typedef Documentation

Definition at line 103 of file MuonRoadTrajectoryBuilder.h.

Definition at line 104 of file MuonRoadTrajectoryBuilder.h.


Constructor & Destructor Documentation

MuonRoadTrajectoryBuilder::MuonRoadTrajectoryBuilder ( const edm::ParameterSet par,
const MeasurementTracker mt,
const MagneticField f,
const Propagator p 
)

constructor from PSet and things from record

Definition at line 40 of file MuonRoadTrajectoryBuilder.cc.

References Chi2MeasurementEstimatorESProducer_cfi::Chi2MeasurementEstimator, f, edm::ParameterSet::getParameter(), KFTrajectorySmootherESProducer_cfi::KFTrajectorySmoother, L1TEmulatorMonitor_cff::p, mathSSE::sqrt(), theBranchonfirstlayer, theCarriedIPatfirstlayer, theCarriedIPatfirstlayerModule, theCategory, theDynamicMaxNumberOfHitPerModule, theField, theHitEstimator, theMaxTrajectories, theMaxTrajectoriesThreshold, theMeasurementTracker, theMinNumberOfHitOnCandidate, theNumberOfHitPerModuleDefault, theNumberOfHitPerModuleThreshold, theOutputAllTraj, thePropagator, theRoadEstimator, theSmoother, and theUpdator.

                                                                      :
  theRoadEstimator(0),theHitEstimator(0),theUpdator(0),theSmoother(0)
{

  theMeasurementTracker = mt;
  theField = f;
  thePropagator = p;

  theCategory = "Muon|RecoMuon|MuonRoadTrajectoryBuilder";

  //get parameters from ParameterSet

  //propagator name to get from muon service
  //  thePropagatorName = iConfig.getParameter<std::string>("propagatorName");

  //chi2 estimator
  double chi2R=iConfig.getParameter<double>("maxChi2Road");
  theRoadEstimator = new Chi2MeasurementEstimator(chi2R,sqrt(chi2R));
  double chi2H=iConfig.getParameter<double>("maxChi2Hit");
  theHitEstimator = new Chi2MeasurementEstimator(chi2H,sqrt(chi2H));

  //trajectory updator
  theUpdator= new KFUpdator();

  //limit the total number of possible trajectories taken into account for a single seed
  theMaxTrajectories = iConfig.getParameter<unsigned int>("maxTrajectories");

  //limit the type of module considered to gather rechits
  theDynamicMaxNumberOfHitPerModule = iConfig.getParameter<bool>("dynamicMaxNumberOfHitPerModule");
  theNumberOfHitPerModuleDefault = iConfig.getParameter<unsigned int>("numberOfHitPerModule");
  theMaxTrajectoriesThreshold = iConfig.getParameter<std::vector<unsigned int> >("maxTrajectoriesThreshold");
  theNumberOfHitPerModuleThreshold = iConfig.getParameter<std::vector<unsigned int> >("numberOfHitPerModuleThreshold");

  //could be configurable, but not
  theBranchonfirstlayer=true;
  theCarriedIPatfirstlayer=false;
  theCarriedIPatfirstlayerModule =true;

  //output track candidate selection
  theMinNumberOfHitOnCandidate = iConfig.getParameter<unsigned int>("minNumberOfHitOnCandidate");
  
  //single or multiple trajectories per muon
  theOutputAllTraj = iConfig.getParameter<bool>("outputAllTraj");

  if (!theSmoother)
    theSmoother = new KFTrajectorySmoother(thePropagator,theUpdator,theRoadEstimator);
}
MuonRoadTrajectoryBuilder::~MuonRoadTrajectoryBuilder ( )

Definition at line 92 of file MuonRoadTrajectoryBuilder.cc.

References theCategory, theHitEstimator, theRoadEstimator, theSmoother, and theUpdator.

{
  edm::LogInfo(theCategory)<<"cleaning the object";
  if (theRoadEstimator) delete theRoadEstimator;
  if (theHitEstimator) delete theHitEstimator;
  if (theUpdator) delete theUpdator;
  if (theSmoother) delete theSmoother;
}

Member Function Documentation

void MuonRoadTrajectoryBuilder::checkDuplicate ( TrajectoryCollection collection) const [private]

Definition at line 688 of file MuonRoadTrajectoryBuilder.cc.

References LogDebug, and theCategory.

Referenced by GatherHits().

{
  LogDebug(theCategory)<<"checking for duplicates here. list size: "<<collection.size();

  //loop over the trajectory list
  TrajectoryCollection::iterator traj1 = collection.begin();
  while(traj1 != collection.end())
    {
      LogDebug(theCategory)<<"";

      //reloop over the trajectory list from traj1
      TrajectoryCollection::iterator traj2 = traj1;
      traj2++;//advance one more
      bool traj1_removed=false;
      while( traj2 != collection.end())
        {
          if (traj2 == traj1 ) continue; //skip itself of course

          //need to start from the back of the list of measurment
          std::list <TrajectoryMeasurement >::reverse_iterator h1 = traj1->measurements.rbegin();
          std::list <TrajectoryMeasurement >::reverse_iterator h2 = traj2->measurements.rbegin();
                  
          bool break_different = false;
          while (h1 != traj1->measurements.rend() && h2!=traj2->measurements.rend())
            {
              TransientTrackingRecHit::ConstRecHitPointer hit1 = h1->recHit();
              TransientTrackingRecHit::ConstRecHitPointer hit2 = h2->recHit();

              LogDebug(theCategory)<<"(hit1->geographicalId().rawId()) "<<(hit1->geographicalId().rawId())<<"(hit1->globalPosition()) "<<(hit1->globalPosition())
                                   <<"(hit2->geographicalId().rawId()) "<<(hit2->geographicalId().rawId())<<"(hit2->globalPosition()) "<<(hit2->globalPosition());

              if (hit1 == hit2){/*those are common hits, everything's alright so far*/ h1++;h2++; continue;}
              else{break_different =true;
                        
                LogDebug(theCategory)<<"list of hits are different";
                        
                break;}
            }
          if (!break_different) 
            //meaning one of the list has been exhausted
            //one list if the subset of the other
            {
              LogDebug(theCategory)<<"list of hits are identical/subset. remove one of them.";
              //there is a common subset to the two list of rechits.
              //remove the one with the fewer hits (iterator that reached the end first)
              //in case they are exactly identical (both iterator reached the end at the same time), traj2 will be removed by default
              if (h1 != traj1->measurements.rend())
                {
                  LogDebug(theCategory)<<"I removed traj2";
                  //traj2 has been exhausted first. remove it and place the iterator on next item
                  traj2=collection.erase(traj2);
                }
              else
                {
                  LogDebug(theCategory)<<"I removed traj1. and decrement traj2";
                  //traj1 has been exhausted first. remove it
                  traj1=collection.erase(traj1);
                  //and set the iterator traj1 so that next++ will set it to the correct place in the list  
                  traj1_removed=true;
                  break; // break the traj2 loop, advance to next traj1 item
                }
            }
          else
            {traj2++;   }
        }//loop on traj2
      if (!traj1_removed)
        {//increment only if you have remove the item at traj1
          traj1++;}
    }//loop on traj1
}
bool MuonRoadTrajectoryBuilder::checkStep ( TrajectoryCollection collection) const [private]

Definition at line 665 of file MuonRoadTrajectoryBuilder.cc.

References theCategory, theDynamicMaxNumberOfHitPerModule, theMaxTrajectories, theMaxTrajectoriesThreshold, theNumberOfHitPerModule, theNumberOfHitPerModuleThreshold, and trajectoryOrder().

Referenced by makeTrajectories_0().

{
  //dynamic cut on the max number of rechit allowed on a single module
  if (theDynamicMaxNumberOfHitPerModule) {
    for (unsigned int vit = 0; vit!= theMaxTrajectoriesThreshold.size() ; vit++){
      if (collection.size() >theMaxTrajectoriesThreshold[vit]){
        theNumberOfHitPerModule= theNumberOfHitPerModuleThreshold[vit];}
      else
        break;}}
  
  //reduce the number of possible trajectories if too many
  if ( collection.size() > theMaxTrajectories) {
    //order with most hits or best chi2
    collection.sort(trajectoryOrder);
    unsigned int prevSize=collection.size();
    collection.resize(theMaxTrajectories);
    edm::LogInfo(theCategory)
      <<" too many possible trajectories ("<<prevSize
      <<"), reduce the possibilities to "<<theMaxTrajectories<<" bests.";
  }
  return true;
}
void MuonRoadTrajectoryBuilder::cleanTrajectory ( Trajectory traj) const [private]

Definition at line 394 of file MuonRoadTrajectoryBuilder.cc.

References GeometricSearchTracker::detLayer(), Trajectory::direction(), MeasurementTracker::geometricSearchTracker(), LogDebug, Trajectory::measurements(), Trajectory::push(), Trajectory::seed(), theCategory, and theMeasurementTracker.

Referenced by smooth().

                                                                        {
   //remove the overlapping recHits since the smoother will chock on it
   Trajectory::DataContainer meas =traj.measurements();
   
   const DetLayer* lastDetLayer =0;
   Trajectory::DataContainer::iterator mit = meas.begin();
   while (mit!=meas.end()){
     {
       const DetLayer * detLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(mit->recHit()->det()->geographicalId());
       LogDebug(theCategory)<<"(mit->recHit()->det()->geographicalId().rawId()) "<<(mit->recHit()->det()->geographicalId().rawId())
                            <<"(detLayer) "<<(detLayer)<<"(lastDetLayer) "<<(lastDetLayer);

       if (detLayer==lastDetLayer)
         {
           mit=meas.erase(mit); //serve as mit++ too
         }
       else {mit++;}
       lastDetLayer=detLayer;
     }
     Trajectory newTraj(traj.seed(),traj.direction());
     for (Trajectory::DataContainer::iterator mit=meas.begin();mit!=meas.end();++mit)
       {newTraj.push(*mit);}
     traj=newTraj;
   }
 }
int MuonRoadTrajectoryBuilder::GatherHits ( const TrajectoryStateOnSurface step,
const DetLayer thislayer,
TrajectoryCollectionFPair Trajectories 
) const [private]

Definition at line 459 of file MuonRoadTrajectoryBuilder.cc.

References checkDuplicate(), MuonRoadTrajectoryBuilder::trajectory::chi2, runregparse::combined, GeometricSearchDet::compatibleDets(), MuonRoadTrajectoryBuilder::trajectory::duplicate, Chi2MeasurementEstimator::estimate(), MuonRoadTrajectoryBuilder::flippingPair< A >::flip(), Trajectory::foundHits(), TrajectoryStateOnSurface::freeState(), TrajectoryStateOnSurface::globalPosition(), MuonRoadTrajectoryBuilder::flippingPair< A >::head(), MeasurementTracker::idToDet(), TrajectoryStateOnSurface::isValid(), MuonRoadTrajectoryBuilder::trajectory::lastmissed, LogDebug, PV3DBase< T, PVType, FrameType >::mag(), MuonRoadTrajectoryBuilder::trajectory::measurements, MuonRoadTrajectoryBuilder::trajectory::missed, MuonRoadTrajectoryBuilder::trajectory::missedinarow, PV3DBase< T, PVType, FrameType >::perp(), Propagator::propagate(), Trajectory::push(), DetId::rawId(), MeasurementDet::recHits(), MuonRoadTrajectoryBuilder::flippingPair< A >::tail(), theBranchonfirstlayer, theCarriedIPatfirstlayer, theCarriedIPatfirstlayerModule, theCategory, theFirstlayer, theHitEstimator, theMeasurementTracker, theNumberOfHitPerModule, thePropagator, theRoadEstimator, theUpdator, MuonRoadTrajectoryBuilder::trajectory::traj, MuonRoadTrajectoryBuilder::trajectory::TSOS, and TrajectoryStateUpdator::update().

Referenced by makeTrajectories_0().

{
  TrajectoryStateOnSurface restep;
  bool atleastoneadded=false;
  int Nhits=0;

  //find compatible modules
  std::vector<DetLayer::DetWithState>  compatible =thislayer->compatibleDets( step, *thePropagator , *theRoadEstimator);
  
  //loop over compatible modules
  for (std::vector< DetLayer::DetWithState > ::iterator dws = compatible.begin();dws != compatible.end();dws++)
    {
      const DetId presentdetid = dws->first->geographicalId();//get the det Id
      restep = dws->second;


        LogDebug(theCategory)<<((dws->first->components ().size()!=0) ? /*stereo layer*/"double sided layer":/*single sided*/"single sided layer")
                             <<"(presentdetid.rawId()) "<<(presentdetid.rawId());

      //get the rechits on this module
      TransientTrackingRecHit::ConstRecHitContainer  thoseHits = theMeasurementTracker->idToDet(presentdetid)->recHits(restep); 
      int additionalHits =0;

      LogDebug(theCategory)<<"(thoseHits.size()) "<<(thoseHits.size());

      if (thoseHits.size()>theNumberOfHitPerModule)
        { edm::LogInfo(theCategory)<<"more than "<<theNumberOfHitPerModule 
                                           <<" compatible hits ("<<thoseHits.size()<<")on module "<<presentdetid.rawId()
                                           <<", skip it";continue; }

      //loop over the rechit on the module
      for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iTThit = thoseHits.begin();iTThit != thoseHits.end(); iTThit++)
        {
          if (!(*iTThit)->isValid())/*skip*/{edm::LogInfo(theCategory)<<"rec hit not valid on module "<<presentdetid.rawId() <<". I skip it"; continue;}

          LogDebug(theCategory)<<((dynamic_cast<const SiStripMatchedRecHit2D *>((*iTThit)->hit())!=0) ? /*matched rechit*/ "matched rechit" : " r-phi rechit");

          //get the surface of the module Id
          const BoundPlane & surface = (*iTThit)->det()->surface();

          //estimate the consistency
          MeasurementEstimator::HitReturnType est_road = theRoadEstimator->estimate(restep,**iTThit);

          LogDebug(theCategory)<<"(restep.globalPosition().perp()) "<<(restep.globalPosition().perp())<<"(restep.globalPosition().mag()) "<<(restep.globalPosition().mag())
                             <<"road step parameters at module: \n"<<restep
                               << "(est_road.first) "<<(est_road.first)<<"(est_road.second) "<<(est_road.second);

          //check consistency
          if (est_road.first)
            { //hit and propagation are consistent : gather the hit
              Nhits++;
              additionalHits++;
              atleastoneadded=true;
              
              LogDebug(theCategory)<<"hit is consistent with road"<<"(presentdetid.rawId()) "<<(presentdetid.rawId())
                                 <<"loop over previous trajectories\n"
                                   <<"(Trajectories.tail().size()) "<<(Trajectories.tail().size())<<"(Trajectories.head().size()) "<<(Trajectories.head().size());
                
              //update the list of trajectory that we have with the consistent rechit.
              //loop over the existing list of trajectories and add the hit if necessary
              for ( TrajectoryCollection::iterator traj =Trajectories.head().begin();traj != Trajectories.head().end(); traj++)
                {
                  //what is the previous state for this trajectory
                  const TrajectoryStateOnSurface & previousState = traj->TSOS;
                  if (!previousState.isValid())
                    {edm::LogError(theCategory)<<"previous free state is invalid at trajectory update, this is WRONG"; continue;}
                  const FreeTrajectoryState & previousFreestate = *previousState.freeState();
                  
                  //propagate it to the current surface
                  TrajectoryStateOnSurface predictedState = thePropagator->propagate(previousFreestate,surface);
                  if (!predictedState.isValid())/*skip*/{edm::LogError(theCategory)
                        <<"predicted state is not valid at trajectory update, rechit surface cannot be reached by the previous updated state.";continue;}

                  MeasurementEstimator::HitReturnType est ;
                  est= theHitEstimator->estimate(predictedState,**iTThit);

                  LogDebug(theCategory)<<"(est.first) "<<(est.first)<<"(est.second) "<<(est.second);

                  //is the current hit consistent with the trajectory
                  if (est.first )
                    {
                      //update the trajectory state with the rechit
                      const TrajectoryStateOnSurface & updatedState = theUpdator->update(predictedState,**iTThit);
                      if (!updatedState.isValid())/*skip*/{edm::LogError(theCategory)<<"updated state is not valid, this is really wrong";continue;}

                      LogDebug(theCategory)<<"updating a previous state with a rechit";

                      //add a combined trajectory to the new list of trajectory, starting from the existing trajecotry
                      Trajectories.tail().push_back(*traj); //from existing
                      trajectory & combined = (*Trajectories.tail().rbegin());
                      combined.duplicate =false; //this is important
                      //increment the chi2
                      //combined.chi2 += est.second;
                      combined.chi2 += est_road.second;//better to add the road chi2 too unbias the chi2 sum towards first hits
                      //some history about the trajectory
                      combined.lastmissed=false;
                      combined.missedinarow=0;
                      //add a new hits to the measurements
                      combined.measurements.push_back(TrajectoryMeasurement(updatedState,*iTThit));
                      TrajectoryMeasurement & trajMeasurement = (*combined.measurements.rbegin());
                      //assigne updated state
                      combined.TSOS = updatedState;
                      //add trajectory measurement 
                      combined.traj.push(trajMeasurement,est.second);

                      LogDebug(theCategory)<<"(combined.traj.foundHits()) "<<(combined.traj.foundHits())
                                         <<"muon measurement on previous module: \n"<<traj->TSOS
                                         <<"muon measurement before update: \n"<<predictedState
                                         <<"muon measurement after update: \n"<<updatedState;

                    }//hit consistent with trajectory
                  else
                    {
                      LogDebug(theCategory)<<"hit failed chi2 test for trajectories update\n"<<"(traj->duplicate) "<<(traj->duplicate);

                      if (!traj->duplicate){
                        Trajectories.tail().push_back(*traj);
                        trajectory & copied = (*Trajectories.tail().rbegin());
                        copied.missed++;
                        copied.lastmissed=true;
                        copied.missedinarow++;
                        traj->duplicate =true; //set this trajectory to already have been copied into next step
                       }//not a duplicate trajectory
                    }//hit not consistent with trajectory
                }//existing trajectories loop
            }//hit is consistent with muon road
        }//rechits on module loop
      
      
      //discard previous list of trajectories
      //if something as been done of course
      if (Trajectories.tail().size()!=0) 
        {
          //if this is not the first "layer" of detector, set updated trajectories as new seed trajectories
          //will branch on every single rechit uncountered for first "layer"
          if (!theFirstlayer || theBranchonfirstlayer)
            {
              LogDebug(theCategory)<<"swapping trajectory list index";
              
              //always carry on the <IP> alone state in the next list of trajectories to avoid bias from the first rechits
              if (theFirstlayer && theCarriedIPatfirstlayerModule)
                {
                  LogDebug(theCategory)<<"push front <IP> to next list of trajectories";

                  Trajectories.tail().push_front(*Trajectories.head().begin()); //[0] is <IP> always
                  trajectory & pushed = *Trajectories.tail().begin();
                  pushed.missed+=additionalHits;
                  pushed.lastmissed = ( additionalHits!=0);
                  pushed.missedinarow++;
                }
              //FIXME, there is a candidate leak at this point
              //if the IP state is carried on at first layer, without update, then it will be duplicated. this is very unlikely though

              Trajectories.head().clear();

              //swap the lists
              Trajectories.flip();
            }
        }//discard previous list of trajectories
    }//module loop

  

  //do some cleaning of the list of trajectories
  if(theFirstlayer && atleastoneadded )
    {
      theFirstlayer =false; //we are not in the first layer if something has been added
      if (!theBranchonfirstlayer)
        {
          LogDebug(theCategory)<<"swapping trajectory list index (end of first layer)";

          //and for consistency, you have to swap index here, because it could ahve not been done in the module loop above
          //always carry on the <IP> alone state in the next list of trajectories to avoid bias from the first rechits
          if (theCarriedIPatfirstlayer)
            {Trajectories.tail().push_front(*Trajectories.head().begin());} //[0] is <IP> always at this stage
          //FIXME, there is a candidate leak at this point
          //if the IP state is carried on at first layer, without update, then it will be duplicated. this is very unlikely though
          
          Trajectories.head().clear();
          //swap the switch
          Trajectories.flip();
        }
      else
        {
          //actually remove the <IP> from first postion of the next source of trajectories
          //since the swaping has already been done. pop up from theTrajectorysource, not !theTrajectorysource
          //only if it has been done at the first layer module though
          if (!theCarriedIPatfirstlayer && theCarriedIPatfirstlayerModule){

            LogDebug(theCategory)<<"pop up <IP> from trajectories";

            Trajectories.head().pop_front(); }
        }


      //check an remove trajectories that are subset of the other in next source
      //after the first layer only
      if (Trajectories.head().size()>=2)
        {checkDuplicate(Trajectories.head());}

    }  //do some cleaning of the list of trajectories
  return Nhits;

}
void MuonRoadTrajectoryBuilder::makeTrajectories ( const TrajectorySeed seed,
std::vector< Trajectory > &  result,
int  version = 0 
) const [private]

Definition at line 143 of file MuonRoadTrajectoryBuilder.cc.

References makeTrajectories_0(), and makeTrajectories_1().

Referenced by trajectories().

{  if (version==0) { makeTrajectories_0(seed,result);}
  else if (version==1) { makeTrajectories_1(seed,result);}}
void MuonRoadTrajectoryBuilder::makeTrajectories_0 ( const TrajectorySeed seed,
std::vector< Trajectory > &  result 
) const [private]

Definition at line 150 of file MuonRoadTrajectoryBuilder.cc.

References alongMomentum, checkStep(), PTrajectoryStateOnDet::detId(), TrajectoryStateOnSurface::freeState(), GatherHits(), MeasurementTracker::geometricSearchTracker(), MeasurementTracker::geomTracker(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), MuonRoadTrajectoryBuilder::flippingPair< A >::head(), TrackingGeometry::idToDet(), SimpleDiskBounds::innerRadius(), TrajectoryStateOnSurface::isValid(), LogDebug, MuonRoadTrajectoryBuilder::trajectory::measurements, GeometricSearchTracker::negTecLayers(), GeometricSearchTracker::negTidLayers(), SimpleDiskBounds::outerRadius(), PV3DBase< T, PVType, FrameType >::perp(), position, GeometricSearchTracker::posTecLayers(), GeometricSearchTracker::posTidLayers(), Propagator::propagate(), csvReporter::r, smooth(), TrajectorySeed::startingState(), GeomDet::surface(), sistripsummary::TEC, theCategory, theField, theFirstlayer, theMeasurementTracker, theMinNumberOfHitOnCandidate, theNumberOfHitPerModule, theNumberOfHitPerModuleDefault, theOutputAllTraj, thePropagator, theTransformer, sistripsummary::TIB, GeometricSearchTracker::tibLayers(), sistripsummary::TID, sistripsummary::TOB, GeometricSearchTracker::tobLayers(), MuonRoadTrajectoryBuilder::trajectory::traj, trajectoryOrder(), TrajectoryStateTransform::transientState(), MuonRoadTrajectoryBuilder::trajectory::TSOS, PV3DBase< T, PVType, FrameType >::z(), and z.

Referenced by makeTrajectories().

{
  Trajectory basicTrajectory(seed,alongMomentum);
  //add the muon system measurement to the basicTrajectory
  //...no yet implemented...

  //get the initial state
  PTrajectoryStateOnDet   PTstart=seed.startingState();
  DetId detId(PTstart.detId());
  const BoundPlane * surface =(&theMeasurementTracker->geomTracker()->idToDet(detId)->surface());
  //start from this point
  TrajectoryStateOnSurface TSOS = theTransformer.transientState(PTstart,surface,theField);
  if (!TSOS.isValid())/*/abort*/{ edm::LogError(theCategory)<<"TSOS from PTSOD is not valid.";return ;}
  
  LogDebug(theCategory) <<"(detId.rawId()) "<<(detId.rawId())<<"(detId.subdetId()) "<<(detId.subdetId())
                        <<"(TSOS.globalPosition()) "<<(TSOS.globalPosition())
                        <<"(TSOS.globalMomentum()) "<<(TSOS.globalMomentum());
 
  //initialization
  //----------------------
  //  WARNING. this is mutable 
  //  you should never remove this line
  theFirstlayer=true;
  theNumberOfHitPerModule = theNumberOfHitPerModuleDefault;
  //----------------------
  int Nhits=0;
  GlobalPoint position;
  double z=0,r=0;
  //flipping pair of list of trajectories
  TrajectoryCollectionFPair Trajectories; 
  //initialize the head() of pair of list of trajectories
  Trajectories.head().push_back(trajectory());
  trajectory & initial = *Trajectories.head().rbegin();
  initial.TSOS=TSOS;  
  initial.traj=basicTrajectory;      
  


  //get the DetLayer in the tracker
  std::vector<BarrelDetLayer*> tiblc = theMeasurementTracker->geometricSearchTracker()->tibLayers();
  std::vector<BarrelDetLayer*> toblc = theMeasurementTracker->geometricSearchTracker()->tobLayers();
  std::vector<ForwardDetLayer*> tidlc[2];
  tidlc[0]=theMeasurementTracker->geometricSearchTracker()->negTidLayers();
  tidlc[1]=theMeasurementTracker->geometricSearchTracker()->posTidLayers();
  std::vector<ForwardDetLayer*> teclc[2];
  teclc[0]=theMeasurementTracker->geometricSearchTracker()->negTecLayers();
  teclc[1]=theMeasurementTracker->geometricSearchTracker()->posTecLayers();
  const int Ntib=tiblc.size();
  const int Ntob=toblc.size();
  const int Ntid=tidlc[0].size();


  const int Ntec=teclc[0].size();

  LogDebug(theCategory)<<"(Ntib) "<<(Ntib)<<"(Ntob) "<<(Ntob)<<"(Ntid) "<<(Ntid)<<"(Ntec) "<<(Ntec);

  const SimpleDiskBounds * sdbounds=0;

  position = TSOS.globalPosition();
  z = position.z();
  r = position.perp();

  //select which part we are in
  enum PART { fault , PXB, PXF, TIB , TID , TOB , TEC};
  PART whichpart = fault;
  switch(detId.subdetId()){
  case 1: {whichpart=PXB;break;}
  case 2: {whichpart=PXF;break;}

  case 3: {whichpart=TIB;break;}
  case 4: {whichpart=TID;break;}
  case 5: {whichpart=TOB;break;}
  case 6: {whichpart=TEC;break;}}

  bool inapart=true;
  bool firstRound =true;
  int indexinpart=0;
  //loop while in a valid part of the tracker
  while(inapart){

      //check on the trajectory list and stuff
      if (!checkStep(Trajectories.head())) break;
      //................

    switch (whichpart){

    case fault: /*abort*/ {edm::LogError(theCategory)<<"something's wrong with the seed";return;}
    case PXB: /*abort*/ {edm::LogError(theCategory)<<"PXB no yet implemented";return;}
    case PXF: /*abort*/ {edm::LogError(theCategory)<<"PXF no yet implemented";return;}
        
      //-------------- into TIB----------------
    case TIB:{
      if (indexinpart==Ntib){/*you have reach the last layer of the TIB.let's go to the TOB*/whichpart = TOB; indexinpart=-1;break; }

      LogDebug(theCategory)<<"within TIB "<<indexinpart+1;
        
      //propagate to corresponding surface
      if (!firstRound) TSOS = thePropagator->propagate(*TSOS.freeState(),tiblc[indexinpart]->surface());
      if (!TSOS.isValid()) {;break;} //go to the next one 
            
      z=TSOS.globalPosition().z();

      //have we reached a boundary
      if (fabs(z) > fabs(tidlc[(z>0)][0]->surface().position().z())) 
        {/*z bigger than the TID min z: go to TID*/
          LogDebug(theCategory)<<"|z| ("<<z<<") bigger than the TID min z("<<tidlc[(z>0)][0]->surface().position().z()<<"): go to TID";
          whichpart = TID; indexinpart=-1; break;}
      else  {/*gather hits in the corresponding TIB layer*/ Nhits+=GatherHits(TSOS,tiblc[indexinpart],Trajectories);}
      break;}

      //-------------- into TID----------------
    case TID: {
      if (indexinpart==Ntid){/*you have reach the last layer of the TID. let's go to the TEC */ whichpart = TEC; indexinpart=-1; break;}

      LogDebug(theCategory)<<"within TID "<<indexinpart+1;

      //propagate to corresponding surface
      if (!firstRound) TSOS = thePropagator->propagate(*TSOS.freeState(),tidlc[(z>0)][indexinpart]->surface());
      if (!TSOS.isValid()){break;}//go to the next one 
          
      position =  TSOS.globalPosition();
      z = position.z();
      r = position.perp();

      sdbounds = dynamic_cast<const SimpleDiskBounds *>(&tidlc[(z>0)][indexinpart]->surface().bounds());
      if (!sdbounds)/*abort*/{edm::LogError(theCategory)<<" detlayer bounds are not SimpleDiskBounds in tid geometry";return;}
        
      //have we reached a boundary
      if (r < sdbounds->innerRadius())
        {/*radius smaller than the TID disk inner radius: next disk please*/
          LogDebug(theCategory)<<"radius ("<<r<<") smaller than the TID disk inner radius ("<<sdbounds->innerRadius()<<"): next disk please";
          break;}
      else if (r >sdbounds->outerRadius())
        {/*radius bigger than the TID disk outer radius: go to TOB*/
          LogDebug(theCategory)<<"radius ("<<r<<") bigger than the TID disk outer radius("<<sdbounds->outerRadius()<<"): go to TOB";
          whichpart = TOB; indexinpart=-1;break;}
      else  {/*gather hits in the corresponding TIB layer*/
        LogDebug(theCategory)<<"collecting hits";
        Nhits+=GatherHits(TSOS,tidlc[(z>0)][indexinpart],Trajectories);}
      break;}

      //-------------- into TOB----------------
    case TOB: {
      if (indexinpart==Ntob){/*you have reach the last layer of the TOB. this is an end*/ inapart=false;break;}

      LogDebug(theCategory)<<"within TOB "<<indexinpart+1;
        
      //propagate to corresponding surface
      if (!firstRound) TSOS = thePropagator->propagate(*TSOS.freeState(),toblc[indexinpart]->surface());
      if (!TSOS.isValid())  {break;} //go to the next one 
            
      z =  TSOS.globalPosition().z();
        
      //have we reached a boundary
      if (fabs(z) > fabs(teclc[(z>0)][0]->surface().position().z()))
        {/*z bigger than the TOB layer max z: go to TEC*/
          LogDebug(theCategory)<<"|z| ("<<z<<") bigger than the TOB layer max z ("<< teclc[(z>0)][0]->surface().position().z()<<"): go to TEC";
          whichpart = TEC; indexinpart=-1;break;}
      else {/*gather hits in the corresponding TOB layer*/Nhits+=GatherHits(TSOS,toblc[indexinpart],Trajectories);}
      break;}

      //-------------- into TEC----------------
    case TEC:   {
      if (indexinpart==Ntec){/*you have reach the last layer of the TEC. let's end here*/inapart=false;break;}
          
      LogDebug(theCategory)<<"within TEC "<<indexinpart+1;
          
      //propagate to corresponding TEC disk
      if (!firstRound) TSOS = thePropagator->propagate(*TSOS.freeState(),teclc[(z>0)][indexinpart]->surface());
      if (!TSOS.isValid()) {break;} //go to the next one
          
      position = TSOS.globalPosition();
      z = position.z();
      r = position.perp();
          
      sdbounds = dynamic_cast<const SimpleDiskBounds *>(&teclc[(z>0)][indexinpart]->surface().bounds());
      if (!sdbounds)/*abort*/ {edm::LogError(theCategory)<<" detlayer bounds are not SimpleDiskBounds in tec geometry";return;}
          
      //have we reached a boundary ?
      if (r < sdbounds->innerRadius())
        {/*radius smaller than the TEC disk inner radius: next disk please*/
          LogDebug(theCategory)<<"radius ("<<r<<") smaller than the TEC disk inner radius ("<<sdbounds->innerRadius()<<"): next disk please";
          break;}
      else if (r > sdbounds->outerRadius())
        {/*radius bigger than the TEC disk outer radius: I can stop here*/
          LogDebug(theCategory)<<"radius ("<<r<<") bigger than the TEC disk outer radius ("<<sdbounds->outerRadius()<<"): I can stop here";
          inapart=false;break;}
      else {/*gather hits in the corresponding TEC layer*/Nhits+=GatherHits(TSOS,teclc[(z>0)][indexinpart],Trajectories);}
          
      break;}

    }//switch
    indexinpart++;
    firstRound=false;
  }//while inapart


  LogDebug(theCategory)<<"propagating through layers is done";
  
  //--------------------------------------- SWIMMING DONE --------------------------------------
  //the list of trajectory has been build. let's find out which one is the best to use.
  
  edm::LogInfo(theCategory)<<"found: "
                                      <<Nhits<<" hits in the road for this seed \n"
                                      << Trajectories.head().size()<<" possible trajectory(ies)";
  
  if (Trajectories.head().size() == 0) /*abort*/{edm::LogError(theCategory)<<" no possible trajectory found"; return;}
  
  //order the list to put the best in front
  Trajectories.head().sort(trajectoryOrder);
  
  //------------------------------------OUTPUT FILL---------------------------------------------
  if (theOutputAllTraj) {
    //output all the possible trajectories found, if they have enough hits on them
    //the best is still the first
    for (TrajectoryCollection::iterator tit = Trajectories.head().begin(); tit!=Trajectories.head().end();tit++)
      {if (tit->measurements.size()>= theMinNumberOfHitOnCandidate){result.push_back(smooth(tit->traj));}}
  }
  else{
    //output only the best one
    trajectory & best = (*Trajectories.head().begin());
    if (best.measurements.size()< theMinNumberOfHitOnCandidate)/*abort*/{edm::LogError(theCategory)<<"best trajectory does not have enough ("<<best.measurements.size()<<") hits on it (<"<<theMinNumberOfHitOnCandidate<<")"; return;}
    //output only the best trajectory
    result.push_back(smooth(best.traj));}
  
  edm::LogInfo(theCategory)<<result.size()<<" trajectory(ies) output";
}//makeTrajectories_0
void MuonRoadTrajectoryBuilder::makeTrajectories_1 ( const TrajectorySeed seed,
std::vector< Trajectory > &  result 
) const [private]

Definition at line 762 of file MuonRoadTrajectoryBuilder.cc.

References alongMomentum.

Referenced by makeTrajectories().

{
  Trajectory basicTrajectory(seed,alongMomentum);
  //add the muon system measurement to the basicTrajectory
  //...

//   //build the trajectories
//   std::vector<Trajectory> unsmoothed = theCkfbuilder->trajectories(seed);
  
//   //smoothed them
//   if (theOutputAllTraj) {
//     for (std::vector<Trajectory>::iterator tit = unsmoothed.begin(); tit!=unsmoothed.end();tit++)
//       {result.push_back(smooth(*tit));}
//   }
//   else{
//     //output only the first one
//     result.push_back(smooth(unsmoothed.front()));} 
    
}
void MuonRoadTrajectoryBuilder::setEvent ( const edm::Event iEvent) const [virtual]
Trajectory MuonRoadTrajectoryBuilder::smooth ( Trajectory traj) const [private]

Definition at line 440 of file MuonRoadTrajectoryBuilder.cc.

References cleanTrajectory(), runTheMatrix::ret, sortTrajectoryMeasurements(), theCategory, theSmoother, and KFTrajectorySmoother::trajectories().

Referenced by makeTrajectories_0().

                                                                    {
  
  //need to order the list of measurements on the trajectory first
  sortTrajectoryMeasurements(traj);

  std::vector<Trajectory> ret=theSmoother->trajectories(traj);

  if (ret.empty()){
    edm::LogError(theCategory)<<"smoother returns an empty vector of trajectories: try cleaning first and resmooth";
    cleanTrajectory(traj);
    ret=theSmoother->trajectories(traj);
  }
  
  if (ret.empty()){edm::LogError(theCategory)<<"smoother returns an empty vector of trajectories.";
    return traj;}
  else{return (ret.front());}
}
void MuonRoadTrajectoryBuilder::trajectories ( const TrajectorySeed seed,
TrajectoryContainer ret 
) const [virtual]

process the seed, in a faster manner.

Reimplemented from TrajectoryBuilder.

Definition at line 119 of file MuonRoadTrajectoryBuilder.cc.

References LogDebug, makeTrajectories(), and theCategory.

{
  LogDebug(theCategory)<<"makeTrajectories START";
  //process the seed through the tracker
  makeTrajectories(seed,ret);
}
std::vector< Trajectory > MuonRoadTrajectoryBuilder::trajectories ( const TrajectorySeed seed) const [virtual]

Implements TrajectoryBuilder.

Definition at line 127 of file MuonRoadTrajectoryBuilder.cc.

References LogDebug, makeTrajectories(), query::result, and theCategory.

{
  LogDebug(theCategory)<<"makeTrajectories START";
  
  //default output
  std::vector<Trajectory> result;
  
  //process the seed through the tracker
  makeTrajectories(seed,result);
  
  //output the result of regional tracking
  return result;
}

Member Data Documentation

Definition at line 142 of file MuonRoadTrajectoryBuilder.h.

Referenced by GatherHits(), and MuonRoadTrajectoryBuilder().

Definition at line 143 of file MuonRoadTrajectoryBuilder.h.

Referenced by GatherHits(), and MuonRoadTrajectoryBuilder().

Definition at line 144 of file MuonRoadTrajectoryBuilder.h.

Referenced by GatherHits(), and MuonRoadTrajectoryBuilder().

Info/Debug category "Muon|RecoMuon|MuonRoadTrajectoryBuilder".

Definition at line 75 of file MuonRoadTrajectoryBuilder.h.

Referenced by checkDuplicate(), checkStep(), cleanTrajectory(), GatherHits(), makeTrajectories_0(), MuonRoadTrajectoryBuilder(), smooth(), trajectories(), and ~MuonRoadTrajectoryBuilder().

Definition at line 135 of file MuonRoadTrajectoryBuilder.h.

Referenced by checkStep(), and MuonRoadTrajectoryBuilder().

Definition at line 171 of file MuonRoadTrajectoryBuilder.h.

Referenced by makeTrajectories_0(), and MuonRoadTrajectoryBuilder().

Definition at line 126 of file MuonRoadTrajectoryBuilder.h.

Referenced by GatherHits(), and makeTrajectories_0().

Definition at line 129 of file MuonRoadTrajectoryBuilder.h.

Referenced by checkStep(), and MuonRoadTrajectoryBuilder().

std::vector<unsigned int> MuonRoadTrajectoryBuilder::theMaxTrajectoriesThreshold [private]

Definition at line 138 of file MuonRoadTrajectoryBuilder.h.

Referenced by checkStep(), and MuonRoadTrajectoryBuilder().

Definition at line 148 of file MuonRoadTrajectoryBuilder.h.

Referenced by makeTrajectories_0(), and MuonRoadTrajectoryBuilder().

unsigned int MuonRoadTrajectoryBuilder::theNumberOfHitPerModule [mutable, private]

Definition at line 137 of file MuonRoadTrajectoryBuilder.h.

Referenced by checkStep(), GatherHits(), and makeTrajectories_0().

Definition at line 136 of file MuonRoadTrajectoryBuilder.h.

Referenced by makeTrajectories_0(), and MuonRoadTrajectoryBuilder().

Definition at line 139 of file MuonRoadTrajectoryBuilder.h.

Referenced by checkStep(), and MuonRoadTrajectoryBuilder().

Definition at line 150 of file MuonRoadTrajectoryBuilder.h.

Referenced by makeTrajectories_0(), and MuonRoadTrajectoryBuilder().

Definition at line 168 of file MuonRoadTrajectoryBuilder.h.

Referenced by makeTrajectories_0().