EventGenerator/interface/BasicHepMCValidation.h

0001 #ifndef BASICHEPMCVALIDATION_H
0002 #define BASICHEPMCVALIDATION_H
0003
0004 /*class BasicHepMCValidation
0005  *
0006  *  Class to fill Event Generator dqm monitor elements; works on HepMCProduct
0007  *
0008  *
0009  */
0010
0011 // framework & common header files
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/Run.h"
0015
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020
0021 //DQM services
0022 #include "DQMServices/Core/interface/DQMStore.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0025 #include "Validation/EventGenerator/interface/DQMHelper.h"
0026
0027 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0028
0029 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0030
0031 #include "Validation/EventGenerator/interface/WeightManager.h"
0032 #include "TVector3.h"
0033
0034 class BasicHepMCValidation : public DQMEDAnalyzer {
0035 public:
0036   explicit BasicHepMCValidation(const edm::ParameterSet &);
0037   ~BasicHepMCValidation() override;
0038
0039   void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override;
0040   void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override;
0041   void analyze(edm::Event const &, edm::EventSetup const &) override;
0042
0043 private:
0044   WeightManager wmanager_;
0045   edm::InputTag hepmcCollection_;
0046
0047   /// PDT table
0048   edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable;
0049   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> fPDGTableToken;
0050
0051   class ParticleMonitor {
0052   public:
0053     ParticleMonitor(std::string name_, int pdgid_, DQMStore::IBooker &i, bool nlog_ = false)
0054         : name(name_), pdgid(pdgid_), count(0), nlog(nlog_) {
0055       DQMHelper dqm(&i);
0056       // Number of analyzed events
0057       if (!nlog) {
0058         numberPerEvent = dqm.book1dHisto(
0059             name + "Number", "Number of  " + name + "'s per event", 20, 0, 20, "No. of " + name, "Number of Events");
0060       } else {
0061         numberPerEvent = dqm.book1dHisto(name + "Number",
0062                                          "Number of  " + name + "'s per event",
0063                                          20,
0064                                          0,
0065                                          20,
0066                                          "log_{10}(No. of " + name + ")",
0067                                          "Number of Events");
0068       }
0069       p_init = dqm.book1dHisto(name + "Momentum",
0070                                "log_{10}(P) of the " + name + "s",
0071                                60,
0072                                -2,
0073                                4,
0074                                "log_{10}(P) (log_{10}(GeV))",
0075                                "Number of " + name);
0076
0077       eta_init = dqm.book1dHisto(name + "Eta", "#eta of the " + name + "s", 100, -5., 5., "#eta", "Number of " + name);
0078
0079       lifetime_init = dqm.book1dHisto(name + "LifeTime",
0080                                       "#phi of the " + name + "s",
0081                                       100,
0082                                       -15,
0083                                       -5,
0084                                       "Log_{10}(life-time^{final}) (log_{10}(s))",
0085                                       "Number of " + name);
0086
0087       p_final = dqm.book1dHisto(name + "MomentumFinal",
0088                                 "log_{10}(P^{final}) of " + name + "s at end of decay chain",
0089                                 60,
0090                                 -2,
0091                                 4,
0092                                 "log_{10}(P^{final}) (log_{10}(GeV))",
0093                                 "Number of " + name);
0094
0095       lifetime_final = dqm.book1dHisto(name + "LifeTimeFinal",
0096                                        "Log_{10}(life-time^{final}) of " + name + "s at end of decay chain",
0097                                        100,
0098                                        -15,
0099                                        -5,
0100                                        "Log_{10}(life-time^{final}) (log_{10}(s))",
0101                                        "Number of " + name);
0102     }
0103
0104     ~ParticleMonitor() {}
0105
0106     bool Fill(const HepMC::GenParticle *p, double weight) {
0107       if (p->pdg_id() == pdgid) {
0108         if (isFirst(p)) {
0109           p_init->Fill(log10(p->momentum().rho()), weight);
0110           eta_init->Fill(p->momentum().eta(), weight);
0111           const HepMC::GenParticle *pf = GetFinal(p);  // inlcude mixing
0112           p_final->Fill(log10(pf->momentum().rho()), weight);
0113           // compute lifetime...
0114           if (p->production_vertex() && p->end_vertex()) {
0115             TVector3 PV(p->production_vertex()->point3d().x(),
0116                         p->production_vertex()->point3d().y(),
0117                         p->production_vertex()->point3d().z());
0118             TVector3 SV(p->end_vertex()->point3d().x(), p->end_vertex()->point3d().y(), p->end_vertex()->point3d().z());
0119             TVector3 DL = SV - PV;
0120             double c(2.99792458E8), Ltau(DL.Mag() / 100) /*cm->m*/, beta(p->momentum().rho() / p->momentum().m());
0121             double lt = Ltau / (c * beta);
0122             if (lt > 1E-16)
0123               lifetime_init->Fill(log10(lt), weight);
0124             if (pf->end_vertex()) {
0125               TVector3 SVf(
0126                   pf->end_vertex()->point3d().x(), pf->end_vertex()->point3d().y(), pf->end_vertex()->point3d().z());
0127               DL = SVf - PV;
0128               Ltau = DL.Mag() / 100;
0129               lt = Ltau / (c * beta);
0130               if (lt > 1E-16)
0131                 lifetime_final->Fill(log10(lt), weight);
0132             }
0133           }
0134           count++;
0135         }
0136         return true;
0137       }
0138       return false;
0139     }
0140
0141     void FillCount(double weight) {
0142       if (nlog)
0143         numberPerEvent->Fill(log10(count), weight);
0144       else
0145         numberPerEvent->Fill(count, weight);
0146       count = 0;
0147     }
0148
0149     int PDGID() { return pdgid; }
0150
0151   private:
0152     bool isFirst(const HepMC::GenParticle *p) {
0153       if (p->production_vertex()) {
0154         for (HepMC::GenVertex::particles_in_const_iterator m = p->production_vertex()->particles_in_const_begin();
0155              m != p->production_vertex()->particles_in_const_end();
0156              m++) {
0157           if (abs((*m)->pdg_id()) == abs(p->pdg_id()))
0158             return false;
0159         }
0160       }
0161       return true;
0162     }
0163
0164     // includes mixing (assuming mixing is not occurring more than 5 times back and forth)
0165     const HepMC::GenParticle *GetFinal(const HepMC::GenParticle *p) {
0166       const HepMC::GenParticle *aPart = p;
0167       for (unsigned int iMix = 0; iMix < 10; iMix++) {
0168         bool foundSimilar = false;
0169         if (aPart->end_vertex()) {
0170           if (aPart->end_vertex()->particles_out_size() != 0) {
0171             for (HepMC::GenVertex::particles_out_const_iterator d = aPart->end_vertex()->particles_out_const_begin();
0172                  d != aPart->end_vertex()->particles_out_const_end();
0173                  d++) {
0174               if (abs((*d)->pdg_id()) == abs(aPart->pdg_id())) {
0175                 aPart = *d;
0176                 foundSimilar = true;
0177                 break;
0178               }
0179             }
0180           }
0181           if (!foundSimilar)
0182             break;
0183         }
0184       }
0185       return aPart;
0186     }
0187
0188     std::string name;
0189     int pdgid;
0190     unsigned int count;
0191     bool nlog;
0192     MonitorElement *p_init, *p_final, *eta_init, *lifetime_init, *lifetime_final, *numberPerEvent;
0193   };
0194
0195   MonitorElement *nEvt;
0196   std::vector<ParticleMonitor> particles;
0197
0198   ///other ME's
0199   MonitorElement *otherPtclNumber;
0200   MonitorElement *otherPtclMomentum;
0201   MonitorElement *genPtclNumber;
0202   MonitorElement *genVrtxNumber;
0203   MonitorElement *unknownPDTNumber;
0204   MonitorElement *outVrtxPtclNumber;
0205   MonitorElement *genPtclStatus;
0206   //
0207   MonitorElement *stablePtclNumber;
0208   MonitorElement *stableChaNumber;
0209   MonitorElement *stablePtclPhi;
0210   MonitorElement *stablePtclEta;
0211   MonitorElement *stablePtclCharge;
0212   MonitorElement *stablePtclp;
0213   MonitorElement *stablePtclpT;
0214   MonitorElement *partonNumber;
0215   MonitorElement *partonpT;
0216   MonitorElement *outVrtxStablePtclNumber;
0217   //
0218   MonitorElement *vrtxZ;
0219   MonitorElement *vrtxRadius;
0220   //
0221   MonitorElement *Bjorken_x;
0222   MonitorElement *pdf_u;
0223   MonitorElement *pdf_ubar;
0224   MonitorElement *pdf_d;
0225   MonitorElement *pdf_dbar;
0226   MonitorElement *pdf_ssbar;
0227   MonitorElement *pdf_ccbar;
0228   MonitorElement *pdf_bbbar;
0229   MonitorElement *pdf_g;
0230   MonitorElement *scalePDF;
0231   MonitorElement *parton1Id;
0232   MonitorElement *parton2Id;
0233
0234   MonitorElement *status1ShortLived;
0235
0236   MonitorElement *log10DeltaEcms;
0237   MonitorElement *DeltaEcms;
0238   MonitorElement *DeltaPx;
0239   MonitorElement *DeltaPy;
0240   MonitorElement *DeltaPz;
0241
0242   edm::EDGetTokenT<edm::HepMCProduct> hepmcCollectionToken_;
0243 };
0244
0245 #endif