Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:10 

0001 /*class BasicHepMCValidation
0002  *  
0003  *  Class to fill dqm monitor elements from existing EDM file
0004  *
0005  */
0006 
0007 #include "Validation/EventGenerator/interface/BasicHepMCValidation.h"
0008 
0009 #include "CLHEP/Units/defs.h"
0010 #include "CLHEP/Units/PhysicalConstants.h"
0011 #include "Validation/EventGenerator/interface/DQMHelper.h"
0012 using namespace edm;
0013 
0014 BasicHepMCValidation::BasicHepMCValidation(const edm::ParameterSet &iPSet)
0015     : wmanager_(iPSet, consumesCollector()), hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")) {
0016   hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
0017   fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0018 }
0019 
0020 BasicHepMCValidation::~BasicHepMCValidation() {}
0021 
0022 void BasicHepMCValidation::dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) {
0023   fPDGTable = c.getHandle(fPDGTableToken);
0024 }
0025 
0026 namespace {
0027   // Set upper bound & lower bound for PDF & Scale related histograms
0028   constexpr double logPdfMax = 3.0;
0029   constexpr double logPdfMin = -3.0;
0030   constexpr int logPdfNbin = 150;
0031   constexpr double logPdfBinsize = (logPdfMax - logPdfMin) / (double)logPdfNbin;
0032   constexpr double logQScaleMax = 4.0;
0033   constexpr double logQScaleMin = -1.0;
0034   constexpr int logQScaleNbin = 500;
0035   constexpr double logQScaleBinsize = (logQScaleMax - logQScaleMin) / (double)logQScaleNbin;
0036 }  // namespace
0037 
0038 void BasicHepMCValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) {
0039   ///Setting the DQM top directories
0040   DQMHelper dqm(&i);
0041   i.setCurrentFolder("Generator/Particles");
0042 
0043   // Number of analyzed events
0044   nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1.);
0045 
0046   ///Booking the ME's
0047   ///multiplicity
0048   // quarks
0049   particles.push_back(ParticleMonitor("u", 1, i));
0050   particles.push_back(ParticleMonitor("ubar", -1, i));
0051   particles.push_back(ParticleMonitor("d", 2, i));
0052   particles.push_back(ParticleMonitor("dbar", -2, i));
0053   particles.push_back(ParticleMonitor("s", 3, i));
0054   particles.push_back(ParticleMonitor("sbar", -3, i));
0055   particles.push_back(ParticleMonitor("c", 4, i));
0056   particles.push_back(ParticleMonitor("cbar", -4, i));
0057   particles.push_back(ParticleMonitor("b", 5, i));
0058   particles.push_back(ParticleMonitor("bbar", -5, i));
0059   particles.push_back(ParticleMonitor("t", 6, i));
0060   particles.push_back(ParticleMonitor("tbar", -6, i));
0061 
0062   //leptons
0063   particles.push_back(ParticleMonitor("eminus", 11, i));
0064   particles.push_back(ParticleMonitor("eplus", -11, i));
0065   particles.push_back(ParticleMonitor("nue", 12, i));
0066   particles.push_back(ParticleMonitor("nuebar", -12, i));
0067   particles.push_back(ParticleMonitor("muminus", 13, i));
0068   particles.push_back(ParticleMonitor("muplus", -13, i));
0069   particles.push_back(ParticleMonitor("numu", 14, i));
0070   particles.push_back(ParticleMonitor("numubar", -14, i));
0071   particles.push_back(ParticleMonitor("tauminus", 15, i));
0072   particles.push_back(ParticleMonitor("tauplus", -15, i));
0073   particles.push_back(ParticleMonitor("nutau", 16, i));
0074   particles.push_back(ParticleMonitor("nutaubar", -16, i));
0075 
0076   //bosons
0077   particles.push_back(ParticleMonitor("Wplus", 24, i));
0078   particles.push_back(ParticleMonitor("Wminus", -24, i));
0079   particles.push_back(ParticleMonitor("Z", 23, i));
0080   particles.push_back(ParticleMonitor("gamma", 22, i));
0081   particles.push_back(ParticleMonitor("gluon", 21, i));
0082 
0083   //mesons
0084   particles.push_back(ParticleMonitor("piplus", 211, i, true));    //log
0085   particles.push_back(ParticleMonitor("piminus", -211, i, true));  //log
0086   particles.push_back(ParticleMonitor("pizero", 111, i, true));    //log
0087   particles.push_back(ParticleMonitor("Kplus", 321, i));
0088   particles.push_back(ParticleMonitor("Kminus", -321, i));
0089   particles.push_back(ParticleMonitor("Klzero", 130, i));
0090   particles.push_back(ParticleMonitor("Kszero", 310, i));
0091 
0092   //baryons
0093   particles.push_back(ParticleMonitor("p", 2212, i, true));      //log
0094   particles.push_back(ParticleMonitor("pbar", -2212, i, true));  //log
0095   particles.push_back(ParticleMonitor("n", 2112, i, true));      //log
0096   particles.push_back(ParticleMonitor("nbar", -2112, i, true));  //log
0097   particles.push_back(ParticleMonitor("lambda0", 3122, i));
0098   particles.push_back(ParticleMonitor("lambda0bar", -3122, i));
0099 
0100   //D mesons
0101   particles.push_back(ParticleMonitor("Dplus", 411, i));
0102   particles.push_back(ParticleMonitor("Dminus", -411, i));
0103   particles.push_back(ParticleMonitor("Dzero", 421, i));
0104   particles.push_back(ParticleMonitor("Dzerobar", -421, i));
0105 
0106   //B mesons
0107   particles.push_back(ParticleMonitor("Bplus", 521, i));
0108   particles.push_back(ParticleMonitor("Bminus", -521, i));
0109   particles.push_back(ParticleMonitor("Bzero", 511, i));
0110   particles.push_back(ParticleMonitor("Bzerobar", -511, i));
0111   particles.push_back(ParticleMonitor("Bszero", 531, i));
0112   particles.push_back(ParticleMonitor("Bszerobar", -531, i));
0113 
0114   //
0115   otherPtclNumber = dqm.book1dHisto(
0116       "otherPtclNumber", "Log10(No. other ptcls)", 60, -1, 5, "log_{10}(No. other ptcls)", "Number of Events");  //Log
0117   otherPtclMomentum = dqm.book1dHisto("otherPtclMomentum",
0118                                       "Log10(p) other ptcls",
0119                                       60,
0120                                       -2,
0121                                       4,
0122                                       "log10(P^{other ptcls}) (log_{10}(GeV))",
0123                                       "Number of Events");
0124 
0125   ///other
0126   genPtclNumber = dqm.book1dHisto(
0127       "genPtclNumber", "Log10(No. all particles)", 60, -1, 5, "log10(No. all particles)", "Number of Events");  //Log
0128   genVrtxNumber = dqm.book1dHisto(
0129       "genVrtxNumber", "Log10(No. all vertexs)", 60, -1, 5, "log10(No. all vertexs)", "Number of Events");  //Log
0130   //
0131   stablePtclNumber = dqm.book1dHisto("stablePtclNumber",
0132                                      "Log10(No. stable particles)",
0133                                      50,
0134                                      0,
0135                                      5,
0136                                      "log10(No. stable particles)",
0137                                      "Number of Events");  //Log
0138   stablePtclPhi = dqm.book1dHisto(
0139       "stablePtclPhi", "stable Ptcl Phi", 360, -180, 180, "#phi^{stable Ptcl} (rad)", "Number of Events");
0140   stablePtclEta = dqm.book1dHisto(
0141       "stablePtclEta", "stable Ptcl Eta (pseudo rapidity)", 220, -11, 11, "#eta^{stable Ptcl}", "Number of Events");
0142   stablePtclCharge =
0143       dqm.book1dHisto("stablePtclCharge", "stablePtclCharge", 5, -2, 2, "Charge^{stable ptcls}", "Number of Events");
0144   stableChaNumber = dqm.book1dHisto("stableChaNumber",
0145                                     "Log10(No. stable charged particles)",
0146                                     50,
0147                                     0,
0148                                     5,
0149                                     "log_{10}(No. stable charged particles)",
0150                                     "Number of Events");  //Log
0151   stablePtclp = dqm.book1dHisto("stablePtclp",
0152                                 "Log10(p) stable ptcl p",
0153                                 80,
0154                                 -4,
0155                                 4,
0156                                 "log_{10}(P^{stable ptcl}) (log_{10}(GeV))",
0157                                 "Number of Events");  //Log
0158   stablePtclpT = dqm.book1dHisto("stablePtclpT",
0159                                  "Log10(pT) stable ptcl pT",
0160                                  80,
0161                                  -4,
0162                                  4,
0163                                  "log_{10}(P_{t}^{stable ptcl}) (log_{10}(GeV))",
0164                                  "Number of Events");  //Log
0165   partonNumber =
0166       dqm.book1dHisto("partonNumber", "number of partons", 100, 0, 100, "number of partons", "Number of Events");
0167   partonpT =
0168       dqm.book1dHisto("partonpT", "Log10(pT) parton pT", 80, -4, 4, "Log10(P_{t}^{parton})", "Number of Events");  //Log
0169   outVrtxStablePtclNumber = dqm.book1dHisto("outVrtxStablePtclNumber",
0170                                             "No. outgoing stable ptcls from vrtx",
0171                                             10,
0172                                             0,
0173                                             10,
0174                                             "No. outgoing stable ptcls from vrtx",
0175                                             "Number of Events");
0176   //
0177   outVrtxPtclNumber = dqm.book1dHisto("outVrtxPtclNumber",
0178                                       "No. outgoing ptcls from vrtx",
0179                                       30,
0180                                       0,
0181                                       30,
0182                                       "No. outgoing ptcls from vrtx",
0183                                       "Number of Events");
0184   vrtxZ = dqm.book1dHisto("VrtxZ", "VrtxZ", 50, -250, 250, "Z_{Vtx}", "Number of Events");
0185   vrtxRadius = dqm.book1dHisto("vrtxRadius", "vrtxRadius", 50, 0, 50, "R_{vtx}", "Number of Events");
0186   //
0187   unknownPDTNumber = dqm.book1dHisto("unknownPDTNumber",
0188                                      "Log10(No. unknown ptcls PDT)",
0189                                      60,
0190                                      -1,
0191                                      5,
0192                                      "log_{10}(No. unknown ptcls PDT)",
0193                                      "Number of Events");  //Log
0194   genPtclStatus = dqm.book1dHisto("genPtclStatus", "Status of genParticle", 200, 0, 200., "", "Number of Events");
0195   //
0196   Bjorken_x = dqm.book1dHisto("Bjorken_x", "Bjorken_x", 1000, 0.0, 1.0, "Bjorken_{x}", "Number of Events");
0197   pdf_u = dqm.book1dHisto(
0198       "pdf_u", "Log10(PDF(u,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events");  //Log
0199   pdf_ubar = dqm.book1dHisto("pdf_ubar",
0200                              "Log10(PDF(ubar,x,Q))",
0201                              logPdfNbin,
0202                              logPdfMin,
0203                              logPdfMax,
0204                              "log_{10}(x*f(x))",
0205                              "Number of Events");  //Log
0206   pdf_d = dqm.book1dHisto(
0207       "pdf_d", "Log10(PDF(d,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events");  //Log
0208   pdf_dbar = dqm.book1dHisto("pdf_dbar",
0209                              "Log10(PDF(dbar,x,Q))",
0210                              logPdfNbin,
0211                              logPdfMin,
0212                              logPdfMax,
0213                              "log_{10}(x*f(x))",
0214                              "Number of Events");  //Log
0215   pdf_ssbar = dqm.book1dHisto("pdf_ssbar",
0216                               "Log10(PDF(ssbar,x,Q))",
0217                               logPdfNbin,
0218                               logPdfMin,
0219                               logPdfMax,
0220                               "log_{10}(x*f(x))",
0221                               "Number of Events");  //Log
0222   pdf_ccbar = dqm.book1dHisto("pdf_ccbar",
0223                               "Log10(PDF(ccbar,x,Q))",
0224                               logPdfNbin,
0225                               logPdfMin,
0226                               logPdfMax,
0227                               "log_{10}(x*f(x))",
0228                               "Number of Events");  //Log
0229   pdf_bbbar = dqm.book1dHisto("pdf_bbbar",
0230                               "Log10(PDF(bbbar,x,Q))",
0231                               logPdfNbin,
0232                               logPdfMin,
0233                               logPdfMax,
0234                               "log_{10}(x*f(x))",
0235                               "Number of Events");  //Log
0236   pdf_g = dqm.book1dHisto(
0237       "pdf_g", "Log10(PDF(g,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events");  //Log
0238   scalePDF = dqm.book1dHisto("scalePDF",
0239                              "Log10(Q-scale(GeV))",
0240                              logQScaleNbin,
0241                              logQScaleMin,
0242                              logQScaleMax,
0243                              "log_{10}(Q-scale(GeV))",
0244                              "Number of Events");
0245   parton1Id = dqm.book1dHisto("parton1Id", "ID of parton 1", 45, -14.5, 30.5, "ID", "Number of Events");
0246   parton2Id = dqm.book1dHisto("parton2Id", "ID of parton 2", 45, -14.5, 30.5, "ID", "Number of Events");
0247   //
0248   status1ShortLived = dqm.book1dHisto("status1ShortLived", "Status 1 short lived", 11, 0, 11, "", "Number of Events");
0249   status1ShortLived->setBinLabel(1, "d/dbar");
0250   status1ShortLived->setBinLabel(2, "u/ubar");
0251   status1ShortLived->setBinLabel(3, "s/sbar");
0252   status1ShortLived->setBinLabel(4, "c/cbar");
0253   status1ShortLived->setBinLabel(5, "b/bbar");
0254   status1ShortLived->setBinLabel(6, "t/tbar");
0255   status1ShortLived->setBinLabel(7, "g");
0256   status1ShortLived->setBinLabel(8, "tau-/tau+");
0257   status1ShortLived->setBinLabel(9, "Z0");
0258   status1ShortLived->setBinLabel(10, "W-/W+");
0259   status1ShortLived->setBinLabel(11, "PDG = 7,8,17,25-99");
0260 
0261   log10DeltaEcms = dqm.book1dHisto("DeltaEcms1log10",
0262                                    "log_{10} of deviation from nominal Ecms",
0263                                    200,
0264                                    -1.,
0265                                    5.,
0266                                    "log_{10}(#DeltaE) (log_{10}(GeV))",
0267                                    "Number of Events");
0268   DeltaEcms =
0269       dqm.book1dHisto("DeltaEcms1", "deviation from nominal Ecms", 200, -1., 1., "#DeltaE (GeV)", "Number of Events");
0270   DeltaPx =
0271       dqm.book1dHisto("DeltaPx1", "deviation from nominal Px", 200, -1., 1., "#DeltaP_{x} (GeV)", "Number of Events");
0272   DeltaPy =
0273       dqm.book1dHisto("DeltaPy1", "deviation from nominal Py", 200, -1., 1., "#DeltaP_{y} (GeV)", "Number of Events");
0274   DeltaPz =
0275       dqm.book1dHisto("DeltaPz1", "deviation from nominal Pz", 200, -1., 1., "#DeltaP_{z} (GeV)", "Number of Events");
0276   return;
0277 }
0278 
0279 void BasicHepMCValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0280   ///counters to zero for every event
0281   int partonNum = 0;
0282   //
0283   int outVrtxStablePtclNum = 0;
0284   int stablePtclNum = 0;
0285   int otherPtclNum = 0;
0286   int unknownPDTNum = 0;
0287   int stableChaNum = 0;
0288   //
0289   double bjorken = 0.;
0290   double logPdf1 = 0.;
0291   double logPdf2 = 0.;
0292   double logQScale = 0.;
0293   //
0294   double etotal = 0.;
0295   double pxtotal = 0.;
0296   double pytotal = 0.;
0297   double pztotal = 0.;
0298 
0299   ///Gathering the HepMCProduct information
0300   edm::Handle<HepMCProduct> evt;
0301   iEvent.getByToken(hepmcCollectionToken_, evt);
0302 
0303   //Get EVENT
0304   HepMC::GenEvent const *myGenEvent = evt->GetEvent();
0305 
0306   double weight = wmanager_.weight(iEvent);
0307 
0308   nEvt->Fill(0.5, weight);
0309 
0310   genPtclNumber->Fill(log10(myGenEvent->particles_size()), weight);
0311   genVrtxNumber->Fill(log10(myGenEvent->vertices_size()), weight);
0312 
0313   ///PDF informations
0314   HepMC::PdfInfo const *pdf = myGenEvent->pdf_info();
0315   if (pdf) {
0316     bjorken = ((pdf->x1()) / ((pdf->x1()) + (pdf->x2())));
0317     logQScale = log10(pdf->scalePDF());
0318     if (logQScale > logQScaleMax)
0319       logQScale = logQScaleMax - 0.5 * logQScaleBinsize;  // visualize overflow & underflow in the histograms
0320     if (logQScale < logQScaleMin)
0321       logQScale = logQScaleMin + 0.5 * logQScaleBinsize;
0322     logPdf1 = log10(pdf->pdf1());
0323     if (logPdf1 > logPdfMax)
0324       logPdf1 = logPdfMax - 0.5 * logPdfBinsize;
0325     if (logPdf1 < logPdfMin)
0326       logPdf1 = logPdfMin + 0.5 * logPdfBinsize;
0327     logPdf2 = log10(pdf->pdf2());
0328     if (logPdf2 > logPdfMax)
0329       logPdf2 = logPdfMax - 0.5 * logPdfBinsize;
0330     if (logPdf2 < logPdfMin)
0331       logPdf2 = logPdfMin + 0.5 * logPdfBinsize;
0332     Bjorken_x->Fill(bjorken, weight);
0333     scalePDF->Fill(logQScale, weight);
0334     parton1Id->Fill((double)pdf->id1(), weight);
0335     parton2Id->Fill((double)pdf->id2(), weight);
0336     if (pdf->id1() == 2)
0337       pdf_u->Fill(logPdf1, weight);
0338     if (pdf->id2() == 2)
0339       pdf_u->Fill(logPdf2, weight);
0340     if (pdf->id1() == -2)
0341       pdf_ubar->Fill(logPdf1, weight);
0342     if (pdf->id2() == -2)
0343       pdf_ubar->Fill(logPdf2, weight);
0344     if (pdf->id1() == 1)
0345       pdf_d->Fill(logPdf1, weight);
0346     if (pdf->id2() == 1)
0347       pdf_d->Fill(logPdf2, weight);
0348     if (pdf->id1() == -1)
0349       pdf_dbar->Fill(logPdf1, weight);
0350     if (pdf->id2() == -1)
0351       pdf_dbar->Fill(logPdf2, weight);
0352     if (std::abs(pdf->id1()) == 3)
0353       pdf_ssbar->Fill(logPdf1, weight);
0354     if (std::abs(pdf->id2()) == 3)
0355       pdf_ssbar->Fill(logPdf2, weight);
0356     if (std::abs(pdf->id1()) == 4)
0357       pdf_ccbar->Fill(logPdf1, weight);
0358     if (std::abs(pdf->id2()) == 4)
0359       pdf_ccbar->Fill(logPdf2, weight);
0360     if (std::abs(pdf->id1()) == 5)
0361       pdf_bbbar->Fill(logPdf1, weight);
0362     if (std::abs(pdf->id2()) == 5)
0363       pdf_bbbar->Fill(logPdf2, weight);
0364     if (std::abs(pdf->id1()) == 21)
0365       pdf_g->Fill(logPdf1, weight);
0366     if (std::abs(pdf->id2()) == 21)
0367       pdf_g->Fill(logPdf2, weight);
0368   }
0369 
0370   //Looping through the VERTICES in the event
0371   HepMC::GenEvent::vertex_const_iterator vrtxBegin = myGenEvent->vertices_begin();
0372   HepMC::GenEvent::vertex_const_iterator vrtxEnd = myGenEvent->vertices_end();
0373   unsigned int nvtx(0);
0374   for (HepMC::GenEvent::vertex_const_iterator vrtxIt = vrtxBegin; vrtxIt != vrtxEnd; ++vrtxIt) {
0375     ///Vertices
0376     HepMC::GenVertex const *vrtx = *vrtxIt;
0377     outVrtxPtclNumber->Fill(vrtx->particles_out_size(), weight);
0378     //std::cout << "all " << vrtx->particles_out_size() << '\n';
0379 
0380     if (nvtx == 0) {
0381       vrtxZ->Fill(vrtx->point3d().z(), weight);
0382       vrtxRadius->Fill(vrtx->point3d().perp(), weight);
0383     }
0384     ///loop on vertex particles
0385     HepMC::GenVertex::particles_out_const_iterator vrtxPtclBegin = vrtx->particles_out_const_begin();
0386     HepMC::GenVertex::particles_out_const_iterator vrtxPtclEnd = vrtx->particles_out_const_end();
0387     outVrtxStablePtclNum = 0;
0388     for (HepMC::GenVertex::particles_out_const_iterator vrtxPtclIt = vrtxPtclBegin; vrtxPtclIt != vrtxPtclEnd;
0389          ++vrtxPtclIt) {
0390       HepMC::GenParticle const *vrtxPtcl = *vrtxPtclIt;
0391       if (vrtxPtcl->status() == 1) {
0392         ++outVrtxStablePtclNum;
0393         //std::cout << "stable " << outVrtxStablePtclNum << '\n';
0394       }
0395     }
0396     outVrtxStablePtclNumber->Fill(outVrtxStablePtclNum, weight);
0397     nvtx++;
0398   }  //vertices
0399 
0400   ///Looping through the PARTICLES in the event
0401   HepMC::GenEvent::particle_const_iterator ptclBegin = myGenEvent->particles_begin();
0402   HepMC::GenEvent::particle_const_iterator ptclEnd = myGenEvent->particles_end();
0403   for (HepMC::GenEvent::particle_const_iterator ptclIt = ptclBegin; ptclIt != ptclEnd; ++ptclIt) {
0404     ///Particles
0405     HepMC::GenParticle const *ptcl = *ptclIt;
0406     int Id = ptcl->pdg_id();  // std::cout << Id << '\n';
0407     float Log_p = log10(ptcl->momentum().rho());
0408     double charge = 999.;  // for the charge it's needed a HepPDT method
0409     int status = ptcl->status();
0410     const HepPDT::ParticleData *PData = fPDGTable->particle(HepPDT::ParticleID(Id));
0411     if (PData == nullptr) {
0412       //        std::cout << "Unknown id = " << Id << '\n';
0413       ++unknownPDTNum;
0414     } else
0415       charge = PData->charge();
0416 
0417     ///Status statistics
0418     genPtclStatus->Fill((float)status, weight);
0419 
0420     ///Stable particles
0421     if (ptcl->status() == 1) {
0422       ++stablePtclNum;
0423       stablePtclPhi->Fill(ptcl->momentum().phi() / CLHEP::degree,
0424                           weight);  //std::cout << ptcl->polarization().phi() << '\n';
0425       stablePtclEta->Fill(ptcl->momentum().pseudoRapidity(), weight);
0426       stablePtclCharge->Fill(charge, weight);  // std::cout << ptclData.charge() << '\n';
0427       stablePtclp->Fill(Log_p, weight);
0428       stablePtclpT->Fill(log10(ptcl->momentum().perp()), weight);
0429       if (charge != 0. && charge != 999.)
0430         ++stableChaNum;
0431       if (std::abs(Id) == 1)
0432         status1ShortLived->Fill(1, weight);
0433       if (std::abs(Id) == 2)
0434         status1ShortLived->Fill(2, weight);
0435       if (std::abs(Id) == 3)
0436         status1ShortLived->Fill(3, weight);
0437       if (std::abs(Id) == 4)
0438         status1ShortLived->Fill(4, weight);
0439       if (std::abs(Id) == 5)
0440         status1ShortLived->Fill(5, weight);
0441       if (std::abs(Id) == 6)
0442         status1ShortLived->Fill(6, weight);
0443       if (Id == 21)
0444         status1ShortLived->Fill(7, weight);
0445       if (std::abs(Id) == 15)
0446         status1ShortLived->Fill(8, weight);
0447       if (Id == 23)
0448         status1ShortLived->Fill(9, weight);
0449       if (std::abs(Id) == 24)
0450         status1ShortLived->Fill(10, weight);
0451       if (std::abs(Id) == 7 || std::abs(Id) == 8 || std::abs(Id) == 17 || (std::abs(Id) >= 25 && std::abs(Id) <= 99))
0452         status1ShortLived->Fill(11, weight);
0453       etotal += ptcl->momentum().e();
0454       pxtotal += ptcl->momentum().px();
0455       pytotal += ptcl->momentum().py();
0456       pztotal += ptcl->momentum().pz();
0457     }
0458 
0459     if (abs(Id) < 6 || abs(Id) == 22) {
0460       ++partonNum;
0461       partonpT->Fill(Log_p, weight);
0462     }
0463 
0464     bool indentified = false;
0465     for (unsigned int i = 0; i < particles.size(); i++) {
0466       if (particles.at(i).Fill(ptcl, weight)) {
0467         indentified = true;
0468         break;
0469       }
0470     }
0471     if (!indentified) {
0472       ++otherPtclNum;
0473       otherPtclMomentum->Fill(Log_p, weight);
0474     }
0475   }  //event particles
0476 
0477   // set a default sqrt(s) and then check in the event
0478   double ecms = 13000.;
0479   if (myGenEvent->valid_beam_particles()) {
0480     ecms = myGenEvent->beam_particles().first->momentum().e() + myGenEvent->beam_particles().second->momentum().e();
0481   }
0482   log10DeltaEcms->Fill(log10(fabs(etotal - ecms)), weight);
0483   DeltaEcms->Fill(etotal - ecms, weight);
0484   DeltaPx->Fill(pxtotal, weight);
0485   DeltaPy->Fill(pytotal, weight);
0486   DeltaPz->Fill(pztotal, weight);
0487 
0488   ///filling multiplicity ME's
0489   stablePtclNumber->Fill(log10(stablePtclNum + 0.1), weight);
0490   stableChaNumber->Fill(log10(stableChaNum + 0.1), weight);
0491   otherPtclNumber->Fill(log10(otherPtclNum + 0.1), weight);
0492   unknownPDTNumber->Fill(log10(unknownPDTNum + 0.1), weight);
0493   //
0494   partonNumber->Fill(partonNum, weight);
0495   for (unsigned int i = 0; i < particles.size(); i++) {
0496     particles.at(i).FillCount(weight);
0497   };
0498 
0499 }  //analyze
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!