Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "Validation/EventGenerator/interface/GenWeightValidation.h"
0002 
0003 #include "DataFormats/Math/interface/LorentzVector.h"
0004 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
0005 
0006 using namespace edm;
0007 
0008 GenWeightValidation::GenWeightValidation(const edm::ParameterSet& iPSet)
0009     : wmanager_(iPSet, consumesCollector()),
0010       genParticleToken_(consumes<reco::GenParticleCollection>(iPSet.getParameter<edm::InputTag>("genParticles"))),
0011       genJetToken_(consumes<reco::GenJetCollection>(iPSet.getParameter<edm::InputTag>("genJets"))),
0012       idxGenEvtInfo_(iPSet.getParameter<int>("whichGenEventInfo")),
0013       idxFSRup_(iPSet.getParameter<int>("idxFSRup")),
0014       idxFSRdown_(iPSet.getParameter<int>("idxFSRdown")),
0015       idxISRup_(iPSet.getParameter<int>("idxISRup")),
0016       idxISRdown_(iPSet.getParameter<int>("idxISRdown")),
0017       leadLepPtNbin_(iPSet.getParameter<int>("leadLepPtNbin")),
0018       rapidityNbin_(iPSet.getParameter<int>("rapidityNbin")),
0019       leadLepPtRange_(iPSet.getParameter<double>("leadLepPtRange")),
0020       leadLepPtCut_(iPSet.getParameter<double>("leadLepPtCut")),
0021       lepEtaCut_(iPSet.getParameter<double>("lepEtaCut")),
0022       rapidityRange_(iPSet.getParameter<double>("rapidityRange")),
0023       nJetsNbin_(iPSet.getParameter<int>("nJetsNbin")),
0024       jetPtNbin_(iPSet.getParameter<int>("jetPtNbin")),
0025       jetPtCut_(iPSet.getParameter<double>("jetPtCut")),
0026       jetEtaCut_(iPSet.getParameter<double>("jetEtaCut")),
0027       jetPtRange_(iPSet.getParameter<double>("jetPtRange")) {
0028   std::vector<int> idxs = {idxFSRup_, idxFSRdown_, idxISRup_, idxISRdown_};
0029   std::sort(idxs.begin(), idxs.end(), std::greater<int>());
0030   idxMax_ = idxs.at(0);
0031 }
0032 
0033 void GenWeightValidation::bookHistograms(DQMStore::IBooker& iBook, edm::Run const&, edm::EventSetup const&) {
0034   ///Setting the DQM top directories
0035   std::string folderName = "Generator/GenWeight";
0036   DQMHelper aDqmHelper(&iBook);
0037   iBook.setCurrentFolder(folderName);
0038 
0039   // Number of analyzed events
0040   nEvt_ = aDqmHelper.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., "bin", "Number of Events");
0041   nlogWgt_ = aDqmHelper.book1dHisto("nlogWgt", "Log10(n weights)", 100, 0., 3., "log_{10}(nWgts)", "Number of Events");
0042   wgtVal_ = aDqmHelper.book1dHisto("wgtVal", "weights", 100, -1.5, 3., "weight", "Number of Weights");
0043   bookTemplates(aDqmHelper,
0044                 leadLepPtTemp_,
0045                 "leadLepPt",
0046                 "leading lepton Pt",
0047                 leadLepPtNbin_,
0048                 0.,
0049                 leadLepPtRange_,
0050                 "Pt_{l} (GeV)",
0051                 "Number of Events");
0052   bookTemplates(aDqmHelper,
0053                 leadLepEtaTemp_,
0054                 "leadLepEta",
0055                 "leading lepton #eta",
0056                 rapidityNbin_,
0057                 -rapidityRange_,
0058                 rapidityRange_,
0059                 "#eta_{l}",
0060                 "Number of Events");
0061   bookTemplates(aDqmHelper,
0062                 jetMultTemp_,
0063                 "JetMultiplicity",
0064                 "Gen jet multiplicity",
0065                 nJetsNbin_,
0066                 0,
0067                 nJetsNbin_,
0068                 "n",
0069                 "Number of Events");
0070   bookTemplates(aDqmHelper,
0071                 leadJetPtTemp_,
0072                 "leadJetPt",
0073                 "leading Gen jet Pt",
0074                 jetPtNbin_,
0075                 0.,
0076                 jetPtRange_,
0077                 "Pt_{j} (GeV)",
0078                 "Number of Events");
0079   bookTemplates(aDqmHelper,
0080                 leadJetEtaTemp_,
0081                 "leadJetEta",
0082                 "leading Gen jet #eta",
0083                 rapidityNbin_,
0084                 -rapidityRange_,
0085                 rapidityRange_,
0086                 "#eta_{j}",
0087                 "Number of Events");
0088 
0089   return;
0090 }
0091 
0092 void GenWeightValidation::bookTemplates(DQMHelper& aDqmHelper,
0093                                         std::vector<MonitorElement*>& tmps,
0094                                         const std::string& name,
0095                                         const std::string& title,
0096                                         int nbin,
0097                                         float low,
0098                                         float high,
0099                                         const std::string& xtitle,
0100                                         const std::string& ytitle) {
0101   tmps.push_back(aDqmHelper.book1dHisto(name, title, nbin, low, high, xtitle, ytitle));
0102   tmps.push_back(aDqmHelper.book1dHisto(name + "FSRup", title + " FSR up", nbin, low, high, xtitle, ytitle));
0103   tmps.push_back(aDqmHelper.book1dHisto(name + "FSRdn", title + " FSR down", nbin, low, high, xtitle, ytitle));
0104   tmps.push_back(aDqmHelper.book1dHisto(
0105       name + "FSRup_ratio", "Ratio of " + title + " FSR up / Nominal", nbin, low, high, xtitle, ytitle));
0106   tmps.at(3)->setEfficiencyFlag();
0107   tmps.push_back(aDqmHelper.book1dHisto(
0108       name + "FSRdn_ratio", "Ratio of " + title + " FSR down / Nominal", nbin, low, high, xtitle, ytitle));
0109   tmps.at(4)->setEfficiencyFlag();
0110   tmps.push_back(aDqmHelper.book1dHisto(name + "ISRup", title + " ISR up", nbin, low, high, xtitle, ytitle));
0111   tmps.push_back(aDqmHelper.book1dHisto(name + "ISRdn", title + " ISR down", nbin, low, high, xtitle, ytitle));
0112   tmps.push_back(aDqmHelper.book1dHisto(
0113       name + "ISRup_ratio", "Ratio of " + title + " ISR up / Nominal", nbin, low, high, xtitle, ytitle));
0114   tmps.at(7)->setEfficiencyFlag();
0115   tmps.push_back(aDqmHelper.book1dHisto(
0116       name + "ISRdn_ratio", "Ratio of " + title + " ISR down / Nominal", nbin, low, high, xtitle, ytitle));
0117   tmps.at(8)->setEfficiencyFlag();
0118 }  // to get ratio plots correctly - need to modify PostProcessor_cff.py as well!
0119 
0120 void GenWeightValidation::dqmBeginRun(const edm::Run&, const edm::EventSetup&) {}
0121 
0122 void GenWeightValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0123   weights_ = wmanager_.weightsCollection(iEvent);
0124 
0125   unsigned weightsSize = weights_.at(idxGenEvtInfo_).size();
0126   if (weightsSize < 2)
0127     return;  // no baseline weight in GenEventInfo
0128 
0129   weight_ = weights_.at(idxGenEvtInfo_)[0] / std::abs(weights_.at(idxGenEvtInfo_)[0]);
0130   nEvt_->Fill(0.5, weight_);
0131   nlogWgt_->Fill(std::log10(weightsSize), weight_);
0132 
0133   for (unsigned idx = 0; idx < weightsSize; idx++)
0134     wgtVal_->Fill(weights_.at(idxGenEvtInfo_)[idx] / weights_.at(idxGenEvtInfo_)[1], weight_);
0135 
0136   if ((int)weightsSize <= idxMax_)
0137     return;  // no PS weights in GenEventInfo
0138 
0139   edm::Handle<reco::GenParticleCollection> ptcls;
0140   iEvent.getByToken(genParticleToken_, ptcls);
0141   edm::Handle<reco::GenJetCollection> genjets;
0142   iEvent.getByToken(genJetToken_, genjets);
0143 
0144   std::vector<reco::GenParticleRef> leptons;
0145 
0146   for (unsigned iptc = 0; iptc < ptcls->size(); iptc++) {
0147     reco::GenParticleRef ptc(ptcls, iptc);
0148     if (ptc->status() == 1 && (std::abs(ptc->pdgId()) == 11 || std::abs(ptc->pdgId()) == 13)) {
0149       if (ptc->pt() > leadLepPtCut_ && std::abs(ptc->eta()) < lepEtaCut_)
0150         leptons.push_back(ptc);
0151     }
0152   }
0153 
0154   std::sort(leptons.begin(), leptons.end(), HepMCValidationHelper::sortByPtRef<reco::GenParticleRef>);
0155 
0156   if (!leptons.empty()) {
0157     reco::GenParticleRef leadLep = leptons.at(0);
0158     fillTemplates(leadLepPtTemp_, leadLep->pt());
0159     fillTemplates(leadLepEtaTemp_, leadLep->eta());
0160   }
0161 
0162   std::vector<reco::GenJetRef> genjetVec;
0163 
0164   for (unsigned igj = 0; igj < genjets->size(); igj++) {
0165     reco::GenJetRef genjet(genjets, igj);
0166 
0167     if (genjet->pt() > jetPtCut_ && std::abs(genjet->eta()) < jetEtaCut_)
0168       genjetVec.push_back(genjet);
0169   }
0170 
0171   fillTemplates(jetMultTemp_, (float)genjetVec.size());
0172 
0173   if (!genjetVec.empty()) {
0174     std::sort(genjetVec.begin(), genjetVec.end(), HepMCValidationHelper::sortByPtRef<reco::GenJetRef>);
0175 
0176     auto leadJet = genjetVec.at(0);
0177     fillTemplates(leadJetPtTemp_, leadJet->pt());
0178     fillTemplates(leadJetEtaTemp_, leadJet->eta());
0179   }
0180 }  //analyze
0181 
0182 void GenWeightValidation::fillTemplates(std::vector<MonitorElement*>& tmps, float obs) {
0183   tmps.at(0)->Fill(obs, weight_);
0184   tmps.at(1)->Fill(obs, weights_.at(idxGenEvtInfo_)[idxFSRup_] / weights_.at(idxGenEvtInfo_)[1]);
0185   tmps.at(2)->Fill(obs, weights_.at(idxGenEvtInfo_)[idxFSRdown_] / weights_.at(idxGenEvtInfo_)[1]);
0186   tmps.at(5)->Fill(obs, weights_.at(idxGenEvtInfo_)[idxISRup_] / weights_.at(idxGenEvtInfo_)[1]);
0187   tmps.at(6)->Fill(obs, weights_.at(idxGenEvtInfo_)[idxISRdown_] / weights_.at(idxGenEvtInfo_)[1]);
0188 }