Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-02-05 03:15:09 

0001 #include "FWCore/Framework/interface/global/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/Run.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0006 #include "DataFormats/NanoAOD/interface/MergeableCounterTable.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009 #include "FWCore/Utilities/interface/transform.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0011 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0012 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0013 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "boost/algorithm/string.hpp"
0016 
0017 #include <array>
0018 #include <memory>
0019 
0020 #include <vector>
0021 #include <unordered_map>
0022 #include <iostream>
0023 #include <regex>
0024 
0025 namespace {
0026   ///  ---- Cache object for running sums of weights ----
0027   struct Counter {
0028     Counter() : num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumRwgt(), sumNamed(), sumPS() {}
0029 
0030     // the counters
0031     long long num;
0032     long double sumw;
0033     long double sumw2;
0034     std::vector<long double> sumPDF, sumScale, sumRwgt, sumNamed, sumPS;
0035 
0036     void clear() {
0037       num = 0;
0038       sumw = 0;
0039       sumw2 = 0;
0040       sumPDF.clear();
0041       sumScale.clear();
0042       sumRwgt.clear();
0043       sumNamed.clear(), sumPS.clear();
0044     }
0045 
0046     // inc the counters
0047     void incGenOnly(double w) {
0048       num++;
0049       sumw += w;
0050       sumw2 += (w * w);
0051     }
0052 
0053     void incPSOnly(double w0, const std::vector<double>& wPS) {
0054       if (!wPS.empty()) {
0055         if (sumPS.empty())
0056           sumPS.resize(wPS.size(), 0);
0057         for (unsigned int i = 0, n = wPS.size(); i < n; ++i)
0058           sumPS[i] += (w0 * wPS[i]);
0059       }
0060     }
0061 
0062     void incLHE(double w0,
0063                 const std::vector<double>& wScale,
0064                 const std::vector<double>& wPDF,
0065                 const std::vector<double>& wRwgt,
0066                 const std::vector<double>& wNamed,
0067                 const std::vector<double>& wPS) {
0068       // add up weights
0069       incGenOnly(w0);
0070       // then add up variations
0071       if (!wScale.empty()) {
0072         if (sumScale.empty())
0073           sumScale.resize(wScale.size(), 0);
0074         for (unsigned int i = 0, n = wScale.size(); i < n; ++i)
0075           sumScale[i] += (w0 * wScale[i]);
0076       }
0077       if (!wPDF.empty()) {
0078         if (sumPDF.empty())
0079           sumPDF.resize(wPDF.size(), 0);
0080         for (unsigned int i = 0, n = wPDF.size(); i < n; ++i)
0081           sumPDF[i] += (w0 * wPDF[i]);
0082       }
0083       if (!wRwgt.empty()) {
0084         if (sumRwgt.empty())
0085           sumRwgt.resize(wRwgt.size(), 0);
0086         for (unsigned int i = 0, n = wRwgt.size(); i < n; ++i)
0087           sumRwgt[i] += (w0 * wRwgt[i]);
0088       }
0089       if (!wNamed.empty()) {
0090         if (sumNamed.empty())
0091           sumNamed.resize(wNamed.size(), 0);
0092         for (unsigned int i = 0, n = wNamed.size(); i < n; ++i)
0093           sumNamed[i] += (w0 * wNamed[i]);
0094       }
0095       incPSOnly(w0, wPS);
0096     }
0097 
0098     void merge(const Counter& other) {
0099       num += other.num;
0100       sumw += other.sumw;
0101       sumw2 += other.sumw2;
0102       if (sumScale.empty() && !other.sumScale.empty())
0103         sumScale.resize(other.sumScale.size(), 0);
0104       if (sumPDF.empty() && !other.sumPDF.empty())
0105         sumPDF.resize(other.sumPDF.size(), 0);
0106       if (sumRwgt.empty() && !other.sumRwgt.empty())
0107         sumRwgt.resize(other.sumRwgt.size(), 0);
0108       if (sumNamed.empty() && !other.sumNamed.empty())
0109         sumNamed.resize(other.sumNamed.size(), 0);
0110       if (sumPS.empty() && !other.sumPS.empty())
0111         sumPS.resize(other.sumPS.size(), 0);
0112       if (!other.sumScale.empty())
0113         for (unsigned int i = 0, n = sumScale.size(); i < n; ++i)
0114           sumScale[i] += other.sumScale[i];
0115       if (!other.sumPDF.empty())
0116         for (unsigned int i = 0, n = sumPDF.size(); i < n; ++i)
0117           sumPDF[i] += other.sumPDF[i];
0118       if (!other.sumRwgt.empty())
0119         for (unsigned int i = 0, n = sumRwgt.size(); i < n; ++i)
0120           sumRwgt[i] += other.sumRwgt[i];
0121       if (!other.sumNamed.empty())
0122         for (unsigned int i = 0, n = sumNamed.size(); i < n; ++i)
0123           sumNamed[i] += other.sumNamed[i];
0124       if (!other.sumPS.empty())
0125         for (unsigned int i = 0, n = sumPS.size(); i < n; ++i)
0126           sumPS[i] += other.sumPS[i];
0127     }
0128   };
0129 
0130   struct CounterMap {
0131     std::map<std::string, Counter> countermap;
0132     Counter* active_el = nullptr;
0133     std::string active_label = "";
0134     void merge(const CounterMap& other) {
0135       for (const auto& y : other.countermap)
0136         countermap[y.first].merge(y.second);
0137       active_el = nullptr;
0138     }
0139     void clear() {
0140       for (auto x : countermap)
0141         x.second.clear();
0142       active_el = nullptr;
0143       active_label = "";
0144     }
0145     void setLabel(std::string label) {
0146       active_el = &(countermap[label]);
0147       active_label = label;
0148     }
0149     void checkLabelSet() {
0150       if (!active_el)
0151         throw cms::Exception("LogicError", "Called CounterMap::get() before setting the active label\n");
0152     }
0153     Counter* get() {
0154       checkLabelSet();
0155       return active_el;
0156     }
0157     std::string& getLabel() {
0158       checkLabelSet();
0159       return active_label;
0160     }
0161   };
0162 
0163   ///  ---- RunCache object for dynamic choice of LHE IDs ----
0164   struct DynamicWeightChoice {
0165     // choice of LHE weights
0166     // ---- scale ----
0167     std::vector<std::string> scaleWeightIDs;
0168     std::string scaleWeightsDoc;
0169     // ---- pdf ----
0170     std::vector<std::string> pdfWeightIDs;
0171     std::string pdfWeightsDoc;
0172     // ---- rwgt ----
0173     std::vector<std::string> rwgtIDs;
0174     std::string rwgtWeightDoc;
0175   };
0176 
0177   constexpr std::array<unsigned int, 4> defPSWeightIDs = {{6, 7, 8, 9}};
0178   constexpr std::array<unsigned int, 4> defPSWeightIDs_alt = {{27, 5, 26, 4}};
0179 
0180   struct DynamicWeightChoiceGenInfo {
0181     // choice of LHE weights
0182     // ---- scale ----
0183     std::vector<unsigned int> scaleWeightIDs;
0184     std::string scaleWeightsDoc;
0185     // ---- pdf ----
0186     std::vector<unsigned int> pdfWeightIDs;
0187     std::string pdfWeightsDoc;
0188     // ---- ps ----
0189     bool matchPS_alt = false;
0190     std::vector<unsigned int> psWeightIDs;
0191     unsigned int psBaselineID = 1;
0192     std::string psWeightsDoc;
0193 
0194     void setMissingWeight(int idx) { psWeightIDs[idx] = (matchPS_alt) ? defPSWeightIDs_alt[idx] : defPSWeightIDs[idx]; }
0195 
0196     bool empty() const { return scaleWeightIDs.empty() && pdfWeightIDs.empty() && psWeightIDs.empty(); }
0197   };
0198 
0199   struct LumiCacheInfoHolder {
0200     CounterMap countermap;
0201     void clear() { countermap.clear(); }
0202   };
0203 
0204   float stof_fortrancomp(const std::string& str) {
0205     std::string::size_type match = str.find('d');
0206     if (match != std::string::npos) {
0207       std::string pre = str.substr(0, match);
0208       std::string post = str.substr(match + 1);
0209       return std::stof(pre) * std::pow(10.0f, std::stof(post));
0210     } else {
0211       return std::stof(str);
0212     }
0213   }
0214   ///  -------------- temporary objects --------------
0215   struct ScaleVarWeight {
0216     std::string wid, label;
0217     std::pair<float, float> scales;
0218     ScaleVarWeight(const std::string& id, const std::string& text, const std::string& muR, const std::string& muF)
0219         : wid(id), label(text), scales(stof_fortrancomp(muR), stof_fortrancomp(muF)) {}
0220     bool operator<(const ScaleVarWeight& other) {
0221       return (scales == other.scales ? wid < other.wid : scales < other.scales);
0222     }
0223   };
0224   struct PDFSetWeights {
0225     std::vector<std::string> wids;
0226     std::pair<unsigned int, unsigned int> lhaIDs;
0227     PDFSetWeights(const std::string& wid, unsigned int lhaID) : wids(1, wid), lhaIDs(lhaID, lhaID) {}
0228     bool operator<(const PDFSetWeights& other) const { return lhaIDs < other.lhaIDs; }
0229     void add(const std::string& wid, unsigned int lhaID) {
0230       wids.push_back(wid);
0231       lhaIDs.second = lhaID;
0232     }
0233     bool maybe_add(const std::string& wid, unsigned int lhaID) {
0234       if (lhaID == lhaIDs.second + 1) {
0235         lhaIDs.second++;
0236         wids.push_back(wid);
0237         return true;
0238       } else {
0239         return false;
0240       }
0241     }
0242   };
0243 }  // namespace
0244 
0245 class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<LumiCacheInfoHolder>,
0246                                                                edm::RunCache<DynamicWeightChoice>,
0247                                                                edm::LuminosityBlockCache<DynamicWeightChoiceGenInfo>,
0248                                                                edm::RunSummaryCache<CounterMap>,
0249                                                                edm::EndRunProducer> {
0250 public:
0251   GenWeightsTableProducer(edm::ParameterSet const& params)
0252       : genTag_(consumes<GenEventInfoProduct>(params.getParameter<edm::InputTag>("genEvent"))),
0253         lheLabel_(params.getParameter<std::vector<edm::InputTag>>("lheInfo")),
0254         lheTag_(edm::vector_transform(lheLabel_,
0255                                       [this](const edm::InputTag& tag) { return mayConsume<LHEEventProduct>(tag); })),
0256         lheRunTag_(edm::vector_transform(
0257             lheLabel_, [this](const edm::InputTag& tag) { return mayConsume<LHERunInfoProduct, edm::InRun>(tag); })),
0258         genLumiInfoHeadTag_(
0259             mayConsume<GenLumiInfoHeader, edm::InLumi>(params.getParameter<edm::InputTag>("genLumiInfoHeader"))),
0260         namedWeightIDs_(params.getParameter<std::vector<std::string>>("namedWeightIDs")),
0261         namedWeightLabels_(params.getParameter<std::vector<std::string>>("namedWeightLabels")),
0262         lheWeightPrecision_(params.getParameter<int32_t>("lheWeightPrecision")),
0263         maxPdfWeights_(params.getParameter<uint32_t>("maxPdfWeights")),
0264         keepAllPSWeights_(params.getParameter<bool>("keepAllPSWeights")),
0265         allowedNumScaleWeights_(params.getParameter<std::vector<uint32_t>>("allowedNumScaleWeights")),
0266         debug_(params.getUntrackedParameter<bool>("debug", false)),
0267         debugRun_(debug_.load()),
0268         hasIssuedWarning_(false),
0269         psWeightWarning_(false) {
0270     produces<nanoaod::FlatTable>();
0271     produces<std::string>("genModel");
0272     produces<nanoaod::FlatTable>("LHEScale");
0273     produces<nanoaod::FlatTable>("LHEPdf");
0274     produces<nanoaod::FlatTable>("LHEReweighting");
0275     produces<nanoaod::FlatTable>("LHENamed");
0276     produces<nanoaod::FlatTable>("PS");
0277     produces<nanoaod::MergeableCounterTable, edm::Transition::EndRun>();
0278     if (namedWeightIDs_.size() != namedWeightLabels_.size()) {
0279       throw cms::Exception("Configuration", "Size mismatch between namedWeightIDs & namedWeightLabels");
0280     }
0281     for (const edm::ParameterSet& pdfps : params.getParameter<std::vector<edm::ParameterSet>>("preferredPDFs")) {
0282       const std::string& name = pdfps.getParameter<std::string>("name");
0283       uint32_t lhaid = pdfps.getParameter<uint32_t>("lhaid");
0284       preferredPDFLHAIDs_.push_back(lhaid);
0285       lhaNameToID_[name] = lhaid;
0286       lhaNameToID_[name + ".LHgrid"] = lhaid;
0287     }
0288   }
0289 
0290   ~GenWeightsTableProducer() override {}
0291 
0292   void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
0293     // get my counter for weights
0294     Counter* counter = streamCache(id)->countermap.get();
0295 
0296     // generator information (always available)
0297     edm::Handle<GenEventInfoProduct> genInfo;
0298     iEvent.getByToken(genTag_, genInfo);
0299     double weight = genInfo->weight();
0300 
0301     // table for gen info, always available
0302     auto out = std::make_unique<nanoaod::FlatTable>(1, "genWeight", true);
0303     out->setDoc("generator weight");
0304     out->addColumnValue<float>("", weight, "generator weight");
0305     iEvent.put(std::move(out));
0306 
0307     std::string model_label = streamCache(id)->countermap.getLabel();
0308     auto outM = std::make_unique<std::string>((!model_label.empty()) ? std::string("GenModel_") + model_label : "");
0309     iEvent.put(std::move(outM), "genModel");
0310     bool getLHEweightsFromGenInfo = !model_label.empty();
0311 
0312     // tables for LHE weights, may not be filled
0313     std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab;
0314     std::unique_ptr<nanoaod::FlatTable> genPSTab;
0315 
0316     edm::Handle<LHEEventProduct> lheInfo;
0317     for (const auto& lheTag : lheTag_) {
0318       iEvent.getByToken(lheTag, lheInfo);
0319       if (lheInfo.isValid()) {
0320         break;
0321       }
0322     }
0323 
0324     const auto genWeightChoice = luminosityBlockCache(iEvent.getLuminosityBlock().index());
0325     if (lheInfo.isValid()) {
0326       if (getLHEweightsFromGenInfo && !hasIssuedWarning_.exchange(true))
0327         edm::LogWarning("LHETablesProducer")
0328             << "Found both a LHEEventProduct and a GenLumiInfoHeader: will only save weights from LHEEventProduct.\n";
0329       // get the dynamic choice of weights
0330       const DynamicWeightChoice* weightChoice = runCache(iEvent.getRun().index());
0331       // go fill tables
0332       fillLHEWeightTables(counter,
0333                           weightChoice,
0334                           genWeightChoice,
0335                           weight,
0336                           *lheInfo,
0337                           *genInfo,
0338                           lheScaleTab,
0339                           lhePdfTab,
0340                           lheRwgtTab,
0341                           lheNamedTab,
0342                           genPSTab);
0343     } else if (getLHEweightsFromGenInfo) {
0344       fillLHEPdfWeightTablesFromGenInfo(
0345           counter, genWeightChoice, weight, *genInfo, lheScaleTab, lhePdfTab, lheNamedTab, genPSTab);
0346       lheRwgtTab = std::make_unique<nanoaod::FlatTable>(1, "LHEReweightingWeights", true);
0347       //lheNamedTab.reset(new nanoaod::FlatTable(1, "LHENamedWeights", true));
0348       //genPSTab.reset(new nanoaod::FlatTable(1, "PSWeight", true));
0349     } else {
0350       // Still try to add the PS weights
0351       fillOnlyPSWeightTable(counter, genWeightChoice, weight, *genInfo, genPSTab);
0352       // make dummy values
0353       lheScaleTab = std::make_unique<nanoaod::FlatTable>(1, "LHEScaleWeights", true);
0354       lhePdfTab = std::make_unique<nanoaod::FlatTable>(1, "LHEPdfWeights", true);
0355       lheRwgtTab = std::make_unique<nanoaod::FlatTable>(1, "LHEReweightingWeights", true);
0356       lheNamedTab = std::make_unique<nanoaod::FlatTable>(1, "LHENamedWeights", true);
0357       if (!hasIssuedWarning_.exchange(true)) {
0358         edm::LogWarning("LHETablesProducer") << "No LHEEventProduct, so there will be no LHE Tables\n";
0359       }
0360     }
0361 
0362     iEvent.put(std::move(lheScaleTab), "LHEScale");
0363     iEvent.put(std::move(lhePdfTab), "LHEPdf");
0364     iEvent.put(std::move(lheRwgtTab), "LHEReweighting");
0365     iEvent.put(std::move(lheNamedTab), "LHENamed");
0366     iEvent.put(std::move(genPSTab), "PS");
0367   }
0368 
0369   void fillLHEWeightTables(Counter* counter,
0370                            const DynamicWeightChoice* weightChoice,
0371                            const DynamicWeightChoiceGenInfo* genWeightChoice,
0372                            double genWeight,
0373                            const LHEEventProduct& lheProd,
0374                            const GenEventInfoProduct& genProd,
0375                            std::unique_ptr<nanoaod::FlatTable>& outScale,
0376                            std::unique_ptr<nanoaod::FlatTable>& outPdf,
0377                            std::unique_ptr<nanoaod::FlatTable>& outRwgt,
0378                            std::unique_ptr<nanoaod::FlatTable>& outNamed,
0379                            std::unique_ptr<nanoaod::FlatTable>& outPS) const {
0380     bool lheDebug = debug_.exchange(
0381         false);  // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
0382 
0383     const std::vector<std::string>& scaleWeightIDs = weightChoice->scaleWeightIDs;
0384     const std::vector<std::string>& pdfWeightIDs = weightChoice->pdfWeightIDs;
0385     const std::vector<std::string>& rwgtWeightIDs = weightChoice->rwgtIDs;
0386 
0387     double w0 = lheProd.originalXWGTUP();
0388 
0389     std::vector<double> wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wRwgt(rwgtWeightIDs.size(), 1),
0390         wNamed(namedWeightIDs_.size(), 1);
0391     for (auto& weight : lheProd.weights()) {
0392       if (lheDebug)
0393         printf("Weight  %+9.5f   rel %+9.5f   for id %s\n", weight.wgt, weight.wgt / w0, weight.id.c_str());
0394       // now we do it slowly, can be optimized
0395       auto mScale = std::find(scaleWeightIDs.begin(), scaleWeightIDs.end(), weight.id);
0396       if (mScale != scaleWeightIDs.end())
0397         wScale[mScale - scaleWeightIDs.begin()] = weight.wgt / w0;
0398 
0399       auto mPDF = std::find(pdfWeightIDs.begin(), pdfWeightIDs.end(), weight.id);
0400       if (mPDF != pdfWeightIDs.end())
0401         wPDF[mPDF - pdfWeightIDs.begin()] = weight.wgt / w0;
0402 
0403       auto mRwgt = std::find(rwgtWeightIDs.begin(), rwgtWeightIDs.end(), weight.id);
0404       if (mRwgt != rwgtWeightIDs.end())
0405         wRwgt[mRwgt - rwgtWeightIDs.begin()] = weight.wgt / w0;
0406 
0407       auto mNamed = std::find(namedWeightIDs_.begin(), namedWeightIDs_.end(), weight.id);
0408       if (mNamed != namedWeightIDs_.end())
0409         wNamed[mNamed - namedWeightIDs_.begin()] = weight.wgt / w0;
0410     }
0411 
0412     std::vector<double> wPS;
0413     std::string psWeightDocStr;
0414     setPSWeightInfo(genProd.weights(), genWeightChoice, wPS, psWeightDocStr);
0415 
0416     outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
0417     outPS->addColumn<float>("", wPS, psWeightDocStr, lheWeightPrecision_);
0418 
0419     outScale = std::make_unique<nanoaod::FlatTable>(wScale.size(), "LHEScaleWeight", false);
0420     outScale->addColumn<float>("", wScale, weightChoice->scaleWeightsDoc, lheWeightPrecision_);
0421 
0422     outPdf = std::make_unique<nanoaod::FlatTable>(wPDF.size(), "LHEPdfWeight", false);
0423     outPdf->addColumn<float>("", wPDF, weightChoice->pdfWeightsDoc, lheWeightPrecision_);
0424 
0425     outRwgt = std::make_unique<nanoaod::FlatTable>(wRwgt.size(), "LHEReweightingWeight", false);
0426     outRwgt->addColumn<float>("", wRwgt, weightChoice->rwgtWeightDoc, lheWeightPrecision_);
0427 
0428     outNamed = std::make_unique<nanoaod::FlatTable>(1, "LHEWeight", true);
0429     outNamed->addColumnValue<float>("originalXWGTUP", lheProd.originalXWGTUP(), "Nominal event weight in the LHE file");
0430     for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) {
0431       outNamed->addColumnValue<float>(namedWeightLabels_[i],
0432                                       wNamed[i],
0433                                       "LHE weight for id " + namedWeightIDs_[i] + ", relative to nominal",
0434                                       lheWeightPrecision_);
0435     }
0436 
0437     counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS);
0438   }
0439 
0440   void fillLHEPdfWeightTablesFromGenInfo(Counter* counter,
0441                                          const DynamicWeightChoiceGenInfo* weightChoice,
0442                                          double genWeight,
0443                                          const GenEventInfoProduct& genProd,
0444                                          std::unique_ptr<nanoaod::FlatTable>& outScale,
0445                                          std::unique_ptr<nanoaod::FlatTable>& outPdf,
0446                                          std::unique_ptr<nanoaod::FlatTable>& outNamed,
0447                                          std::unique_ptr<nanoaod::FlatTable>& outPS) const {
0448     const std::vector<unsigned int>& scaleWeightIDs = weightChoice->scaleWeightIDs;
0449     const std::vector<unsigned int>& pdfWeightIDs = weightChoice->pdfWeightIDs;
0450 
0451     auto weights = genProd.weights();
0452     double w0 = (weights.size() > 1) ? weights.at(1) : 1.;
0453     double originalXWGTUP = (weights.size() > 1) ? weights.at(1) : 1.;
0454 
0455     std::vector<double> wScale, wPDF, wPS;
0456     for (auto id : scaleWeightIDs)
0457       wScale.push_back(weights.at(id) / w0);
0458     for (auto id : pdfWeightIDs) {
0459       wPDF.push_back(weights.at(id) / w0);
0460     }
0461 
0462     std::string psWeightsDocStr;
0463     setPSWeightInfo(genProd.weights(), weightChoice, wPS, psWeightsDocStr);
0464 
0465     outScale = std::make_unique<nanoaod::FlatTable>(wScale.size(), "LHEScaleWeight", false);
0466     outScale->addColumn<float>("", wScale, weightChoice->scaleWeightsDoc, lheWeightPrecision_);
0467 
0468     outPdf = std::make_unique<nanoaod::FlatTable>(wPDF.size(), "LHEPdfWeight", false);
0469     outPdf->addColumn<float>("", wPDF, weightChoice->pdfWeightsDoc, lheWeightPrecision_);
0470 
0471     outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
0472     outPS->addColumn<float>("", wPS, psWeightsDocStr, lheWeightPrecision_);
0473 
0474     outNamed = std::make_unique<nanoaod::FlatTable>(1, "LHEWeight", true);
0475     outNamed->addColumnValue<float>("originalXWGTUP", originalXWGTUP, "Nominal event weight in the LHE file");
0476     /*for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) {
0477       outNamed->addColumnValue<float>(namedWeightLabels_[i], wNamed[i], "LHE weight for id "+namedWeightIDs_[i]+", relative to nominal", lheWeightPrecision_);
0478       }*/
0479 
0480     counter->incLHE(genWeight, wScale, wPDF, std::vector<double>(), std::vector<double>(), wPS);
0481   }
0482 
0483   void fillOnlyPSWeightTable(Counter* counter,
0484                              const DynamicWeightChoiceGenInfo* genWeightChoice,
0485                              double genWeight,
0486                              const GenEventInfoProduct& genProd,
0487                              std::unique_ptr<nanoaod::FlatTable>& outPS) const {
0488     std::vector<double> wPS;
0489     std::string psWeightDocStr;
0490     setPSWeightInfo(genProd.weights(), genWeightChoice, wPS, psWeightDocStr);
0491     outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
0492     outPS->addColumn<float>("", wPS, psWeightDocStr, lheWeightPrecision_);
0493 
0494     counter->incGenOnly(genWeight);
0495     counter->incPSOnly(genWeight, wPS);
0496   }
0497 
0498   void setPSWeightInfo(const std::vector<double>& genWeights,
0499                        const DynamicWeightChoiceGenInfo* genWeightChoice,
0500                        std::vector<double>& wPS,
0501                        std::string& psWeightDocStr) const {
0502     wPS.clear();
0503     // isRegularPSSet = keeping all weights and the weights are a usual size, ie
0504     //                  all weights are PS weights (don't use header incase missing names)
0505     bool isRegularPSSet = keepAllPSWeights_ && (genWeights.size() == 14 || genWeights.size() == 46);
0506     if (!genWeightChoice->psWeightIDs.empty() && !isRegularPSSet) {
0507       psWeightDocStr = genWeightChoice->psWeightsDoc;
0508       double psNom = genWeights.at(genWeightChoice->psBaselineID);
0509       for (auto wgtidx : genWeightChoice->psWeightIDs) {
0510         wPS.push_back(genWeights.at(wgtidx) / psNom);
0511       }
0512     } else {
0513       int vectorSize =
0514           keepAllPSWeights_ ? (genWeights.size() - 2) : ((genWeights.size() == 14 || genWeights.size() == 46) ? 4 : 1);
0515 
0516       if (vectorSize > 1) {
0517         double nominal = genWeights.at(1);  // Called 'Baseline' in GenLumiInfoHeader
0518         if (keepAllPSWeights_) {
0519           for (int i = 0; i < vectorSize; i++) {
0520             wPS.push_back(genWeights.at(i + 2) / nominal);
0521           }
0522           psWeightDocStr = ((int)genWeightChoice->psWeightIDs.size() == vectorSize)
0523                                ? genWeightChoice->psWeightsDoc
0524                                : "All PS weights (w_var / w_nominal)";
0525         } else {
0526           if (!psWeightWarning_.exchange(true))
0527             edm::LogWarning("LHETablesProducer")
0528                 << "GenLumiInfoHeader not found: Central PartonShower weights will fill with the 6-10th entries \n"
0529                 << "    This may incorrect for some mcs (madgraph 2.6.1 with its `isr:murfact=0.5` have a differnt "
0530                    "order )";
0531           for (std::size_t i = 6; i < 10; i++) {
0532             wPS.push_back(genWeights.at(i) / nominal);
0533           }
0534           psWeightDocStr =
0535               "PS weights (w_var / w_nominal); default (maybe incorrect) order is: "
0536               "[0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2; [2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;";
0537         }
0538       } else {
0539         wPS.push_back(1.0);
0540         psWeightDocStr = "dummy PS weight (1.0) ";
0541       }
0542     }
0543   }
0544 
0545   // create an empty counter
0546   std::shared_ptr<DynamicWeightChoice> globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override {
0547     edm::Handle<LHERunInfoProduct> lheInfo;
0548 
0549     bool lheDebug = debugRun_.exchange(
0550         false);  // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
0551     auto weightChoice = std::make_shared<DynamicWeightChoice>();
0552 
0553     // getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
0554     //if (iRun.getByToken(lheRunTag_, lheInfo)) {
0555     for (const auto& lheLabel : lheLabel_) {
0556       iRun.getByLabel(lheLabel, lheInfo);
0557       if (lheInfo.isValid()) {
0558         break;
0559       }
0560     }
0561     if (lheInfo.isValid()) {
0562       std::vector<ScaleVarWeight> scaleVariationIDs;
0563       std::vector<PDFSetWeights> pdfSetWeightIDs;
0564       std::vector<std::string> lheReweighingIDs;
0565       bool preScaleVariationGroup = true;
0566 
0567       std::regex weightgroupmg26x("<weightgroup\\s+(?:name|type)=\"(.*)\"\\s+combine=\"(.*)\"\\s*>");
0568       std::regex weightgroup("<weightgroup\\s+combine=\"(.*)\"\\s+(?:name|type)=\"(.*)\"\\s*>");
0569       std::regex weightgroupRwgt("<weightgroup\\s+(?:name|type)=\"(.*)\"\\s*>");
0570       std::regex endweightgroup("</weightgroup>");
0571       std::regex scalewmg26x(
0572           "<weight\\s+(?:.*\\s+)?id=\"(\\d+)\"\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:[mM][uU][rR]|renscfact)=\"("
0573           "\\S+)\"\\s+(?:[mM][uU][Ff]|facscfact)=\"(\\S+)\")(\\s+.*)?</weight>");
0574       std::regex scalewmg26xNew(
0575           "<weight\\s*((?:[mM][uU][fF]|facscfact)=\"(\\S+)\"\\s+(?:[mM][uU][Rr]|renscfact)=\"(\\S+)\").+id=\"(\\d+)\"(."
0576           "*)?</weight>");
0577 
0578       //<weight MUF="1.0" MUR="2.0" PDF="306000" id="1006"> MUR=2.0  </weight>
0579       std::regex scalew(
0580           "<weight\\s+(?:.*\\s+)?id=\"(\\d+|\\d+-NNLOPS)\">\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:mu[rR]|renscfact)"
0581           "=(\\S+)\\s+(?:mu[Ff]|facscfact)=(\\S+)(\\s+.*)?)</weight>");
0582       std::regex pdfw(
0583           "<weight\\s+id=\"(\\d+)\">\\s*(?:PDF set|lhapdf|PDF|pdfset)\\s*=\\s*(\\d+)\\s*(?:\\s.*)?</weight>");
0584       std::regex pdfwOld("<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*Member \\s*(\\d+)\\s*(?:.*)</weight>");
0585       std::regex pdfwmg26x(
0586           "<weight\\s+id=\"(\\d+)\"\\s*MUR=\"(?:\\S+)\"\\s*MUF=\"(?:\\S+)\"\\s*(?:PDF "
0587           "set|lhapdf|PDF|pdfset)\\s*=\\s*\"(\\d+)\"\\s*>\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</"
0588           "weight>");
0589       //<weightgroup combine="symmhessian+as" name="NNPDF31_nnlo_as_0118_mc_hessian_pdfas">
0590 
0591       //<weight MUF="1.0" MUR="1.0" PDF="325300" id="1048"> PDF=325300 MemberID=0 </weight>
0592       std::regex pdfwmg26xNew(
0593           "<weight\\s+MUF=\"(?:\\S+)\"\\s*MUR=\"(?:\\S+)\"\\s*PDF=\"(?:\\S+)\"\\s*id=\"(\\S+)\"\\s*>"
0594           "\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</"
0595           "weight>");
0596 
0597       std::regex mgVerRegex(R"(VERSION\s+(\d+)\.(\d+)\.(\d+))");
0598       bool isMGVer2x = false;
0599 
0600       std::regex rwgt("<weight\\s+id=\"(.+)\">(.+)?(</weight>)?");
0601       std::smatch groups;
0602       for (auto iter = lheInfo->headers_begin(), end = lheInfo->headers_end(); iter != end; ++iter) {
0603         if (iter->tag() == "MG5ProcCard") {
0604           for (const auto& line : iter->lines()) {
0605             if (std::regex_search(line, groups, mgVerRegex)) {
0606               isMGVer2x = (groups[1].str() == "2");
0607               break;
0608             }
0609           }
0610         }
0611         if (iter->tag() != "initrwgt") {
0612           if (lheDebug)
0613             std::cout << "Skipping LHE header with tag" << iter->tag() << std::endl;
0614           continue;
0615         }
0616         if (lheDebug)
0617           std::cout << "Found LHE header with tag" << iter->tag() << std::endl;
0618         std::vector<std::string> lines = iter->lines();
0619         bool missed_weightgroup =
0620             false;  //Needed because in some of the samples ( produced with MG26X ) a small part of the header info is ordered incorrectly
0621         bool ismg26x = false;
0622         bool ismg26xNew = false;
0623         for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines;
0624              ++iLine) {  //First start looping through the lines to see which weightgroup pattern is matched
0625           boost::replace_all(lines[iLine], "&lt;", "<");
0626           boost::replace_all(lines[iLine], "&gt;", ">");
0627           if (std::regex_search(lines[iLine], groups, weightgroupmg26x)) {
0628             ismg26x = true;
0629           } else if (std::regex_search(lines[iLine], groups, scalewmg26xNew) ||
0630                      std::regex_search(lines[iLine], groups, pdfwmg26xNew)) {
0631             ismg26xNew = true;
0632           }
0633         }
0634         for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) {
0635           if (lheDebug)
0636             std::cout << lines[iLine];
0637           auto foundWeightGroup = std::regex_search(lines[iLine], groups, ismg26x ? weightgroupmg26x : weightgroup);
0638           if (foundWeightGroup || preScaleVariationGroup) {
0639             std::string groupname;
0640             if (foundWeightGroup) {
0641               groupname = ismg26x ? groups.str(1) : groups.str(2);
0642             } else {
0643               // rewind by one line and check later in the inner loop
0644               --iLine;
0645             }
0646             if (lheDebug)
0647               std::cout << ">>> Looks like the beginning of a weight group for '" << groupname << "'" << std::endl;
0648             if (groupname.find("scale_variation") == 0 || groupname == "Central scale variation" ||
0649                 preScaleVariationGroup) {
0650               if (lheDebug && groupname.find("scale_variation") != 0 && groupname != "Central scale variation")
0651                 std::cout << ">>> First weight is not scale variation, but assuming is the Central Weight" << std::endl;
0652               else if (lheDebug)
0653                 std::cout << ">>> Looks like scale variation for theory uncertainties" << std::endl;
0654               if (groupname.find("scale_variation") == 0 || groupname == "Central scale variation") {
0655                 preScaleVariationGroup = false;
0656               }
0657               for (++iLine; iLine < nLines; ++iLine) {
0658                 if (lheDebug) {
0659                   std::cout << "    " << lines[iLine];
0660                 }
0661                 if (std::regex_search(
0662                         lines[iLine], groups, ismg26x ? scalewmg26x : (ismg26xNew ? scalewmg26xNew : scalew))) {
0663                   if (lheDebug)
0664                     std::cout << "    >>> Scale weight " << groups[1].str() << " for " << groups[3].str() << " , "
0665                               << groups[4].str() << " , " << groups[5].str() << std::endl;
0666                   if (ismg26xNew) {
0667                     scaleVariationIDs.emplace_back(groups.str(4), groups.str(1), groups.str(3), groups.str(2));
0668                   } else {
0669                     scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
0670                   }
0671                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0672                   if (lheDebug)
0673                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0674                   if (!missed_weightgroup) {
0675                     break;
0676                   } else
0677                     missed_weightgroup = false;
0678                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0679                   if (lheDebug)
0680                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0681                                  "of the group."
0682                               << std::endl;
0683                   if (ismg26x || ismg26xNew)
0684                     missed_weightgroup = true;
0685                   --iLine;  // rewind by one, and go back to the outer loop
0686                   break;
0687                 }
0688               }
0689             } else if (groupname == "PDF_variation" || groupname.find("PDF_variation ") == 0) {
0690               if (lheDebug)
0691                 std::cout << ">>> Looks like a new-style block of PDF weights for one or more pdfs" << std::endl;
0692               for (++iLine; iLine < nLines; ++iLine) {
0693                 if (lheDebug)
0694                   std::cout << "    " << lines[iLine];
0695                 if (std::regex_search(lines[iLine], groups, pdfw)) {
0696                   unsigned int lhaID = std::stoi(groups.str(2));
0697                   if (lheDebug)
0698                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
0699                               << std::endl;
0700                   if (pdfSetWeightIDs.empty() || !pdfSetWeightIDs.back().maybe_add(groups.str(1), lhaID)) {
0701                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0702                   }
0703                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0704                   if (lheDebug)
0705                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0706                   if (!missed_weightgroup) {
0707                     break;
0708                   } else
0709                     missed_weightgroup = false;
0710                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0711                   if (lheDebug)
0712                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0713                                  "of the group."
0714                               << std::endl;
0715                   if (ismg26x || ismg26xNew)
0716                     missed_weightgroup = true;
0717                   --iLine;  // rewind by one, and go back to the outer loop
0718                   break;
0719                 }
0720               }
0721             } else if (groupname == "PDF_variation1" || groupname == "PDF_variation2") {
0722               if (lheDebug)
0723                 std::cout << ">>> Looks like a new-style block of PDF weights for multiple pdfs" << std::endl;
0724               unsigned int lastid = 0;
0725               for (++iLine; iLine < nLines; ++iLine) {
0726                 if (lheDebug)
0727                   std::cout << "    " << lines[iLine];
0728                 if (std::regex_search(lines[iLine], groups, pdfw)) {
0729                   unsigned int id = std::stoi(groups.str(1));
0730                   unsigned int lhaID = std::stoi(groups.str(2));
0731                   if (lheDebug)
0732                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
0733                               << std::endl;
0734                   if (id != (lastid + 1) || pdfSetWeightIDs.empty()) {
0735                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0736                   } else {
0737                     pdfSetWeightIDs.back().add(groups.str(1), lhaID);
0738                   }
0739                   lastid = id;
0740                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0741                   if (lheDebug)
0742                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0743                   if (!missed_weightgroup) {
0744                     break;
0745                   } else
0746                     missed_weightgroup = false;
0747                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0748                   if (lheDebug)
0749                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0750                                  "of the group."
0751                               << std::endl;
0752                   if (ismg26x || ismg26xNew)
0753                     missed_weightgroup = true;
0754                   --iLine;  // rewind by one, and go back to the outer loop
0755                   break;
0756                 }
0757               }
0758             } else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) {
0759               if (lheDebug)
0760                 std::cout << ">>> Looks like an old-style PDF weight for an individual pdf" << std::endl;
0761               unsigned int firstLhaID = lhaNameToID_.find(groupname)->second;
0762               bool first = true;
0763               for (++iLine; iLine < nLines; ++iLine) {
0764                 if (lheDebug)
0765                   std::cout << "    " << lines[iLine];
0766                 if (std::regex_search(
0767                         lines[iLine], groups, ismg26x ? pdfwmg26x : (ismg26xNew ? pdfwmg26xNew : pdfwOld))) {
0768                   unsigned int member = 0;
0769                   if (!ismg26x && !ismg26xNew) {
0770                     member = std::stoi(groups.str(2));
0771                   } else if (ismg26xNew) {
0772                     if (!groups.str(3).empty()) {
0773                       member = std::stoi(groups.str(3));
0774                     }
0775                   } else {
0776                     if (!groups.str(4).empty()) {
0777                       member = std::stoi(groups.str(4));
0778                     }
0779                   }
0780                   unsigned int lhaID = member + firstLhaID;
0781                   if (lheDebug)
0782                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << member << " = " << lhaID
0783                               << std::endl;
0784                   //if (member == 0) continue; // let's keep also the central value for now
0785                   if (first) {
0786                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0787                     first = false;
0788                   } else {
0789                     pdfSetWeightIDs.back().add(groups.str(1), lhaID);
0790                   }
0791                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0792                   if (lheDebug)
0793                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0794                   if (!missed_weightgroup) {
0795                     break;
0796                   } else
0797                     missed_weightgroup = false;
0798                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0799                   if (lheDebug)
0800                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0801                                  "of the group."
0802                               << std::endl;
0803                   if (ismg26x || ismg26xNew)
0804                     missed_weightgroup = true;
0805                   --iLine;  // rewind by one, and go back to the outer loop
0806                   break;
0807                 }
0808               }
0809             } else if (groupname == "mass_variation" || groupname == "sthw2_variation" ||
0810                        groupname == "width_variation") {
0811               if (lheDebug)
0812                 std::cout << ">>> Looks like an EW parameter weight" << std::endl;
0813               for (++iLine; iLine < nLines; ++iLine) {
0814                 if (lheDebug)
0815                   std::cout << "    " << lines[iLine];
0816                 if (std::regex_search(lines[iLine], groups, rwgt)) {
0817                   std::string rwgtID = groups.str(1);
0818                   if (lheDebug)
0819                     std::cout << "    >>> LHE reweighting weight: " << rwgtID << std::endl;
0820                   if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
0821                     // we're only interested in the beggining of the block
0822                     lheReweighingIDs.emplace_back(rwgtID);
0823                   }
0824                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0825                   if (lheDebug)
0826                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0827                 }
0828               }
0829             } else {
0830               for (++iLine; iLine < nLines; ++iLine) {
0831                 if (lheDebug)
0832                   std::cout << "    " << lines[iLine];
0833                 if (std::regex_search(lines[iLine], groups, endweightgroup)) {
0834                   if (lheDebug)
0835                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0836                   if (!missed_weightgroup) {
0837                     break;
0838                   } else
0839                     missed_weightgroup = false;
0840                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0841                   if (lheDebug)
0842                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0843                                  "of the group."
0844                               << std::endl;
0845                   if (ismg26x || ismg26xNew)
0846                     missed_weightgroup = true;
0847                   --iLine;  // rewind by one, and go back to the outer loop
0848                   break;
0849                 }
0850               }
0851             }
0852           } else if (std::regex_search(lines[iLine], groups, weightgroupRwgt)) {
0853             std::string groupname = groups.str(1);
0854             if (groupname.find("mg_reweighting") != std::string::npos) {
0855               if (lheDebug)
0856                 std::cout << ">>> Looks like a LHE weights for reweighting" << std::endl;
0857               for (++iLine; iLine < nLines; ++iLine) {
0858                 if (lheDebug)
0859                   std::cout << "    " << lines[iLine];
0860                 if (std::regex_search(lines[iLine], groups, rwgt)) {
0861                   std::string rwgtID = groups.str(1);
0862                   if (lheDebug)
0863                     std::cout << "    >>> LHE reweighting weight: " << rwgtID << std::endl;
0864                   if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
0865                     // we're only interested in the beggining of the block
0866                     lheReweighingIDs.emplace_back(rwgtID);
0867                   }
0868                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0869                   if (lheDebug)
0870                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0871                   if (!missed_weightgroup) {
0872                     break;
0873                   } else
0874                     missed_weightgroup = false;
0875                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0876                   if (lheDebug)
0877                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0878                                  "of the group."
0879                               << std::endl;
0880                   if (ismg26x)
0881                     missed_weightgroup = true;
0882                   --iLine;  // rewind by one, and go back to the outer loop
0883                   break;
0884                 }
0885               }
0886             }
0887           }
0888         }
0889         //std::cout << "============= END [ " << iter->tag() << " ] ============ \n\n" << std::endl;
0890 
0891         // ----- SCALE VARIATIONS -----
0892         std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
0893         if (lheDebug)
0894           std::cout << "Found " << scaleVariationIDs.size() << " scale variations: " << std::endl;
0895         std::stringstream scaleDoc;
0896         scaleDoc << "LHE scale variation weights (w_var / w_nominal); ";
0897         for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
0898           const auto& sw = scaleVariationIDs[isw];
0899           if (isw)
0900             scaleDoc << "; ";
0901           scaleDoc << "[" << isw << "] is " << sw.label;
0902           weightChoice->scaleWeightIDs.push_back(sw.wid);
0903           if (lheDebug)
0904             printf("    id %s: scales ren = % .2f  fact = % .2f  text = %s\n",
0905                    sw.wid.c_str(),
0906                    sw.scales.first,
0907                    sw.scales.second,
0908                    sw.label.c_str());
0909         }
0910         if (!scaleVariationIDs.empty())
0911           weightChoice->scaleWeightsDoc = scaleDoc.str();
0912 
0913         // ------ PDF VARIATIONS (take the preferred one) -----
0914         if (lheDebug) {
0915           std::cout << "Found " << pdfSetWeightIDs.size() << " PDF set errors: " << std::endl;
0916           for (const auto& pw : pdfSetWeightIDs)
0917             printf("lhaIDs %6d - %6d (%3lu weights: %s, ... )\n",
0918                    pw.lhaIDs.first,
0919                    pw.lhaIDs.second,
0920                    pw.wids.size(),
0921                    pw.wids.front().c_str());
0922         }
0923 
0924         // ------ LHE REWEIGHTING -------
0925         if (lheDebug) {
0926           std::cout << "Found " << lheReweighingIDs.size() << " reweighting weights" << std::endl;
0927         }
0928         std::copy(lheReweighingIDs.begin(), lheReweighingIDs.end(), std::back_inserter(weightChoice->rwgtIDs));
0929 
0930         std::stringstream pdfDoc;
0931         pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA IDs ";
0932         bool found = false;
0933         for (const auto& pw : pdfSetWeightIDs) {
0934           for (uint32_t lhaid : preferredPDFLHAIDs_) {
0935             if (pw.lhaIDs.first != lhaid && pw.lhaIDs.first != (lhaid + 1))
0936               continue;  // sometimes the first weight is not saved if that PDF is the nominal one for the sample
0937             if (pw.wids.size() == 1)
0938               continue;  // only consider error sets
0939             pdfDoc << pw.lhaIDs.first << " - " << pw.lhaIDs.second;
0940             weightChoice->pdfWeightIDs = pw.wids;
0941             if (maxPdfWeights_ < pw.wids.size()) {
0942               weightChoice->pdfWeightIDs.resize(maxPdfWeights_);  // drop some replicas
0943               pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
0944             }
0945             weightChoice->pdfWeightsDoc = pdfDoc.str();
0946             found = true;
0947             break;
0948           }
0949           if (found)
0950             break;
0951         }
0952       }
0953       // check the number of scale variations
0954       if (isMGVer2x && !allowedNumScaleWeights_.empty()) {
0955         auto it = std::find(allowedNumScaleWeights_.begin(), allowedNumScaleWeights_.end(), scaleVariationIDs.size());
0956         if (it == allowedNumScaleWeights_.end()) {
0957           throw cms::Exception("LogicError")
0958               << "Number of scale variations found (" << scaleVariationIDs.size() << ") is invalid.";
0959         }
0960       }
0961     }
0962 
0963     return weightChoice;
0964   }
0965 
0966   // create an empty counter
0967   std::unique_ptr<LumiCacheInfoHolder> beginStream(edm::StreamID) const override {
0968     return std::make_unique<LumiCacheInfoHolder>();
0969   }
0970   // inizialize to zero at begin run
0971   void streamBeginRun(edm::StreamID id, edm::Run const&, edm::EventSetup const&) const override {
0972     streamCache(id)->clear();
0973   }
0974 
0975   std::shared_ptr<DynamicWeightChoiceGenInfo> globalBeginLuminosityBlock(edm::LuminosityBlock const& lumiBlock,
0976                                                                          edm::EventSetup const&) const override {
0977     auto dynamicWeightChoiceGenInfo = std::make_shared<DynamicWeightChoiceGenInfo>();
0978 
0979     edm::Handle<GenLumiInfoHeader> genLumiInfoHead;
0980     lumiBlock.getByToken(genLumiInfoHeadTag_, genLumiInfoHead);
0981     if (!genLumiInfoHead.isValid())
0982       edm::LogWarning("LHETablesProducer")
0983           << "No GenLumiInfoHeader product found, will not fill generator model string.\n";
0984 
0985     if (genLumiInfoHead.isValid()) {
0986       auto weightChoice = dynamicWeightChoiceGenInfo.get();
0987 
0988       std::vector<ScaleVarWeight> scaleVariationIDs;
0989       std::vector<PDFSetWeights> pdfSetWeightIDs;
0990       weightChoice->psWeightIDs.clear();
0991 
0992       std::regex scalew("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+)\\,\\s+mur=(\\S+)\\smuf=(\\S+)");
0993       std::regex pdfw("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+),\\s+Member\\s+(\\d+)\\s+of\\ssets\\s+(\\w+\\b)");
0994       std::regex mainPSw("sr(Def|:murfac=|\\.murfac=)(Hi|Lo|_dn|_up|0\\.5|2\\.0)");
0995       std::smatch groups;
0996       auto weightNames = genLumiInfoHead->weightNames();
0997       std::unordered_map<std::string, uint32_t> knownPDFSetsFromGenInfo_;
0998       unsigned int weightIter = 0;
0999       for (const auto& line : weightNames) {
1000         if (std::regex_search(line, groups, scalew)) {  // scale variation
1001           auto id = groups.str(1);
1002           auto group = groups.str(2);
1003           auto mur = groups.str(3);
1004           auto muf = groups.str(4);
1005           if (group.find("Central scale variation") != std::string::npos)
1006             scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
1007         } else if (std::regex_search(line, groups, pdfw)) {  // PDF variation
1008           auto id = groups.str(1);
1009           auto group = groups.str(2);
1010           auto memberid = groups.str(3);
1011           auto pdfset = groups.str(4);
1012           if (group.find(pdfset) != std::string::npos) {
1013             if (knownPDFSetsFromGenInfo_.find(pdfset) == knownPDFSetsFromGenInfo_.end()) {
1014               knownPDFSetsFromGenInfo_[pdfset] = std::atoi(id.c_str());
1015               pdfSetWeightIDs.emplace_back(id, std::atoi(id.c_str()));
1016             } else
1017               pdfSetWeightIDs.back().add(id, std::atoi(id.c_str()));
1018           }
1019         } else if (line == "Baseline") {
1020           weightChoice->psBaselineID = weightIter;
1021         } else if (line.find("isr") != std::string::npos || line.find("fsr") != std::string::npos) {
1022           weightChoice->matchPS_alt = line.find("sr:") != std::string::npos ||
1023                                       line.find("sr.") != std::string::npos;  // (f/i)sr: for new weights
1024           if (keepAllPSWeights_) {
1025             weightChoice->psWeightIDs.push_back(weightIter);  // PS variations
1026           } else if (std::regex_search(line, groups, mainPSw)) {
1027             if (weightChoice->psWeightIDs.empty())
1028               weightChoice->psWeightIDs = std::vector<unsigned int>(4, -1);
1029             int psIdx = (line.find("fsr") != std::string::npos) ? 1 : 0;
1030             psIdx += (groups.str(2) == "Hi" || groups.str(2) == "_up" || groups.str(2) == "2.0") ? 0 : 2;
1031             weightChoice->psWeightIDs[psIdx] = weightIter;
1032           }
1033         }
1034         weightIter++;
1035       }
1036       if (keepAllPSWeights_) {
1037         weightChoice->psWeightsDoc = "All PS weights (w_var / w_nominal) ";
1038       } else if (weightChoice->psWeightIDs.size() == 4) {
1039         weightChoice->psWeightsDoc = "PS weights (w_var / w_nominal) ";
1040         for (int i = 0; i < 4; i++) {
1041           if (static_cast<int>(weightChoice->psWeightIDs[i]) == -1)
1042             weightChoice->setMissingWeight(i);
1043         }
1044       } else {
1045         weightChoice->psWeightsDoc = "dummy PS weight (1.0) ";
1046       }
1047       for (unsigned i = 0; i < weightChoice->psWeightIDs.size(); ++i) {
1048         weightChoice->psWeightsDoc +=
1049             "[" + std::to_string(i) + "] " + weightNames.at(weightChoice->psWeightIDs.at(i)) + "; ";
1050       }
1051 
1052       weightChoice->scaleWeightIDs.clear();
1053       weightChoice->pdfWeightIDs.clear();
1054 
1055       std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
1056       std::stringstream scaleDoc;
1057       scaleDoc << "LHE scale variation weights (w_var / w_nominal); ";
1058       for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
1059         const auto& sw = scaleVariationIDs[isw];
1060         if (isw)
1061           scaleDoc << "; ";
1062         scaleDoc << "[" << isw << "] is " << sw.label;
1063         weightChoice->scaleWeightIDs.push_back(std::atoi(sw.wid.c_str()));
1064       }
1065       if (!scaleVariationIDs.empty())
1066         weightChoice->scaleWeightsDoc = scaleDoc.str();
1067       std::stringstream pdfDoc;
1068       pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA names ";
1069       bool found = false;
1070       for (const auto& pw : pdfSetWeightIDs) {
1071         if (pw.wids.size() == 1)
1072           continue;  // only consider error sets
1073         for (const auto& wantedpdf : lhaNameToID_) {
1074           auto pdfname = wantedpdf.first;
1075           if (knownPDFSetsFromGenInfo_.find(pdfname) == knownPDFSetsFromGenInfo_.end())
1076             continue;
1077           uint32_t lhaid = knownPDFSetsFromGenInfo_.at(pdfname);
1078           if (pw.lhaIDs.first != lhaid)
1079             continue;
1080           pdfDoc << pdfname;
1081           for (const auto& x : pw.wids)
1082             weightChoice->pdfWeightIDs.push_back(std::atoi(x.c_str()));
1083           if (maxPdfWeights_ < pw.wids.size()) {
1084             weightChoice->pdfWeightIDs.resize(maxPdfWeights_);  // drop some replicas
1085             pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
1086           }
1087           weightChoice->pdfWeightsDoc = pdfDoc.str();
1088           found = true;
1089           break;
1090         }
1091         if (found)
1092           break;
1093       }
1094     }
1095     return dynamicWeightChoiceGenInfo;
1096   }
1097 
1098   void globalEndLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override {}
1099 
1100   void streamBeginLuminosityBlock(edm::StreamID id,
1101                                   edm::LuminosityBlock const& lumiBlock,
1102                                   edm::EventSetup const&) const override {
1103     auto counterMap = &(streamCache(id)->countermap);
1104     edm::Handle<GenLumiInfoHeader> genLumiInfoHead;
1105     lumiBlock.getByToken(genLumiInfoHeadTag_, genLumiInfoHead);
1106 
1107     std::string label;
1108     if (genLumiInfoHead.isValid()) {
1109       label = genLumiInfoHead->configDescription();
1110       boost::replace_all(label, "-", "_");
1111       boost::replace_all(label, "/", "_");
1112     }
1113     counterMap->setLabel(label);
1114   }
1115   // create an empty counter
1116   std::shared_ptr<CounterMap> globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override {
1117     return std::make_shared<CounterMap>();
1118   }
1119   // add this stream to the summary
1120   void streamEndRunSummary(edm::StreamID id,
1121                            edm::Run const&,
1122                            edm::EventSetup const&,
1123                            CounterMap* runCounterMap) const override {
1124     runCounterMap->merge(streamCache(id)->countermap);
1125   }
1126   // nothing to do per se
1127   void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {}
1128   // write the total to the run
1129   void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, CounterMap const* runCounterMap) const override {
1130     auto out = std::make_unique<nanoaod::MergeableCounterTable>();
1131 
1132     for (const auto& x : runCounterMap->countermap) {
1133       auto runCounter = &(x.second);
1134       std::string label = (!x.first.empty()) ? (std::string("_") + x.first) : "";
1135       std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : "";
1136 
1137       out->addInt("genEventCount" + label, "event count" + doclabel, runCounter->num);
1138       out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter->sumw);
1139       out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter->sumw2);
1140 
1141       double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1;
1142       auto sumScales = runCounter->sumScale;
1143       for (auto& val : sumScales)
1144         val *= norm;
1145       out->addVFloatWithNorm("LHEScaleSumw" + label,
1146                              "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel,
1147                              sumScales,
1148                              runCounter->sumw);
1149       auto sumPDFs = runCounter->sumPDF;
1150       for (auto& val : sumPDFs)
1151         val *= norm;
1152       out->addVFloatWithNorm("LHEPdfSumw" + label,
1153                              "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel,
1154                              sumPDFs,
1155                              runCounter->sumw);
1156       auto sumPS = runCounter->sumPS;
1157       for (auto& val : sumPS)
1158         val *= norm;
1159       out->addVFloatWithNorm("PSSumw" + label,
1160                              "Sum of genEventWeight * PSWeight[i], divided by genEventSumw" + doclabel,
1161                              sumPS,
1162                              runCounter->sumw);
1163       if (!runCounter->sumRwgt.empty()) {
1164         auto sumRwgts = runCounter->sumRwgt;
1165         for (auto& val : sumRwgts)
1166           val *= norm;
1167         out->addVFloatWithNorm("LHEReweightingSumw" + label,
1168                                "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel,
1169                                sumRwgts,
1170                                runCounter->sumw);
1171       }
1172       if (!runCounter->sumNamed.empty()) {  // it could be empty if there's no LHE info in the sample
1173         for (unsigned int i = 0, n = namedWeightLabels_.size(); i < n; ++i) {
1174           out->addFloatWithNorm(
1175               "LHESumw_" + namedWeightLabels_[i] + label,
1176               "Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[i] + ", divided by genEventSumw" + doclabel,
1177               runCounter->sumNamed[i] * norm,
1178               runCounter->sumw);
1179         }
1180       }
1181     }
1182     iRun.put(std::move(out));
1183   }
1184   // nothing to do here
1185   void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {}
1186 
1187   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1188     edm::ParameterSetDescription desc;
1189     desc.add<edm::InputTag>("genEvent", edm::InputTag("generator"))
1190         ->setComment("tag for the GenEventInfoProduct, to get the main weight");
1191     desc.add<edm::InputTag>("genLumiInfoHeader", edm::InputTag("generator"))
1192         ->setComment("tag for the GenLumiInfoProduct, to get the model string");
1193     desc.add<std::vector<edm::InputTag>>("lheInfo", std::vector<edm::InputTag>{{"externalLHEProducer"}, {"source"}})
1194         ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
1195 
1196     edm::ParameterSetDescription prefpdf;
1197     prefpdf.add<std::string>("name");
1198     prefpdf.add<uint32_t>("lhaid");
1199     desc.addVPSet("preferredPDFs", prefpdf, std::vector<edm::ParameterSet>())
1200         ->setComment(
1201             "LHA PDF Ids of the preferred PDF sets, in order of preference (the first matching one will be used)");
1202     desc.add<std::vector<std::string>>("namedWeightIDs")->setComment("set of LHA weight IDs for named LHE weights");
1203     desc.add<std::vector<std::string>>("namedWeightLabels")
1204         ->setComment("output names for the namedWeightIDs (in the same order)");
1205     desc.add<int32_t>("lheWeightPrecision")->setComment("Number of bits in the mantissa for LHE weights");
1206     desc.add<uint32_t>("maxPdfWeights")->setComment("Maximum number of PDF weights to save (to crop NN replicas)");
1207     desc.add<bool>("keepAllPSWeights")->setComment("Store all PS weights found");
1208     desc.add<std::vector<uint32_t>>("allowedNumScaleWeights")
1209         ->setComment(
1210             "Allowed numbers of scale weights parsed from the header. Empty list means any number is allowed.");
1211     desc.addOptionalUntracked<bool>("debug")->setComment("dump out all LHE information for one event");
1212     descriptions.add("genWeightsTable", desc);
1213   }
1214 
1215 protected:
1216   const edm::EDGetTokenT<GenEventInfoProduct> genTag_;
1217   const std::vector<edm::InputTag> lheLabel_;
1218   const std::vector<edm::EDGetTokenT<LHEEventProduct>> lheTag_;
1219   const std::vector<edm::EDGetTokenT<LHERunInfoProduct>> lheRunTag_;
1220   const edm::EDGetTokenT<GenLumiInfoHeader> genLumiInfoHeadTag_;
1221 
1222   std::vector<uint32_t> preferredPDFLHAIDs_;
1223   std::unordered_map<std::string, uint32_t> lhaNameToID_;
1224   std::vector<std::string> namedWeightIDs_;
1225   std::vector<std::string> namedWeightLabels_;
1226   int lheWeightPrecision_;
1227   unsigned int maxPdfWeights_;
1228   bool keepAllPSWeights_;
1229   std::vector<uint32_t> allowedNumScaleWeights_;
1230 
1231   mutable std::atomic<bool> debug_, debugRun_, hasIssuedWarning_, psWeightWarning_;
1232 };
1233 
1234 #include "FWCore/Framework/interface/MakerMacros.h"
1235 DEFINE_FWK_MODULE(GenWeightsTableProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!