9 #include "CLHEP/Units/defs.h"
10 #include "CLHEP/Units/PhysicalConstants.h"
15 : wmanager_(iPSet, consumesCollector()), hepmcCollection_(iPSet.getParameter<
edm::
InputTag>(
"hepmcCollection")) {
25 constexpr
double logPdfMax = 3.0;
26 constexpr
double logPdfMin = -3.0;
27 constexpr
int logPdfNbin = 150;
28 constexpr
double logPdfBinsize = (logPdfMax - logPdfMin) / (
double)logPdfNbin;
29 constexpr
double logQScaleMax = 4.0;
30 constexpr
double logQScaleMin = -1.0;
31 constexpr
int logQScaleNbin = 500;
32 constexpr
double logQScaleBinsize = (logQScaleMax - logQScaleMin) / (
double)logQScaleNbin;
38 i.setCurrentFolder(
"Generator/Particles");
41 nEvt =
dqm.book1dHisto(
"nEvt",
"n analyzed Events", 1, 0., 1.);
113 "otherPtclNumber",
"Log10(No. other ptcls)", 60, -1, 5,
"log_{10}(No. other ptcls)",
"Number of Events");
115 "Log10(p) other ptcls",
119 "log10(P^{other ptcls}) (log_{10}(GeV))",
124 "genPtclNumber",
"Log10(No. all particles)", 60, -1, 5,
"log10(No. all particles)",
"Number of Events");
126 "genVrtxNumber",
"Log10(No. all vertexs)", 60, -1, 5,
"log10(No. all vertexs)",
"Number of Events");
129 "Log10(No. stable particles)",
133 "log10(No. stable particles)",
136 "stablePtclPhi",
"stable Ptcl Phi", 360, -180, 180,
"#phi^{stable Ptcl} (rad)",
"Number of Events");
138 "stablePtclEta",
"stable Ptcl Eta (pseudo rapidity)", 220, -11, 11,
"#eta^{stable Ptcl}",
"Number of Events");
140 dqm.book1dHisto(
"stablePtclCharge",
"stablePtclCharge", 5, -2, 2,
"Charge^{stable ptcls}",
"Number of Events");
142 "Log10(No. stable charged particles)",
146 "log_{10}(No. stable charged particles)",
149 "Log10(p) stable ptcl p",
153 "log_{10}(P^{stable ptcl}) (log_{10}(GeV))",
156 "Log10(pT) stable ptcl pT",
160 "log_{10}(P_{t}^{stable ptcl}) (log_{10}(GeV))",
163 dqm.book1dHisto(
"partonNumber",
"number of partons", 100, 0, 100,
"number of partons",
"Number of Events");
165 dqm.book1dHisto(
"partonpT",
"Log10(pT) parton pT", 80, -4, 4,
"Log10(P_{t}^{parton})",
"Number of Events");
167 "No. outgoing stable ptcls from vrtx",
171 "No. outgoing stable ptcls from vrtx",
175 "No. outgoing ptcls from vrtx",
179 "No. outgoing ptcls from vrtx",
181 vrtxZ =
dqm.book1dHisto(
"VrtxZ",
"VrtxZ", 50, -250, 250,
"Z_{Vtx}",
"Number of Events");
182 vrtxRadius =
dqm.book1dHisto(
"vrtxRadius",
"vrtxRadius", 50, 0, 50,
"R_{vtx}",
"Number of Events");
185 "Log10(No. unknown ptcls PDT)",
189 "log_{10}(No. unknown ptcls PDT)",
191 genPtclStatus =
dqm.book1dHisto(
"genPtclStatus",
"Status of genParticle", 200, 0, 200.,
"",
"Number of Events");
193 Bjorken_x =
dqm.book1dHisto(
"Bjorken_x",
"Bjorken_x", 1000, 0.0, 1.0,
"Bjorken_{x}",
"Number of Events");
195 "pdf_u",
"Log10(PDF(u,x,Q))", logPdfNbin, logPdfMin, logPdfMax,
"log_{10}(x*f(x))",
"Number of Events");
197 "Log10(PDF(ubar,x,Q))",
204 "pdf_d",
"Log10(PDF(d,x,Q))", logPdfNbin, logPdfMin, logPdfMax,
"log_{10}(x*f(x))",
"Number of Events");
206 "Log10(PDF(dbar,x,Q))",
213 "Log10(PDF(ssbar,x,Q))",
220 "Log10(PDF(ccbar,x,Q))",
227 "Log10(PDF(bbbar,x,Q))",
234 "pdf_g",
"Log10(PDF(g,x,Q))", logPdfNbin, logPdfMin, logPdfMax,
"log_{10}(x*f(x))",
"Number of Events");
236 "Log10(Q-scale(GeV))",
240 "log_{10}(Q-scale(GeV))",
242 parton1Id =
dqm.book1dHisto(
"parton1Id",
"ID of parton 1", 45, -14.5, 30.5,
"ID",
"Number of Events");
243 parton2Id =
dqm.book1dHisto(
"parton2Id",
"ID of parton 2", 45, -14.5, 30.5,
"ID",
"Number of Events");
245 status1ShortLived =
dqm.book1dHisto(
"status1ShortLived",
"Status 1 short lived", 11, 0, 11,
"",
"Number of Events");
259 "log_{10} of deviation from nominal Ecms",
263 "log_{10}(#DeltaE) (log_{10}(GeV))",
266 dqm.book1dHisto(
"DeltaEcms1",
"deviation from nominal Ecms", 200, -1., 1.,
"#DeltaE (GeV)",
"Number of Events");
268 dqm.book1dHisto(
"DeltaPx1",
"deviation from nominal Px", 200, -1., 1.,
"#DeltaP_{x} (GeV)",
"Number of Events");
270 dqm.book1dHisto(
"DeltaPy1",
"deviation from nominal Py", 200, -1., 1.,
"#DeltaP_{y} (GeV)",
"Number of Events");
272 dqm.book1dHisto(
"DeltaPz1",
"deviation from nominal Pz", 200, -1., 1.,
"#DeltaP_{z} (GeV)",
"Number of Events");
280 int outVrtxStablePtclNum = 0;
281 int stablePtclNum = 0;
282 int otherPtclNum = 0;
283 int unknownPDTNum = 0;
284 int stableChaNum = 0;
289 double logQScale = 0.;
311 HepMC::PdfInfo
const *pdf = myGenEvent->pdf_info();
313 bjorken = ((pdf->x1()) / ((pdf->x1()) + (pdf->x2())));
314 logQScale = log10(pdf->scalePDF());
315 if (logQScale > logQScaleMax)
316 logQScale = logQScaleMax - 0.5 * logQScaleBinsize;
317 if (logQScale < logQScaleMin)
318 logQScale = logQScaleMin + 0.5 * logQScaleBinsize;
319 logPdf1 = log10(pdf->pdf1());
320 if (logPdf1 > logPdfMax)
321 logPdf1 = logPdfMax - 0.5 * logPdfBinsize;
322 if (logPdf1 < logPdfMin)
323 logPdf1 = logPdfMin + 0.5 * logPdfBinsize;
324 logPdf2 = log10(pdf->pdf2());
325 if (logPdf2 > logPdfMax)
326 logPdf2 = logPdfMax - 0.5 * logPdfBinsize;
327 if (logPdf2 < logPdfMin)
328 logPdf2 = logPdfMin + 0.5 * logPdfBinsize;
337 if (pdf->id1() == -2)
339 if (pdf->id2() == -2)
345 if (pdf->id1() == -1)
347 if (pdf->id2() == -1)
368 HepMC::GenEvent::vertex_const_iterator vrtxBegin = myGenEvent->vertices_begin();
369 HepMC::GenEvent::vertex_const_iterator vrtxEnd = myGenEvent->vertices_end();
370 unsigned int nvtx(0);
371 for (HepMC::GenEvent::vertex_const_iterator vrtxIt = vrtxBegin; vrtxIt != vrtxEnd; ++vrtxIt) {
373 HepMC::GenVertex
const *vrtx = *vrtxIt;
382 HepMC::GenVertex::particles_out_const_iterator vrtxPtclBegin = vrtx->particles_out_const_begin();
383 HepMC::GenVertex::particles_out_const_iterator vrtxPtclEnd = vrtx->particles_out_const_end();
384 outVrtxStablePtclNum = 0;
385 for (HepMC::GenVertex::particles_out_const_iterator vrtxPtclIt = vrtxPtclBegin; vrtxPtclIt != vrtxPtclEnd;
388 if (vrtxPtcl->status() == 1) {
389 ++outVrtxStablePtclNum;
398 HepMC::GenEvent::particle_const_iterator ptclBegin = myGenEvent->particles_begin();
399 HepMC::GenEvent::particle_const_iterator ptclEnd = myGenEvent->particles_end();
400 for (HepMC::GenEvent::particle_const_iterator ptclIt = ptclBegin; ptclIt != ptclEnd; ++ptclIt) {
403 int Id = ptcl->pdg_id();
404 float Log_p = log10(ptcl->momentum().rho());
406 int status = ptcl->status();
408 if (PData ==
nullptr) {
418 if (ptcl->status() == 1) {
450 etotal += ptcl->momentum().e();
451 pxtotal += ptcl->momentum().px();
452 pytotal += ptcl->momentum().py();
453 pztotal += ptcl->momentum().pz();
456 if (
abs(Id) < 6 ||
abs(Id) == 22) {
461 bool indentified =
false;
475 double ecms = 13000.;
476 if (myGenEvent->valid_beam_particles()) {
477 ecms = myGenEvent->beam_particles().first->momentum().e() + myGenEvent->beam_particles().second->momentum().e();