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
0035 std::string folderName = "Generator/GenWeight";
0036 DQMHelper aDqmHelper(&iBook);
0037 iBook.setCurrentFolder(folderName);
0038
0039
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 }
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;
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;
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 }
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 }