RecoTrack/plugins/MultiTrackValidatorGenPs.cc

0001 #include "Validation/RecoTrack/interface/MultiTrackValidatorGenPs.h"
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "DataFormats/TrackReco/interface/Track.h"
0007 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0008 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0009 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0010 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0011 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0012 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0013 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0015 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0016 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0017 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0018 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0019 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0020
0021 #include "DataFormats/TrackReco/interface/DeDxData.h"
0022 #include "DataFormats/Common/interface/ValueMap.h"
0023 #include "DataFormats/Common/interface/Ref.h"
0024
0025 #include "TMath.h"
0026 #include <TF1.h>
0027
0028 //#include <iostream>
0029
0030 using namespace std;
0031 using namespace edm;
0032
0033 static const std::string kTrackAssociatorByChi2("trackAssociatorByChi2");
0034
0035 MultiTrackValidatorGenPs::MultiTrackValidatorGenPs(const edm::ParameterSet& pset) : MultiTrackValidator(pset) {
0036   gpSelector = GenParticleCustomSelector(pset.getParameter<double>("ptMinGP"),
0037                                          pset.getParameter<double>("minRapidityGP"),
0038                                          pset.getParameter<double>("maxRapidityGP"),
0039                                          pset.getParameter<double>("tipGP"),
0040                                          pset.getParameter<double>("lipGP"),
0041                                          pset.getParameter<bool>("chargedOnlyGP"),
0042                                          pset.getParameter<int>("statusGP"),
0043                                          pset.getParameter<std::vector<int> >("pdgIdGP"));
0044
0045   if (useAssociators_) {
0046     for (auto const& src : associators) {
0047       if (src.label() == kTrackAssociatorByChi2) {
0048         label_gen_associator = consumes<reco::TrackToGenParticleAssociator>(src);
0049         break;
0050       }
0051     }
0052   } else {
0053     for (auto const& src : associators) {
0054       associatormapGtR = consumes<reco::GenToRecoCollection>(src);
0055       associatormapRtG = consumes<reco::RecoToGenCollection>(src);
0056       break;
0057     }
0058   }
0059 }
0060
0061 MultiTrackValidatorGenPs::~MultiTrackValidatorGenPs() {}
0062
0063 void MultiTrackValidatorGenPs::dqmAnalyze(const edm::Event& event,
0064                                           const edm::EventSetup& setup,
0065                                           const Histograms& histograms) const {
0066   using namespace reco;
0067
0068   edm::LogInfo("TrackValidator") << "\n===================================================="
0069                                  << "\n"
0070                                  << "Analyzing new event"
0071                                  << "\n"
0072                                  << "====================================================\n"
0073                                  << "\n";
0074
0075   const TrackerTopology& ttopo = setup.getData(tTopoEsToken);
0076
0077   edm::Handle<GenParticleCollection> TPCollectionHeff;
0078   event.getByToken(label_tp_effic, TPCollectionHeff);
0079   const GenParticleCollection tPCeff = *(TPCollectionHeff.product());
0080
0081   edm::Handle<GenParticleCollection> TPCollectionHfake;
0082   event.getByToken(label_tp_fake, TPCollectionHfake);
0083   const GenParticleCollection tPCfake = *(TPCollectionHfake.product());
0084
0085   //if (tPCeff.size()==0) {edm::LogInfo("TrackValidator")
0086   //<< "TP Collection for efficiency studies has size = 0! Skipping Event." ; return;}
0087   //if (tPCfake.size()==0) {edm::LogInfo("TrackValidator")
0088   //<< "TP Collection for fake rate studies has size = 0! Skipping Event." ; return;}
0089
0090   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0091   event.getByToken(bsSrc, recoBeamSpotHandle);
0092   reco::BeamSpot bs = *recoBeamSpotHandle;
0093
0094   edm::Handle<vector<PileupSummaryInfo> > puinfoH;
0095   event.getByToken(label_pileupinfo, puinfoH);
0096   PileupSummaryInfo puinfo;
0097
0098   for (unsigned int puinfo_ite = 0; puinfo_ite < (*puinfoH).size(); ++puinfo_ite) {
0099     if ((*puinfoH)[puinfo_ite].getBunchCrossing() == 0) {
0100       puinfo = (*puinfoH)[puinfo_ite];
0101       break;
0102     }
0103   }
0104
0105   const reco::TrackToGenParticleAssociator* trackGenAssociator = nullptr;
0106   if (useAssociators_) {
0107     if (label_gen_associator.isUninitialized()) {
0108       return;
0109     } else {
0110       edm::Handle<reco::TrackToGenParticleAssociator> trackGenAssociatorH;
0111       event.getByToken(label_gen_associator, trackGenAssociatorH);
0112       trackGenAssociator = trackGenAssociatorH.product();
0113     }
0114   } else if (associatormapGtR.isUninitialized()) {
0115     return;
0116   }
0117
0118   // dE/dx
0119   // at some point this could be generalized, with a vector of tags and a corresponding vector of Handles
0120   // I'm writing the interface such to take vectors of ValueMaps
0121   std::vector<const edm::ValueMap<reco::DeDxData>*> v_dEdx;
0122   if (dodEdxPlots_) {
0123     edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx1Handle;
0124     edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx2Handle;
0125     event.getByToken(m_dEdx1Tag, dEdx1Handle);
0126     event.getByToken(m_dEdx2Tag, dEdx2Handle);
0127     v_dEdx.push_back(dEdx1Handle.product());
0128     v_dEdx.push_back(dEdx2Handle.product());
0129   }
0130
0131   std::vector<float> mvaDummy;
0132
0133   int w = 0;  //counter counting the number of sets of histograms
0134   for (unsigned int www = 0; www < label.size(); www++) {
0135     //
0136     //get collections from the event
0137     //
0138     edm::Handle<View<Track> > trackCollection;
0139     if (!event.getByToken(labelToken[www], trackCollection) && ignoremissingtkcollection_)
0140       continue;
0141     //if (trackCollection->size()==0)
0142     //edm::LogInfo("TrackValidator") << "TrackCollection size = 0!" ;
0143     //continue;
0144     //}
0145     reco::RecoToGenCollection recGenColl;
0146     reco::GenToRecoCollection genRecColl;
0147     //associate tracks
0148     if (useAssociators_) {
0149       edm::LogVerbatim("TrackValidator") << "Analyzing " << label[www].process() << ":" << label[www].label() << ":"
0150                                          << label[www].instance() << " with " << kTrackAssociatorByChi2 << "\n";
0151
0152       LogTrace("TrackValidator") << "Calling associateRecoToGen method"
0153                                  << "\n";
0154       recGenColl = trackGenAssociator->associateRecoToGen(trackCollection, TPCollectionHfake);
0155       LogTrace("TrackValidator") << "Calling associateGenToReco method"
0156                                  << "\n";
0157       genRecColl = trackGenAssociator->associateGenToReco(trackCollection, TPCollectionHeff);
0158     } else {
0159       edm::LogVerbatim("TrackValidator") << "Analyzing " << label[www].process() << ":" << label[www].label() << ":"
0160                                          << label[www].instance() << " with " << associators[0] << "\n";
0161
0162       Handle<reco::GenToRecoCollection> gentorecoCollectionH;
0163       event.getByToken(associatormapGtR, gentorecoCollectionH);
0164       genRecColl = *(gentorecoCollectionH.product());
0165
0166       Handle<reco::RecoToGenCollection> recotogenCollectionH;
0167       event.getByToken(associatormapRtG, recotogenCollectionH);
0168       recGenColl = *(recotogenCollectionH.product());
0169     }
0170
0171     // ########################################################
0172     // fill simulation histograms (LOOP OVER TRACKINGPARTICLES)
0173     // ########################################################
0174
0175     //compute number of tracks per eta interval
0176     //
0177     edm::LogVerbatim("TrackValidator") << "\n# of GenParticles: " << tPCeff.size() << "\n";
0178     int ats(0);  //This counter counts the number of simTracks that are "associated" to recoTracks
0179     int st(0);   //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
0180     for (GenParticleCollection::size_type i = 0; i < tPCeff.size();
0181          i++) {  //loop over TPs collection for tracking efficiency
0182       GenParticleRef tpr(TPCollectionHeff, i);
0183       GenParticle* tp = const_cast<GenParticle*>(tpr.get());
0184       TrackingParticle::Vector momentumTP;
0185       TrackingParticle::Point vertexTP;
0186       double dxyGen(0);
0187       double dzGen(0);
0188
0189       //---------- THIS PART HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
0190       //If the GenParticle is collison like, get the momentum and vertex at production state
0191       if (!parametersDefinerIsCosmic_) {
0192         //fixme this one shold be implemented
0193         if (!gpSelector(*tp))
0194           continue;
0195         momentumTP = tp->momentum();
0196         vertexTP = tp->vertex();
0197         //Calcualte the impact parameters w.r.t. PCA
0198         TrackingParticle::Vector momentum = parametersDefinerTP_->momentum(event, setup, *tp);
0199         TrackingParticle::Point vertex = parametersDefinerTP_->vertex(event, setup, *tp);
0200         dxyGen = (-vertex.x() * sin(momentum.phi()) + vertex.y() * cos(momentum.phi()));
0201         dzGen = vertex.z() - (vertex.x() * momentum.x() + vertex.y() * momentum.y()) / sqrt(momentum.perp2()) *
0202                                  momentum.z() / sqrt(momentum.perp2());
0203       }
0204       //If the GenParticle is comics, get the momentum and vertex at PCA
0205       else {
0206         //if(! cosmictpSelector(*tp,&bs,event,setup)) continue;
0207         momentumTP = parametersDefinerTP_->momentum(event, setup, *tp);
0208         vertexTP = parametersDefinerTP_->vertex(event, setup, *tp);
0209         dxyGen = (-vertexTP.x() * sin(momentumTP.phi()) + vertexTP.y() * cos(momentumTP.phi()));
0210         dzGen = vertexTP.z() - (vertexTP.x() * momentumTP.x() + vertexTP.y() * momentumTP.y()) /
0211                                    sqrt(momentumTP.perp2()) * momentumTP.z() / sqrt(momentumTP.perp2());
0212       }
0213       //---------- THE PART ABOVE HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
0214
0215       st++;  //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
0216
0217       // in the coming lines, histos are filled using as input
0218       // - momentumTP
0219       // - vertexTP
0220       // - dxyGen
0221       // - dzGen
0222
0223       if (doSimPlots_ && w == 0) {
0224         histoProducerAlgo_->fill_generic_simTrack_histos(histograms.histoProducerAlgo,
0225                                                          momentumTP,
0226                                                          vertexTP,
0227                                                          tp->collisionId());  //fixme: check meaning of collisionId
0228       }
0229       if (!doSimTrackPlots_)
0230         continue;
0231
0232       // ##############################################
0233       // fill RecoAssociated GenTracks' histograms
0234       // ##############################################
0235       // bool isRecoMatched(false); // UNUSED
0236       const reco::Track* matchedTrackPointer = nullptr;
0237       std::vector<std::pair<RefToBase<Track>, double> > rt;
0238       if (genRecColl.find(tpr) != genRecColl.end()) {
0239         rt = (std::vector<std::pair<RefToBase<Track>, double> >)genRecColl[tpr];
0240         if (!rt.empty()) {
0241           ats++;  //This counter counts the number of simTracks that have a recoTrack associated
0242           // isRecoMatched = true; // UNUSED
0243           matchedTrackPointer = rt.begin()->first.get();
0244           edm::LogVerbatim("TrackValidator") << "GenParticle #" << st << " with pt=" << sqrt(momentumTP.perp2())
0245                                              << " associated with quality:" << rt.begin()->second << "\n";
0246         }
0247       } else {
0248         edm::LogVerbatim("TrackValidator") << "GenParticle #" << st << " with pt,eta,phi: " << sqrt(momentumTP.perp2())
0249                                            << " , " << momentumTP.eta() << " , " << momentumTP.phi() << " , "
0250                                            << " NOT associated to any reco::Track"
0251                                            << "\n";
0252       }
0253
0254       int nSimHits = 0;
0255       histoProducerAlgo_->fill_recoAssociated_simTrack_histos(histograms.histoProducerAlgo,
0256                                                               w,
0257                                                               *tp,
0258                                                               momentumTP,
0259                                                               vertexTP,
0260                                                               dxyGen,
0261                                                               dzGen,
0262                                                               nSimHits,
0263                                                               matchedTrackPointer,
0264                                                               puinfo.getPU_NumInteractions());
0265
0266     }  // End  for (GenParticleCollection::size_type i=0; i<tPCeff.size(); i++){
0267
0268     if (doSimPlots_ && w == 0) {
0269       histoProducerAlgo_->fill_simTrackBased_histos(histograms.histoProducerAlgo, st);
0270     }
0271
0272     // ##############################################
0273     // fill recoTracks histograms (LOOP OVER TRACKS)
0274     // ##############################################
0275     if (!doRecoTrackPlots_)
0276       continue;
0277     edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with " << label[www].process() << ":"
0278                                        << label[www].label() << ":" << label[www].instance() << ": "
0279                                        << trackCollection->size() << "\n";
0280
0281     //int sat(0); //This counter counts the number of recoTracks that are associated to GenTracks from Signal only
0282     int at(0);  //This counter counts the number of recoTracks that are associated to GenTracks
0283     int rT(0);  //This counter counts the number of recoTracks in general
0284
0285     for (View<Track>::size_type i = 0; i < trackCollection->size(); ++i) {
0286       RefToBase<Track> track(trackCollection, i);
0287       rT++;
0288
0289       bool isSigGenMatched(false);
0290       bool isGenMatched(false);
0291       bool isChargeMatched(true);
0292       int numAssocRecoTracks = 0;
0293       int nSimHits = 0;
0294       double sharedFraction = 0.;
0295       std::vector<std::pair<GenParticleRef, double> > tp;
0296       if (recGenColl.find(track) != recGenColl.end()) {
0297         tp = recGenColl[track];
0298         if (!tp.empty()) {
0299           /*
0300         std::vector<PSimHit> simhits=tp[0].first->trackPSimHit(DetId::Tracker);
0301             nSimHits = simhits.end()-simhits.begin();
0302           */
0303           sharedFraction = tp[0].second;
0304           isGenMatched = true;
0305           if (tp[0].first->charge() != track->charge())
0306             isChargeMatched = false;
0307           if (genRecColl.find(tp[0].first) != genRecColl.end())
0308             numAssocRecoTracks = genRecColl[tp[0].first].size();
0309           //std::cout << numAssocRecoTracks << std::endl;
0310           at++;
0311           for (unsigned int tp_ite = 0; tp_ite < tp.size(); ++tp_ite) {
0312             GenParticle trackpart = *(tp[tp_ite].first);
0313             /*
0314           if ((trackpart.eventId().event() == 0) && (trackpart.eventId().bunchCrossing() == 0)){
0315               isSigGenMatched = true;
0316               sat++;
0317               break;
0318           }
0319             */
0320           }
0321           edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
0322                                              << " associated with quality:" << tp.begin()->second << "\n";
0323         }
0324       } else {
0325         edm::LogVerbatim("TrackValidator")
0326             << "reco::Track #" << rT << " with pt=" << track->pt() << " NOT associated to any GenParticle"
0327             << "\n";
0328       }
0329
0330       double dR = 0;  //fixme: plots vs dR and vs dRjet not implemented for now
0331       histoProducerAlgo_->fill_generic_recoTrack_histos(histograms.histoProducerAlgo,
0332                                                         w,
0333                                                         *track,
0334                                                         ttopo,
0335                                                         bs.position(),
0336                                                         nullptr,
0337                                                         nullptr,
0338                                                         isGenMatched,
0339                                                         isSigGenMatched,
0340                                                         isChargeMatched,
0341                                                         numAssocRecoTracks,
0342                                                         puinfo.getPU_NumInteractions(),
0343                                                         nSimHits,
0344                                                         sharedFraction,
0345                                                         dR,
0346                                                         dR,
0347                                                         mvaDummy,
0348                                                         0,
0349                                                         0);
0350
0351       // dE/dx
0352       if (dodEdxPlots_)
0353         histoProducerAlgo_->fill_dedx_recoTrack_histos(histograms.histoProducerAlgo, w, track, v_dEdx);
0354
0355       //Fill other histos
0356       //try{ //Is this really necessary ????
0357
0358       if (tp.empty())
0359         continue;
0360
0361       histoProducerAlgo_->fill_simAssociated_recoTrack_histos(histograms.histoProducerAlgo, w, *track);
0362
0363       GenParticleRef tpr = tp.begin()->first;
0364
0365       /* TO BE FIXED LATER
0366          if (associators[ww]=="TrackAssociatorByChi2"){
0367          //association chi2
0368          double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
0369          h_assochi2[www]->Fill(assocChi2);
0370          h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
0371          }
0372          else if (associators[ww]=="quickTrackAssociatorByHits"){
0373          double fraction = tp.begin()->second;
0374          h_assocFraction[www]->Fill(fraction);
0375          h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
0376     }
0377       */
0378
0379       //Get tracking particle parameters at point of closest approach to the beamline
0380       TrackingParticle::Vector momentumTP = parametersDefinerTP_->momentum(event, setup, *(tpr.get()));
0381       TrackingParticle::Point vertexTP = parametersDefinerTP_->vertex(event, setup, *(tpr.get()));
0382       int chargeTP = tpr->charge();
0383
0384       histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(
0385           histograms.histoProducerAlgo, w, momentumTP, vertexTP, chargeTP, *track, bs.position());
0386
0387       //TO BE FIXED
0388       //std::vector<PSimHit> simhits=tpr.get()->trackPSimHit(DetId::Tracker);
0389       //nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
0390
0391       /*
0392         } // End of try{
0393         catch (cms::Exception e){
0394         LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
0395         }
0396       */
0397
0398     }  // End of for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
0399
0400     histoProducerAlgo_->fill_trackBased_histos(histograms.histoProducerAlgo, w, at, rT, rT, st);
0401
0402     edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
0403                                        << "Total Associated (genToReco): " << ats << "\n"
0404                                        << "Total Reconstructed: " << rT << "\n"
0405                                        << "Total Associated (recoToGen): " << at << "\n"
0406                                        << "Total Fakes: " << rT - at << "\n";
0407
0408     w++;
0409   }  // End of  for (unsigned int www=0;www<label.size();www++){
0410 }