Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-18 22:30:56

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     unsigned sts(0);  //This counter counts the number of simTracks surviving the bunchcrossing cut
0181     unsigned asts(
0182         0);  //This counter counts the number of simTracks that are "associated" to recoTracks surviving the bunchcrossing cut
0183     for (GenParticleCollection::size_type i = 0; i < tPCeff.size();
0184          i++) {  //loop over TPs collection for tracking efficiency
0185       GenParticleRef tpr(TPCollectionHeff, i);
0186       GenParticle* tp = const_cast<GenParticle*>(tpr.get());
0187       TrackingParticle::Vector momentumTP;
0188       TrackingParticle::Point vertexTP;
0189       double dxyGen(0);
0190       double dzGen(0);
0191 
0192       //---------- THIS PART HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
0193       //If the GenParticle is collison like, get the momentum and vertex at production state
0194       if (!parametersDefinerIsCosmic_) {
0195         //fixme this one shold be implemented
0196         if (!gpSelector(*tp))
0197           continue;
0198         momentumTP = tp->momentum();
0199         vertexTP = tp->vertex();
0200         //Calcualte the impact parameters w.r.t. PCA
0201         TrackingParticle::Vector momentum = parametersDefinerTP_->momentum(event, setup, *tp);
0202         TrackingParticle::Point vertex = parametersDefinerTP_->vertex(event, setup, *tp);
0203         dxyGen = (-vertex.x() * sin(momentum.phi()) + vertex.y() * cos(momentum.phi()));
0204         dzGen = vertex.z() - (vertex.x() * momentum.x() + vertex.y() * momentum.y()) / sqrt(momentum.perp2()) *
0205                                  momentum.z() / sqrt(momentum.perp2());
0206       }
0207       //If the GenParticle is comics, get the momentum and vertex at PCA
0208       else {
0209         //if(! cosmictpSelector(*tp,&bs,event,setup)) continue;
0210         momentumTP = parametersDefinerTP_->momentum(event, setup, *tp);
0211         vertexTP = parametersDefinerTP_->vertex(event, setup, *tp);
0212         dxyGen = (-vertexTP.x() * sin(momentumTP.phi()) + vertexTP.y() * cos(momentumTP.phi()));
0213         dzGen = vertexTP.z() - (vertexTP.x() * momentumTP.x() + vertexTP.y() * momentumTP.y()) /
0214                                    sqrt(momentumTP.perp2()) * momentumTP.z() / sqrt(momentumTP.perp2());
0215       }
0216       //---------- THE PART ABOVE HAS TO BE CLEANED UP. THE PARAMETER DEFINER WAS NOT MEANT TO BE USED IN THIS WAY ----------
0217 
0218       st++;  //This counter counts the number of simulated tracks passing the MTV selection (i.e. tpSelector(tp) )
0219 
0220       // in the coming lines, histos are filled using as input
0221       // - momentumTP
0222       // - vertexTP
0223       // - dxyGen
0224       // - dzGen
0225 
0226       if (doSimPlots_ && w == 0) {
0227         histoProducerAlgo_->fill_generic_simTrack_histos(histograms.histoProducerAlgo,
0228                                                          momentumTP,
0229                                                          vertexTP,
0230                                                          tp->collisionId());  //fixme: check meaning of collisionId
0231       }
0232       if (!doSimTrackPlots_)
0233         continue;
0234 
0235       // ##############################################
0236       // fill RecoAssociated GenTracks' histograms
0237       // ##############################################
0238       // bool isRecoMatched(false); // UNUSED
0239       const reco::Track* matchedTrackPointer = nullptr;
0240       std::vector<std::pair<RefToBase<Track>, double> > rt;
0241       if (genRecColl.find(tpr) != genRecColl.end()) {
0242         rt = (std::vector<std::pair<RefToBase<Track>, double> >)genRecColl[tpr];
0243         if (!rt.empty()) {
0244           ats++;  //This counter counts the number of simTracks that have a recoTrack associated
0245           // isRecoMatched = true; // UNUSED
0246           matchedTrackPointer = rt.begin()->first.get();
0247           edm::LogVerbatim("TrackValidator") << "GenParticle #" << st << " with pt=" << sqrt(momentumTP.perp2())
0248                                              << " associated with quality:" << rt.begin()->second << "\n";
0249         }
0250       } else {
0251         edm::LogVerbatim("TrackValidator") << "GenParticle #" << st << " with pt,eta,phi: " << sqrt(momentumTP.perp2())
0252                                            << " , " << momentumTP.eta() << " , " << momentumTP.phi() << " , "
0253                                            << " NOT associated to any reco::Track"
0254                                            << "\n";
0255       }
0256 
0257       int nSimHits = 0;
0258       histoProducerAlgo_->fill_recoAssociated_simTrack_histos(histograms.histoProducerAlgo,
0259                                                               w,
0260                                                               *tp,
0261                                                               momentumTP,
0262                                                               vertexTP,
0263                                                               dxyGen,
0264                                                               dzGen,
0265                                                               nSimHits,
0266                                                               matchedTrackPointer,
0267                                                               puinfo.getPU_NumInteractions());
0268 
0269       sts++;
0270       if (matchedTrackPointer)
0271         asts++;
0272 
0273     }  // End  for (GenParticleCollection::size_type i=0; i<tPCeff.size(); i++){
0274 
0275     if (doSimPlots_ && w == 0) {
0276       histoProducerAlgo_->fill_simTrackBased_histos(histograms.histoProducerAlgo, st);
0277     }
0278 
0279     // ##############################################
0280     // fill recoTracks histograms (LOOP OVER TRACKS)
0281     // ##############################################
0282     if (!doRecoTrackPlots_)
0283       continue;
0284     edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with " << label[www].process() << ":"
0285                                        << label[www].label() << ":" << label[www].instance() << ": "
0286                                        << trackCollection->size() << "\n";
0287 
0288     //int sat(0); //This counter counts the number of recoTracks that are associated to GenTracks from Signal only
0289     int at(0);  //This counter counts the number of recoTracks that are associated to GenTracks
0290     int rT(0);  //This counter counts the number of recoTracks in general
0291 
0292     for (View<Track>::size_type i = 0; i < trackCollection->size(); ++i) {
0293       RefToBase<Track> track(trackCollection, i);
0294       rT++;
0295 
0296       bool isSigGenMatched(false);
0297       bool isGenMatched(false);
0298       bool isChargeMatched(true);
0299       int numAssocRecoTracks = 0;
0300       int nSimHits = 0;
0301       double sharedFraction = 0.;
0302       std::vector<std::pair<GenParticleRef, double> > tp;
0303       if (recGenColl.find(track) != recGenColl.end()) {
0304         tp = recGenColl[track];
0305         if (!tp.empty()) {
0306           /*
0307         std::vector<PSimHit> simhits=tp[0].first->trackPSimHit(DetId::Tracker);
0308             nSimHits = simhits.end()-simhits.begin();
0309           */
0310           sharedFraction = tp[0].second;
0311           isGenMatched = true;
0312           if (tp[0].first->charge() != track->charge())
0313             isChargeMatched = false;
0314           if (genRecColl.find(tp[0].first) != genRecColl.end())
0315             numAssocRecoTracks = genRecColl[tp[0].first].size();
0316           //std::cout << numAssocRecoTracks << std::endl;
0317           at++;
0318           for (unsigned int tp_ite = 0; tp_ite < tp.size(); ++tp_ite) {
0319             GenParticle trackpart = *(tp[tp_ite].first);
0320             /*
0321           if ((trackpart.eventId().event() == 0) && (trackpart.eventId().bunchCrossing() == 0)){
0322               isSigGenMatched = true;
0323               sat++;
0324               break;
0325           }
0326             */
0327           }
0328           edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
0329                                              << " associated with quality:" << tp.begin()->second << "\n";
0330         }
0331       } else {
0332         edm::LogVerbatim("TrackValidator")
0333             << "reco::Track #" << rT << " with pt=" << track->pt() << " NOT associated to any GenParticle"
0334             << "\n";
0335       }
0336 
0337       double dR = 0;  //fixme: plots vs dR and vs dRjet not implemented for now
0338       histoProducerAlgo_->fill_generic_recoTrack_histos(histograms.histoProducerAlgo,
0339                                                         w,
0340                                                         *track,
0341                                                         ttopo,
0342                                                         bs.position(),
0343                                                         nullptr,
0344                                                         nullptr,
0345                                                         isGenMatched,
0346                                                         isSigGenMatched,
0347                                                         isChargeMatched,
0348                                                         numAssocRecoTracks,
0349                                                         puinfo.getPU_NumInteractions(),
0350                                                         nSimHits,
0351                                                         sharedFraction,
0352                                                         dR,
0353                                                         dR,
0354                                                         mvaDummy,
0355                                                         0,
0356                                                         0);
0357 
0358       // dE/dx
0359       if (dodEdxPlots_)
0360         histoProducerAlgo_->fill_dedx_recoTrack_histos(histograms.histoProducerAlgo, w, track, v_dEdx);
0361 
0362       //Fill other histos
0363       //try{ //Is this really necessary ????
0364 
0365       if (tp.empty())
0366         continue;
0367 
0368       histoProducerAlgo_->fill_simAssociated_recoTrack_histos(histograms.histoProducerAlgo, w, *track);
0369 
0370       GenParticleRef tpr = tp.begin()->first;
0371 
0372       /* TO BE FIXED LATER
0373          if (associators[ww]=="TrackAssociatorByChi2"){
0374          //association chi2
0375          double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
0376          h_assochi2[www]->Fill(assocChi2);
0377          h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
0378          }
0379          else if (associators[ww]=="quickTrackAssociatorByHits"){
0380          double fraction = tp.begin()->second;
0381          h_assocFraction[www]->Fill(fraction);
0382          h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
0383     }
0384       */
0385 
0386       //Get tracking particle parameters at point of closest approach to the beamline
0387       TrackingParticle::Vector momentumTP = parametersDefinerTP_->momentum(event, setup, *(tpr.get()));
0388       TrackingParticle::Point vertexTP = parametersDefinerTP_->vertex(event, setup, *(tpr.get()));
0389       int chargeTP = tpr->charge();
0390 
0391       histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(
0392           histograms.histoProducerAlgo, w, momentumTP, vertexTP, chargeTP, *track, bs.position());
0393 
0394       //TO BE FIXED
0395       //std::vector<PSimHit> simhits=tpr.get()->trackPSimHit(DetId::Tracker);
0396       //nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
0397 
0398       /*
0399         } // End of try{
0400         catch (cms::Exception e){
0401         LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
0402         }
0403       */
0404 
0405     }  // End of for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
0406 
0407     histoProducerAlgo_->fill_trackBased_histos(histograms.histoProducerAlgo, w, at, rT, rT, st);
0408 
0409     edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
0410                                        << "Total Associated (genToReco): " << ats << "\n"
0411                                        << "Total Reconstructed: " << rT << "\n"
0412                                        << "Total Associated (recoToGen): " << at << "\n"
0413                                        << "Total Fakes: " << rT - at << "\n";
0414 
0415     w++;
0416   }  // End of  for (unsigned int www=0;www<label.size();www++){
0417 }