CMSSW Reference Manual

00001 #include "RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h"
00002 #include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h"
00003 #include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h"
00004 
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 #include "MagneticField/Engine/interface/MagneticField.h"
00008 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00009 
00010 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
00011 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00012 #include "DataFormats/TrackReco/interface/Track.h"
00013 
00014 #include "TrackingTools/PatternTools/interface/TrajectoryFitter.h"
00015 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00016 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00017 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00018 #include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h"
00019 //#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00020 #include "TrackingTools/PatternTools/interface/TrajectoryStateClosestToBeamLineBuilder.h"
00021 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
00022 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
00023 
00024 #include "Utilities/General/interface/CMSexception.h"
00025 
00026 void DAFTrackProducerAlgorithm::runWithCandidate(const TrackingGeometry * theG,
00027                                                  const MagneticField * theMF,
00028                                                  const TrackCandidateCollection& theTCCollection,
00029                                                  const TrajectoryFitter * theFitter,
00030                                                  const TransientTrackingRecHitBuilder* builder,
00031                                                  const MultiRecHitCollector* measurementCollector,
00032                                                  const SiTrackerMultiRecHitUpdator* updator,
00033                                                  const reco::BeamSpot& bs,
00034                                                  AlgoProductCollection& algoResults) const
00035 {
00036   edm::LogInfo("TrackProducer") << "Number of TrackCandidates: " << theTCCollection.size() << "\n";
00037 
00038   int cont = 0;
00039   for (TrackCandidateCollection::const_iterator i=theTCCollection.begin(); i!=theTCCollection.end();i++){
00040     
00041       
00042       const TrackCandidate * theTC = &(*i);
00043       PTrajectoryStateOnDet state = theTC->trajectoryStateOnDet();
00044       const TrackCandidate::range& recHitVec=theTC->recHits();
00045       const TrajectorySeed& seed = theTC->seed();
00046 
00047       TrajectoryStateTransform transformer;
00048   
00049       DetId  detId(state.detId());
00050       TrajectoryStateOnSurface theTSOS = transformer.transientState( state,
00051                                                                      &(theG->idToDet(detId)->surface()), 
00052                                                                      theMF);
00053 
00054       LogDebug("TrackProducer") << "Initial TSOS\n" << theTSOS << "\n";
00055       
00056       //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
00057       TransientTrackingRecHit::RecHitContainer hits;
00058       
00059       float ndof=0;
00060       LogDebug("DAFTrackProducerAlgorithm") << "about to build the rechits for the first round....";        
00061       for (edm::OwnVector<TrackingRecHit>::const_iterator i=recHitVec.first;
00062            i!=recHitVec.second; i++){
00063         hits.push_back(builder->build(&(*i) ));
00064       }
00065       LogDebug("DAFTrackProducerAlgorithm") << "done" << std::endl;
00066       //first of all do a fit-smooth round to get the trajectory
00067       LogDebug("DAFTrackProducerAlgorithm") << "About to build tha trajectory for the first round....";         
00068       //theMRHChi2Estimator->setAnnealingFactor(1);
00069       std::vector<Trajectory> vtraj = theFitter->fit(seed, hits, theTSOS);
00070       LogDebug("DAFTrackProducerAlgorithm") << "done" << std::endl;     
00071       LogDebug("DAFTrackProducerAlgorithm") << "after the first round found " << vtraj.size() << " trajectories"  << std::endl; 
00072 
00073       if (vtraj.size() != 0){
00074         //bool isFirstIteration=true;
00075         std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> hits = collectHits(vtraj, measurementCollector);
00076         fit(hits, theFitter, vtraj);
00077 
00078         for (std::vector<double>::const_iterator ian = updator->getAnnealingProgram().begin(); ian != updator->getAnnealingProgram().end(); ian++){
00079                 if (vtraj.size()){
00080                         LogDebug("DAFTrackProducerAlgorithm") << "Seed direction is " << vtraj.front().seed().direction() 
00081                                                               << "  Traj direction is " << vtraj.front().direction(); 
00082                         std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> curiterationhits = updateHits(vtraj,updator,*ian);
00083                         fit(curiterationhits, theFitter, vtraj);
00084                                 
00085                         LogDebug("DAFTrackProducerAlgorithm") << "done annealing value "  <<  (*ian) << " with " << vtraj.size() << " trajectories";
00086                 } else break;
00087                 LogDebug("DAFTrackProducerAlgorithm") << "number of trajectories after annealing value " << (*ian) << " annealing step " << vtraj.size();
00088         }
00089         LogDebug("DAFTrackProducerAlgorithm") << "Ended annealing program with " << vtraj.size() << " trajectories" << std::endl;
00090 
00091         //check the trajectory to see that the number of valid hits with 
00092         //reasonable weight (>1e-6) is higher than the minimum set in the DAFTrackProducer.
00093         //This is a kind of filter to remove tracks with too many outliers.
00094         //Outliers are mrhits with all components with weight less than 1e-6 
00095   
00096         std::vector<Trajectory> filtered;
00097         filter(theFitter, vtraj, conf_.getParameter<int>("MinHits"), filtered);                         
00098         ndof = calculateNdof(filtered); 
00099         bool ok = buildTrack(filtered, algoResults, ndof, bs);
00100         if(ok) cont++;
00101 
00102       }
00103   }
00104   edm::LogInfo("TrackProducer") << "Number of Tracks found: " << cont << "\n";
00105 }
00106 
00107 std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> 
00108 DAFTrackProducerAlgorithm::collectHits(const std::vector<Trajectory>& vtraj, 
00109                                        const MultiRecHitCollector* measurementCollector) const{
00110         TransientTrackingRecHit::RecHitContainer hits;
00111         std::vector<TrajectoryMeasurement> collectedmeas = measurementCollector->recHits(vtraj.front());
00112         if (collectedmeas.empty()) return std::make_pair(TransientTrackingRecHit::RecHitContainer(), TrajectoryStateOnSurface());
00113         for (std::vector<TrajectoryMeasurement>::const_iterator iter = collectedmeas.begin(); iter!=collectedmeas.end(); iter++){
00114                 hits.push_back(iter->recHit());
00115         }
00116         return std::make_pair(hits,TrajectoryStateWithArbitraryError()(collectedmeas.front().predictedState()));        
00117 }
00118 
00119 std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface>
00120 DAFTrackProducerAlgorithm::updateHits(const std::vector<Trajectory>& vtraj,
00121                                       const SiTrackerMultiRecHitUpdator* updator,
00122                                       double annealing) const {
00123         TransientTrackingRecHit::RecHitContainer hits;
00124         std::vector<TrajectoryMeasurement> vmeas = vtraj.front().measurements();
00125         std::vector<TrajectoryMeasurement>::reverse_iterator imeas;
00126         for (imeas = vmeas.rbegin(); imeas != vmeas.rend(); imeas++){
00127                 TransientTrackingRecHit::RecHitPointer updated = updator->update(imeas->recHit(), imeas->updatedState(), annealing);
00128                 hits.push_back(updated);
00129         }
00130         return std::make_pair(hits,TrajectoryStateWithArbitraryError()(vmeas.back().updatedState()));
00131 }
00132 
00133 void DAFTrackProducerAlgorithm::fit(const std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface>& hits, 
00134                                     const TrajectoryFitter * theFitter,
00135                                     std::vector<Trajectory>& vtraj) const {
00136         std::vector<Trajectory> newVec = theFitter->fit(TrajectorySeed(PTrajectoryStateOnDet(),
00137                                                                        BasicTrajectorySeed::recHitContainer(),
00138                                                                        vtraj.front().seed().direction()),
00139                                                         hits.first,
00140                                                         hits.second);
00141         vtraj.reserve(newVec.size());
00142         vtraj.swap(newVec);
00143         LogDebug("DAFTrackProducerAlgorithm") << "swapped!" << std::endl;
00144 }
00145 
00146 
00147 bool DAFTrackProducerAlgorithm::buildTrack (const std::vector<Trajectory>& vtraj,
00148                                             AlgoProductCollection& algoResults,
00149                                             float ndof,
00150                                             const reco::BeamSpot& bs) const{
00151   //variable declarations
00152   reco::Track * theTrack;
00153   Trajectory * theTraj; 
00154       
00155   //perform the fit: the result's size is 1 if it succeded, 0 if fails
00156   
00157   
00158   //LogDebug("TrackProducer") <<" FITTER FOUND "<< vtraj.size() << " TRAJECTORIES" <<"\n";
00159   LogDebug("DAFTrackProducerAlgorithm") <<" FITTER FOUND "<< vtraj.size() << " TRAJECTORIES" << std::endl;;
00160   TrajectoryStateOnSurface innertsos;
00161   
00162   if (vtraj.size() != 0){
00163 
00164     theTraj = new Trajectory( vtraj.front() );
00165     
00166     if (theTraj->direction() == alongMomentum) {
00167     //if (theTraj->direction() == oppositeToMomentum) {
00168       innertsos = theTraj->firstMeasurement().updatedState();
00169       //std::cout << "Inner momentum " << innertsos.globalParameters().momentum().mag() << std::endl;   
00170     } else { 
00171       innertsos = theTraj->lastMeasurement().updatedState();
00172     }
00173     
00174     
00175     TrajectoryStateClosestToBeamLineBuilder tscblBuilder;
00176     TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*(innertsos.freeState()),bs);
00177 
00178     if (tscbl.isValid()==false) return false;
00179 
00180     GlobalPoint v = tscbl.trackStateAtPCA().position();
00181     math::XYZPoint  pos( v.x(), v.y(), v.z() );
00182     GlobalVector p = tscbl.trackStateAtPCA().momentum();
00183     math::XYZVector mom( p.x(), p.y(), p.z() );
00184 
00185     LogDebug("TrackProducer") <<v<<p<<std::endl;
00186 
00187     theTrack = new reco::Track(theTraj->chiSquared(),
00188                                ndof, //in the DAF the ndof is not-integer
00189                                pos, mom, tscbl.trackStateAtPCA().charge(), tscbl.trackStateAtPCA().curvilinearError());
00190 
00191 
00192     LogDebug("TrackProducer") <<"track done\n";
00193 
00194     AlgoProduct aProduct(theTraj,std::make_pair(theTrack,theTraj->direction()));
00195     LogDebug("TrackProducer") <<"track done1\n";
00196     algoResults.push_back(aProduct);
00197     LogDebug("TrackProducer") <<"track done2\n";
00198     
00199     return true;
00200   } 
00201   else  return false;
00202 }
00203 
00204 void  DAFTrackProducerAlgorithm::filter(const TrajectoryFitter* fitter, std::vector<Trajectory>& input, int minhits, std::vector<Trajectory>& output) const {
00205         if (input.empty()) return;
00206 
00207         int ngoodhits = 0;
00208 
00209         std::vector<TrajectoryMeasurement> vtm = input[0].measurements();       
00210 
00211         TransientTrackingRecHit::RecHitContainer hits;
00212 
00213         //count the number of non-outlier and non-invalid hits  
00214         for (std::vector<TrajectoryMeasurement>::reverse_iterator tm=vtm.rbegin(); tm!=vtm.rend();tm++){
00215                 //if the rechit is valid
00216                 if (tm->recHit()->isValid()) {
00217                         TransientTrackingRecHit::ConstRecHitContainer components = tm->recHit()->transientHits();
00218                         bool isGood = false;
00219                         for (TransientTrackingRecHit::ConstRecHitContainer::iterator rechit = components.begin(); rechit != components.end(); rechit++){
00220                                 //if there is at least one component with weight higher than 1e-6 then the hit is not an outlier
00221                                 if ((*rechit)->weight()>1e-6) {ngoodhits++; isGood = true; break;}
00222                         }
00223                         if (isGood) hits.push_back(tm->recHit()->clone(tm->updatedState()));
00224                         else hits.push_back(InvalidTransientRecHit::build(tm->recHit()->det(), TrackingRecHit::missing));
00225                 } else {
00226                         hits.push_back(tm->recHit()->clone(tm->updatedState()));
00227                 }
00228         }
00229 
00230 
00231         LogDebug("DAFTrackProducerAlgorithm") << "Original number of valid hits " << input[0].foundHits() << "; after filtering " << ngoodhits;
00232         //debug
00233         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();
00234         
00235         if (ngoodhits < minhits) return;        
00236 
00237         TrajectoryStateOnSurface curstartingTSOS = input.front().lastMeasurement().updatedState();
00238         LogDebug("DAFTrackProducerAlgorithm") << "starting tsos for final refitting " << curstartingTSOS ;
00239         //curstartingTSOS.rescaleError(100);
00240 
00241         output = fitter->fit(TrajectorySeed(PTrajectoryStateOnDet(),
00242                                                 BasicTrajectorySeed::recHitContainer(),
00243                                                 input.front().seed().direction()),
00244                                 hits,
00245                                 TrajectoryStateWithArbitraryError()(curstartingTSOS));
00246         LogDebug("DAFTrackProducerAlgorithm") << "After filtering " << output.size() << " trajectories";
00247 
00248 }
00249 
00250 float DAFTrackProducerAlgorithm::calculateNdof(const std::vector<Trajectory>& vtraj) const {
00251         if (vtraj.empty()) return 0;
00252         float ndof = 0;
00253         const std::vector<TrajectoryMeasurement>& meas = vtraj.front().measurements();
00254         for (std::vector<TrajectoryMeasurement>::const_iterator iter = meas.begin(); iter != meas.end(); iter++){
00255                 if (iter->recHit()->isValid()){
00256                         TransientTrackingRecHit::ConstRecHitContainer components = iter->recHit()->transientHits();
00257                         TransientTrackingRecHit::ConstRecHitContainer::const_iterator iter2;
00258                         for (iter2 = components.begin(); iter2 != components.end(); iter2++){
00259                                 if ((*iter2)->isValid())ndof+=((*iter2)->dimension())*(*iter2)->weight();
00260                         }
00261                 }
00262         }
00263         return ndof-5;
00264 }
DAFTrackProducerAlgorithm.cc
