EventGenerator/plugins/BasicGenParticleValidation.cc

0001 /*class BasicGenParticleValidation
0002  *
0003  *  Class to fill dqm monitor elements from existing EDM file
0004  *
0005  */
0006
0007 #include "Validation/EventGenerator/interface/BasicGenParticleValidation.h"
0008
0009 #include "CLHEP/Units/defs.h"
0010 #include "CLHEP/Units/PhysicalConstants.h"
0011 #include "Validation/EventGenerator/interface/DQMHelper.h"
0012
0013 using namespace edm;
0014
0015 BasicGenParticleValidation::BasicGenParticleValidation(const edm::ParameterSet& iPSet)
0016     : wmanager_(iPSet, consumesCollector()),
0017       hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
0018       genparticleCollection_(iPSet.getParameter<edm::InputTag>("genparticleCollection")),
0019       genjetCollection_(iPSet.getParameter<edm::InputTag>("genjetsCollection")),
0020       matchPr_(iPSet.getParameter<double>("matchingPrecision")),
0021       verbosity_(iPSet.getUntrackedParameter<unsigned int>("verbosity", 0)),
0022       signalParticlesOnly_(iPSet.getParameter<bool>("signalParticlesOnly")) {
0023   hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
0024   genparticleCollectionToken_ = consumes<reco::GenParticleCollection>(genparticleCollection_);
0025   genjetCollectionToken_ = consumes<reco::GenJetCollection>(genjetCollection_);
0026 }
0027
0028 BasicGenParticleValidation::~BasicGenParticleValidation() {}
0029
0030 void BasicGenParticleValidation::bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) {
0031   ///Setting the DQM top directories
0032   DQMHelper dqm(&i);
0033   i.setCurrentFolder("Generator/GenParticles");
0034   ///Booking the ME's
0035
0036   // Number of analyzed events
0037   nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., "", "Number of Events");
0038
0039   ///multiplicity
0040   genPMultiplicity = dqm.book1dHisto("genPMultiplicty",
0041                                      "Log(No. all GenParticles)",
0042                                      50,
0043                                      -1,
0044                                      5,
0045                                      "log_{10}(N_{All GenParticles})",
0046                                      "Number of Events");  //Log
0047   //difference in HepMC and reco multiplicity
0048   genMatched = dqm.book1dHisto(
0049       "genMatched", "Difference reco - matched", 50, -25, 25, "N_{All GenParticles}-N_{Matched}", "Number of Events");
0050   //multiple matching
0051   multipleMatching = dqm.book1dHisto("multipleMatching",
0052                                      "multiple reco HepMC matching",
0053                                      50,
0054                                      0,
0055                                      50,
0056                                      "N_{multiple reco HepMC matching}",
0057                                      "Number of Events");
0058   //momentum difference of matched particles
0059   matchedResolution = dqm.book1dHisto("matchedResolution",
0060                                       "log10(momentum difference of matched particles)",
0061                                       70,
0062                                       -10.,
0063                                       -3.,
0064                                       "log_{10}(#DeltaP_{matched Particles})",
0065                                       "Number of Events");
0066
0067   // GenJet general distributions
0068   genJetMult = dqm.book1dHisto("genJetMult", "GenJet multiplicity", 50, 0, 50, "N_{gen-jets}", "Number of Events");
0069   genJetEnergy = dqm.book1dHisto(
0070       "genJetEnergy", "Log10(GenJet energy)", 60, -1, 5, "log_{10}(E^{gen-jets}) (log_{10}(GeV))", "Number of Events");
0071   genJetPt = dqm.book1dHisto(
0072       "genJetPt", "Log10(GenJet pt)", 60, -1, 5, "log_{10}(P_{t}^{gen-jets}) (log_{10}(GeV))", "Number of Events");
0073   genJetEta = dqm.book1dHisto("genJetEta", "GenJet eta", 220, -11, 11, "#eta^{gen-jets}", "Number of Events");
0074   genJetPhi = dqm.book1dHisto("genJetPhi", "GenJet phi", 360, -180, 180, "#phi^{gen-jets} (rad)", "Number of Events");
0075   genJetDeltaEtaMin = dqm.book1dHisto(
0076       "genJetDeltaEtaMin", "GenJet minimum rapidity gap", 30, 0, 30, "#delta#eta_{min}^{gen-jets}", "Number of Events");
0077
0078   genJetPto1 = dqm.book1dHisto(
0079       "genJetPto1", "GenJet multiplicity above 1 GeV", 50, 0, 50, "N_{gen-jets P_{t}>1GeV}", "Number of Events");
0080   genJetPto10 = dqm.book1dHisto(
0081       "genJetPto10", "GenJet multiplicity above 10 GeV", 50, 0, 50, "N_{gen-jets P_{t}>10GeV}", "Number of Events");
0082   genJetPto100 = dqm.book1dHisto(
0083       "genJetPto100", "GenJet multiplicity above 100 GeV", 50, 0, 50, "N_{gen-jets P_{t}>100GeV}", "Number of Events");
0084   genJetCentral = dqm.book1dHisto(
0085       "genJetCentral", "GenJet multiplicity |eta|.lt.2.5", 50, 0, 50, "N_{gen-jets |#eta|#leq2.5}", "Number of Events");
0086
0087   genJetTotPt = dqm.book1dHisto("genJetTotPt",
0088                                 "Log10(GenJet total pt)",
0089                                 100,
0090                                 -5,
0091                                 5,
0092                                 "log_{10}(#SigmaP_{t}^{gen-jets}) (log_{10}(GeV))",
0093                                 "Number of Events");
0094
0095   return;
0096 }
0097
0098 void BasicGenParticleValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0099   unsigned int initSize = 1000;
0100
0101   ///Gathering the HepMCProduct information
0102   edm::Handle<HepMCProduct> evt;
0103   iEvent.getByToken(hepmcCollectionToken_, evt);
0104
0105   //Get HepMC EVENT
0106   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0107
0108   double weight = wmanager_.weight(iEvent);
0109
0110   nEvt->Fill(0.5, weight);
0111
0112   std::vector<const HepMC::GenParticle*> hepmcGPCollection;
0113   std::vector<int> barcodeList;
0114   hepmcGPCollection.reserve(initSize);
0115   barcodeList.reserve(initSize);
0116
0117   //Looping through HepMC::GenParticle collection to search for status 1 particles
0118   for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
0119        iter != myGenEvent->particles_end();
0120        ++iter) {
0121     if ((*iter)->status() == 1) {
0122       hepmcGPCollection.push_back(*iter);
0123       barcodeList.push_back((*iter)->barcode());
0124       if (verbosity_ > 0) {
0125         std::cout << "HepMC " << std::setw(14) << std::fixed << (*iter)->pdg_id() << std::setw(14) << std::fixed
0126                   << (*iter)->momentum().px() << std::setw(14) << std::fixed << (*iter)->momentum().py()
0127                   << std::setw(14) << std::fixed << (*iter)->momentum().pz() << std::endl;
0128       }
0129     }
0130   }
0131
0132   // Gather information on the reco::GenParticle collection
0133   edm::Handle<reco::GenParticleCollection> genParticles;
0134   iEvent.getByToken(genparticleCollectionToken_, genParticles);
0135
0136   std::vector<const reco::GenParticle*> particles;
0137   particles.reserve(initSize);
0138   for (reco::GenParticleCollection::const_iterator iter = genParticles->begin(); iter != genParticles->end(); ++iter) {
0139     if ((*iter).status() == 1) {
0140       particles.push_back(&*iter);
0141       if (verbosity_ > 0) {
0142         std::cout << "reco  " << std::setw(14) << std::fixed << (*iter).pdgId() << std::setw(14) << std::fixed
0143                   << (*iter).px() << std::setw(14) << std::fixed << (*iter).py() << std::setw(14) << std::fixed
0144                   << (*iter).pz() << std::endl;
0145       }
0146     }
0147   }
0148
0149   unsigned int nReco = particles.size();
0150   unsigned int nHepMC = hepmcGPCollection.size();
0151
0152   if (signalParticlesOnly_) {
0153     for (unsigned int i = 0; i < particles.size(); ++i) {
0154       if (particles[i]->collisionId() != 0)
0155         nReco -= 1;
0156     }
0157   }
0158
0159   genPMultiplicity->Fill(std::log10(nReco), weight);
0160
0161   // Define vector containing index of hepmc corresponding to the reco::GenParticle
0162   std::vector<int> hepmcMatchIndex;
0163   hepmcMatchIndex.reserve(initSize);
0164
0165   // Matching procedure
0166
0167   // Check array size consistency
0168
0169   if (nReco != nHepMC) {
0170     edm::LogWarning("CollectionSizeInconsistency")
0171         << "Collection size inconsistency: HepMC::GenParticle = " << nHepMC << " reco::GenParticle = " << nReco;
0172   }
0173
0174   // Match each HepMC with a reco
0175
0176   for (unsigned int i = 0; i < nReco; ++i) {
0177     for (unsigned int j = 0; j < nHepMC; ++j) {
0178       if (matchParticles(hepmcGPCollection[j], particles[i])) {
0179         hepmcMatchIndex.push_back((int)j);
0180         if (hepmcGPCollection[j]->momentum().rho() != 0.) {
0181           double reso = 1. - particles[i]->p() / hepmcGPCollection[j]->momentum().rho();
0182           if (verbosity_ > 0) {
0183             std::cout << "Matching momentum: reco = " << particles[i]->p()
0184                       << " HepMC = " << hepmcGPCollection[j]->momentum().rho() << " resoultion = " << reso << std::endl;
0185           }
0186           matchedResolution->Fill(std::log10(std::fabs(reso)), weight);
0187         }
0188         continue;
0189       }
0190     }
0191   }
0192
0193   // Check unicity of matching
0194
0195   unsigned int nMatched = hepmcMatchIndex.size();
0196
0197   if (nMatched != nReco) {
0198     edm::LogWarning("IncorrectMatching") << "Incorrect number of matched indexes: GenParticle = " << nReco
0199                                          << " matched indexes = " << nMatched;
0200   }
0201   genMatched->Fill(int(nReco - nMatched), weight);
0202
0203   unsigned int nWrMatch = 0;
0204
0205   for (unsigned int i = 0; i < nMatched; ++i) {
0206     for (unsigned int j = i + 1; j < nMatched; ++j) {
0207       if (hepmcMatchIndex[i] == hepmcMatchIndex[j]) {
0208         int theIndex = hepmcMatchIndex[i];
0209         edm::LogWarning("DuplicatedMatching")
0210             << "Multiple matching occurencies for GenParticle barcode = " << barcodeList[theIndex];
0211         nWrMatch++;
0212       }
0213     }
0214   }
0215   multipleMatching->Fill(int(nWrMatch), weight);
0216
0217   // Gather information in the GenJet collection
0218   edm::Handle<reco::GenJetCollection> genJets;
0219   iEvent.getByToken(genjetCollectionToken_, genJets);
0220
0221   int nJets = 0;
0222   int nJetso1 = 0;
0223   int nJetso10 = 0;
0224   int nJetso100 = 0;
0225   int nJetsCentral = 0;
0226   double totPt = 0.;
0227
0228   std::vector<double> jetEta;
0229   jetEta.reserve(initSize);
0230
0231   for (reco::GenJetCollection::const_iterator iter = genJets->begin(); iter != genJets->end(); ++iter) {
0232     nJets++;
0233     double pt = (*iter).pt();
0234     totPt += pt;
0235     if (pt > 1.)
0236       nJetso1++;
0237     if (pt > 10.)
0238       nJetso10++;
0239     if (pt > 100.)
0240       nJetso100++;
0241     double eta = (*iter).eta();
0242     if (std::fabs(eta) < 2.5)
0243       nJetsCentral++;
0244     jetEta.push_back(eta);
0245
0246     genJetEnergy->Fill(std::log10((*iter).energy()), weight);
0247     genJetPt->Fill(std::log10(pt), weight);
0248     genJetEta->Fill(eta, weight);
0249     genJetPhi->Fill((*iter).phi() / CLHEP::degree, weight);
0250   }
0251
0252   genJetMult->Fill(nJets, weight);
0253   genJetPto1->Fill(nJetso1, weight);
0254   genJetPto10->Fill(nJetso10, weight);
0255   genJetPto100->Fill(nJetso100, weight);
0256   genJetCentral->Fill(nJetsCentral, weight);
0257
0258   genJetTotPt->Fill(std::log10(totPt), weight);
0259
0260   double deltaEta = 999.;
0261   if (jetEta.size() > 1) {
0262     for (unsigned int i = 0; i < jetEta.size(); i++) {
0263       for (unsigned int j = i + 1; j < jetEta.size(); j++) {
0264         deltaEta = std::min(deltaEta, std::fabs(jetEta[i] - jetEta[j]));
0265       }
0266     }
0267   }
0268
0269   genJetDeltaEtaMin->Fill(deltaEta, weight);
0270
0271   delete myGenEvent;
0272 }  //analyze
0273
0274 bool BasicGenParticleValidation::matchParticles(const HepMC::GenParticle*& hepmcP, const reco::GenParticle*& recoP) {
0275   bool state = false;
0276
0277   if (hepmcP->pdg_id() != recoP->pdgId())
0278     return state;
0279   if (std::fabs(hepmcP->momentum().px() - recoP->px()) <= std::fabs(matchPr_ * hepmcP->momentum().px()) &&
0280       std::fabs(hepmcP->momentum().py() - recoP->py()) <= std::fabs(matchPr_ * hepmcP->momentum().py()) &&
0281       std::fabs(hepmcP->momentum().pz() - recoP->pz()) <= std::fabs(matchPr_ * hepmcP->momentum().pz())) {
0282     state = true;
0283   }
0284
0285   return state;
0286 }