Reference Manual

Go to the documentation of this file.
00001 #include "Validation/RecoTrack/interface/MultiTrackValidator.h"
00002 #include "DQMServices/ClientConfig/interface/FitSlicesYTool.h"
00003 
00004 #include "FWCore/Framework/interface/MakerMacros.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00009 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00010 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00011 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00012 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00013 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByChi2.h"
00014 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"
00015 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00016 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00018 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00019 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
00020 #include "SimTracker/TrackAssociation/plugins/ParametersDefinerForTPESProducer.h"
00021 #include "SimTracker/TrackAssociation/plugins/CosmicParametersDefinerForTPESProducer.h"
00022 #include "Validation/RecoTrack/interface/MTVHistoProducerAlgoFactory.h"
00023 
00024 #include "DataFormats/TrackReco/interface/DeDxData.h"
00025 #include "DataFormats/Common/interface/ValueMap.h"
00026 #include "DataFormats/Common/interface/Ref.h"
00027 
00028 #include "TMath.h"
00029 #include <TF1.h>
00030 
00031 using namespace std;
00032 using namespace edm;
00033 
00034 MultiTrackValidator::MultiTrackValidator(const edm::ParameterSet& pset):MultiTrackValidatorBase(pset){
00035   //theExtractor = IsoDepositExtractorFactory::get()->create( extractorName, extractorPSet);
00036 
00037   ParameterSet psetForHistoProducerAlgo = pset.getParameter<ParameterSet>("histoProducerAlgoBlock");
00038   string histoProducerAlgoName = psetForHistoProducerAlgo.getParameter<string>("ComponentName");
00039   histoProducerAlgo_ = MTVHistoProducerAlgoFactory::get()->create(histoProducerAlgoName ,psetForHistoProducerAlgo);
00040   histoProducerAlgo_->setDQMStore(dbe_);
00041 
00042   dirName_ = pset.getParameter<std::string>("dirName");
00043   associatormap = pset.getParameter< edm::InputTag >("associatormap");
00044   UseAssociators = pset.getParameter< bool >("UseAssociators");
00045 
00046   m_dEdx1Tag = pset.getParameter< edm::InputTag >("dEdx1Tag");
00047   m_dEdx2Tag = pset.getParameter< edm::InputTag >("dEdx2Tag");
00048 
00049   tpSelector = TrackingParticleSelector(pset.getParameter<double>("ptMinTP"),
00050                                         pset.getParameter<double>("minRapidityTP"),
00051                                         pset.getParameter<double>("maxRapidityTP"),
00052                                         pset.getParameter<double>("tipTP"),
00053                                         pset.getParameter<double>("lipTP"),
00054                                         pset.getParameter<int>("minHitTP"),
00055                                         pset.getParameter<bool>("signalOnlyTP"),
00056                                         pset.getParameter<bool>("chargedOnlyTP"),
00057                                         pset.getParameter<bool>("stableOnlyTP"),
00058                                         pset.getParameter<std::vector<int> >("pdgIdTP"));
00059 
00060   cosmictpSelector = CosmicTrackingParticleSelector(pset.getParameter<double>("ptMinTP"),
00061                                                     pset.getParameter<double>("minRapidityTP"),
00062                                                     pset.getParameter<double>("maxRapidityTP"),
00063                                                     pset.getParameter<double>("tipTP"),
00064                                                     pset.getParameter<double>("lipTP"),
00065                                                     pset.getParameter<int>("minHitTP"),
00066                                                     pset.getParameter<bool>("chargedOnlyTP"),
00067                                                     pset.getParameter<std::vector<int> >("pdgIdTP"));
00068 
00069   useGsf = pset.getParameter<bool>("useGsf");
00070   runStandalone = pset.getParameter<bool>("runStandalone");
00071 
00072 
00073     
00074   if (!UseAssociators) {
00075     associators.clear();
00076     associators.push_back(associatormap.label());
00077   }
00078 }
00079 
00080 
00081 MultiTrackValidator::~MultiTrackValidator(){delete histoProducerAlgo_;}
00082 
00083 void MultiTrackValidator::beginRun(Run const&, EventSetup const& setup) {
00084   //  dbe_->showDirStructure();
00085 
00086   //int j=0;  //is This Necessary ???
00087   for (unsigned int ww=0;ww<associators.size();ww++){
00088     for (unsigned int www=0;www<label.size();www++){
00089       dbe_->cd();
00090       InputTag algo = label[www];
00091       string dirName=dirName_;
00092       if (algo.process()!="")
00093         dirName+=algo.process()+"_";
00094       if(algo.label()!="")
00095         dirName+=algo.label()+"_";
00096       if(algo.instance()!="")
00097         dirName+=algo.instance()+"_";      
00098       if (dirName.find("Tracks")<dirName.length()){
00099         dirName.replace(dirName.find("Tracks"),6,"");
00100       }
00101       string assoc= associators[ww];
00102       if (assoc.find("Track")<assoc.length()){
00103         assoc.replace(assoc.find("Track"),5,"");
00104       }
00105       dirName+=assoc;
00106       std::replace(dirName.begin(), dirName.end(), ':', '_');
00107 
00108       dbe_->setCurrentFolder(dirName.c_str());
00109       
00110       // vector of vector initialization
00111       histoProducerAlgo_->initialize(); //TO BE FIXED. I'D LIKE TO AVOID THIS CALL
00112 
00113       dbe_->goUp(); //Is this really necessary ???
00114       string subDirName = dirName + "/simulation";
00115       dbe_->setCurrentFolder(subDirName.c_str());
00116 
00117       //Booking histograms concerning with simulated tracks
00118       histoProducerAlgo_->bookSimHistos();
00119 
00120       dbe_->cd();
00121       dbe_->setCurrentFolder(dirName.c_str());
00122 
00123       //Booking histograms concerning with reconstructed tracks
00124       histoProducerAlgo_->bookRecoHistos();
00125       if (runStandalone) histoProducerAlgo_->bookRecoHistosForStandaloneRunning();
00126 
00127       if (UseAssociators) {
00128         edm::ESHandle<TrackAssociatorBase> theAssociator;
00129         for (unsigned int w=0;w<associators.size();w++) {
00130           setup.get<TrackAssociatorRecord>().get(associators[w],theAssociator);
00131           associator.push_back( theAssociator.product() );
00132         }//end loop w
00133       }
00134     }//end loop www
00135   }// end loop ww
00136 }
00137 
00138 void MultiTrackValidator::analyze(const edm::Event& event, const edm::EventSetup& setup){
00139   using namespace reco;
00140   
00141   edm::LogInfo("TrackValidator") << "\n====================================================" << "\n"
00142                                  << "Analyzing new event" << "\n"
00143                                  << "====================================================\n" << "\n";
00144   edm::ESHandle<ParametersDefinerForTP> parametersDefinerTP; 
00145   setup.get<TrackAssociatorRecord>().get(parametersDefiner,parametersDefinerTP);    
00146   
00147   edm::Handle<TrackingParticleCollection>  TPCollectionHeff ;
00148   event.getByLabel(label_tp_effic,TPCollectionHeff);
00149   const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
00150   
00151   edm::Handle<TrackingParticleCollection>  TPCollectionHfake ;
00152   event.getByLabel(label_tp_fake,TPCollectionHfake);
00153   const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
00154   
00155   //if (tPCeff.size()==0) {edm::LogInfo("TrackValidator") 
00156   //<< "TP Collection for efficiency studies has size = 0! Skipping Event." ; return;}
00157   //if (tPCfake.size()==0) {edm::LogInfo("TrackValidator") 
00158   //<< "TP Collection for fake rate studies has size = 0! Skipping Event." ; return;}
00159   
00160   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00161   event.getByLabel(bsSrc,recoBeamSpotHandle);
00162   reco::BeamSpot bs = *recoBeamSpotHandle;      
00163   
00164   int w=0; //counter counting the number of sets of histograms
00165   for (unsigned int ww=0;ww<associators.size();ww++){
00166     for (unsigned int www=0;www<label.size();www++){
00167       //
00168       //get collections from the event
00169       //
00170       edm::Handle<View<Track> >  trackCollection;
00171       if(!event.getByLabel(label[www], trackCollection)&&ignoremissingtkcollection_)continue;
00172       //if (trackCollection->size()==0) {
00173       //edm::LogInfo("TrackValidator") << "TrackCollection size = 0!" ; 
00174       //continue;
00175       //}
00176       reco::RecoToSimCollection recSimColl;
00177       reco::SimToRecoCollection simRecColl;
00178       //associate tracks
00179       if(UseAssociators){
00180         edm::LogVerbatim("TrackValidator") << "Analyzing " 
00181                                            << label[www].process()<<":"
00182                                            << label[www].label()<<":"
00183                                            << label[www].instance()<<" with "
00184                                            << associators[ww].c_str() <<"\n";
00185         
00186         LogTrace("TrackValidator") << "Calling associateRecoToSim method" << "\n";
00187         recSimColl=associator[ww]->associateRecoToSim(trackCollection,
00188                                                       TPCollectionHfake,
00189                                                       &event);
00190         LogTrace("TrackValidator") << "Calling associateSimToReco method" << "\n";
00191         simRecColl=associator[ww]->associateSimToReco(trackCollection,
00192                                                       TPCollectionHeff, 
00193                                                       &event);
00194       }
00195       else{
00196         edm::LogVerbatim("TrackValidator") << "Analyzing " 
00197                                            << label[www].process()<<":"
00198                                            << label[www].label()<<":"
00199                                            << label[www].instance()<<" with "
00200                                            << associatormap.process()<<":"
00201                                            << associatormap.label()<<":"
00202                                            << associatormap.instance()<<"\n";
00203         
00204         Handle<reco::SimToRecoCollection > simtorecoCollectionH;
00205         event.getByLabel(associatormap,simtorecoCollectionH);
00206         simRecColl= *(simtorecoCollectionH.product()); 
00207         
00208         Handle<reco::RecoToSimCollection > recotosimCollectionH;
00209         event.getByLabel(associatormap,recotosimCollectionH);
00210         recSimColl= *(recotosimCollectionH.product()); 
00211       }
00212 
00213 
00214 
00215       // ########################################################
00216       // fill simulation histograms (LOOP OVER TRACKINGPARTICLES)
00217       // ########################################################
00218 
00219       //compute number of tracks per eta interval
00220       //
00221       edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles: " << tPCeff.size() << "\n";
00222       int ats(0);  //This counter counts the number of simTracks that are "associated" to recoTracks
00223       int st(0);   //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
00224       for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){ //loop over TPs collection for tracking efficiency
00225         TrackingParticleRef tpr(TPCollectionHeff, i);
00226         TrackingParticle* tp=const_cast<TrackingParticle*>(tpr.get());
00227         ParticleBase::Vector momentumTP; 
00228         ParticleBase::Point vertexTP;
00229         double dxySim(0);
00230         double dzSim(0); 
00231 
00232         
00233         //---------- THIS PART HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
00234         //If the TrackingParticle is collison like, get the momentum and vertex at production state
00235         if(parametersDefiner=="LhcParametersDefinerForTP")
00236           {
00237             if(! tpSelector(*tp)) continue;
00238             momentumTP = tp->momentum();
00239             vertexTP = tp->vertex();
00240             //Calcualte the impact parameters w.r.t. PCA
00241             ParticleBase::Vector momentum = parametersDefinerTP->momentum(event,setup,*tp);
00242             ParticleBase::Point vertex = parametersDefinerTP->vertex(event,setup,*tp);
00243             dxySim = (-vertex.x()*sin(momentum.phi())+vertex.y()*cos(momentum.phi()));
00244             dzSim = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y())/sqrt(momentum.perp2()) 
00245               * momentum.z()/sqrt(momentum.perp2());
00246           }
00247         //If the TrackingParticle is comics, get the momentum and vertex at PCA
00248         if(parametersDefiner=="CosmicParametersDefinerForTP")
00249           {
00250             if(! cosmictpSelector(*tp,&bs,event,setup)) continue;       
00251             momentumTP = parametersDefinerTP->momentum(event,setup,*tp);
00252             vertexTP = parametersDefinerTP->vertex(event,setup,*tp);
00253             dxySim = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
00254             dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y())/sqrt(momentumTP.perp2()) 
00255               * momentumTP.z()/sqrt(momentumTP.perp2());
00256           }
00257         //---------- THE PART ABOVE HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
00258 
00259         st++;   //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
00260 
00261         // in the coming lines, histos are filled using as input
00262         // - momentumTP 
00263         // - vertexTP 
00264         // - dxySim
00265         // - dzSim
00266 
00267         histoProducerAlgo_->fill_generic_simTrack_histos(w,momentumTP,vertexTP);
00268 
00269 
00270         // ##############################################
00271         // fill RecoAssociated SimTracks' histograms
00272         // ##############################################
00273         bool isRecoMatched(false);
00274         const reco::Track* matchedTrackPointer=0;
00275         std::vector<std::pair<RefToBase<Track>, double> > rt;
00276         if(simRecColl.find(tpr) != simRecColl.end()){
00277           rt = (std::vector<std::pair<RefToBase<Track>, double> >) simRecColl[tpr];
00278           if (rt.size()!=0) {
00279             ats++; //This counter counts the number of simTracks that have a recoTrack associated
00280             isRecoMatched = true; matchedTrackPointer = rt.begin()->first.get();
00281             edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st 
00282                                                << " with pt=" << sqrt(momentumTP.perp2()) 
00283                                                << " associated with quality:" << rt.begin()->second <<"\n";
00284           }
00285         }else{
00286           edm::LogVerbatim("TrackValidator") 
00287             << "TrackingParticle #" << st
00288             << " with pt,eta,phi: " 
00289             << sqrt(momentumTP.perp2()) << " , "
00290             << momentumTP.eta() << " , "
00291             << momentumTP.phi() << " , "
00292             << " NOT associated to any reco::Track" << "\n";
00293         }
00294         
00295 
00296         
00297 
00298         std::vector<PSimHit> simhits=tp->trackPSimHit(DetId::Tracker);
00299         int nSimHits = simhits.end()-simhits.begin();
00300 
00301         histoProducerAlgo_->fill_recoAssociated_simTrack_histos(w,*tp,momentumTP,vertexTP,dxySim,dzSim,nSimHits,matchedTrackPointer);
00302 
00303 
00304 
00305 
00306       } // End  for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00307 
00308       //if (st!=0) h_tracksSIM[w]->Fill(st);  // TO BE FIXED
00309       
00310       
00311       // ##############################################
00312       // fill recoTracks histograms (LOOP OVER TRACKS)
00313       // ##############################################
00314       edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with "
00315                                          << label[www].process()<<":"
00316                                          << label[www].label()<<":"
00317                                          << label[www].instance()
00318                                          << ": " << trackCollection->size() << "\n";
00319 
00320       int at(0); //This counter counts the number of recoTracks that are associated to SimTracks
00321       int rT(0); //This counter counts the number of recoTracks in general
00322 
00323 
00324       // dE/dx
00325       // at some point this could be generalized, with a vector of tags and a corresponding vector of Handles
00326       // I'm writing the interface such to take vectors of ValueMaps
00327       edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx1Handle;
00328       edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx2Handle;
00329       std::vector<edm::ValueMap<reco::DeDxData> > v_dEdx;
00330       v_dEdx.clear();
00331       //std::cout << "PIPPO: label is " << label[www] << std::endl;
00332       if (label[www].label()=="generalTracks") {
00333         try {
00334           event.getByLabel(m_dEdx1Tag, dEdx1Handle);
00335           const edm::ValueMap<reco::DeDxData> dEdx1 = *dEdx1Handle.product();
00336           event.getByLabel(m_dEdx2Tag, dEdx2Handle);
00337           const edm::ValueMap<reco::DeDxData> dEdx2 = *dEdx2Handle.product();
00338           v_dEdx.push_back(dEdx1);
00339           v_dEdx.push_back(dEdx2);
00340         } catch (cms::Exception e){
00341           LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
00342         }
00343       }
00344       //end dE/dx
00345 
00346       for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
00347         RefToBase<Track> track(trackCollection, i);
00348         rT++;
00349         
00350         bool isSimMatched(false);
00351         std::vector<std::pair<TrackingParticleRef, double> > tp;
00352         if(recSimColl.find(track) != recSimColl.end()){
00353           tp = recSimColl[track];
00354           if (tp.size()!=0) {
00355             isSimMatched = true;
00356             at++;
00357             edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt() 
00358                                                << " associated with quality:" << tp.begin()->second <<"\n";
00359           }
00360         } else {
00361           edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00362                                              << " NOT associated to any TrackingParticle" << "\n";                
00363         }
00364         
00365 
00366         histoProducerAlgo_->fill_generic_recoTrack_histos(w,*track,bs.position(),isSimMatched);
00367 
00368         // dE/dx
00369         //      reco::TrackRef track2  = reco::TrackRef( trackCollection, i );
00370         if (v_dEdx.size() > 0) histoProducerAlgo_->fill_dedx_recoTrack_histos(w,track, v_dEdx);
00371         //if (v_dEdx.size() > 0) histoProducerAlgo_->fill_dedx_recoTrack_histos(track2, v_dEdx);
00372 
00373 
00374         //Fill other histos
00375         //try{ //Is this really necessary ????
00376 
00377         if (tp.size()==0) continue;     
00378 
00379         histoProducerAlgo_->fill_simAssociated_recoTrack_histos(w,*track);
00380 
00381         TrackingParticleRef tpr = tp.begin()->first;
00382         
00383         /* TO BE FIXED LATER
00384         if (associators[ww]=="TrackAssociatorByChi2"){
00385           //association chi2
00386           double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
00387           h_assochi2[www]->Fill(assocChi2);
00388           h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
00389         }
00390         else if (associators[ww]=="TrackAssociatorByHits"){
00391           double fraction = tp.begin()->second;
00392           h_assocFraction[www]->Fill(fraction);
00393           h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
00394         }
00395         */
00396 
00397           
00398         //Get tracking particle parameters at point of closest approach to the beamline
00399         ParticleBase::Vector momentumTP = parametersDefinerTP->momentum(event,setup,*(tpr.get()));
00400         ParticleBase::Point vertexTP = parametersDefinerTP->vertex(event,setup,*(tpr.get()));                    
00401         int chargeTP = tpr->charge();
00402 
00403         histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(w,momentumTP,vertexTP,chargeTP,
00404                                                              *track,bs.position());
00405         
00406         
00407         //TO BE FIXED
00408         //std::vector<PSimHit> simhits=tpr.get()->trackPSimHit(DetId::Tracker);
00409         //nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
00410         
00411         /*
00412           } // End of try{
00413           catch (cms::Exception e){
00414           LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
00415           }
00416         */
00417         
00418       } // End of for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
00419 
00420       //TO BE FIXED
00421       //if (at!=0) h_tracks[w]->Fill(at);
00422       //h_fakes[w]->Fill(rT-at);
00423       //nrec_vs_nsim[w]->Fill(rT,st);
00424 
00425       edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
00426                                          << "Total Associated (simToReco): " << ats << "\n"
00427                                          << "Total Reconstructed: " << rT << "\n"
00428                                          << "Total Associated (recoToSim): " << at << "\n"
00429                                          << "Total Fakes: " << rT-at << "\n";
00430 
00431       w++;
00432     } // End of  for (unsigned int www=0;www<label.size();www++){
00433   } //END of for (unsigned int ww=0;ww<associators.size();ww++){
00434 }
00435 
00436 void MultiTrackValidator::endRun(Run const&, EventSetup const&) {
00437   int w=0;
00438   for (unsigned int ww=0;ww<associators.size();ww++){
00439     for (unsigned int www=0;www<label.size();www++){
00440       if(!skipHistoFit && runStandalone)        histoProducerAlgo_->finalHistoFits(w);
00441       if (runStandalone) histoProducerAlgo_->fillProfileHistosFromVectors(w);
00442       histoProducerAlgo_->fillHistosFromVectors(w);
00443       w++;
00444     }    
00445   }
00446   if ( out.size() != 0 && dbe_ ) dbe_->save(out);
00447 }
00448 
00449 
00450
CMSSW_4_4_3_patch1/src/Validation/RecoTrack/plugins/MultiTrackValidator.cc
