Reference Manual

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