Reference Manual

Go to the documentation of this file.
00001 #include "Validation/RecoTrack/interface/MultiTrackValidatorGenPs.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/QuickTrackAssociatorByHits.h"
00015 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00016 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00018 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
00019 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
00020 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
00021 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
00022 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00023 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
00024 #include "SimTracker/TrackAssociation/plugins/ParametersDefinerForTPESProducer.h"
00025 #include "SimTracker/TrackAssociation/plugins/CosmicParametersDefinerForTPESProducer.h"
00026 #include "Validation/RecoTrack/interface/MTVHistoProducerAlgoFactory.h"
00027 
00028 #include "DataFormats/TrackReco/interface/DeDxData.h"
00029 #include "DataFormats/Common/interface/ValueMap.h"
00030 #include "DataFormats/Common/interface/Ref.h"
00031 
00032 #include "TMath.h"
00033 #include <TF1.h>
00034 
00035 //#include <iostream>
00036 
00037 using namespace std;
00038 using namespace edm;
00039 
00040 MultiTrackValidatorGenPs::MultiTrackValidatorGenPs(const edm::ParameterSet& pset):MultiTrackValidator(pset){
00041 
00042   gpSelector = GenParticleCustomSelector(pset.getParameter<double>("ptMinGP"),
00043                                          pset.getParameter<double>("minRapidityGP"),
00044                                          pset.getParameter<double>("maxRapidityGP"),
00045                                          pset.getParameter<double>("tipGP"),
00046                                          pset.getParameter<double>("lipGP"),
00047                                          pset.getParameter<bool>("chargedOnlyGP"),
00048                                          pset.getParameter<int>("statusGP"),
00049                                          pset.getParameter<std::vector<int> >("pdgIdGP"));
00050 
00051 }
00052 
00053 MultiTrackValidatorGenPs::~MultiTrackValidatorGenPs(){}
00054 
00055 
00056 void MultiTrackValidatorGenPs::analyze(const edm::Event& event, const edm::EventSetup& setup){
00057   using namespace reco;
00058   
00059   edm::LogInfo("TrackValidator") << "\n====================================================" << "\n"
00060                                  << "Analyzing new event" << "\n"
00061                                  << "====================================================\n" << "\n";
00062 
00063   edm::ESHandle<ParametersDefinerForTP> parametersDefinerTP; 
00064   setup.get<TrackAssociatorRecord>().get(parametersDefiner,parametersDefinerTP);    
00065   
00066   edm::Handle<GenParticleCollection>  TPCollectionHeff ;
00067   event.getByLabel(label_tp_effic,TPCollectionHeff);
00068   const GenParticleCollection tPCeff = *(TPCollectionHeff.product());
00069   
00070   edm::Handle<GenParticleCollection>  TPCollectionHfake ;
00071   event.getByLabel(label_tp_fake,TPCollectionHfake);
00072   const GenParticleCollection tPCfake = *(TPCollectionHfake.product());
00073   
00074   //if (tPCeff.size()==0) {edm::LogInfo("TrackValidator") 
00075   //<< "TP Collection for efficiency studies has size = 0! Skipping Event." ; return;}
00076   //if (tPCfake.size()==0) {edm::LogInfo("TrackValidator") 
00077   //<< "TP Collection for fake rate studies has size = 0! Skipping Event." ; return;}
00078   
00079   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00080   event.getByLabel(bsSrc,recoBeamSpotHandle);
00081   reco::BeamSpot bs = *recoBeamSpotHandle;      
00082   
00083   edm::Handle< vector<PileupSummaryInfo> > puinfoH;
00084   event.getByLabel(label_pileupinfo,puinfoH);
00085   PileupSummaryInfo puinfo;      
00086   
00087   for (unsigned int puinfo_ite=0;puinfo_ite<(*puinfoH).size();++puinfo_ite){ 
00088     if ((*puinfoH)[puinfo_ite].getBunchCrossing()==0){
00089       puinfo=(*puinfoH)[puinfo_ite];
00090       break;
00091     }
00092   }
00093   
00094   int w=0; //counter counting the number of sets of histograms
00095   for (unsigned int ww=0;ww<associators.size();ww++){
00096 
00097     associatorByChi2 = dynamic_cast<const TrackAssociatorByChi2*>(associator[ww]);
00098     if (associatorByChi2==0) continue;
00099 
00100     for (unsigned int www=0;www<label.size();www++){
00101       //
00102       //get collections from the event
00103       //
00104       edm::Handle<View<Track> >  trackCollection;
00105       if(!event.getByLabel(label[www], trackCollection)&&ignoremissingtkcollection_)continue;
00106       //if (trackCollection->size()==0) 
00107       //edm::LogInfo("TrackValidator") << "TrackCollection size = 0!" ; 
00108       //continue;
00109       //}
00110       reco::RecoToGenCollection recGenColl;
00111       reco::GenToRecoCollection genRecColl;
00112       //associate tracks
00113       if(UseAssociators){
00114         edm::LogVerbatim("TrackValidator") << "Analyzing " 
00115                                            << label[www].process()<<":"
00116                                            << label[www].label()<<":"
00117                                            << label[www].instance()<<" with "
00118                                            << associators[ww].c_str() <<"\n";
00119         
00120         LogTrace("TrackValidator") << "Calling associateRecoToGen method" << "\n";
00121         recGenColl=associatorByChi2->associateRecoToGen(trackCollection,
00122                                                         TPCollectionHfake,
00123                                                         &event);
00124         LogTrace("TrackValidator") << "Calling associateGenToReco method" << "\n";
00125         genRecColl=associatorByChi2->associateGenToReco(trackCollection,
00126                                                         TPCollectionHeff, 
00127                                                         &event);
00128       }
00129       else{
00130         edm::LogVerbatim("TrackValidator") << "Analyzing " 
00131                                            << label[www].process()<<":"
00132                                            << label[www].label()<<":"
00133                                            << label[www].instance()<<" with "
00134                                            << associatormap.process()<<":"
00135                                            << associatormap.label()<<":"
00136                                            << associatormap.instance()<<"\n";
00137         
00138         Handle<reco::GenToRecoCollection > gentorecoCollectionH;
00139         event.getByLabel(associatormap,gentorecoCollectionH);
00140         genRecColl= *(gentorecoCollectionH.product()); 
00141         
00142         Handle<reco::RecoToGenCollection > recotogenCollectionH;
00143         event.getByLabel(associatormap,recotogenCollectionH);
00144         recGenColl= *(recotogenCollectionH.product()); 
00145       }
00146 
00147 
00148 
00149       // ########################################################
00150       // fill simulation histograms (LOOP OVER TRACKINGPARTICLES)
00151       // ########################################################
00152 
00153       //compute number of tracks per eta interval
00154       //
00155       edm::LogVerbatim("TrackValidator") << "\n# of GenParticles: " << tPCeff.size() << "\n";
00156       int ats(0);         //This counter counts the number of simTracks that are "associated" to recoTracks
00157       int st(0);          //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
00158       unsigned sts(0);   //This counter counts the number of simTracks surviving the bunchcrossing cut 
00159       unsigned asts(0);  //This counter counts the number of simTracks that are "associated" to recoTracks surviving the bunchcrossing cut
00160       for (GenParticleCollection::size_type i=0; i<tPCeff.size(); i++){ //loop over TPs collection for tracking efficiency
00161         GenParticleRef tpr(TPCollectionHeff, i);
00162         GenParticle* tp=const_cast<GenParticle*>(tpr.get());
00163         TrackingParticle::Vector momentumTP; 
00164         TrackingParticle::Point vertexTP;
00165         double dxyGen(0);
00166         double dzGen(0);
00167         
00168 
00169         //---------- THIS PART HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
00170         //If the GenParticle is collison like, get the momentum and vertex at production state
00171         if(parametersDefiner=="LhcParametersDefinerForTP")
00172           {
00173             //fixme this one shold be implemented
00174             if(! gpSelector(*tp)) continue;
00175             momentumTP = tp->momentum();
00176             vertexTP = tp->vertex();
00177             //Calcualte the impact parameters w.r.t. PCA
00178             TrackingParticle::Vector momentum = parametersDefinerTP->momentum(event,setup,*tp);
00179             TrackingParticle::Point vertex = parametersDefinerTP->vertex(event,setup,*tp);
00180             dxyGen = (-vertex.x()*sin(momentum.phi())+vertex.y()*cos(momentum.phi()));
00181             dzGen = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y())/sqrt(momentum.perp2()) 
00182               * momentum.z()/sqrt(momentum.perp2());
00183           }
00184         //If the GenParticle is comics, get the momentum and vertex at PCA
00185         if(parametersDefiner=="CosmicParametersDefinerForTP")
00186           {
00187             //if(! cosmictpSelector(*tp,&bs,event,setup)) continue;     
00188             momentumTP = parametersDefinerTP->momentum(event,setup,*tp);
00189             vertexTP = parametersDefinerTP->vertex(event,setup,*tp);
00190             dxyGen = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
00191             dzGen = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y())/sqrt(momentumTP.perp2()) 
00192               * momentumTP.z()/sqrt(momentumTP.perp2());
00193           }
00194         //---------- THE PART ABOVE HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
00195 
00196         st++;   //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
00197 
00198         // in the coming lines, histos are filled using as input
00199         // - momentumTP 
00200         // - vertexTP 
00201         // - dxyGen
00202         // - dzGen
00203 
00204         histoProducerAlgo_->fill_generic_simTrack_histos(w,momentumTP,vertexTP, tp->collisionId());//fixme: check meaning of collisionId
00205 
00206 
00207         // ##############################################
00208         // fill RecoAssociated GenTracks' histograms
00209         // ##############################################
00210         // bool isRecoMatched(false); // UNUSED
00211         const reco::Track* matchedTrackPointer=0;
00212         std::vector<std::pair<RefToBase<Track>, double> > rt;
00213         if(genRecColl.find(tpr) != genRecColl.end()){
00214           rt = (std::vector<std::pair<RefToBase<Track>, double> >) genRecColl[tpr];
00215           if (rt.size()!=0) {
00216             ats++; //This counter counts the number of simTracks that have a recoTrack associated
00217             // isRecoMatched = true; // UNUSED
00218             matchedTrackPointer = rt.begin()->first.get();
00219             edm::LogVerbatim("TrackValidator") << "GenParticle #" << st 
00220                                                << " with pt=" << sqrt(momentumTP.perp2()) 
00221                                                << " associated with quality:" << rt.begin()->second <<"\n";
00222           }
00223         }else{
00224           edm::LogVerbatim("TrackValidator") 
00225             << "GenParticle #" << st
00226             << " with pt,eta,phi: " 
00227             << sqrt(momentumTP.perp2()) << " , "
00228             << momentumTP.eta() << " , "
00229             << momentumTP.phi() << " , "
00230             << " NOT associated to any reco::Track" << "\n";
00231         }
00232         
00233 
00234         
00235         /*
00236         std::vector<PSimHit> simhits=tp->trackPSimHit(DetId::Tracker);
00237         int nSimHits = simhits.end()-simhits.begin();
00238         */
00239         int nSimHits = 0;
00240 
00241         double vtx_z_PU = vertexTP.z();
00242         /*
00243         for (size_t j = 0; j < tv.size(); j++) {
00244             if (tp->eventId().event() == tv[j].eventId().event()) {
00245                 vtx_z_PU = tv[j].position().z();
00246                 break;
00247             }
00248         }
00249         */
00250 
00251         histoProducerAlgo_->fill_recoAssociated_simTrack_histos(w,*tp,momentumTP,vertexTP,dxyGen,dzGen,nSimHits,matchedTrackPointer,puinfo.getPU_NumInteractions(), vtx_z_PU);
00252 
00253         sts++;
00254         if (matchedTrackPointer) asts++;
00255 
00256 
00257 
00258 
00259       } // End  for (GenParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00260 
00261       //if (st!=0) h_tracksSIM[w]->Fill(st);  // TO BE FIXED
00262       
00263       
00264       // ##############################################
00265       // fill recoTracks histograms (LOOP OVER TRACKS)
00266       // ##############################################
00267       edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with "
00268                                          << label[www].process()<<":"
00269                                          << label[www].label()<<":"
00270                                          << label[www].instance()
00271                                          << ": " << trackCollection->size() << "\n";
00272 
00273       //int sat(0); //This counter counts the number of recoTracks that are associated to GenTracks from Signal only
00274       int at(0); //This counter counts the number of recoTracks that are associated to GenTracks
00275       int rT(0); //This counter counts the number of recoTracks in general
00276 
00277 
00278       // dE/dx
00279       // at some point this could be generalized, with a vector of tags and a corresponding vector of Handles
00280       // I'm writing the interface such to take vectors of ValueMaps
00281       edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx1Handle;
00282       edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx2Handle;
00283       std::vector<edm::ValueMap<reco::DeDxData> > v_dEdx;
00284       v_dEdx.clear();
00285       //std::cout << "PIPPO: label is " << label[www] << std::endl;
00286       if (label[www].label()=="generalTracks") {
00287         try {
00288           event.getByLabel(m_dEdx1Tag, dEdx1Handle);
00289           const edm::ValueMap<reco::DeDxData> dEdx1 = *dEdx1Handle.product();
00290           event.getByLabel(m_dEdx2Tag, dEdx2Handle);
00291           const edm::ValueMap<reco::DeDxData> dEdx2 = *dEdx2Handle.product();
00292           v_dEdx.push_back(dEdx1);
00293           v_dEdx.push_back(dEdx2);
00294         } catch (cms::Exception e){
00295           LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
00296         }
00297       }
00298       //end dE/dx
00299 
00300       for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
00301 
00302         RefToBase<Track> track(trackCollection, i);
00303         rT++;
00304         
00305         bool isSigGenMatched(false);
00306         bool isGenMatched(false);
00307         bool isChargeMatched(true);
00308         int numAssocRecoTracks = 0;
00309         int tpbx = 0;
00310         int nSimHits = 0;
00311         double sharedFraction = 0.;
00312         std::vector<std::pair<GenParticleRef, double> > tp;
00313         if(recGenColl.find(track) != recGenColl.end()){
00314           tp = recGenColl[track];
00315           if (tp.size()!=0) {
00316             /*
00317             std::vector<PSimHit> simhits=tp[0].first->trackPSimHit(DetId::Tracker);
00318             nSimHits = simhits.end()-simhits.begin();
00319             */
00320             sharedFraction = tp[0].second;
00321             isGenMatched = true;
00322             if (tp[0].first->charge() != track->charge()) isChargeMatched = false;
00323             if(genRecColl.find(tp[0].first) != genRecColl.end()) numAssocRecoTracks = genRecColl[tp[0].first].size();
00324             //std::cout << numAssocRecoTracks << std::endl;
00325             /*
00326             tpbx = tp[0].first->eventId().bunchCrossing();
00327             */
00328             at++;
00329             for (unsigned int tp_ite=0;tp_ite<tp.size();++tp_ite){ 
00330               GenParticle trackpart = *(tp[tp_ite].first);
00331               /*
00332               if ((trackpart.eventId().event() == 0) && (trackpart.eventId().bunchCrossing() == 0)){
00333                 isSigGenMatched = true;
00334                 sat++;
00335                 break;
00336               }
00337               */
00338             }
00339             edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt() 
00340                                                << " associated with quality:" << tp.begin()->second <<"\n";
00341           }
00342         } else {
00343           edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00344                                              << " NOT associated to any GenParticle" << "\n";             
00345         }
00346         
00347 
00348         histoProducerAlgo_->fill_generic_recoTrack_histos(w,*track,bs.position(),isGenMatched,isSigGenMatched, isChargeMatched, numAssocRecoTracks, puinfo.getPU_NumInteractions(), tpbx, nSimHits, sharedFraction);
00349 
00350         // dE/dx
00351         //      reco::TrackRef track2  = reco::TrackRef( trackCollection, i );
00352         if (v_dEdx.size() > 0) histoProducerAlgo_->fill_dedx_recoTrack_histos(w,track, v_dEdx);
00353         //if (v_dEdx.size() > 0) histoProducerAlgo_->fill_dedx_recoTrack_histos(track2, v_dEdx);
00354 
00355 
00356         //Fill other histos
00357         //try{ //Is this really necessary ????
00358 
00359         if (tp.size()==0) continue;     
00360 
00361         histoProducerAlgo_->fill_simAssociated_recoTrack_histos(w,*track);
00362 
00363         GenParticleRef tpr = tp.begin()->first;
00364         
00365         /* TO BE FIXED LATER
00366         if (associators[ww]=="TrackAssociatorByChi2"){
00367           //association chi2
00368           double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
00369           h_assochi2[www]->Fill(assocChi2);
00370           h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
00371         }
00372         else if (associators[ww]=="quickTrackAssociatorByHits"){
00373           double fraction = tp.begin()->second;
00374           h_assocFraction[www]->Fill(fraction);
00375           h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
00376         }
00377         */
00378 
00379           
00380         //Get tracking particle parameters at point of closest approach to the beamline
00381         TrackingParticle::Vector momentumTP = parametersDefinerTP->momentum(event,setup,*(tpr.get()));
00382         TrackingParticle::Point vertexTP = parametersDefinerTP->vertex(event,setup,*(tpr.get()));                        
00383         int chargeTP = tpr->charge();
00384 
00385         histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(w,momentumTP,vertexTP,chargeTP,
00386                                                              *track,bs.position());
00387         
00388         
00389         //TO BE FIXED
00390         //std::vector<PSimHit> simhits=tpr.get()->trackPSimHit(DetId::Tracker);
00391         //nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
00392         
00393         /*
00394           } // End of try{
00395           catch (cms::Exception e){
00396           LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
00397           }
00398         */
00399         
00400       } // End of for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
00401 
00402       histoProducerAlgo_->fill_trackBased_histos(w,at,rT,st);
00403 
00404       edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
00405                                          << "Total Associated (genToReco): " << ats << "\n"
00406                                          << "Total Reconstructed: " << rT << "\n"
00407                                          << "Total Associated (recoToGen): " << at << "\n"
00408                                          << "Total Fakes: " << rT-at << "\n";
00409 
00410       w++;
00411     } // End of  for (unsigned int www=0;www<label.size();www++){
00412   } //END of for (unsigned int ww=0;ww<associators.size();ww++){
00413 
00414 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Validation/RecoTrack/plugins/MultiTrackValidatorGenPs.cc
