/
githubmirror
/
cmssw
Обзор
Документация
Войти
/
githubmirror
/
cmssw
Код
Запросы
0
Пакеты
0
Релизы
0
Аналитика
Безопасность
master
Validation/EventGenerator/plugins/BasicHepMCValidation.cc
499 строк
21 KB
Christopher Jones
Use esConsumes for modules in Validation/EventGenerator
18 окт 2021, 18:28
18 окт 2021, 18:28
e850268
Код
Авторство
О чём код?
/*class BasicHepMCValidation * * Class to fill dqm monitor elements from existing EDM file * */ #include "Validation/EventGenerator/interface/BasicHepMCValidation.h" #include "CLHEP/Units/defs.h" #include "CLHEP/Units/PhysicalConstants.h" #include "Validation/EventGenerator/interface/DQMHelper.h" using namespace edm; BasicHepMCValidation::BasicHepMCValidation(const edm::ParameterSet &iPSet) : wmanager_(iPSet, consumesCollector()), hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")) { hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_); fPDGTableToken = esConsumes<edm::Transition::BeginRun>(); } BasicHepMCValidation::~BasicHepMCValidation() {} void BasicHepMCValidation::dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) { fPDGTable = c.getHandle(fPDGTableToken); } namespace { // Set upper bound & lower bound for PDF & Scale related histograms constexpr double logPdfMax = 3.0; constexpr double logPdfMin = -3.0; constexpr int logPdfNbin = 150; constexpr double logPdfBinsize = (logPdfMax - logPdfMin) / (double)logPdfNbin; constexpr double logQScaleMax = 4.0; constexpr double logQScaleMin = -1.0; constexpr int logQScaleNbin = 500; constexpr double logQScaleBinsize = (logQScaleMax - logQScaleMin) / (double)logQScaleNbin; } // namespace void BasicHepMCValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) { ///Setting the DQM top directories DQMHelper dqm(&i); i.setCurrentFolder("Generator/Particles"); // Number of analyzed events nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1.); ///Booking the ME's ///multiplicity // quarks particles.push_back(ParticleMonitor("u", 1, i)); particles.push_back(ParticleMonitor("ubar", -1, i)); particles.push_back(ParticleMonitor("d", 2, i)); particles.push_back(ParticleMonitor("dbar", -2, i)); particles.push_back(ParticleMonitor("s", 3, i)); particles.push_back(ParticleMonitor("sbar", -3, i)); particles.push_back(ParticleMonitor("c", 4, i)); particles.push_back(ParticleMonitor("cbar", -4, i)); particles.push_back(ParticleMonitor("b", 5, i)); particles.push_back(ParticleMonitor("bbar", -5, i)); particles.push_back(ParticleMonitor("t", 6, i)); particles.push_back(ParticleMonitor("tbar", -6, i)); //leptons particles.push_back(ParticleMonitor("eminus", 11, i)); particles.push_back(ParticleMonitor("eplus", -11, i)); particles.push_back(ParticleMonitor("nue", 12, i)); particles.push_back(ParticleMonitor("nuebar", -12, i)); particles.push_back(ParticleMonitor("muminus", 13, i)); particles.push_back(ParticleMonitor("muplus", -13, i)); particles.push_back(ParticleMonitor("numu", 14, i)); particles.push_back(ParticleMonitor("numubar", -14, i)); particles.push_back(ParticleMonitor("tauminus", 15, i)); particles.push_back(ParticleMonitor("tauplus", -15, i)); particles.push_back(ParticleMonitor("nutau", 16, i)); particles.push_back(ParticleMonitor("nutaubar", -16, i)); //bosons particles.push_back(ParticleMonitor("Wplus", 24, i)); particles.push_back(ParticleMonitor("Wminus", -24, i)); particles.push_back(ParticleMonitor("Z", 23, i)); particles.push_back(ParticleMonitor("gamma", 22, i)); particles.push_back(ParticleMonitor("gluon", 21, i)); //mesons particles.push_back(ParticleMonitor("piplus", 211, i, true)); //log particles.push_back(ParticleMonitor("piminus", -211, i, true)); //log particles.push_back(ParticleMonitor("pizero", 111, i, true)); //log particles.push_back(ParticleMonitor("Kplus", 321, i)); particles.push_back(ParticleMonitor("Kminus", -321, i)); particles.push_back(ParticleMonitor("Klzero", 130, i)); particles.push_back(ParticleMonitor("Kszero", 310, i)); //baryons particles.push_back(ParticleMonitor("p", 2212, i, true)); //log particles.push_back(ParticleMonitor("pbar", -2212, i, true)); //log particles.push_back(ParticleMonitor("n", 2112, i, true)); //log particles.push_back(ParticleMonitor("nbar", -2112, i, true)); //log particles.push_back(ParticleMonitor("lambda0", 3122, i)); particles.push_back(ParticleMonitor("lambda0bar", -3122, i)); //D mesons particles.push_back(ParticleMonitor("Dplus", 411, i)); particles.push_back(ParticleMonitor("Dminus", -411, i)); particles.push_back(ParticleMonitor("Dzero", 421, i)); particles.push_back(ParticleMonitor("Dzerobar", -421, i)); //B mesons particles.push_back(ParticleMonitor("Bplus", 521, i)); particles.push_back(ParticleMonitor("Bminus", -521, i)); particles.push_back(ParticleMonitor("Bzero", 511, i)); particles.push_back(ParticleMonitor("Bzerobar", -511, i)); particles.push_back(ParticleMonitor("Bszero", 531, i)); particles.push_back(ParticleMonitor("Bszerobar", -531, i)); // otherPtclNumber = dqm.book1dHisto( "otherPtclNumber", "Log10(No. other ptcls)", 60, -1, 5, "log_{10}(No. other ptcls)", "Number of Events"); //Log otherPtclMomentum = dqm.book1dHisto("otherPtclMomentum", "Log10(p) other ptcls", 60, -2, 4, "log10(P^{other ptcls}) (log_{10}(GeV))", "Number of Events"); ///other genPtclNumber = dqm.book1dHisto( "genPtclNumber", "Log10(No. all particles)", 60, -1, 5, "log10(No. all particles)", "Number of Events"); //Log genVrtxNumber = dqm.book1dHisto( "genVrtxNumber", "Log10(No. all vertexs)", 60, -1, 5, "log10(No. all vertexs)", "Number of Events"); //Log // stablePtclNumber = dqm.book1dHisto("stablePtclNumber", "Log10(No. stable particles)", 50, 0, 5, "log10(No. stable particles)", "Number of Events"); //Log stablePtclPhi = dqm.book1dHisto( "stablePtclPhi", "stable Ptcl Phi", 360, -180, 180, "#phi^{stable Ptcl} (rad)", "Number of Events"); stablePtclEta = dqm.book1dHisto( "stablePtclEta", "stable Ptcl Eta (pseudo rapidity)", 220, -11, 11, "#eta^{stable Ptcl}", "Number of Events"); stablePtclCharge = dqm.book1dHisto("stablePtclCharge", "stablePtclCharge", 5, -2, 2, "Charge^{stable ptcls}", "Number of Events"); stableChaNumber = dqm.book1dHisto("stableChaNumber", "Log10(No. stable charged particles)", 50, 0, 5, "log_{10}(No. stable charged particles)", "Number of Events"); //Log stablePtclp = dqm.book1dHisto("stablePtclp", "Log10(p) stable ptcl p", 80, -4, 4, "log_{10}(P^{stable ptcl}) (log_{10}(GeV))", "Number of Events"); //Log stablePtclpT = dqm.book1dHisto("stablePtclpT", "Log10(pT) stable ptcl pT", 80, -4, 4, "log_{10}(P_{t}^{stable ptcl}) (log_{10}(GeV))", "Number of Events"); //Log partonNumber = dqm.book1dHisto("partonNumber", "number of partons", 100, 0, 100, "number of partons", "Number of Events"); partonpT = dqm.book1dHisto("partonpT", "Log10(pT) parton pT", 80, -4, 4, "Log10(P_{t}^{parton})", "Number of Events"); //Log outVrtxStablePtclNumber = dqm.book1dHisto("outVrtxStablePtclNumber", "No. outgoing stable ptcls from vrtx", 10, 0, 10, "No. outgoing stable ptcls from vrtx", "Number of Events"); // outVrtxPtclNumber = dqm.book1dHisto("outVrtxPtclNumber", "No. outgoing ptcls from vrtx", 30, 0, 30, "No. outgoing ptcls from vrtx", "Number of Events"); vrtxZ = dqm.book1dHisto("VrtxZ", "VrtxZ", 50, -250, 250, "Z_{Vtx}", "Number of Events"); vrtxRadius = dqm.book1dHisto("vrtxRadius", "vrtxRadius", 50, 0, 50, "R_{vtx}", "Number of Events"); // unknownPDTNumber = dqm.book1dHisto("unknownPDTNumber", "Log10(No. unknown ptcls PDT)", 60, -1, 5, "log_{10}(No. unknown ptcls PDT)", "Number of Events"); //Log genPtclStatus = dqm.book1dHisto("genPtclStatus", "Status of genParticle", 200, 0, 200., "", "Number of Events"); // Bjorken_x = dqm.book1dHisto("Bjorken_x", "Bjorken_x", 1000, 0.0, 1.0, "Bjorken_{x}", "Number of Events"); pdf_u = dqm.book1dHisto( "pdf_u", "Log10(PDF(u,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_ubar = dqm.book1dHisto("pdf_ubar", "Log10(PDF(ubar,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_d = dqm.book1dHisto( "pdf_d", "Log10(PDF(d,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_dbar = dqm.book1dHisto("pdf_dbar", "Log10(PDF(dbar,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_ssbar = dqm.book1dHisto("pdf_ssbar", "Log10(PDF(ssbar,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_ccbar = dqm.book1dHisto("pdf_ccbar", "Log10(PDF(ccbar,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_bbbar = dqm.book1dHisto("pdf_bbbar", "Log10(PDF(bbbar,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log pdf_g = dqm.book1dHisto( "pdf_g", "Log10(PDF(g,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events"); //Log scalePDF = dqm.book1dHisto("scalePDF", "Log10(Q-scale(GeV))", logQScaleNbin, logQScaleMin, logQScaleMax, "log_{10}(Q-scale(GeV))", "Number of Events"); parton1Id = dqm.book1dHisto("parton1Id", "ID of parton 1", 45, -14.5, 30.5, "ID", "Number of Events"); parton2Id = dqm.book1dHisto("parton2Id", "ID of parton 2", 45, -14.5, 30.5, "ID", "Number of Events"); // status1ShortLived = dqm.book1dHisto("status1ShortLived", "Status 1 short lived", 11, 0, 11, "", "Number of Events"); status1ShortLived->setBinLabel(1, "d/dbar"); status1ShortLived->setBinLabel(2, "u/ubar"); status1ShortLived->setBinLabel(3, "s/sbar"); status1ShortLived->setBinLabel(4, "c/cbar"); status1ShortLived->setBinLabel(5, "b/bbar"); status1ShortLived->setBinLabel(6, "t/tbar"); status1ShortLived->setBinLabel(7, "g"); status1ShortLived->setBinLabel(8, "tau-/tau+"); status1ShortLived->setBinLabel(9, "Z0"); status1ShortLived->setBinLabel(10, "W-/W+"); status1ShortLived->setBinLabel(11, "PDG = 7,8,17,25-99"); log10DeltaEcms = dqm.book1dHisto("DeltaEcms1log10", "log_{10} of deviation from nominal Ecms", 200, -1., 5., "log_{10}(#DeltaE) (log_{10}(GeV))", "Number of Events"); DeltaEcms = dqm.book1dHisto("DeltaEcms1", "deviation from nominal Ecms", 200, -1., 1., "#DeltaE (GeV)", "Number of Events"); DeltaPx = dqm.book1dHisto("DeltaPx1", "deviation from nominal Px", 200, -1., 1., "#DeltaP_{x} (GeV)", "Number of Events"); DeltaPy = dqm.book1dHisto("DeltaPy1", "deviation from nominal Py", 200, -1., 1., "#DeltaP_{y} (GeV)", "Number of Events"); DeltaPz = dqm.book1dHisto("DeltaPz1", "deviation from nominal Pz", 200, -1., 1., "#DeltaP_{z} (GeV)", "Number of Events"); return; } void BasicHepMCValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) { ///counters to zero for every event int partonNum = 0; // int outVrtxStablePtclNum = 0; int stablePtclNum = 0; int otherPtclNum = 0; int unknownPDTNum = 0; int stableChaNum = 0; // double bjorken = 0.; double logPdf1 = 0.; double logPdf2 = 0.; double logQScale = 0.; // double etotal = 0.; double pxtotal = 0.; double pytotal = 0.; double pztotal = 0.; ///Gathering the HepMCProduct information edm::Handle<HepMCProduct> evt; iEvent.getByToken(hepmcCollectionToken_, evt); //Get EVENT HepMC::GenEvent const *myGenEvent = evt->GetEvent(); double weight = wmanager_.weight(iEvent); nEvt->Fill(0.5, weight); genPtclNumber->Fill(log10(myGenEvent->particles_size()), weight); genVrtxNumber->Fill(log10(myGenEvent->vertices_size()), weight); ///PDF informations HepMC::PdfInfo const *pdf = myGenEvent->pdf_info(); if (pdf) { bjorken = ((pdf->x1()) / ((pdf->x1()) + (pdf->x2()))); logQScale = log10(pdf->scalePDF()); if (logQScale > logQScaleMax) logQScale = logQScaleMax - 0.5 * logQScaleBinsize; // visualize overflow & underflow in the histograms if (logQScale < logQScaleMin) logQScale = logQScaleMin + 0.5 * logQScaleBinsize; logPdf1 = log10(pdf->pdf1()); if (logPdf1 > logPdfMax) logPdf1 = logPdfMax - 0.5 * logPdfBinsize; if (logPdf1 < logPdfMin) logPdf1 = logPdfMin + 0.5 * logPdfBinsize; logPdf2 = log10(pdf->pdf2()); if (logPdf2 > logPdfMax) logPdf2 = logPdfMax - 0.5 * logPdfBinsize; if (logPdf2 < logPdfMin) logPdf2 = logPdfMin + 0.5 * logPdfBinsize; Bjorken_x->Fill(bjorken, weight); scalePDF->Fill(logQScale, weight); parton1Id->Fill((double)pdf->id1(), weight); parton2Id->Fill((double)pdf->id2(), weight); if (pdf->id1() == 2) pdf_u->Fill(logPdf1, weight); if (pdf->id2() == 2) pdf_u->Fill(logPdf2, weight); if (pdf->id1() == -2) pdf_ubar->Fill(logPdf1, weight); if (pdf->id2() == -2) pdf_ubar->Fill(logPdf2, weight); if (pdf->id1() == 1) pdf_d->Fill(logPdf1, weight); if (pdf->id2() == 1) pdf_d->Fill(logPdf2, weight); if (pdf->id1() == -1) pdf_dbar->Fill(logPdf1, weight); if (pdf->id2() == -1) pdf_dbar->Fill(logPdf2, weight); if (std::abs(pdf->id1()) == 3) pdf_ssbar->Fill(logPdf1, weight); if (std::abs(pdf->id2()) == 3) pdf_ssbar->Fill(logPdf2, weight); if (std::abs(pdf->id1()) == 4) pdf_ccbar->Fill(logPdf1, weight); if (std::abs(pdf->id2()) == 4) pdf_ccbar->Fill(logPdf2, weight); if (std::abs(pdf->id1()) == 5) pdf_bbbar->Fill(logPdf1, weight); if (std::abs(pdf->id2()) == 5) pdf_bbbar->Fill(logPdf2, weight); if (std::abs(pdf->id1()) == 21) pdf_g->Fill(logPdf1, weight); if (std::abs(pdf->id2()) == 21) pdf_g->Fill(logPdf2, weight); } //Looping through the VERTICES in the event HepMC::GenEvent::vertex_const_iterator vrtxBegin = myGenEvent->vertices_begin(); HepMC::GenEvent::vertex_const_iterator vrtxEnd = myGenEvent->vertices_end(); unsigned int nvtx(0); for (HepMC::GenEvent::vertex_const_iterator vrtxIt = vrtxBegin; vrtxIt != vrtxEnd; ++vrtxIt) { ///Vertices HepMC::GenVertex const *vrtx = *vrtxIt; outVrtxPtclNumber->Fill(vrtx->particles_out_size(), weight); //std::cout << "all " << vrtx->particles_out_size() << '\n'; if (nvtx == 0) { vrtxZ->Fill(vrtx->point3d().z(), weight); vrtxRadius->Fill(vrtx->point3d().perp(), weight); } ///loop on vertex particles HepMC::GenVertex::particles_out_const_iterator vrtxPtclBegin = vrtx->particles_out_const_begin(); HepMC::GenVertex::particles_out_const_iterator vrtxPtclEnd = vrtx->particles_out_const_end(); outVrtxStablePtclNum = 0; for (HepMC::GenVertex::particles_out_const_iterator vrtxPtclIt = vrtxPtclBegin; vrtxPtclIt != vrtxPtclEnd; ++vrtxPtclIt) { HepMC::GenParticle const *vrtxPtcl = *vrtxPtclIt; if (vrtxPtcl->status() == 1) { ++outVrtxStablePtclNum; //std::cout << "stable " << outVrtxStablePtclNum << '\n'; } } outVrtxStablePtclNumber->Fill(outVrtxStablePtclNum, weight); nvtx++; } //vertices ///Looping through the PARTICLES in the event HepMC::GenEvent::particle_const_iterator ptclBegin = myGenEvent->particles_begin(); HepMC::GenEvent::particle_const_iterator ptclEnd = myGenEvent->particles_end(); for (HepMC::GenEvent::particle_const_iterator ptclIt = ptclBegin; ptclIt != ptclEnd; ++ptclIt) { ///Particles HepMC::GenParticle const *ptcl = *ptclIt; int Id = ptcl->pdg_id(); // std::cout << Id << '\n'; float Log_p = log10(ptcl->momentum().rho()); double charge = 999.; // for the charge it's needed a HepPDT method int status = ptcl->status(); const HepPDT::ParticleData *PData = fPDGTable->particle(HepPDT::ParticleID(Id)); if (PData == nullptr) { // std::cout << "Unknown id = " << Id << '\n'; ++unknownPDTNum; } else charge = PData->charge(); ///Status statistics genPtclStatus->Fill((float)status, weight); ///Stable particles if (ptcl->status() == 1) { ++stablePtclNum; stablePtclPhi->Fill(ptcl->momentum().phi() / CLHEP::degree, weight); //std::cout << ptcl->polarization().phi() << '\n'; stablePtclEta->Fill(ptcl->momentum().pseudoRapidity(), weight); stablePtclCharge->Fill(charge, weight); // std::cout << ptclData.charge() << '\n'; stablePtclp->Fill(Log_p, weight); stablePtclpT->Fill(log10(ptcl->momentum().perp()), weight); if (charge != 0. && charge != 999.) ++stableChaNum; if (std::abs(Id) == 1) status1ShortLived->Fill(1, weight); if (std::abs(Id) == 2) status1ShortLived->Fill(2, weight); if (std::abs(Id) == 3) status1ShortLived->Fill(3, weight); if (std::abs(Id) == 4) status1ShortLived->Fill(4, weight); if (std::abs(Id) == 5) status1ShortLived->Fill(5, weight); if (std::abs(Id) == 6) status1ShortLived->Fill(6, weight); if (Id == 21) status1ShortLived->Fill(7, weight); if (std::abs(Id) == 15) status1ShortLived->Fill(8, weight); if (Id == 23) status1ShortLived->Fill(9, weight); if (std::abs(Id) == 24) status1ShortLived->Fill(10, weight); if (std::abs(Id) == 7 || std::abs(Id) == 8 || std::abs(Id) == 17 || (std::abs(Id) >= 25 && std::abs(Id) <= 99)) status1ShortLived->Fill(11, weight); etotal += ptcl->momentum().e(); pxtotal += ptcl->momentum().px(); pytotal += ptcl->momentum().py(); pztotal += ptcl->momentum().pz(); } if (abs(Id) < 6 || abs(Id) == 22) { ++partonNum; partonpT->Fill(Log_p, weight); } bool indentified = false; for (unsigned int i = 0; i < particles.size(); i++) { if (particles.at(i).Fill(ptcl, weight)) { indentified = true; break; } } if (!indentified) { ++otherPtclNum; otherPtclMomentum->Fill(Log_p, weight); } } //event particles // set a default sqrt(s) and then check in the event double ecms = 13000.; if (myGenEvent->valid_beam_particles()) { ecms = myGenEvent->beam_particles().first->momentum().e() + myGenEvent->beam_particles().second->momentum().e(); } log10DeltaEcms->Fill(log10(fabs(etotal - ecms)), weight); DeltaEcms->Fill(etotal - ecms, weight); DeltaPx->Fill(pxtotal, weight); DeltaPy->Fill(pytotal, weight); DeltaPz->Fill(pztotal, weight); ///filling multiplicity ME's stablePtclNumber->Fill(log10(stablePtclNum + 0.1), weight); stableChaNumber->Fill(log10(stableChaNum + 0.1), weight); otherPtclNumber->Fill(log10(otherPtclNum + 0.1), weight); unknownPDTNumber->Fill(log10(unknownPDTNum + 0.1), weight); // partonNumber->Fill(partonNum, weight); for (unsigned int i = 0; i < particles.size(); i++) { particles.at(i).FillCount(weight); }; } //analyze