EventGenerator/plugins/LheWeightValidation.cc

0001 #include "Validation/EventGenerator/interface/LheWeightValidation.h"
0002
0003 #include "DataFormats/Math/interface/LorentzVector.h"
0004 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
0005
0006 using namespace edm;
0007
0008 LheWeightValidation::LheWeightValidation(const edm::ParameterSet& iPSet)
0009     : lheLabel_(iPSet.getParameter<edm::InputTag>("lheProduct")),
0010       genParticleToken_(consumes<reco::GenParticleCollection>(iPSet.getParameter<edm::InputTag>("genParticles"))),
0011       lheEvtToken_(consumes<LHEEventProduct>(lheLabel_)),
0012       lheRunToken_(consumes<LHERunInfoProduct, edm::InRun>(lheLabel_)),
0013       genJetToken_(consumes<reco::GenJetCollection>(iPSet.getParameter<edm::InputTag>("genJets"))),
0014       dumpLHEheader_(iPSet.getParameter<bool>("dumpLHEheader")),
0015       leadLepPtNbin_(iPSet.getParameter<int>("leadLepPtNbin")),
0016       rapidityNbin_(iPSet.getParameter<int>("rapidityNbin")),
0017       leadLepPtRange_(iPSet.getParameter<double>("leadLepPtRange")),
0018       leadLepPtCut_(iPSet.getParameter<double>("leadLepPtCut")),
0019       lepEtaCut_(iPSet.getParameter<double>("lepEtaCut")),
0020       rapidityRange_(iPSet.getParameter<double>("rapidityRange")),
0021       nJetsNbin_(iPSet.getParameter<int>("nJetsNbin")),
0022       jetPtNbin_(iPSet.getParameter<int>("jetPtNbin")),
0023       jetPtCut_(iPSet.getParameter<double>("jetPtCut")),
0024       jetEtaCut_(iPSet.getParameter<double>("jetEtaCut")),
0025       jetPtRange_(iPSet.getParameter<double>("jetPtRange")),
0026       nScaleVar_(iPSet.getParameter<int>("nScaleVar")),
0027       idxPdfStart_(iPSet.getParameter<int>("idxPdfStart")),
0028       idxPdfEnd_(iPSet.getParameter<int>("idxPdfEnd")),
0029       nPdfVar_(idxPdfEnd_ - idxPdfStart_ + 1) {}
0030
0031 void LheWeightValidation::bookHistograms(DQMStore::IBooker& iBook, edm::Run const& iRun, edm::EventSetup const&) {
0032   // check LHE product exists
0033   edm::Handle<LHERunInfoProduct> lheInfo;
0034   // getByToken throws an exception unless we're in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
0035   iRun.getByLabel(lheLabel_, lheInfo);
0036
0037   if (!lheInfo.isValid())
0038     return;
0039
0040   ///Setting the DQM top directories
0041   std::string folderName = "Generator/LHEWeight";
0042   DQMHelper aDqmHelper(&iBook);
0043   iBook.setCurrentFolder(folderName);
0044
0045   // Number of analyzed events
0046   nEvt_ = aDqmHelper.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., "bin", "Number of Events");
0047   nlogWgt_ = aDqmHelper.book1dHisto("nlogWgt", "Log10(n weights)", 100, 0., 5., "log_{10}(nWgts)", "Number of Events");
0048   wgtVal_ = aDqmHelper.book1dHisto("wgtVal", "weights", 100, -1.5, 3., "weight", "Number of Weights");
0049
0050   bookTemplates(aDqmHelper,
0051                 leadLepPtScaleVar_,
0052                 leadLepPtPdfVar_,
0053                 leadLepPtTemp_,
0054                 "leadLepPt",
0055                 "leading lepton Pt",
0056                 leadLepPtNbin_,
0057                 0.,
0058                 leadLepPtRange_,
0059                 "Pt_{l} (GeV)",
0060                 "Number of Events");
0061   bookTemplates(aDqmHelper,
0062                 leadLepEtaScaleVar_,
0063                 leadLepEtaPdfVar_,
0064                 leadLepEtaTemp_,
0065                 "leadLepEta",
0066                 "leading lepton #eta",
0067                 rapidityNbin_,
0068                 -rapidityRange_,
0069                 rapidityRange_,
0070                 "#eta_{l}",
0071                 "Number of Events");
0072   bookTemplates(aDqmHelper,
0073                 jetMultScaleVar_,
0074                 jetMultPdfVar_,
0075                 jetMultTemp_,
0076                 "JetMultiplicity",
0077                 "Gen jet multiplicity",
0078                 nJetsNbin_,
0079                 0,
0080                 nJetsNbin_,
0081                 "n",
0082                 "Number of Events");
0083   bookTemplates(aDqmHelper,
0084                 leadJetPtScaleVar_,
0085                 leadJetPtPdfVar_,
0086                 leadJetPtTemp_,
0087                 "leadJetPt",
0088                 "leading Gen jet Pt",
0089                 jetPtNbin_,
0090                 0.,
0091                 jetPtRange_,
0092                 "Pt_{j} (GeV)",
0093                 "Number of Events");
0094   bookTemplates(aDqmHelper,
0095                 leadJetEtaScaleVar_,
0096                 leadJetEtaPdfVar_,
0097                 leadJetEtaTemp_,
0098                 "leadJetEta",
0099                 "leading Gen jet #eta",
0100                 rapidityNbin_,
0101                 -rapidityRange_,
0102                 rapidityRange_,
0103                 "#eta_{j}",
0104                 "Number of Events");
0105
0106   return;
0107 }
0108
0109 void LheWeightValidation::bookTemplates(DQMHelper& aDqmHelper,
0110                                         std::vector<std::unique_ptr<TH1F>>& scaleVar,
0111                                         std::vector<std::unique_ptr<TH1F>>& pdfVar,
0112                                         std::vector<MonitorElement*>& tmps,
0113                                         const std::string& name,
0114                                         const std::string& title,
0115                                         int nbin,
0116                                         float low,
0117                                         float high,
0118                                         const std::string& xtitle,
0119                                         const std::string& ytitle) {
0120   tmps.push_back(aDqmHelper.book1dHisto(name, title, nbin, low, high, xtitle, ytitle));
0121   tmps.push_back(aDqmHelper.book1dHisto(name + "ScaleUp", title + " scale up", nbin, low, high, xtitle, ytitle));
0122   tmps.at(1)->getTH1()->Sumw2(false);
0123   tmps.push_back(aDqmHelper.book1dHisto(name + "ScaleDn", title + " scale down", nbin, low, high, xtitle, ytitle));
0124   tmps.at(2)->getTH1()->Sumw2(false);
0125   tmps.push_back(aDqmHelper.book1dHisto(
0126       name + "ScaleUp_ratio", "Ratio of " + title + " scale upper envelop / Nominal", nbin, low, high, xtitle, ytitle));
0127   tmps.at(3)->setEfficiencyFlag();
0128   tmps.push_back(aDqmHelper.book1dHisto(
0129       name + "ScaleDn_ratio", "Ratio of " + title + " scale lower envelop / Nominal", nbin, low, high, xtitle, ytitle));
0130   tmps.at(4)->setEfficiencyFlag();
0131   tmps.push_back(aDqmHelper.book1dHisto(name + "PdfUp", title + " PDF upper RMS", nbin, low, high, xtitle, ytitle));
0132   tmps.at(5)->getTH1()->Sumw2(false);
0133   tmps.push_back(aDqmHelper.book1dHisto(name + "PdfDn", title + " PDF lower RMS", nbin, low, high, xtitle, ytitle));
0134   tmps.at(6)->getTH1()->Sumw2(false);
0135   tmps.push_back(aDqmHelper.book1dHisto(
0136       name + "PdfUp_ratio", "Ratio of " + title + " PDF upper RMS / Nominal", nbin, low, high, xtitle, ytitle));
0137   tmps.at(7)->setEfficiencyFlag();
0138   tmps.push_back(aDqmHelper.book1dHisto(
0139       name + "PdfDn_ratio", "Ratio of " + title + " PDF lower RMS / Nominal", nbin, low, high, xtitle, ytitle));
0140   tmps.at(8)->setEfficiencyFlag();
0141
0142   for (int idx = 0; idx < nScaleVar_; idx++) {
0143     scaleVar.push_back(
0144         std::make_unique<TH1F>(std::string(name + "Scale" + std::to_string(idx)).c_str(),
0145                                std::string(";" + std::string(xtitle) + ";" + std::string(ytitle)).c_str(),
0146                                nbin,
0147                                low,
0148                                high));
0149     scaleVar.at(idx)->Sumw2();
0150   }
0151
0152   for (int idx = 0; idx < nPdfVar_; idx++) {
0153     pdfVar.push_back(std::make_unique<TH1F>(std::string(name + "Pdf" + std::to_string(idx)).c_str(),
0154                                             std::string(";" + std::string(xtitle) + ";" + std::string(ytitle)).c_str(),
0155                                             nbin,
0156                                             low,
0157                                             high));
0158     pdfVar.at(idx)->Sumw2();
0159   }
0160 }  // to get ratio plots correctly - need to modify PostProcessor_cff.py as well!
0161
0162 void LheWeightValidation::dqmBeginRun(const edm::Run&, const edm::EventSetup&) {}
0163
0164 void LheWeightValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0165   edm::Handle<LHEEventProduct> lheEvt;
0166
0167   if (!lheEvtToken_.isUninitialized())
0168     iEvent.getByToken(lheEvtToken_, lheEvt);
0169
0170   if (!lheEvt.isValid())
0171     return;  // do nothing if there is no LHE product
0172
0173   orgWgt_ = lheEvt->originalXWGTUP();
0174   weight_ = orgWgt_ / std::abs(orgWgt_);
0175   weights_ = lheEvt->weights();
0176
0177   nEvt_->Fill(0.5, weight_);
0178   nlogWgt_->Fill(std::log10(lheEvt->weights().size()), weight_);
0179
0180   for (unsigned idx = 0; idx < lheEvt->weights().size(); idx++)
0181     wgtVal_->Fill(weights_[idx].wgt / orgWgt_);
0182
0183   edm::Handle<reco::GenParticleCollection> ptcls;
0184   iEvent.getByToken(genParticleToken_, ptcls);
0185   edm::Handle<reco::GenJetCollection> genjets;
0186   iEvent.getByToken(genJetToken_, genjets);
0187
0188   std::vector<reco::GenParticleRef> leptons;
0189
0190   for (unsigned iptc = 0; iptc < ptcls->size(); iptc++) {
0191     reco::GenParticleRef ptc(ptcls, iptc);
0192     if (ptc->status() == 1 && (std::abs(ptc->pdgId()) == 11 || std::abs(ptc->pdgId()) == 13)) {
0193       if (ptc->pt() > leadLepPtCut_ && std::abs(ptc->eta()) < lepEtaCut_)
0194         leptons.push_back(ptc);
0195     }
0196   }
0197
0198   std::sort(leptons.begin(), leptons.end(), HepMCValidationHelper::sortByPtRef<reco::GenParticleRef>);
0199
0200   if (!leptons.empty()) {
0201     reco::GenParticleRef leadLep = leptons.at(0);
0202     fillTemplates(leadLepPtScaleVar_, leadLepPtPdfVar_, leadLepPtTemp_, leadLep->pt());
0203     fillTemplates(leadLepEtaScaleVar_, leadLepEtaPdfVar_, leadLepEtaTemp_, leadLep->eta());
0204   }
0205
0206   std::vector<reco::GenJetRef> genjetVec;
0207
0208   for (unsigned igj = 0; igj < genjets->size(); igj++) {
0209     reco::GenJetRef genjet(genjets, igj);
0210
0211     if (genjet->pt() > jetPtCut_ && std::abs(genjet->eta()) < jetEtaCut_)
0212       genjetVec.push_back(genjet);
0213   }
0214
0215   fillTemplates(jetMultScaleVar_, jetMultPdfVar_, jetMultTemp_, (float)genjetVec.size());
0216
0217   if (!genjetVec.empty()) {
0218     std::sort(genjetVec.begin(), genjetVec.end(), HepMCValidationHelper::sortByPtRef<reco::GenJetRef>);
0219
0220     auto leadJet = genjetVec.at(0);
0221     fillTemplates(leadJetPtScaleVar_, leadJetPtPdfVar_, leadJetPtTemp_, leadJet->pt());
0222     fillTemplates(leadJetEtaScaleVar_, leadJetEtaPdfVar_, leadJetEtaTemp_, leadJet->eta());
0223   }
0224 }  //analyze
0225
0226 void LheWeightValidation::fillTemplates(std::vector<std::unique_ptr<TH1F>>& scaleVar,
0227                                         std::vector<std::unique_ptr<TH1F>>& pdfVar,
0228                                         std::vector<MonitorElement*>& tmps,
0229                                         float obs) {
0230   tmps.at(0)->Fill(obs, weight_);
0231
0232   if (static_cast<int>(weights_.size()) >= nScaleVar_) {
0233     for (int iWgt = 0; iWgt < nScaleVar_; iWgt++)
0234       scaleVar.at(iWgt)->Fill(obs, weights_[iWgt].wgt / orgWgt_);
0235   }
0236
0237   if (static_cast<int>(weights_.size()) >= idxPdfEnd_) {
0238     for (int iWgt = 0; iWgt < nPdfVar_; iWgt++)
0239       pdfVar.at(iWgt)->Fill(obs, weights_[idxPdfStart_ + iWgt].wgt / orgWgt_);
0240   }
0241 }
0242
0243 void LheWeightValidation::dqmEndRun(const edm::Run& iRun, const edm::EventSetup&) {
0244   if (lheRunToken_.isUninitialized())
0245     return;
0246
0247   edm::Handle<LHERunInfoProduct> lheInfo;
0248   iRun.getByToken(lheRunToken_, lheInfo);
0249
0250   if (!lheInfo.isValid())
0251     return;
0252
0253   envelop(leadLepPtScaleVar_, leadLepPtTemp_);
0254   pdfRMS(leadLepPtPdfVar_, leadLepPtTemp_);
0255   envelop(leadLepEtaScaleVar_, leadLepEtaTemp_);
0256   pdfRMS(leadLepEtaPdfVar_, leadLepEtaTemp_);
0257   envelop(jetMultScaleVar_, jetMultTemp_);
0258   pdfRMS(jetMultPdfVar_, jetMultTemp_);
0259   envelop(leadJetPtScaleVar_, leadJetPtTemp_);
0260   pdfRMS(leadJetPtPdfVar_, leadJetPtTemp_);
0261   envelop(leadJetEtaScaleVar_, leadJetEtaTemp_);
0262   pdfRMS(leadJetEtaPdfVar_, leadJetEtaTemp_);
0263
0264   if (dumpLHEheader_) {
0265     for (auto it = lheInfo->headers_begin(); it != lheInfo->headers_end(); it++) {
0266       std::cout << "Header start" << std::endl;
0267       std::cout << "Tag: " << it->tag() << std::endl;
0268       for (const auto& l : it->lines()) {
0269         std::cout << l << std::endl;
0270       }
0271       std::cout << "Header end" << std::endl;
0272     }
0273   }
0274 }
0275
0276 void LheWeightValidation::envelop(const std::vector<std::unique_ptr<TH1F>>& var, std::vector<MonitorElement*>& tmps) {
0277   if (var.empty())
0278     return;
0279
0280   for (int b = 0; b < var.at(0)->GetNbinsX() + 2; b++) {
0281     float valU = var.at(0)->GetBinContent(b);
0282     float valD = valU;
0283
0284     if (valU == 0.)
0285       continue;
0286
0287     for (unsigned v = 1; v < var.size(); v++) {
0288       if (var.at(v)->GetEntries() == 0.)
0289         continue;
0290
0291       valU = std::max(valU, (float)var.at(v)->GetBinContent(b));
0292       valD = std::min(valD, (float)var.at(v)->GetBinContent(b));
0293     }
0294     tmps.at(1)->setBinContent(b, valU);
0295     tmps.at(2)->setBinContent(b, valD);
0296   }
0297
0298   tmps.at(1)->setEntries(var.at(0)->GetEntries());
0299   tmps.at(2)->setEntries(var.at(0)->GetEntries());
0300   tmps.at(1)->getTH1()->Sumw2(true);
0301   tmps.at(2)->getTH1()->Sumw2(true);
0302
0303   return;
0304 }
0305
0306 void LheWeightValidation::pdfRMS(const std::vector<std::unique_ptr<TH1F>>& var, std::vector<MonitorElement*>& tmps) {
0307   if (var.empty())
0308     return;
0309
0310   float denom = var.size();
0311   for (int b = 0; b < tmps.at(0)->getNbinsX() + 2; b++) {
0312     float valNom = tmps.at(0)->getBinContent(b);
0313     float rmsSq = 0.;
0314     if (valNom == 0.)
0315       continue;
0316
0317     for (unsigned v = 0; v < var.size(); v++) {
0318       if (var.at(v)->GetEntries() == 0.)
0319         continue;
0320
0321       float dev = (float)var.at(v)->GetBinContent(b) - valNom;
0322       rmsSq += dev * dev;
0323     }
0324
0325     float rms = std::sqrt(rmsSq / denom);
0326     float rmsup = valNom + rms;
0327     float rmsdn = valNom - rms;
0328     tmps.at(5)->setBinContent(b, rmsup);
0329     tmps.at(6)->setBinContent(b, rmsdn);
0330   }
0331
0332   tmps.at(5)->setEntries(tmps.at(0)->getTH1F()->GetEntries());
0333   tmps.at(6)->setEntries(tmps.at(0)->getTH1F()->GetEntries());
0334   tmps.at(5)->getTH1()->Sumw2(true);
0335   tmps.at(6)->getTH1()->Sumw2(true);
0336
0337   return;
0338 }