EventGenerator/plugins/WValidation.cc

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