EventGenerator/plugins/HiggsValidation.cc

0001 /*class HiggsValidation
0002  *
0003  *  Class to fill dqm monitor elements from existing EDM file
0004  *
0005  */
0006
0007 #include "Validation/EventGenerator/interface/HiggsValidation.h"
0008
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "CLHEP/Units/defs.h"
0012 #include "CLHEP/Units/PhysicalConstants.h"
0013
0014 #include "DataFormats/Math/interface/LorentzVector.h"
0015
0016 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0017 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
0018 #include "Validation/EventGenerator/interface/DQMHelper.h"
0019 using namespace edm;
0020
0021 HiggsValidation::HiggsValidation(const edm::ParameterSet& iPSet)
0022     : wmanager_(iPSet, consumesCollector()),
0023       hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
0024       particle_id(iPSet.getParameter<int>("pdg_id")),
0025       particle_name(iPSet.getParameter<std::string>("particleName")) {
0026   monitoredDecays = new MonitoredDecays(iPSet);
0027   hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
0028   fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0029 }
0030
0031 HiggsValidation::~HiggsValidation() {}
0032
0033 void HiggsValidation::dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) {
0034   fPDGTable = c.getHandle(fPDGTableToken);
0035 }
0036
0037 void HiggsValidation::bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) {
0038   ///Setting the DQM top directories
0039   std::string dir = "Generator/";
0040   dir += particle_name;
0041   DQMHelper dqm(&i);
0042   i.setCurrentFolder(dir);
0043
0044   // Number of analyzed events
0045   nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., "bin", "Number of Events");
0046
0047   //decay type
0048
0049   std::string channel = particle_name + "_DecayChannels";
0050   HiggsDecayChannels = dqm.book1dHisto(channel,
0051                                        particle_name + " decay channels",
0052                                        monitoredDecays->size(),
0053                                        0,
0054                                        monitoredDecays->size(),
0055                                        "Decay Channels",
0056                                        "Number of Events");
0057
0058   for (size_t j = 0; j < monitoredDecays->size(); ++j) {
0059     HiggsDecayChannels->setBinLabel(1 + j, monitoredDecays->channel(j));
0060   }
0061
0062   //Kinematics
0063   Higgs_pt = dqm.book1dHisto((particle_name + "_pt"),
0064                              (particle_name + " p_{t}"),
0065                              50,
0066                              0,
0067                              250,
0068                              "P_{t}^{" + particle_name + "} (GeV)",
0069                              "Number of Events");
0070   Higgs_eta = dqm.book1dHisto((particle_name + "_eta"),
0071                               (particle_name + " #eta"),
0072                               50,
0073                               -5,
0074                               5,
0075                               "#eta^{" + particle_name + "}",
0076                               "Number of Events");
0077   Higgs_mass = dqm.book1dHisto(
0078       (particle_name + "_m"), (particle_name + " M"), 500, 0, 500, "M^{" + particle_name + "}", "Number of Events");
0079
0080   int idx = 0;
0081   for (unsigned int j = 0; j < monitoredDecays->NDecayParticles(); j++) {
0082     HiggsDecayProd_pt.push_back(dqm.book1dHisto((monitoredDecays->ConvertIndex(idx) + "_pt"),
0083                                                 (monitoredDecays->ConvertIndex(idx) + " p_{t}"),
0084                                                 50,
0085                                                 0,
0086                                                 250,
0087                                                 "P_{t}^{" + monitoredDecays->ConvertIndex(idx) + "} (GeV)",
0088                                                 "Number of Events"));
0089     HiggsDecayProd_eta.push_back(dqm.book1dHisto((monitoredDecays->ConvertIndex(idx) + "_eta"),
0090                                                  (monitoredDecays->ConvertIndex(idx) + " #eta"),
0091                                                  50,
0092                                                  -5,
0093                                                  5,
0094                                                  "#eta^{" + monitoredDecays->ConvertIndex(idx) + "}",
0095                                                  "Number of Events"));
0096     idx++;
0097   }
0098
0099   return;
0100 }
0101
0102 void HiggsValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0103   double weight = wmanager_.weight(iEvent);
0104   nEvt->Fill(0.5, weight);
0105
0106   //Gathering the HepMCProduct information
0107   edm::Handle<HepMCProduct> evt;
0108   iEvent.getByToken(hepmcCollectionToken_, evt);
0109
0110   //Get EVENT
0111   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
0112
0113   // loop over all particles
0114   bool filled = false;
0115   for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
0116        iter != myGenEvent->particles_end() && !filled;
0117        ++iter) {
0118     if (particle_id == std::abs((*iter)->pdg_id())) {
0119       std::vector<HepMC::GenParticle*> decayprod;
0120       int channel = findHiggsDecayChannel(*iter, decayprod);
0121       HiggsDecayChannels->Fill(channel, weight);
0122       Higgs_pt->Fill((*iter)->momentum().perp(), weight);
0123       Higgs_eta->Fill((*iter)->momentum().eta(), weight);
0124       Higgs_mass->Fill((*iter)->momentum().m(), weight);
0125       for (unsigned int i = 0; i < decayprod.size(); i++) {
0126         int idx = monitoredDecays->isDecayParticle(decayprod.at(i)->pdg_id());
0127         if (0 <= idx && idx <= (int)HiggsDecayProd_pt.size()) {
0128           HiggsDecayProd_pt.at(idx)->Fill(decayprod.at(i)->momentum().perp(), weight);
0129           HiggsDecayProd_eta.at(idx)->Fill(decayprod.at(i)->momentum().eta(), weight);
0130         }
0131       }
0132       filled = true;
0133     }
0134   }
0135
0136   delete myGenEvent;
0137
0138 }  //analyze
0139
0140 int HiggsValidation::findHiggsDecayChannel(const HepMC::GenParticle* genParticle,
0141                                            std::vector<HepMC::GenParticle*>& decayprod) {
0142   if (genParticle->status() == 1)
0143     return monitoredDecays->stable();
0144   std::vector<int> children;
0145   if (genParticle->end_vertex()) {
0146     HepMC::GenVertex::particle_iterator des;
0147     for (des = genParticle->end_vertex()->particles_begin(HepMC::descendants);
0148          des != genParticle->end_vertex()->particles_end(HepMC::descendants);
0149          ++des) {
0150       if ((*des)->pdg_id() == genParticle->pdg_id())
0151         continue;
0152
0153       HepMC::GenVertex::particle_iterator mother = (*des)->production_vertex()->particles_begin(HepMC::parents);
0154       if ((*mother)->pdg_id() == genParticle->pdg_id()) {
0155         children.push_back((*des)->pdg_id());
0156         decayprod.push_back((*des));
0157       }
0158     }
0159   }
0160
0161   if (children.size() == 2 && children.at(0) != 0 && children.at(1) != 0) {
0162     return monitoredDecays->position(children.at(0), children.at(1));
0163   }
0164   return monitoredDecays->undetermined();
0165 }