BaseCkfTrajectoryBuilder.cc

#include "RecoTracker/CkfPattern/interface/BaseCkfTrajectoryBuilder.h"

#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
#include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"

  
#include "TrackingTools/DetLayers/interface/NavigationSchool.h"
#include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
#include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"

#include "TrackingTools/PatternTools/interface/TempTrajectory.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"

BaseCkfTrajectoryBuilder::BaseCkfTrajectoryBuilder(const edm::ParameterSet& conf,
                                                   TrajectoryFilter *filter,
                                                   TrajectoryFilter *inOutFilter):
  theUpdator(nullptr),
  thePropagatorAlong(nullptr),
  thePropagatorOpposite(nullptr),
  theEstimator(nullptr),
  theTTRHBuilder(nullptr),
  theMeasurementTracker(nullptr),
  theFilter(filter),
  theInOutFilter(inOutFilter),
  theUpdatorName(conf.getParameter<std::string>("updator")),
  thePropagatorAlongName(conf.getParameter<std::string>("propagatorAlong")),
  thePropagatorOppositeName(conf.getParameter<std::string>("propagatorOpposite")),
  theEstimatorName(conf.getParameter<std::string>("estimator")),
  theRecHitBuilderName(conf.getParameter<std::string>("TTRHBuilder"))
{
  if (conf.exists("clustersToSkip")) edm::LogError("BaseCkfTrajectoryBuilder") << "ERROR: " << typeid(*this).name() << " has a clustersToSkip parameter set";
}


BaseCkfTrajectoryBuilder::~BaseCkfTrajectoryBuilder(){
}

TrajectoryFilter *BaseCkfTrajectoryBuilder::createTrajectoryFilter(const edm::ParameterSet& pset, edm::ConsumesCollector& iC) {
  return TrajectoryFilterFactory::get()->create(pset.getParameter<std::string>("ComponentType"), pset, iC);
}

void
BaseCkfTrajectoryBuilder::seedMeasurements(const TrajectorySeed& seed,  TempTrajectory & result) const
{
  

  TrajectorySeed::range hitRange = seed.recHits();

  PTrajectoryStateOnDet pState( seed.startingState());
  const GeomDet* gdet = theMeasurementTracker->geomTracker()->idToDet(pState.detId());
  TSOS outerState = trajectoryStateTransform::transientState(pState, &(gdet->surface()),
							     forwardPropagator(seed)->magneticField());


  for (TrajectorySeed::const_iterator ihit = hitRange.first; ihit != hitRange.second; ihit++) {
 
    TrackingRecHit::RecHitPointer recHit = ihit->cloneSH();
    const GeomDet* hitGeomDet = recHit->det();
 
    const DetLayer* hitLayer = 
      theMeasurementTracker->geometricSearchTracker()->detLayer(ihit->geographicalId());

    TSOS invalidState( hitGeomDet->surface());
    if (ihit == hitRange.second - 1) {
      // the seed trajectory state should correspond to this hit
      if (&gdet->surface() != &hitGeomDet->surface()) {
	edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: the seed state is not on the surface of the detector of the last seed hit";
	return; // FIXME: should throw exception
      }

      //TSOS updatedState = outerstate;
      result.emplace(invalidState, outerState, recHit, 0, hitLayer);
    }
    else {
      TSOS innerState   = backwardPropagator(seed)->propagate(outerState,hitGeomDet->surface());

      // try to recover if propagation failed
      if unlikely(!innerState.isValid())
        innerState =
          trajectoryStateTransform::transientState(pState, &(hitGeomDet->surface()),
                                                            forwardPropagator(seed)->magneticField());

      if(innerState.isValid()) {
	TSOS innerUpdated = theUpdator->update(innerState,*recHit);
	result.emplace(invalidState, innerUpdated, recHit, 0, hitLayer);
      }
    }
  }

  // method for debugging
  // fix somehow
  // fillSeedHistoDebugger(result.begin(),result.end());

}


TempTrajectory BaseCkfTrajectoryBuilder::
createStartingTrajectory( const TrajectorySeed& seed) const
{
  TempTrajectory result(seed.direction(),seed.nHits());
  seedMeasurements(seed, result);

  LogDebug("CkfPattern")
    <<" initial trajectory from the seed: "<<PrintoutHelper::dumpCandidate(result,true);
  
  return result;
}


bool BaseCkfTrajectoryBuilder::toBeContinued (TempTrajectory& traj, bool inOut) const
{
  if unlikely(traj.measurements().size() > 400) {
    edm::LogError("BaseCkfTrajectoryBuilder_InfiniteLoop");
    LogTrace("BaseCkfTrajectoryBuilder_InfiniteLoop") << 
              "Cropping Track After 400 Measurements:\n" <<
              "   Last predicted state: " << traj.lastMeasurement().predictedState() << "\n" <<
              "   Last layer subdetector: " << (traj.lastLayer() ? traj.lastLayer()->subDetector() : -1) << "\n" <<
              "   Found hits: " << traj.foundHits() << ", lost hits: " << traj.lostHits() << "\n\n";
    return false;
  }
  // Called after each new hit is added to the trajectory, to see if it is 
  // worth continuing to build this track candidate.
  if (inOut) {
    // if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
    return theInOutFilter->toBeContinued(traj);
  } else {
    return theFilter->toBeContinued(traj);
  }
}


 bool BaseCkfTrajectoryBuilder::qualityFilter( const TempTrajectory& traj, bool inOut) const
{
  // Called after building a trajectory is completed, to see if it is good enough
  // to keep.
  if (inOut) {
    // if (theInOutFilter == 0) edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: trying to use dedicated filter for in-out tracking phase, when none specified";
    return theInOutFilter->qualityFilter(traj);
  } else {
    return theFilter->qualityFilter(traj);
  }
}


void 
BaseCkfTrajectoryBuilder::addToResult (boost::shared_ptr<const TrajectorySeed> const & seed, TempTrajectory& tmptraj, 
				       TrajectoryContainer& result,
                                       bool inOut) const
{
  // quality check
  if ( !qualityFilter(tmptraj, inOut) )  return;
  Trajectory traj = tmptraj.toTrajectory();
  traj.setSharedSeed(seed);
  // discard latest dummy measurements
  while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid()) traj.pop();
  LogDebug("CkfPattern")<<inOut<<"=inOut option. pushing a Trajectory with: "<<traj.foundHits()<<" found hits. "<<traj.lostHits()
			<<" lost hits. Popped :"<<(tmptraj.measurements().size())-(traj.measurements().size())<<" hits.";
  result.push_back(std::move(traj));
}


void 
BaseCkfTrajectoryBuilder::addToResult (TempTrajectory const & tmptraj, 
				       TempTrajectoryContainer& result,
                                       bool inOut) const
{
  // quality check
  if ( !qualityFilter(tmptraj, inOut) )  return;
  // discard latest dummy measurements
  TempTrajectory traj = tmptraj;
  while (!traj.empty() && !traj.lastMeasurement().recHit()->isValid()) traj.pop();
  LogDebug("CkfPattern")<<inOut<<"=inOut option. pushing a TempTrajectory with: "<<traj.foundHits()<<" found hits. "<<traj.lostHits()
			<<" lost hits. Popped :"<<(tmptraj.measurements().size())-(traj.measurements().size())<<" hits.";
  result.push_back(std::move(traj));
}

void 
BaseCkfTrajectoryBuilder::moveToResult (TempTrajectory&& traj, 
				       TempTrajectoryContainer& result,
                                       bool inOut) const
{
  // quality check
  if ( !qualityFilter(traj, inOut) )  return;
  // discard latest dummy measurements
  while (!traj.empty() && !traj.lastMeasurement().recHitR().isValid()) traj.pop();
  LogDebug("CkfPattern")<<inOut<<"=inOut option. pushing a TempTrajectory with: "<<traj.foundHits()<<" found hits. "<<traj.lostHits();
    //			<<" lost hits. Popped :"<<(ttraj.measurements().size())-(traj.measurements().size())<<" hits.";
  result.push_back(std::move(traj));
}



BaseCkfTrajectoryBuilder::StateAndLayers
BaseCkfTrajectoryBuilder::findStateAndLayers(const TrajectorySeed& seed, const TempTrajectory& traj) const
{
  if (traj.empty())
    {
      //set the currentState to be the one from the trajectory seed starting point
      PTrajectoryStateOnDet const & ptod = seed.startingState();
      DetId id(ptod.detId());
      const GeomDet * g = theMeasurementTracker->geomTracker()->idToDet(id);                    
      const Surface * surface=&g->surface();
      
      
      TSOS currentState(trajectoryStateTransform::transientState(ptod,surface,forwardPropagator(seed)->magneticField()));      
      const DetLayer* lastLayer = theMeasurementTracker->geometricSearchTracker()->detLayer(id);      
      return StateAndLayers(currentState,theNavigationSchool->nextLayers(*lastLayer,*currentState.freeState(), traj.direction()) );
    }
  else
    {  
      TSOS const & currentState = traj.lastMeasurement().updatedState();
      return StateAndLayers(currentState,theNavigationSchool->nextLayers(*traj.lastLayer(), *currentState.freeState(), traj.direction()) );
    }
}

void BaseCkfTrajectoryBuilder::setData(const MeasurementTrackerEvent *data) 
{
    // possibly do some sanity check here
    theMeasurementTracker = data;
}

void BaseCkfTrajectoryBuilder::setEvent(const edm::Event& event) const
{
    std::cerr << "ERROR SetEvent called on " << typeid(*this).name() << ( theMeasurementTracker ? " with valid " : "witout any ") << "MeasurementTrackerEvent" << std::endl;
}

void BaseCkfTrajectoryBuilder::unset() const
{
    std::cerr << "ERROR unSet called on " << typeid(*this).name() << ( theMeasurementTracker ? " with valid " : "witout any ") << "MeasurementTrackerEvent" << std::endl;
}

void BaseCkfTrajectoryBuilder::setEvent(const edm::Event& iEvent, const edm::EventSetup& iSetup, const MeasurementTrackerEvent *data) {
  edm::ESHandle<TrajectoryStateUpdator> updatorHandle;
  edm::ESHandle<Propagator>             propagatorAlongHandle;
  edm::ESHandle<Propagator>             propagatorOppositeHandle;
  edm::ESHandle<Chi2MeasurementEstimatorBase> estimatorHandle;
  edm::ESHandle<TransientTrackingRecHitBuilder> recHitBuilderHandle;

  iSetup.get<TrackingComponentsRecord>().get(theUpdatorName, updatorHandle);
  iSetup.get<TrackingComponentsRecord>().get(thePropagatorAlongName, propagatorAlongHandle);
  iSetup.get<TrackingComponentsRecord>().get(thePropagatorOppositeName, propagatorOppositeHandle);
  iSetup.get<TrackingComponentsRecord>().get(theEstimatorName, estimatorHandle);
  iSetup.get<TransientRecHitRecord>().get(theRecHitBuilderName, recHitBuilderHandle);

  theUpdator = updatorHandle.product();
  thePropagatorAlong = propagatorAlongHandle.product();
  thePropagatorOpposite = propagatorOppositeHandle.product();
  theEstimator = estimatorHandle.product();
  theTTRHBuilder = recHitBuilderHandle.product();

  setData(data);
  if(theFilter) theFilter->setEvent(iEvent, iSetup);
  if(theInOutFilter) theInOutFilter->setEvent(iEvent, iSetup);
  setEvent_(iEvent, iSetup);
}