EventGenerator/plugins/DrellYanValidation.cc

0001 /*class DrellYanValidation
0002  *
0003  *  Class to fill dqm monitor elements from existing EDM file
0004  *
0005  */
0006
0007 #include "Validation/EventGenerator/interface/DrellYanValidation.h"
0008 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010 #include "TLorentzVector.h"
0011
0012 #include "CLHEP/Units/defs.h"
0013 #include "CLHEP/Units/PhysicalConstants.h"
0014
0015 #include "DataFormats/Math/interface/LorentzVector.h"
0016 #include "Validation/EventGenerator/interface/DQMHelper.h"
0017 using namespace edm;
0018
0019 DrellYanValidation::DrellYanValidation(const edm::ParameterSet& iPSet)
0020     : wmanager_(iPSet, consumesCollector()),
0021       hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
0022       _flavor(iPSet.getParameter<int>("decaysTo")),
0023       _name(iPSet.getParameter<std::string>("name")) {
0024   hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
0025   fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0026 }
0027
0028 DrellYanValidation::~DrellYanValidation() {}
0029
0030 void DrellYanValidation::dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) {
0031   fPDGTable = c.getHandle(fPDGTableToken);
0032 }
0033
0034 void DrellYanValidation::bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) {
0035   ///Setting the DQM top directories
0036   std::string folderName = "Generator/DrellYan";
0037   folderName += _name;
0038   DQMHelper dqm(&i);
0039   i.setCurrentFolder(folderName);
0040
0041   // Number of analyzed events
0042   nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., "bin", "Number of Events");
0043
0044   //Kinematics
0045   Zmass = dqm.book1dHisto("Zmass", "inv. Mass Z", 70, 0, 140, "M_{Z} (GeV)", "Number of Events");
0046   ZmassPeak = dqm.book1dHisto("ZmassPeak", "inv. Mass Z", 80, 80, 100, "M_{Z} (GeV)", "Number of Events");
0047   Zpt = dqm.book1dHisto("Zpt", "Z pt", 100, 0, 200, "P_{t}^{Z} (GeV)", "Number of Events");
0048   ZptLog =
0049       dqm.book1dHisto("ZptLog", "log(Z pt)", 100, 0., 5., "log_{10}(P_{t}^{Z}) (log_{10}(GeV))", "Number of Events");
0050   Zrap = dqm.book1dHisto("Zrap", "Z y", 100, -5, 5, "Y_{Z}", "Number of Events");
0051   Zdaughters = dqm.book1dHisto("Zdaughters", "Z daughters", 60, -30, 30, "Z daughters (PDG ID)", "Number of Events");
0052
0053   dilep_mass = dqm.book1dHisto("dilep_mass", "inv. Mass dilepton", 70, 0, 140, "M_{ll} (GeV)", "Number of Events");
0054   dilep_massPeak =
0055       dqm.book1dHisto("dilep_massPeak", "inv. Mass dilepton", 80, 80, 100, "M_{ll} (GeV)", "Number of Events");
0056   dilep_pt = dqm.book1dHisto("dilep_pt", "dilepton pt", 100, 0, 200, "P_{t}^{ll} (GeV)", "Number of Events");
0057   dilep_ptLog = dqm.book1dHisto(
0058       "dilep_ptLog", "log(dilepton pt)", 100, 0., 5., "log_{10}(P_{t}^{ll}) (log_{10}(GeV))", "Number of Events");
0059   dilep_rap = dqm.book1dHisto("dilep_rap", "dilepton y", 100, -5, 5, "Y_{ll}", "Number of Events");
0060
0061   gamma_energy = dqm.book1dHisto("gamma_energy",
0062                                  "photon energy in Z rest frame",
0063                                  200,
0064                                  0.,
0065                                  100.,
0066                                  "E_{#gamma}^{Z rest-frame} (GeV)",
0067                                  "Number of Events");
0068   cos_theta_gamma_lepton = dqm.book1dHisto("cos_theta_gamma_lepton",
0069                                            "cos_theta_gamma_lepton in Z rest frame",
0070                                            200,
0071                                            -1,
0072                                            1,
0073                                            "cos(#theta_{#gamma-lepton}^{Z rest-frame})",
0074                                            "Number of Events");
0075
0076   leadpt = dqm.book1dHisto("leadpt", "leading lepton pt", 200, 0., 200., "P_{t}^{1st-lepton}", "Number of Events");
0077   secpt = dqm.book1dHisto("secpt", "second lepton pt", 200, 0., 200., "P_{t}^{2nd-lepton}", "Number of Events");
0078   leadeta = dqm.book1dHisto("leadeta", "leading lepton eta", 100, -5., 5., "#eta^{1st-lepton}", "Number of Events");
0079   seceta = dqm.book1dHisto("seceta", "second lepton eta", 100, -5., 5., "#eta^{2nd-lepton}", "Number of Events");
0080
0081   return;
0082 }
0083
0084 void DrellYanValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0085   // we *DO NOT* rely on a Z entry in the particle listings!
0086
0087   ///Gathering the HepMCProduct information
0088   edm::Handle<HepMCProduct> evt;
0089   iEvent.getByToken(hepmcCollectionToken_, evt);
0090
0091   //Get EVENT
0092   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0093
0094   double weight = wmanager_.weight(iEvent);
0095
0096   //std::cout << "weight: " << weight << std::endl;
0097
0098   nEvt->Fill(0.5, weight);
0099
0100   std::vector<const HepMC::GenParticle*> allproducts;
0101
0102   //requires status 1 for leptons and neutrinos (except tau)
0103   int requiredstatus = (std::abs(_flavor) == 11 || std::abs(_flavor) == 12 || std::abs(_flavor) == 13 ||
0104                         std::abs(_flavor) == 14 || std::abs(_flavor) == 16)
0105                            ? 1
0106                            : 3;
0107
0108   bool vetotau =
0109       true;  //(std::abs(_flavor) == 11 || std::abs(_flavor) == 12 || std::abs(_flavor) ==13 || std::abs(_flavor) ==14 || std::abs(_flavor) ==16) ? true : false;
0110
0111   for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
0112        iter != myGenEvent->particles_end();
0113        ++iter) {
0114     if (vetotau) {
0115       if ((*iter)->status() == 3 && std::abs((*iter)->pdg_id()) == 15)
0116         return;
0117     }
0118     if ((*iter)->status() == requiredstatus) {
0119       if (std::abs((*iter)->pdg_id()) == _flavor)
0120         allproducts.push_back(*iter);
0121     }
0122   }
0123
0124   //nothing to do if we don't have 2 particles
0125   if (allproducts.size() < 2)
0126     return;
0127
0128   //sort them in pt
0129   std::sort(allproducts.begin(), allproducts.end(), HepMCValidationHelper::sortByPt);
0130
0131   //get the first element and the first following element with opposite charge
0132   std::vector<const HepMC::GenParticle*> products;
0133   products.push_back(allproducts.front());
0134   const HepPDT::ParticleData* PData1 = fPDGTable->particle(HepPDT::ParticleID(allproducts.front()->pdg_id()));
0135   double charge1 = PData1->charge();
0136   for (unsigned int i = 1; i < allproducts.size(); ++i) {
0137     const HepPDT::ParticleData* PData2 = fPDGTable->particle(HepPDT::ParticleID(allproducts[i]->pdg_id()));
0138     double charge2 = PData2->charge();
0139     if (charge1 * charge2 < 0)
0140       products.push_back(allproducts[i]);
0141   }
0142
0143   //if we did not find any opposite charge pair there is nothing to do
0144   if (products.size() < 2)
0145     return;
0146
0147   assert(products[0]->momentum().perp() >= products[1]->momentum().perp());
0148
0149   //leading lepton with pt > 20.
0150   if (products[0]->momentum().perp() < 20.)
0151     return;
0152
0153   //assemble FourMomenta
0154   TLorentzVector lep1(products[0]->momentum().x(),
0155                       products[0]->momentum().y(),
0156                       products[0]->momentum().z(),
0157                       products[0]->momentum().t());
0158   TLorentzVector lep2(products[1]->momentum().x(),
0159                       products[1]->momentum().y(),
0160                       products[1]->momentum().z(),
0161                       products[1]->momentum().t());
0162   TLorentzVector dilepton_mom = lep1 + lep2;
0163   TLorentzVector dilepton_andphoton_mom = dilepton_mom;
0164
0165   //mass > 60.
0166   if (dilepton_mom.M() < 60.)
0167     return;
0168
0169   //find possible qed fsr photons
0170   std::vector<const HepMC::GenParticle*> fsrphotons;
0171   HepMCValidationHelper::findFSRPhotons(products, myGenEvent, 0.1, fsrphotons);
0172
0173   Zdaughters->Fill(products[0]->pdg_id(), weight);
0174   Zdaughters->Fill(products[1]->pdg_id(), weight);
0175
0176   std::vector<TLorentzVector> gammasMomenta;
0177   for (unsigned int ipho = 0; ipho < fsrphotons.size(); ++ipho) {
0178     TLorentzVector phomom(fsrphotons[ipho]->momentum().x(),
0179                           fsrphotons[ipho]->momentum().y(),
0180                           fsrphotons[ipho]->momentum().z(),
0181                           fsrphotons[ipho]->momentum().t());
0182     dilepton_andphoton_mom += phomom;
0183     Zdaughters->Fill(fsrphotons[ipho]->pdg_id(), weight);
0184     gammasMomenta.push_back(phomom);
0185   }
0186   //Fill Z histograms
0187   Zmass->Fill(dilepton_andphoton_mom.M(), weight);
0188   ZmassPeak->Fill(dilepton_andphoton_mom.M(), weight);
0189   Zpt->Fill(dilepton_andphoton_mom.Pt(), weight);
0190   ZptLog->Fill(log10(dilepton_andphoton_mom.Pt()), weight);
0191   Zrap->Fill(dilepton_andphoton_mom.Rapidity(), weight);
0192
0193   //Fill dilepton histograms
0194   dilep_mass->Fill(dilepton_mom.M(), weight);
0195   dilep_massPeak->Fill(dilepton_mom.M(), weight);
0196   dilep_pt->Fill(dilepton_mom.Pt(), weight);
0197   dilep_ptLog->Fill(log10(dilepton_mom.Pt()), weight);
0198   dilep_rap->Fill(dilepton_mom.Rapidity(), weight);
0199
0200   //Fill lepton histograms
0201   leadpt->Fill(lep1.Pt(), weight);
0202   secpt->Fill(lep2.Pt(), weight);
0203   leadeta->Fill(lep1.Eta(), weight);
0204   seceta->Fill(lep2.Eta(), weight);
0205
0206   //boost everything in the Z frame
0207   TVector3 boost = dilepton_andphoton_mom.BoostVector();
0208   boost *= -1.;
0209   lep1.Boost(boost);
0210   lep2.Boost(boost);
0211   for (unsigned int ipho = 0; ipho < gammasMomenta.size(); ++ipho) {
0212     gammasMomenta[ipho].Boost(boost);
0213   }
0214   std::sort(gammasMomenta.begin(), gammasMomenta.end(), HepMCValidationHelper::GreaterByE<TLorentzVector>);
0215
0216   //fill gamma histograms
0217   if (!gammasMomenta.empty() && dilepton_andphoton_mom.M() > 50.) {
0218     gamma_energy->Fill(gammasMomenta.front().E(), weight);
0219     double dphi = lep1.DeltaR(gammasMomenta.front()) < lep2.DeltaR(gammasMomenta.front())
0220                       ? lep1.DeltaPhi(gammasMomenta.front())
0221                       : lep2.DeltaPhi(gammasMomenta.front());
0222     cos_theta_gamma_lepton->Fill(cos(dphi), weight);
0223   }
0224
0225 }  //analyze