DAFTrackProducerAlgorithm Class Reference

#include <DAFTrackProducerAlgorithm.h>

Public Member Functions
	DAFTrackProducerAlgorithm (const edm::ParameterSet &conf)
void	runWithCandidate (const TrackingGeometry *, const MagneticField *, const std::vector< Trajectory > &, const MeasurementTrackerEvent *measTk, const TrajectoryFitter *, const TransientTrackingRecHitBuilder *, const MultiRecHitCollector *measurementTracker, const SiTrackerMultiRecHitUpdator *, const reco::BeamSpot &, AlgoProductCollection &, TrajAnnealingCollection &, bool) const
	Run the Final Fit taking TrackCandidates as input. More...
	~DAFTrackProducerAlgorithm ()

Private Types
typedef std::pair< Trajectory *, std::pair< reco::Track *, PropagationDirection > >	AlgoProduct
typedef std::vector< AlgoProduct >	AlgoProductCollection
typedef std::vector < TrajAnnealing >	TrajAnnealingCollection

Private Member Functions
bool	buildTrack (const Trajectory, AlgoProductCollection &algoResults, float, const reco::BeamSpot &) const
	Construct Tracks to be put in the event. More...
float	calculateNdof (const Trajectory vtraj) const
int	checkHits (Trajectory iInitTraj, const Trajectory iFinalTraj) const
std::pair < TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface >	collectHits (const Trajectory vtraj, const MultiRecHitCollector *measurementCollector, const MeasurementTrackerEvent *measTk) const
void	filter (const TrajectoryFitter *fitter, std::vector< Trajectory > &input, int minhits, std::vector< Trajectory > &output, const TransientTrackingRecHitBuilder *builder) const
Trajectory	fit (const std::pair< TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface > &hits, const TrajectoryFitter *theFitter, Trajectory vtraj) const
	accomplishes the fitting-smoothing step for each annealing value More...
std::pair < TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface >	updateHits (const Trajectory vtraj, const SiTrackerMultiRecHitUpdator *updator, double annealing) const

Private Attributes
edm::ParameterSet	conf_
int	minHits_

Detailed Description

All is needed to run the deterministic annealing algorithm. Ported from ORCA.

Author

tropiano, genta in May 2014 by brondolin

Definition at line 31 of file DAFTrackProducerAlgorithm.h.

Member Typedef Documentation

typedef std::pair<Trajectory*, std::pair<reco::Track*,PropagationDirection> > DAFTrackProducerAlgorithm::AlgoProduct

private

Definition at line 33 of file DAFTrackProducerAlgorithm.h.

typedef std::vector< AlgoProduct > DAFTrackProducerAlgorithm::AlgoProductCollection

private

Definition at line 34 of file DAFTrackProducerAlgorithm.h.

typedef std::vector<TrajAnnealing> DAFTrackProducerAlgorithm::TrajAnnealingCollection

private

Definition at line 35 of file DAFTrackProducerAlgorithm.h.

Constructor & Destructor Documentation

DAFTrackProducerAlgorithm::DAFTrackProducerAlgorithm

(

const edm::ParameterSet &

conf

)

Definition at line 27 of file DAFTrackProducerAlgorithm.cc.

                                                                               :
  conf_(conf),
  minHits_(conf.getParameter<int>("MinHits")){}

DAFTrackProducerAlgorithm::~DAFTrackProducerAlgorithm ( )

inline

Definition at line 40 of file DAFTrackProducerAlgorithm.h.

{}

Member Function Documentation

bool DAFTrackProducerAlgorithm::buildTrack ( const Trajectory  vtraj, AlgoProductCollection &  algoResults, float  ndof, const reco::BeamSpot &  bs  ) const

private

Construct Tracks to be put in the event.

Definition at line 208 of file DAFTrackProducerAlgorithm.cc.

References alongMomentum, FreeTrajectoryState::charge(), Trajectory::chiSquared(), FreeTrajectoryState::curvilinearError(), Trajectory::direction(), Trajectory::firstMeasurement(), TrajectoryStateOnSurface::freeState(), TrajectoryStateClosestToBeamLine::isValid(), Trajectory::isValid(), Trajectory::lastMeasurement(), LogDebug, FreeTrajectoryState::momentum(), AlCaHLTBitMon_ParallelJobs::p, FreeTrajectoryState::position(), TrajectoryStateClosestToBeamLine::trackStateAtPCA(), TrajectoryMeasurement::updatedState(), findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by runWithCandidate().

{
  //variable declarations
  reco::Track * theTrack;
  Trajectory * theTraj; 
      
  LogDebug("DAFTrackProducerAlgorithm") <<" BUILDER " << std::endl;;
  TrajectoryStateOnSurface innertsos;
  
  if ( vtraj.isValid() ){

    theTraj = new Trajectory( vtraj );
    
    if (vtraj.direction() == alongMomentum) {
    //if (theTraj->direction() == oppositeToMomentum) {
      innertsos = vtraj.firstMeasurement().updatedState();
    } else { 
      innertsos = vtraj.lastMeasurement().updatedState();
    }
    
    TSCBLBuilderNoMaterial tscblBuilder;
    TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*(innertsos.freeState()),bs);

    if (tscbl.isValid()==false) return false;

    GlobalPoint v = tscbl.trackStateAtPCA().position();
    math::XYZPoint  pos( v.x(), v.y(), v.z() );
    GlobalVector p = tscbl.trackStateAtPCA().momentum();
    math::XYZVector mom( p.x(), p.y(), p.z() );

    //    LogDebug("TrackProducer") <<v<<p<<std::endl;

    theTrack = new reco::Track(vtraj.chiSquared(),
                   ndof, //in the DAF the ndof is not-integer
                   pos, mom, tscbl.trackStateAtPCA().charge(), 
                   tscbl.trackStateAtPCA().curvilinearError());

    AlgoProduct aProduct(theTraj,std::make_pair(theTrack,vtraj.direction()));
    algoResults.push_back(aProduct);

    return true;
  } 
  else  
    return false;
}

float DAFTrackProducerAlgorithm::calculateNdof ( const Trajectory  vtraj ) const

private

Definition at line 312 of file DAFTrackProducerAlgorithm.cc.

References makeMuonMisalignmentScenario::components, Trajectory::isValid(), getDQMSummary::iter, Trajectory::measurements(), SiTrackerMultiRecHit::recHits(), and SiTrackerMultiRecHit::weight().

Referenced by runWithCandidate().

{

  if (!vtraj.isValid()) return 0;
  float ndof = 0;
  const std::vector<TrajectoryMeasurement>& meas = vtraj.measurements();
  for (std::vector<TrajectoryMeasurement>::const_iterator iter = meas.begin(); iter != meas.end(); iter++){
    if (iter->recHit()->isValid()){
      SiTrackerMultiRecHit const & mHit = dynamic_cast<SiTrackerMultiRecHit const &>(*iter->recHit());
      std::vector<const TrackingRecHit*> components = mHit.recHits();
      int iComp = 0;
      for(std::vector<const TrackingRecHit*>::const_iterator iter2 = components.begin(); iter2 != components.end(); iter2++, iComp++){
        if ((*iter2)->isValid())
      ndof += ((*iter2)->dimension())*mHit.weight(iComp);
      }

    }
  }

  return ndof-5;

}

int DAFTrackProducerAlgorithm::checkHits ( Trajectory  iInitTraj, const Trajectory  iFinalTraj  ) const

private

Definition at line 335 of file DAFTrackProducerAlgorithm.cc.

References makeMuonMisalignmentScenario::components, TrackingRecHit::hit(), TrackingRecHit::isValid(), LogDebug, Trajectory::measurements(), TrajectoryMeasurement::recHit(), SiTrackerMultiRecHit::recHits(), histoStyle::weight, and SiTrackerMultiRecHit::weight().

Referenced by runWithCandidate().

                                                                                                 {

  std::vector<TrajectoryMeasurement> initmeasurements = iInitTraj.measurements();
  std::vector<TrajectoryMeasurement> finalmeasurements = iFinalTraj.measurements();
  std::vector<TrajectoryMeasurement>::iterator imeas;
  std::vector<TrajectoryMeasurement>::iterator jmeas;
  int nSame = 0;
  int ihit = 0;

  if( initmeasurements.empty() || finalmeasurements.empty() ){
    LogDebug("DAFTrackProducerAlgorithm") << "Initial or Final Trajectory empty.";
    return 0;
  }

  if( initmeasurements.size() != finalmeasurements.size() ) {
    LogDebug("DAFTrackProducerAlgorithm") << "Initial and Final Trajectory have different size.";
    return 0;
  }
          
  for(jmeas = initmeasurements.begin(); jmeas != initmeasurements.end(); jmeas++){

    const TrackingRecHit* initHit = jmeas->recHit()->hit();

    if(!initHit->isValid() && ihit == 0 ) continue;

    if(initHit->isValid()){

      TrajectoryMeasurement imeas = finalmeasurements.at(ihit);
      const TrackingRecHit* finalHit = imeas.recHit()->hit();
      const TrackingRecHit* MaxWeightHit=0;
      float maxweight = 0;
  
      const SiTrackerMultiRecHit* mrh = dynamic_cast<const SiTrackerMultiRecHit*>(finalHit);
      if (mrh){
        std::vector<const TrackingRecHit*> components = mrh->recHits();
        std::vector<const TrackingRecHit*>::const_iterator icomp;
        int hitcounter=0;

        for (icomp = components.begin(); icomp != components.end(); icomp++) {
          if((*icomp)->isValid()) {
            double weight = mrh->weight(hitcounter);
            if(weight > maxweight) {
              MaxWeightHit = *icomp;
              maxweight = weight;
            }
          }
          hitcounter++;
        }

        auto myref1 = reinterpret_cast<const BaseTrackerRecHit *>(initHit)->firstClusterRef();
        auto myref2 = reinterpret_cast<const BaseTrackerRecHit *>(MaxWeightHit)->firstClusterRef();

        if( myref1 == myref2 ){
          nSame++;  
        }
      }
    } else {

      TrajectoryMeasurement imeas = finalmeasurements.at(ihit);
      const TrackingRecHit* finalHit = imeas.recHit()->hit();
      if(!finalHit->isValid()){
        nSame++;
      }
    }
    
      ihit++;
  }

  return nSame;
}

std::pair< TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface > DAFTrackProducerAlgorithm::collectHits ( const Trajectory  vtraj, const MultiRecHitCollector *  measurementCollector, const MeasurementTrackerEvent *  measTk  ) const

private

Definition at line 128 of file DAFTrackProducerAlgorithm.cc.

References getDQMSummary::iter, LogDebug, and MultiRecHitCollector::recHits().

Referenced by runWithCandidate().

{

  LogDebug("DAFTrackProducerAlgorithm") << "Calling DAFTrackProducerAlgorithm::collectHits";

  //getting the traj measurements from the MeasurementCollector 
  int nHits = 0;
  TransientTrackingRecHit::RecHitContainer hits;
  std::vector<TrajectoryMeasurement> collectedmeas = measurementCollector->recHits(vtraj, measTk);

  //if the MeasurementCollector is empty, make an "empty" pair 
  //else taking the collected measured hits and building the pair
  if( collectedmeas.empty() ) 
    return std::make_pair(TransientTrackingRecHit::RecHitContainer(), TrajectoryStateOnSurface());

  for(  std::vector<TrajectoryMeasurement>::const_iterator iter = collectedmeas.begin(); 
    iter!=collectedmeas.end(); iter++ ){

    nHits++;
    hits.push_back(iter->recHit());
  
  }
  

  //TrajectoryStateWithArbitraryError() == Creates a TrajectoryState with the same parameters 
  //     as the input one, but with "infinite" errors, i.e. errors so big that they don't
  //     bias a fit starting with this state. 
  //return std::make_pair(hits,TrajectoryStateWithArbitraryError()(collectedmeas.front().predictedState()));

  // I do not have to rescale the error because it is already rescaled in the fit code 
  TrajectoryStateOnSurface initialStateFromTrack = collectedmeas.front().predictedState();

  return std::make_pair(hits, initialStateFromTrack);

}

void DAFTrackProducerAlgorithm::filter ( const TrajectoryFitter *  fitter, std::vector< Trajectory > &  input, int  minhits, std::vector< Trajectory > &  output, const TransientTrackingRecHitBuilder *  builder  ) const

private

Definition at line 257 of file DAFTrackProducerAlgorithm.cc.

References makeMuonMisalignmentScenario::components, alignCSCRings::e, TrajectoryFitter::fit(), edm::OwnVector< T, P >::front(), AnalysisDataFormats_SUSYBSMObjects::hc, getDQMSummary::iter, LogDebug, TkClonerImpl::makeShared(), TrackingRecHit::missing, SiTrackerMultiRecHit::recHits(), and SiTrackerMultiRecHit::weight().

{
  if (input.empty()) return;

  int ngoodhits = 0;
  std::vector<TrajectoryMeasurement> vtm = input[0].measurements(); 
  TransientTrackingRecHit::RecHitContainer hits;

  //count the number of non-outlier and non-invalid hits    
  for (std::vector<TrajectoryMeasurement>::reverse_iterator tm=vtm.rbegin(); tm!=vtm.rend();tm++){
    //if the rechit is valid
    if (tm->recHit()->isValid()) {
      SiTrackerMultiRecHit const & mHit = dynamic_cast<SiTrackerMultiRecHit const &>(*tm->recHit());
      std::vector<const TrackingRecHit*> components = mHit.recHits();
      int iComp = 0;
      bool isGood = false;
      for(std::vector<const TrackingRecHit*>::const_iterator iter = components.begin(); iter != components.end(); iter++, iComp++){        
    //if there is at least one component with weight higher than 1e-6 then the hit is not an outlier
        if (mHit.weight(iComp)>1e-6) {ngoodhits++; iComp++; isGood = true; break;}
      }
      if (isGood) {
        TkClonerImpl hc = static_cast<TkTransientTrackingRecHitBuilder const *>(builder)->cloner();
        auto tempHit = hc.makeShared(tm->recHit(),tm->updatedState());
        hits.push_back(tempHit);
      }
      else hits.push_back(std::make_shared<InvalidTrackingRecHit>(*tm->recHit()->det(), TrackingRecHit::missing));
    } else {
      TkClonerImpl hc = static_cast<TkTransientTrackingRecHitBuilder const *>(builder)->cloner();
      auto tempHit = hc.makeShared(tm->recHit(),tm->updatedState());
      hits.push_back(tempHit);
    }
  }


    LogDebug("DAFTrackProducerAlgorithm") << "Original number of valid hits " << input[0].foundHits() << "; after filtering " << ngoodhits;
    //debug
    if (ngoodhits>input[0].foundHits()) edm::LogError("DAFTrackProducerAlgorithm") << "Something wrong: the number of good hits from DAFTrackProducerAlgorithm::filter " << ngoodhits << " is higher than the original one " << input[0].foundHits();
    
    if (ngoodhits < minhits) return;    

    TrajectoryStateOnSurface curstartingTSOS = input.front().lastMeasurement().updatedState();
    LogDebug("DAFTrackProducerAlgorithm") << "starting tsos for final refitting " << curstartingTSOS ;
        //curstartingTSOS.rescaleError(100);

    output = fitter->fit(TrajectorySeed(PTrajectoryStateOnDet(),
                                                BasicTrajectorySeed::recHitContainer(),
                                                input.front().seed().direction()),
                                hits,
                                TrajectoryStateWithArbitraryError()(curstartingTSOS));
    LogDebug("DAFTrackProducerAlgorithm") << "After filtering " << output.size() << " trajectories";

}

Trajectory DAFTrackProducerAlgorithm::fit ( const std::pair< TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface > &  hits, const TrajectoryFitter *  theFitter, Trajectory  vtraj  ) const

private

accomplishes the fitting-smoothing step for each annealing value

Definition at line 189 of file DAFTrackProducerAlgorithm.cc.

References TrajectorySeed::direction(), TrajectoryFitter::fitOne(), Trajectory::isValid(), LogDebug, and Trajectory::seed().

Referenced by runWithCandidate().

                                                            {

  //creating a new trajectory starting from the direction of the seed of the input one and the hits
  Trajectory newVec = theFitter->fitOne(TrajectorySeed(PTrajectoryStateOnDet(),
                                                                 BasicTrajectorySeed::recHitContainer(),
                                                                 vtraj.seed().direction()),
                                                                 hits.first, hits.second);

  if( newVec.isValid() )  return newVec; 
  else{
    LogDebug("DAFTrackProducerAlgorithm") << "Fit no valid.";
    return Trajectory();
  }

}

void DAFTrackProducerAlgorithm::runWithCandidate	(	const TrackingGeometry *	theG,
		const MagneticField *	theMF,
		const std::vector< Trajectory > &	theTrajectoryCollection,
		const MeasurementTrackerEvent *	measTk,
		const TrajectoryFitter *	theFitter,
		const TransientTrackingRecHitBuilder *	builder,
		const MultiRecHitCollector *	measurementTracker,
		const SiTrackerMultiRecHitUpdator *	updator,
		const reco::BeamSpot &	bs,
		AlgoProductCollection &	algoResults,
		TrajAnnealingCollection &	trajann,
		bool	TrajAnnSaving_
	)		const

Run the Final Fit taking TrackCandidates as input.

Definition at line 32 of file DAFTrackProducerAlgorithm.cc.

References buildTrack(), calculateNdof(), checkHits(), collectHits(), cont, TrajectorySeed::direction(), Trajectory::direction(), fit(), Trajectory::foundHits(), SiTrackerMultiRecHitUpdator::getAnnealingProgram(), Trajectory::isValid(), LogDebug, minHits_, convertSQLiteXML::ok, Trajectory::seed(), groupFilesInBlocks::temp, and updateHits().

Referenced by DAFTrackProducer::produce().

{
  LogDebug("DAFTrackProducerAlgorithm") << "Number of Trajectories: " << theTrajectoryCollection.size() << "\n";
  int cont = 0;

  //running on src trajectory collection
  for (std::vector<Trajectory>::const_iterator ivtraj = theTrajectoryCollection.begin(); 
       ivtraj != theTrajectoryCollection.end(); ivtraj++)
  {

    float ndof = 0;
    Trajectory currentTraj;

    //getting trajectory
    //no need to have std::vector<Trajectory> vtraj !
    if ( (*ivtraj).isValid() ){

      LogDebug("DAFTrackProducerAlgorithm") << "The trajectory is valid. \n";

      //getting the MultiRecHit collection and the trajectory with a first fit-smooth round
      std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> hits = 
                            collectHits(*ivtraj, measurementCollector, &*measTk);
      
      currentTraj = fit(hits, theFitter, *ivtraj);

      //starting the annealing program
      for (std::vector<double>::const_iterator ian = updator->getAnnealingProgram().begin(); 
           ian != updator->getAnnealingProgram().end(); ian++){

        if (currentTraj.isValid()){

      LogDebug("DAFTrackProducerAlgorithm") << "Seed direction is " << currentTraj.seed().direction() 
                                        << ".Traj direction is " << currentTraj.direction();

          //updating MultiRecHits and fit-smooth again 
      std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> curiterationhits = 
                                    updateHits(currentTraj,updator,*ian);
      currentTraj = fit(curiterationhits, theFitter, currentTraj);

      //saving trajectory for each annealing cycle ...
      if(TrajAnnSaving_){
            TrajAnnealing temp(currentTraj, *ian);
        trajann.push_back(temp);
          }

          LogDebug("DAFTrackProducerAlgorithm") << "done annealing value "  <<  (*ian) ;

    } 
        else break;


      } //end of annealing program

      LogDebug("DAFTrackProducerAlgorithm") << "Ended annealing program with " << (1.*checkHits(*ivtraj, currentTraj))/(1.*(*ivtraj).measurements().size())*100. << " unchanged." << std::endl;

      //computing the ndof keeping into account the weights
      ndof = calculateNdof(currentTraj);

      //checking if the trajectory has the minimum number of valid hits ( weight (>1e-6) )
      //in order to remove tracks with too many outliers.

      //std::vector<Trajectory> filtered;
      //filter(theFitter, vtraj, minHits_, filtered, builder);              

      if(currentTraj.foundHits() >= minHits_) {
      
        bool ok = buildTrack(currentTraj, algoResults, ndof, bs) ;
        if(ok) cont++;

      }
      else{
        LogDebug("DAFTrackProducerAlgorithm")  << "Rejecting trajectory with " 
                           << currentTraj.foundHits()<<" hits"; 
      }
    }
    else 
      LogDebug("DAFTrackProducerAlgorithm") << "Rejecting empty trajectory" << std::endl;

  } //end run on track collection

  LogDebug("DAFTrackProducerAlgorithm") << "Number of Tracks found: " << cont << "\n";

}

std::pair< TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface > DAFTrackProducerAlgorithm::updateHits ( const Trajectory  vtraj, const SiTrackerMultiRecHitUpdator *  updator, double  annealing  ) const

private

Definition at line 167 of file DAFTrackProducerAlgorithm.cc.

References Trajectory::measurements(), and SiTrackerMultiRecHitUpdator::update().

Referenced by runWithCandidate().

{
    TransientTrackingRecHit::RecHitContainer hits;
    std::vector<TrajectoryMeasurement> vmeas = vtraj.measurements();
        std::vector<TrajectoryMeasurement>::reverse_iterator imeas;

    //I run inversely on the trajectory obtained and update the state
        for (imeas = vmeas.rbegin(); imeas != vmeas.rend(); imeas++){
              TransientTrackingRecHit::RecHitPointer updated = updator->update(imeas->recHit(), 
                            imeas->updatedState(), annealing);
             hits.push_back(updated);
        }

    TrajectoryStateOnSurface updatedStateFromTrack = vmeas.back().predictedState();
    //updatedStateFromTrack.rescaleError(10);

        //return std::make_pair(hits,TrajectoryStateWithArbitraryError()(vmeas.back().updatedState()));
        return std::make_pair(hits,updatedStateFromTrack);
}

Member Data Documentation

edm::ParameterSet DAFTrackProducerAlgorithm::conf_

private

Definition at line 93 of file DAFTrackProducerAlgorithm.h.

int DAFTrackProducerAlgorithm::minHits_

private

Definition at line 94 of file DAFTrackProducerAlgorithm.h.

Referenced by runWithCandidate().