Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:23:41 

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         debug_(params.getUntrackedParameter<bool>("debug", false)),
0266         debugRun_(debug_.load()),
0267         hasIssuedWarning_(false),
0268         psWeightWarning_(false) {
0269     produces<nanoaod::FlatTable>();
0270     produces<std::string>("genModel");
0271     produces<nanoaod::FlatTable>("LHEScale");
0272     produces<nanoaod::FlatTable>("LHEPdf");
0273     produces<nanoaod::FlatTable>("LHEReweighting");
0274     produces<nanoaod::FlatTable>("LHENamed");
0275     produces<nanoaod::FlatTable>("PS");
0276     produces<nanoaod::MergeableCounterTable, edm::Transition::EndRun>();
0277     if (namedWeightIDs_.size() != namedWeightLabels_.size()) {
0278       throw cms::Exception("Configuration", "Size mismatch between namedWeightIDs & namedWeightLabels");
0279     }
0280     for (const edm::ParameterSet& pdfps : params.getParameter<std::vector<edm::ParameterSet>>("preferredPDFs")) {
0281       const std::string& name = pdfps.getParameter<std::string>("name");
0282       uint32_t lhaid = pdfps.getParameter<uint32_t>("lhaid");
0283       preferredPDFLHAIDs_.push_back(lhaid);
0284       lhaNameToID_[name] = lhaid;
0285       lhaNameToID_[name + ".LHgrid"] = lhaid;
0286     }
0287   }
0288 
0289   ~GenWeightsTableProducer() override {}
0290 
0291   void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
0292     // get my counter for weights
0293     Counter* counter = streamCache(id)->countermap.get();
0294 
0295     // generator information (always available)
0296     edm::Handle<GenEventInfoProduct> genInfo;
0297     iEvent.getByToken(genTag_, genInfo);
0298     double weight = genInfo->weight();
0299 
0300     // table for gen info, always available
0301     auto out = std::make_unique<nanoaod::FlatTable>(1, "genWeight", true);
0302     out->setDoc("generator weight");
0303     out->addColumnValue<float>("", weight, "generator weight");
0304     iEvent.put(std::move(out));
0305 
0306     std::string model_label = streamCache(id)->countermap.getLabel();
0307     auto outM = std::make_unique<std::string>((!model_label.empty()) ? std::string("GenModel_") + model_label : "");
0308     iEvent.put(std::move(outM), "genModel");
0309     bool getLHEweightsFromGenInfo = !model_label.empty();
0310 
0311     // tables for LHE weights, may not be filled
0312     std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab;
0313     std::unique_ptr<nanoaod::FlatTable> genPSTab;
0314 
0315     edm::Handle<LHEEventProduct> lheInfo;
0316     for (const auto& lheTag : lheTag_) {
0317       iEvent.getByToken(lheTag, lheInfo);
0318       if (lheInfo.isValid()) {
0319         break;
0320       }
0321     }
0322 
0323     const auto genWeightChoice = luminosityBlockCache(iEvent.getLuminosityBlock().index());
0324     if (lheInfo.isValid()) {
0325       if (getLHEweightsFromGenInfo && !hasIssuedWarning_.exchange(true))
0326         edm::LogWarning("LHETablesProducer")
0327             << "Found both a LHEEventProduct and a GenLumiInfoHeader: will only save weights from LHEEventProduct.\n";
0328       // get the dynamic choice of weights
0329       const DynamicWeightChoice* weightChoice = runCache(iEvent.getRun().index());
0330       // go fill tables
0331       fillLHEWeightTables(counter,
0332                           weightChoice,
0333                           genWeightChoice,
0334                           weight,
0335                           *lheInfo,
0336                           *genInfo,
0337                           lheScaleTab,
0338                           lhePdfTab,
0339                           lheRwgtTab,
0340                           lheNamedTab,
0341                           genPSTab);
0342     } else if (getLHEweightsFromGenInfo) {
0343       fillLHEPdfWeightTablesFromGenInfo(
0344           counter, genWeightChoice, weight, *genInfo, lheScaleTab, lhePdfTab, lheNamedTab, genPSTab);
0345       lheRwgtTab = std::make_unique<nanoaod::FlatTable>(1, "LHEReweightingWeights", true);
0346       //lheNamedTab.reset(new nanoaod::FlatTable(1, "LHENamedWeights", true));
0347       //genPSTab.reset(new nanoaod::FlatTable(1, "PSWeight", true));
0348     } else {
0349       // Still try to add the PS weights
0350       fillOnlyPSWeightTable(counter, genWeightChoice, weight, *genInfo, genPSTab);
0351       // make dummy values
0352       lheScaleTab = std::make_unique<nanoaod::FlatTable>(1, "LHEScaleWeights", true);
0353       lhePdfTab = std::make_unique<nanoaod::FlatTable>(1, "LHEPdfWeights", true);
0354       lheRwgtTab = std::make_unique<nanoaod::FlatTable>(1, "LHEReweightingWeights", true);
0355       lheNamedTab = std::make_unique<nanoaod::FlatTable>(1, "LHENamedWeights", true);
0356       if (!hasIssuedWarning_.exchange(true)) {
0357         edm::LogWarning("LHETablesProducer") << "No LHEEventProduct, so there will be no LHE Tables\n";
0358       }
0359     }
0360 
0361     iEvent.put(std::move(lheScaleTab), "LHEScale");
0362     iEvent.put(std::move(lhePdfTab), "LHEPdf");
0363     iEvent.put(std::move(lheRwgtTab), "LHEReweighting");
0364     iEvent.put(std::move(lheNamedTab), "LHENamed");
0365     iEvent.put(std::move(genPSTab), "PS");
0366   }
0367 
0368   void fillLHEWeightTables(Counter* counter,
0369                            const DynamicWeightChoice* weightChoice,
0370                            const DynamicWeightChoiceGenInfo* genWeightChoice,
0371                            double genWeight,
0372                            const LHEEventProduct& lheProd,
0373                            const GenEventInfoProduct& genProd,
0374                            std::unique_ptr<nanoaod::FlatTable>& outScale,
0375                            std::unique_ptr<nanoaod::FlatTable>& outPdf,
0376                            std::unique_ptr<nanoaod::FlatTable>& outRwgt,
0377                            std::unique_ptr<nanoaod::FlatTable>& outNamed,
0378                            std::unique_ptr<nanoaod::FlatTable>& outPS) const {
0379     bool lheDebug = debug_.exchange(
0380         false);  // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
0381 
0382     const std::vector<std::string>& scaleWeightIDs = weightChoice->scaleWeightIDs;
0383     const std::vector<std::string>& pdfWeightIDs = weightChoice->pdfWeightIDs;
0384     const std::vector<std::string>& rwgtWeightIDs = weightChoice->rwgtIDs;
0385 
0386     double w0 = lheProd.originalXWGTUP();
0387 
0388     std::vector<double> wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wRwgt(rwgtWeightIDs.size(), 1),
0389         wNamed(namedWeightIDs_.size(), 1);
0390     for (auto& weight : lheProd.weights()) {
0391       if (lheDebug)
0392         printf("Weight  %+9.5f   rel %+9.5f   for id %s\n", weight.wgt, weight.wgt / w0, weight.id.c_str());
0393       // now we do it slowly, can be optimized
0394       auto mScale = std::find(scaleWeightIDs.begin(), scaleWeightIDs.end(), weight.id);
0395       if (mScale != scaleWeightIDs.end())
0396         wScale[mScale - scaleWeightIDs.begin()] = weight.wgt / w0;
0397 
0398       auto mPDF = std::find(pdfWeightIDs.begin(), pdfWeightIDs.end(), weight.id);
0399       if (mPDF != pdfWeightIDs.end())
0400         wPDF[mPDF - pdfWeightIDs.begin()] = weight.wgt / w0;
0401 
0402       auto mRwgt = std::find(rwgtWeightIDs.begin(), rwgtWeightIDs.end(), weight.id);
0403       if (mRwgt != rwgtWeightIDs.end())
0404         wRwgt[mRwgt - rwgtWeightIDs.begin()] = weight.wgt / w0;
0405 
0406       auto mNamed = std::find(namedWeightIDs_.begin(), namedWeightIDs_.end(), weight.id);
0407       if (mNamed != namedWeightIDs_.end())
0408         wNamed[mNamed - namedWeightIDs_.begin()] = weight.wgt / w0;
0409     }
0410 
0411     std::vector<double> wPS;
0412     std::string psWeightDocStr;
0413     setPSWeightInfo(genProd.weights(), genWeightChoice, wPS, psWeightDocStr);
0414 
0415     outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
0416     outPS->addColumn<float>("", wPS, psWeightDocStr, lheWeightPrecision_);
0417 
0418     outScale = std::make_unique<nanoaod::FlatTable>(wScale.size(), "LHEScaleWeight", false);
0419     outScale->addColumn<float>("", wScale, weightChoice->scaleWeightsDoc, lheWeightPrecision_);
0420 
0421     outPdf = std::make_unique<nanoaod::FlatTable>(wPDF.size(), "LHEPdfWeight", false);
0422     outPdf->addColumn<float>("", wPDF, weightChoice->pdfWeightsDoc, lheWeightPrecision_);
0423 
0424     outRwgt = std::make_unique<nanoaod::FlatTable>(wRwgt.size(), "LHEReweightingWeight", false);
0425     outRwgt->addColumn<float>("", wRwgt, weightChoice->rwgtWeightDoc, lheWeightPrecision_);
0426 
0427     outNamed = std::make_unique<nanoaod::FlatTable>(1, "LHEWeight", true);
0428     outNamed->addColumnValue<float>("originalXWGTUP", lheProd.originalXWGTUP(), "Nominal event weight in the LHE file");
0429     for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) {
0430       outNamed->addColumnValue<float>(namedWeightLabels_[i],
0431                                       wNamed[i],
0432                                       "LHE weight for id " + namedWeightIDs_[i] + ", relative to nominal",
0433                                       lheWeightPrecision_);
0434     }
0435 
0436     counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS);
0437   }
0438 
0439   void fillLHEPdfWeightTablesFromGenInfo(Counter* counter,
0440                                          const DynamicWeightChoiceGenInfo* weightChoice,
0441                                          double genWeight,
0442                                          const GenEventInfoProduct& genProd,
0443                                          std::unique_ptr<nanoaod::FlatTable>& outScale,
0444                                          std::unique_ptr<nanoaod::FlatTable>& outPdf,
0445                                          std::unique_ptr<nanoaod::FlatTable>& outNamed,
0446                                          std::unique_ptr<nanoaod::FlatTable>& outPS) const {
0447     const std::vector<unsigned int>& scaleWeightIDs = weightChoice->scaleWeightIDs;
0448     const std::vector<unsigned int>& pdfWeightIDs = weightChoice->pdfWeightIDs;
0449 
0450     auto weights = genProd.weights();
0451     double w0 = (weights.size() > 1) ? weights.at(1) : 1.;
0452     double originalXWGTUP = (weights.size() > 1) ? weights.at(1) : 1.;
0453 
0454     std::vector<double> wScale, wPDF, wPS;
0455     for (auto id : scaleWeightIDs)
0456       wScale.push_back(weights.at(id) / w0);
0457     for (auto id : pdfWeightIDs) {
0458       wPDF.push_back(weights.at(id) / w0);
0459     }
0460 
0461     std::string psWeightsDocStr;
0462     setPSWeightInfo(genProd.weights(), weightChoice, wPS, psWeightsDocStr);
0463 
0464     outScale = std::make_unique<nanoaod::FlatTable>(wScale.size(), "LHEScaleWeight", false);
0465     outScale->addColumn<float>("", wScale, weightChoice->scaleWeightsDoc, lheWeightPrecision_);
0466 
0467     outPdf = std::make_unique<nanoaod::FlatTable>(wPDF.size(), "LHEPdfWeight", false);
0468     outPdf->addColumn<float>("", wPDF, weightChoice->pdfWeightsDoc, lheWeightPrecision_);
0469 
0470     outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
0471     outPS->addColumn<float>("", wPS, psWeightsDocStr, lheWeightPrecision_);
0472 
0473     outNamed = std::make_unique<nanoaod::FlatTable>(1, "LHEWeight", true);
0474     outNamed->addColumnValue<float>("originalXWGTUP", originalXWGTUP, "Nominal event weight in the LHE file");
0475     /*for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) {
0476       outNamed->addColumnValue<float>(namedWeightLabels_[i], wNamed[i], "LHE weight for id "+namedWeightIDs_[i]+", relative to nominal", lheWeightPrecision_);
0477       }*/
0478 
0479     counter->incLHE(genWeight, wScale, wPDF, std::vector<double>(), std::vector<double>(), wPS);
0480   }
0481 
0482   void fillOnlyPSWeightTable(Counter* counter,
0483                              const DynamicWeightChoiceGenInfo* genWeightChoice,
0484                              double genWeight,
0485                              const GenEventInfoProduct& genProd,
0486                              std::unique_ptr<nanoaod::FlatTable>& outPS) const {
0487     std::vector<double> wPS;
0488     std::string psWeightDocStr;
0489     setPSWeightInfo(genProd.weights(), genWeightChoice, wPS, psWeightDocStr);
0490     outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
0491     outPS->addColumn<float>("", wPS, psWeightDocStr, lheWeightPrecision_);
0492 
0493     counter->incGenOnly(genWeight);
0494     counter->incPSOnly(genWeight, wPS);
0495   }
0496 
0497   void setPSWeightInfo(const std::vector<double>& genWeights,
0498                        const DynamicWeightChoiceGenInfo* genWeightChoice,
0499                        std::vector<double>& wPS,
0500                        std::string& psWeightDocStr) const {
0501     wPS.clear();
0502     // isRegularPSSet = keeping all weights and the weights are a usual size, ie
0503     //                  all weights are PS weights (don't use header incase missing names)
0504     bool isRegularPSSet = keepAllPSWeights_ && (genWeights.size() == 14 || genWeights.size() == 46);
0505     if (!genWeightChoice->psWeightIDs.empty() && !isRegularPSSet) {
0506       psWeightDocStr = genWeightChoice->psWeightsDoc;
0507       double psNom = genWeights.at(genWeightChoice->psBaselineID);
0508       for (auto wgtidx : genWeightChoice->psWeightIDs) {
0509         wPS.push_back(genWeights.at(wgtidx) / psNom);
0510       }
0511     } else {
0512       int vectorSize =
0513           keepAllPSWeights_ ? (genWeights.size() - 2) : ((genWeights.size() == 14 || genWeights.size() == 46) ? 4 : 1);
0514 
0515       if (vectorSize > 1) {
0516         double nominal = genWeights.at(1);  // Called 'Baseline' in GenLumiInfoHeader
0517         if (keepAllPSWeights_) {
0518           for (int i = 0; i < vectorSize; i++) {
0519             wPS.push_back(genWeights.at(i + 2) / nominal);
0520           }
0521           psWeightDocStr = "All PS weights (w_var / w_nominal)";
0522         } else {
0523           if (!psWeightWarning_.exchange(true))
0524             edm::LogWarning("LHETablesProducer")
0525                 << "GenLumiInfoHeader not found: Central PartonShower weights will fill with the 6-10th entries \n"
0526                 << "    This may incorrect for some mcs (madgraph 2.6.1 with its `isr:murfact=0.5` have a differnt "
0527                    "order )";
0528           for (std::size_t i = 6; i < 10; i++) {
0529             wPS.push_back(genWeights.at(i) / nominal);
0530           }
0531           psWeightDocStr =
0532               "PS weights (w_var / w_nominal);   [0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2"
0533               "[2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;";
0534         }
0535       } else {
0536         wPS.push_back(1.0);
0537         psWeightDocStr = "dummy PS weight (1.0) ";
0538       }
0539     }
0540   }
0541 
0542   // create an empty counter
0543   std::shared_ptr<DynamicWeightChoice> globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override {
0544     edm::Handle<LHERunInfoProduct> lheInfo;
0545 
0546     bool lheDebug = debugRun_.exchange(
0547         false);  // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
0548     auto weightChoice = std::make_shared<DynamicWeightChoice>();
0549 
0550     // getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
0551     //if (iRun.getByToken(lheRunTag_, lheInfo)) {
0552     for (const auto& lheLabel : lheLabel_) {
0553       iRun.getByLabel(lheLabel, lheInfo);
0554       if (lheInfo.isValid()) {
0555         break;
0556       }
0557     }
0558     if (lheInfo.isValid()) {
0559       std::vector<ScaleVarWeight> scaleVariationIDs;
0560       std::vector<PDFSetWeights> pdfSetWeightIDs;
0561       std::vector<std::string> lheReweighingIDs;
0562       bool isFirstGroup = true;
0563 
0564       std::regex weightgroupmg26x("<weightgroup\\s+(?:name|type)=\"(.*)\"\\s+combine=\"(.*)\"\\s*>");
0565       std::regex weightgroup("<weightgroup\\s+combine=\"(.*)\"\\s+(?:name|type)=\"(.*)\"\\s*>");
0566       std::regex weightgroupRwgt("<weightgroup\\s+(?:name|type)=\"(.*)\"\\s*>");
0567       std::regex endweightgroup("</weightgroup>");
0568       std::regex scalewmg26x(
0569           "<weight\\s+(?:.*\\s+)?id=\"(\\d+)\"\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:[mM][uU][rR]|renscfact)=\"("
0570           "\\S+)\"\\s+(?:[mM][uU][Ff]|facscfact)=\"(\\S+)\")(\\s+.*)?</weight>");
0571       std::regex scalewmg26xNew(
0572           "<weight\\s*((?:[mM][uU][fF]|facscfact)=\"(\\S+)\"\\s+(?:[mM][uU][Rr]|renscfact)=\"(\\S+)\").+id=\"(\\d+)\"(."
0573           "*)?</weight>");
0574 
0575       //<weight MUF="1.0" MUR="2.0" PDF="306000" id="1006"> MUR=2.0  </weight>
0576       std::regex scalew(
0577           "<weight\\s+(?:.*\\s+)?id=\"(\\d+|\\d+-NNLOPS)\">\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:mu[rR]|renscfact)"
0578           "=(\\S+)\\s+(?:mu[Ff]|facscfact)=(\\S+)(\\s+.*)?)</weight>");
0579       std::regex pdfw(
0580           "<weight\\s+id=\"(\\d+)\">\\s*(?:PDF set|lhapdf|PDF|pdfset)\\s*=\\s*(\\d+)\\s*(?:\\s.*)?</weight>");
0581       std::regex pdfwOld("<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*Member \\s*(\\d+)\\s*(?:.*)</weight>");
0582       std::regex pdfwmg26x(
0583           "<weight\\s+id=\"(\\d+)\"\\s*MUR=\"(?:\\S+)\"\\s*MUF=\"(?:\\S+)\"\\s*(?:PDF "
0584           "set|lhapdf|PDF|pdfset)\\s*=\\s*\"(\\d+)\"\\s*>\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</"
0585           "weight>");
0586       //<weightgroup combine="symmhessian+as" name="NNPDF31_nnlo_as_0118_mc_hessian_pdfas">
0587 
0588       //<weight MUF="1.0" MUR="1.0" PDF="325300" id="1048"> PDF=325300 MemberID=0 </weight>
0589       std::regex pdfwmg26xNew(
0590           "<weight\\s+MUF=\"(?:\\S+)\"\\s*MUR=\"(?:\\S+)\"\\s*PDF=\"(?:\\S+)\"\\s*id=\"(\\S+)\"\\s*>"
0591           "\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</"
0592           "weight>");
0593 
0594       std::regex rwgt("<weight\\s+id=\"(.+)\">(.+)?(</weight>)?");
0595       std::smatch groups;
0596       for (auto iter = lheInfo->headers_begin(), end = lheInfo->headers_end(); iter != end; ++iter) {
0597         if (iter->tag() != "initrwgt") {
0598           if (lheDebug)
0599             std::cout << "Skipping LHE header with tag" << iter->tag() << std::endl;
0600           continue;
0601         }
0602         if (lheDebug)
0603           std::cout << "Found LHE header with tag" << iter->tag() << std::endl;
0604         std::vector<std::string> lines = iter->lines();
0605         bool missed_weightgroup =
0606             false;  //Needed because in some of the samples ( produced with MG26X ) a small part of the header info is ordered incorrectly
0607         bool ismg26x = false;
0608         bool ismg26xNew = false;
0609         for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines;
0610              ++iLine) {  //First start looping through the lines to see which weightgroup pattern is matched
0611           boost::replace_all(lines[iLine], "&lt;", "<");
0612           boost::replace_all(lines[iLine], "&gt;", ">");
0613           if (std::regex_search(lines[iLine], groups, weightgroupmg26x)) {
0614             ismg26x = true;
0615           } else if (std::regex_search(lines[iLine], groups, scalewmg26xNew) ||
0616                      std::regex_search(lines[iLine], groups, pdfwmg26xNew)) {
0617             ismg26xNew = true;
0618           }
0619         }
0620         for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) {
0621           if (lheDebug)
0622             std::cout << lines[iLine];
0623           if (std::regex_search(lines[iLine], groups, ismg26x ? weightgroupmg26x : weightgroup)) {
0624             std::string groupname = groups.str(2);
0625             if (ismg26x)
0626               groupname = groups.str(1);
0627             if (lheDebug)
0628               std::cout << ">>> Looks like the beginning of a weight group for '" << groupname << "'" << std::endl;
0629             if (groupname.find("scale_variation") == 0 || groupname == "Central scale variation" || isFirstGroup) {
0630               if (lheDebug && groupname.find("scale_variation") != 0 && groupname != "Central scale variation")
0631                 std::cout << ">>> First weight is not scale variation, but assuming is the Central Weight" << std::endl;
0632               else if (lheDebug)
0633                 std::cout << ">>> Looks like scale variation for theory uncertainties" << std::endl;
0634               isFirstGroup = false;
0635               for (++iLine; iLine < nLines; ++iLine) {
0636                 if (lheDebug) {
0637                   std::cout << "    " << lines[iLine];
0638                 }
0639                 if (std::regex_search(
0640                         lines[iLine], groups, ismg26x ? scalewmg26x : (ismg26xNew ? scalewmg26xNew : scalew))) {
0641                   if (lheDebug)
0642                     std::cout << "    >>> Scale weight " << groups[1].str() << " for " << groups[3].str() << " , "
0643                               << groups[4].str() << " , " << groups[5].str() << std::endl;
0644                   if (ismg26xNew) {
0645                     scaleVariationIDs.emplace_back(groups.str(4), groups.str(1), groups.str(3), groups.str(2));
0646                   } else {
0647                     scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
0648                   }
0649                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0650                   if (lheDebug)
0651                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0652                   if (!missed_weightgroup) {
0653                     break;
0654                   } else
0655                     missed_weightgroup = false;
0656                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0657                   if (lheDebug)
0658                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0659                                  "of the group."
0660                               << std::endl;
0661                   if (ismg26x || ismg26xNew)
0662                     missed_weightgroup = true;
0663                   --iLine;  // rewind by one, and go back to the outer loop
0664                   break;
0665                 }
0666               }
0667             } else if (groupname == "PDF_variation" || groupname.find("PDF_variation ") == 0) {
0668               if (lheDebug)
0669                 std::cout << ">>> Looks like a new-style block of PDF weights for one or more pdfs" << std::endl;
0670               for (++iLine; iLine < nLines; ++iLine) {
0671                 if (lheDebug)
0672                   std::cout << "    " << lines[iLine];
0673                 if (std::regex_search(lines[iLine], groups, pdfw)) {
0674                   unsigned int lhaID = std::stoi(groups.str(2));
0675                   if (lheDebug)
0676                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
0677                               << std::endl;
0678                   if (pdfSetWeightIDs.empty() || !pdfSetWeightIDs.back().maybe_add(groups.str(1), lhaID)) {
0679                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0680                   }
0681                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0682                   if (lheDebug)
0683                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0684                   if (!missed_weightgroup) {
0685                     break;
0686                   } else
0687                     missed_weightgroup = false;
0688                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0689                   if (lheDebug)
0690                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0691                                  "of the group."
0692                               << std::endl;
0693                   if (ismg26x || ismg26xNew)
0694                     missed_weightgroup = true;
0695                   --iLine;  // rewind by one, and go back to the outer loop
0696                   break;
0697                 }
0698               }
0699             } else if (groupname == "PDF_variation1" || groupname == "PDF_variation2") {
0700               if (lheDebug)
0701                 std::cout << ">>> Looks like a new-style block of PDF weights for multiple pdfs" << std::endl;
0702               unsigned int lastid = 0;
0703               for (++iLine; iLine < nLines; ++iLine) {
0704                 if (lheDebug)
0705                   std::cout << "    " << lines[iLine];
0706                 if (std::regex_search(lines[iLine], groups, pdfw)) {
0707                   unsigned int id = std::stoi(groups.str(1));
0708                   unsigned int lhaID = std::stoi(groups.str(2));
0709                   if (lheDebug)
0710                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
0711                               << std::endl;
0712                   if (id != (lastid + 1) || pdfSetWeightIDs.empty()) {
0713                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0714                   } else {
0715                     pdfSetWeightIDs.back().add(groups.str(1), lhaID);
0716                   }
0717                   lastid = id;
0718                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0719                   if (lheDebug)
0720                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0721                   if (!missed_weightgroup) {
0722                     break;
0723                   } else
0724                     missed_weightgroup = false;
0725                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0726                   if (lheDebug)
0727                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0728                                  "of the group."
0729                               << std::endl;
0730                   if (ismg26x || ismg26xNew)
0731                     missed_weightgroup = true;
0732                   --iLine;  // rewind by one, and go back to the outer loop
0733                   break;
0734                 }
0735               }
0736             } else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) {
0737               if (lheDebug)
0738                 std::cout << ">>> Looks like an old-style PDF weight for an individual pdf" << std::endl;
0739               unsigned int firstLhaID = lhaNameToID_.find(groupname)->second;
0740               bool first = true;
0741               for (++iLine; iLine < nLines; ++iLine) {
0742                 if (lheDebug)
0743                   std::cout << "    " << lines[iLine];
0744                 if (std::regex_search(
0745                         lines[iLine], groups, ismg26x ? pdfwmg26x : (ismg26xNew ? pdfwmg26xNew : pdfwOld))) {
0746                   unsigned int member = 0;
0747                   if (!ismg26x && !ismg26xNew) {
0748                     member = std::stoi(groups.str(2));
0749                   } else if (ismg26xNew) {
0750                     if (!groups.str(3).empty()) {
0751                       member = std::stoi(groups.str(3));
0752                     }
0753                   } else {
0754                     if (!groups.str(4).empty()) {
0755                       member = std::stoi(groups.str(4));
0756                     }
0757                   }
0758                   unsigned int lhaID = member + firstLhaID;
0759                   if (lheDebug)
0760                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << member << " = " << lhaID
0761                               << std::endl;
0762                   //if (member == 0) continue; // let's keep also the central value for now
0763                   if (first) {
0764                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0765                     first = false;
0766                   } else {
0767                     pdfSetWeightIDs.back().add(groups.str(1), lhaID);
0768                   }
0769                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0770                   if (lheDebug)
0771                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0772                   if (!missed_weightgroup) {
0773                     break;
0774                   } else
0775                     missed_weightgroup = false;
0776                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0777                   if (lheDebug)
0778                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0779                                  "of the group."
0780                               << std::endl;
0781                   if (ismg26x || ismg26xNew)
0782                     missed_weightgroup = true;
0783                   --iLine;  // rewind by one, and go back to the outer loop
0784                   break;
0785                 }
0786               }
0787             } else if (groupname == "mass_variation" || groupname == "sthw2_variation" ||
0788                        groupname == "width_variation") {
0789               if (lheDebug)
0790                 std::cout << ">>> Looks like an EW parameter weight" << std::endl;
0791               for (++iLine; iLine < nLines; ++iLine) {
0792                 if (lheDebug)
0793                   std::cout << "    " << lines[iLine];
0794                 if (std::regex_search(lines[iLine], groups, rwgt)) {
0795                   std::string rwgtID = groups.str(1);
0796                   if (lheDebug)
0797                     std::cout << "    >>> LHE reweighting weight: " << rwgtID << std::endl;
0798                   if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
0799                     // we're only interested in the beggining of the block
0800                     lheReweighingIDs.emplace_back(rwgtID);
0801                   }
0802                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0803                   if (lheDebug)
0804                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0805                 }
0806               }
0807             } else {
0808               for (++iLine; iLine < nLines; ++iLine) {
0809                 if (lheDebug)
0810                   std::cout << "    " << lines[iLine];
0811                 if (std::regex_search(lines[iLine], groups, endweightgroup)) {
0812                   if (lheDebug)
0813                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0814                   if (!missed_weightgroup) {
0815                     break;
0816                   } else
0817                     missed_weightgroup = false;
0818                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0819                   if (lheDebug)
0820                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0821                                  "of the group."
0822                               << std::endl;
0823                   if (ismg26x || ismg26xNew)
0824                     missed_weightgroup = true;
0825                   --iLine;  // rewind by one, and go back to the outer loop
0826                   break;
0827                 }
0828               }
0829             }
0830           } else if (std::regex_search(lines[iLine], groups, weightgroupRwgt)) {
0831             std::string groupname = groups.str(1);
0832             if (groupname.find("mg_reweighting") != std::string::npos) {
0833               if (lheDebug)
0834                 std::cout << ">>> Looks like a LHE weights for reweighting" << std::endl;
0835               for (++iLine; iLine < nLines; ++iLine) {
0836                 if (lheDebug)
0837                   std::cout << "    " << lines[iLine];
0838                 if (std::regex_search(lines[iLine], groups, rwgt)) {
0839                   std::string rwgtID = groups.str(1);
0840                   if (lheDebug)
0841                     std::cout << "    >>> LHE reweighting weight: " << rwgtID << std::endl;
0842                   if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
0843                     // we're only interested in the beggining of the block
0844                     lheReweighingIDs.emplace_back(rwgtID);
0845                   }
0846                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0847                   if (lheDebug)
0848                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0849                   if (!missed_weightgroup) {
0850                     break;
0851                   } else
0852                     missed_weightgroup = false;
0853                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0854                   if (lheDebug)
0855                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0856                                  "of the group."
0857                               << std::endl;
0858                   if (ismg26x)
0859                     missed_weightgroup = true;
0860                   --iLine;  // rewind by one, and go back to the outer loop
0861                   break;
0862                 }
0863               }
0864             }
0865           }
0866         }
0867         //std::cout << "============= END [ " << iter->tag() << " ] ============ \n\n" << std::endl;
0868 
0869         // ----- SCALE VARIATIONS -----
0870         std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
0871         if (lheDebug)
0872           std::cout << "Found " << scaleVariationIDs.size() << " scale variations: " << std::endl;
0873         std::stringstream scaleDoc;
0874         scaleDoc << "LHE scale variation weights (w_var / w_nominal); ";
0875         for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
0876           const auto& sw = scaleVariationIDs[isw];
0877           if (isw)
0878             scaleDoc << "; ";
0879           scaleDoc << "[" << isw << "] is " << sw.label;
0880           weightChoice->scaleWeightIDs.push_back(sw.wid);
0881           if (lheDebug)
0882             printf("    id %s: scales ren = % .2f  fact = % .2f  text = %s\n",
0883                    sw.wid.c_str(),
0884                    sw.scales.first,
0885                    sw.scales.second,
0886                    sw.label.c_str());
0887         }
0888         if (!scaleVariationIDs.empty())
0889           weightChoice->scaleWeightsDoc = scaleDoc.str();
0890 
0891         // ------ PDF VARIATIONS (take the preferred one) -----
0892         if (lheDebug) {
0893           std::cout << "Found " << pdfSetWeightIDs.size() << " PDF set errors: " << std::endl;
0894           for (const auto& pw : pdfSetWeightIDs)
0895             printf("lhaIDs %6d - %6d (%3lu weights: %s, ... )\n",
0896                    pw.lhaIDs.first,
0897                    pw.lhaIDs.second,
0898                    pw.wids.size(),
0899                    pw.wids.front().c_str());
0900         }
0901 
0902         // ------ LHE REWEIGHTING -------
0903         if (lheDebug) {
0904           std::cout << "Found " << lheReweighingIDs.size() << " reweighting weights" << std::endl;
0905         }
0906         std::copy(lheReweighingIDs.begin(), lheReweighingIDs.end(), std::back_inserter(weightChoice->rwgtIDs));
0907 
0908         std::stringstream pdfDoc;
0909         pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA IDs ";
0910         bool found = false;
0911         for (const auto& pw : pdfSetWeightIDs) {
0912           for (uint32_t lhaid : preferredPDFLHAIDs_) {
0913             if (pw.lhaIDs.first != lhaid && pw.lhaIDs.first != (lhaid + 1))
0914               continue;  // sometimes the first weight is not saved if that PDF is the nominal one for the sample
0915             if (pw.wids.size() == 1)
0916               continue;  // only consider error sets
0917             pdfDoc << pw.lhaIDs.first << " - " << pw.lhaIDs.second;
0918             weightChoice->pdfWeightIDs = pw.wids;
0919             if (maxPdfWeights_ < pw.wids.size()) {
0920               weightChoice->pdfWeightIDs.resize(maxPdfWeights_);  // drop some replicas
0921               pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
0922             }
0923             weightChoice->pdfWeightsDoc = pdfDoc.str();
0924             found = true;
0925             break;
0926           }
0927           if (found)
0928             break;
0929         }
0930       }
0931     }
0932     return weightChoice;
0933   }
0934 
0935   // create an empty counter
0936   std::unique_ptr<LumiCacheInfoHolder> beginStream(edm::StreamID) const override {
0937     return std::make_unique<LumiCacheInfoHolder>();
0938   }
0939   // inizialize to zero at begin run
0940   void streamBeginRun(edm::StreamID id, edm::Run const&, edm::EventSetup const&) const override {
0941     streamCache(id)->clear();
0942   }
0943 
0944   std::shared_ptr<DynamicWeightChoiceGenInfo> globalBeginLuminosityBlock(edm::LuminosityBlock const& lumiBlock,
0945                                                                          edm::EventSetup const&) const override {
0946     auto dynamicWeightChoiceGenInfo = std::make_shared<DynamicWeightChoiceGenInfo>();
0947 
0948     edm::Handle<GenLumiInfoHeader> genLumiInfoHead;
0949     lumiBlock.getByToken(genLumiInfoHeadTag_, genLumiInfoHead);
0950     if (!genLumiInfoHead.isValid())
0951       edm::LogWarning("LHETablesProducer")
0952           << "No GenLumiInfoHeader product found, will not fill generator model string.\n";
0953 
0954     if (genLumiInfoHead.isValid()) {
0955       auto weightChoice = dynamicWeightChoiceGenInfo.get();
0956 
0957       std::vector<ScaleVarWeight> scaleVariationIDs;
0958       std::vector<PDFSetWeights> pdfSetWeightIDs;
0959       weightChoice->psWeightIDs.clear();
0960 
0961       std::regex scalew("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+)\\,\\s+mur=(\\S+)\\smuf=(\\S+)");
0962       std::regex pdfw("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+),\\s+Member\\s+(\\d+)\\s+of\\ssets\\s+(\\w+\\b)");
0963       std::regex mainPSw("sr(Def|:murfac=)(Hi|Lo|_dn|_up|0.5|2.0)");
0964       std::smatch groups;
0965       auto weightNames = genLumiInfoHead->weightNames();
0966       std::unordered_map<std::string, uint32_t> knownPDFSetsFromGenInfo_;
0967       unsigned int weightIter = 0;
0968       for (const auto& line : weightNames) {
0969         if (std::regex_search(line, groups, scalew)) {  // scale variation
0970           auto id = groups.str(1);
0971           auto group = groups.str(2);
0972           auto mur = groups.str(3);
0973           auto muf = groups.str(4);
0974           if (group.find("Central scale variation") != std::string::npos)
0975             scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
0976         } else if (std::regex_search(line, groups, pdfw)) {  // PDF variation
0977           auto id = groups.str(1);
0978           auto group = groups.str(2);
0979           auto memberid = groups.str(3);
0980           auto pdfset = groups.str(4);
0981           if (group.find(pdfset) != std::string::npos) {
0982             if (knownPDFSetsFromGenInfo_.find(pdfset) == knownPDFSetsFromGenInfo_.end()) {
0983               knownPDFSetsFromGenInfo_[pdfset] = std::atoi(id.c_str());
0984               pdfSetWeightIDs.emplace_back(id, std::atoi(id.c_str()));
0985             } else
0986               pdfSetWeightIDs.back().add(id, std::atoi(id.c_str()));
0987           }
0988         } else if (line == "Baseline") {
0989           weightChoice->psBaselineID = weightIter;
0990         } else if (line.find("isr") != std::string::npos || line.find("fsr") != std::string::npos) {
0991           weightChoice->matchPS_alt = line.find("sr:") != std::string::npos;  // (f/i)sr: for new weights
0992           if (keepAllPSWeights_) {
0993             weightChoice->psWeightIDs.push_back(weightIter);  // PS variations
0994           } else if (std::regex_search(line, groups, mainPSw)) {
0995             if (weightChoice->psWeightIDs.empty())
0996               weightChoice->psWeightIDs = std::vector<unsigned int>(4, -1);
0997             int psIdx = (line.find("fsr") != std::string::npos) ? 1 : 0;
0998             psIdx += (groups.str(2) == "Hi" || groups.str(2) == "_up" || groups.str(2) == "2.0") ? 0 : 2;
0999             weightChoice->psWeightIDs[psIdx] = weightIter;
1000           }
1001         }
1002         weightIter++;
1003       }
1004       if (keepAllPSWeights_) {
1005         weightChoice->psWeightsDoc = "All PS weights (w_var / w_nominal)";
1006       } else if (weightChoice->psWeightIDs.size() == 4) {
1007         weightChoice->psWeightsDoc =
1008             "PS weights (w_var / w_nominal);   [0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2"
1009             "[2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;";
1010         for (int i = 0; i < 4; i++) {
1011           if (static_cast<int>(weightChoice->psWeightIDs[i]) == -1)
1012             weightChoice->setMissingWeight(i);
1013         }
1014       } else {
1015         weightChoice->psWeightsDoc = "dummy PS weight (1.0) ";
1016       }
1017 
1018       weightChoice->scaleWeightIDs.clear();
1019       weightChoice->pdfWeightIDs.clear();
1020 
1021       std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
1022       std::stringstream scaleDoc;
1023       scaleDoc << "LHE scale variation weights (w_var / w_nominal); ";
1024       for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
1025         const auto& sw = scaleVariationIDs[isw];
1026         if (isw)
1027           scaleDoc << "; ";
1028         scaleDoc << "[" << isw << "] is " << sw.label;
1029         weightChoice->scaleWeightIDs.push_back(std::atoi(sw.wid.c_str()));
1030       }
1031       if (!scaleVariationIDs.empty())
1032         weightChoice->scaleWeightsDoc = scaleDoc.str();
1033       std::stringstream pdfDoc;
1034       pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA names ";
1035       bool found = false;
1036       for (const auto& pw : pdfSetWeightIDs) {
1037         if (pw.wids.size() == 1)
1038           continue;  // only consider error sets
1039         for (const auto& wantedpdf : lhaNameToID_) {
1040           auto pdfname = wantedpdf.first;
1041           if (knownPDFSetsFromGenInfo_.find(pdfname) == knownPDFSetsFromGenInfo_.end())
1042             continue;
1043           uint32_t lhaid = knownPDFSetsFromGenInfo_.at(pdfname);
1044           if (pw.lhaIDs.first != lhaid)
1045             continue;
1046           pdfDoc << pdfname;
1047           for (const auto& x : pw.wids)
1048             weightChoice->pdfWeightIDs.push_back(std::atoi(x.c_str()));
1049           if (maxPdfWeights_ < pw.wids.size()) {
1050             weightChoice->pdfWeightIDs.resize(maxPdfWeights_);  // drop some replicas
1051             pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
1052           }
1053           weightChoice->pdfWeightsDoc = pdfDoc.str();
1054           found = true;
1055           break;
1056         }
1057         if (found)
1058           break;
1059       }
1060     }
1061     return dynamicWeightChoiceGenInfo;
1062   }
1063 
1064   void globalEndLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override {}
1065 
1066   void streamBeginLuminosityBlock(edm::StreamID id,
1067                                   edm::LuminosityBlock const& lumiBlock,
1068                                   edm::EventSetup const&) const override {
1069     auto counterMap = &(streamCache(id)->countermap);
1070     edm::Handle<GenLumiInfoHeader> genLumiInfoHead;
1071     lumiBlock.getByToken(genLumiInfoHeadTag_, genLumiInfoHead);
1072 
1073     std::string label;
1074     if (genLumiInfoHead.isValid()) {
1075       label = genLumiInfoHead->configDescription();
1076       boost::replace_all(label, "-", "_");
1077       boost::replace_all(label, "/", "_");
1078     }
1079     counterMap->setLabel(label);
1080   }
1081   // create an empty counter
1082   std::shared_ptr<CounterMap> globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override {
1083     return std::make_shared<CounterMap>();
1084   }
1085   // add this stream to the summary
1086   void streamEndRunSummary(edm::StreamID id,
1087                            edm::Run const&,
1088                            edm::EventSetup const&,
1089                            CounterMap* runCounterMap) const override {
1090     runCounterMap->merge(streamCache(id)->countermap);
1091   }
1092   // nothing to do per se
1093   void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {}
1094   // write the total to the run
1095   void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, CounterMap const* runCounterMap) const override {
1096     auto out = std::make_unique<nanoaod::MergeableCounterTable>();
1097 
1098     for (const auto& x : runCounterMap->countermap) {
1099       auto runCounter = &(x.second);
1100       std::string label = (!x.first.empty()) ? (std::string("_") + x.first) : "";
1101       std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : "";
1102 
1103       out->addInt("genEventCount" + label, "event count" + doclabel, runCounter->num);
1104       out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter->sumw);
1105       out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter->sumw2);
1106 
1107       double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1;
1108       auto sumScales = runCounter->sumScale;
1109       for (auto& val : sumScales)
1110         val *= norm;
1111       out->addVFloatWithNorm("LHEScaleSumw" + label,
1112                              "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel,
1113                              sumScales,
1114                              runCounter->sumw);
1115       auto sumPDFs = runCounter->sumPDF;
1116       for (auto& val : sumPDFs)
1117         val *= norm;
1118       out->addVFloatWithNorm("LHEPdfSumw" + label,
1119                              "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel,
1120                              sumPDFs,
1121                              runCounter->sumw);
1122       auto sumPS = runCounter->sumPS;
1123       for (auto& val : sumPS)
1124         val *= norm;
1125       out->addVFloatWithNorm("PSSumw" + label,
1126                              "Sum of genEventWeight * PSWeight[i], divided by genEventSumw" + doclabel,
1127                              sumPS,
1128                              runCounter->sumw);
1129       if (!runCounter->sumRwgt.empty()) {
1130         auto sumRwgts = runCounter->sumRwgt;
1131         for (auto& val : sumRwgts)
1132           val *= norm;
1133         out->addVFloatWithNorm("LHEReweightingSumw" + label,
1134                                "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel,
1135                                sumRwgts,
1136                                runCounter->sumw);
1137       }
1138       if (!runCounter->sumNamed.empty()) {  // it could be empty if there's no LHE info in the sample
1139         for (unsigned int i = 0, n = namedWeightLabels_.size(); i < n; ++i) {
1140           out->addFloatWithNorm(
1141               "LHESumw_" + namedWeightLabels_[i] + label,
1142               "Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[i] + ", divided by genEventSumw" + doclabel,
1143               runCounter->sumNamed[i] * norm,
1144               runCounter->sumw);
1145         }
1146       }
1147     }
1148     iRun.put(std::move(out));
1149   }
1150   // nothing to do here
1151   void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {}
1152 
1153   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1154     edm::ParameterSetDescription desc;
1155     desc.add<edm::InputTag>("genEvent", edm::InputTag("generator"))
1156         ->setComment("tag for the GenEventInfoProduct, to get the main weight");
1157     desc.add<edm::InputTag>("genLumiInfoHeader", edm::InputTag("generator"))
1158         ->setComment("tag for the GenLumiInfoProduct, to get the model string");
1159     desc.add<std::vector<edm::InputTag>>("lheInfo", std::vector<edm::InputTag>{{"externalLHEProducer"}, {"source"}})
1160         ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
1161 
1162     edm::ParameterSetDescription prefpdf;
1163     prefpdf.add<std::string>("name");
1164     prefpdf.add<uint32_t>("lhaid");
1165     desc.addVPSet("preferredPDFs", prefpdf, std::vector<edm::ParameterSet>())
1166         ->setComment(
1167             "LHA PDF Ids of the preferred PDF sets, in order of preference (the first matching one will be used)");
1168     desc.add<std::vector<std::string>>("namedWeightIDs")->setComment("set of LHA weight IDs for named LHE weights");
1169     desc.add<std::vector<std::string>>("namedWeightLabels")
1170         ->setComment("output names for the namedWeightIDs (in the same order)");
1171     desc.add<int32_t>("lheWeightPrecision")->setComment("Number of bits in the mantissa for LHE weights");
1172     desc.add<uint32_t>("maxPdfWeights")->setComment("Maximum number of PDF weights to save (to crop NN replicas)");
1173     desc.add<bool>("keepAllPSWeights")->setComment("Store all PS weights found");
1174     desc.addOptionalUntracked<bool>("debug")->setComment("dump out all LHE information for one event");
1175     descriptions.add("genWeightsTable", desc);
1176   }
1177 
1178 protected:
1179   const edm::EDGetTokenT<GenEventInfoProduct> genTag_;
1180   const std::vector<edm::InputTag> lheLabel_;
1181   const std::vector<edm::EDGetTokenT<LHEEventProduct>> lheTag_;
1182   const std::vector<edm::EDGetTokenT<LHERunInfoProduct>> lheRunTag_;
1183   const edm::EDGetTokenT<GenLumiInfoHeader> genLumiInfoHeadTag_;
1184 
1185   std::vector<uint32_t> preferredPDFLHAIDs_;
1186   std::unordered_map<std::string, uint32_t> lhaNameToID_;
1187   std::vector<std::string> namedWeightIDs_;
1188   std::vector<std::string> namedWeightLabels_;
1189   int lheWeightPrecision_;
1190   unsigned int maxPdfWeights_;
1191   bool keepAllPSWeights_;
1192 
1193   mutable std::atomic<bool> debug_, debugRun_, hasIssuedWarning_, psWeightWarning_;
1194 };
1195 
1196 #include "FWCore/Framework/interface/MakerMacros.h"
1197 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!