Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-24 02:15:00 

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 <memory>
0018 
0019 #include <vector>
0020 #include <unordered_map>
0021 #include <iostream>
0022 #include <regex>
0023 
0024 namespace {
0025   ///  ---- Cache object for running sums of weights ----
0026   struct Counter {
0027     Counter() : num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumRwgt(), sumNamed(), sumPS() {}
0028 
0029     // the counters
0030     long long num;
0031     long double sumw;
0032     long double sumw2;
0033     std::vector<long double> sumPDF, sumScale, sumRwgt, sumNamed, sumPS;
0034 
0035     void clear() {
0036       num = 0;
0037       sumw = 0;
0038       sumw2 = 0;
0039       sumPDF.clear();
0040       sumScale.clear();
0041       sumRwgt.clear();
0042       sumNamed.clear(), sumPS.clear();
0043     }
0044 
0045     // inc the counters
0046     void incGenOnly(double w) {
0047       num++;
0048       sumw += w;
0049       sumw2 += (w * w);
0050     }
0051 
0052     void incPSOnly(double w0, const std::vector<double>& wPS) {
0053       if (!wPS.empty()) {
0054         if (sumPS.empty())
0055           sumPS.resize(wPS.size(), 0);
0056         for (unsigned int i = 0, n = wPS.size(); i < n; ++i)
0057           sumPS[i] += (w0 * wPS[i]);
0058       }
0059     }
0060 
0061     void incLHE(double w0,
0062                 const std::vector<double>& wScale,
0063                 const std::vector<double>& wPDF,
0064                 const std::vector<double>& wRwgt,
0065                 const std::vector<double>& wNamed,
0066                 const std::vector<double>& wPS) {
0067       // add up weights
0068       incGenOnly(w0);
0069       // then add up variations
0070       if (!wScale.empty()) {
0071         if (sumScale.empty())
0072           sumScale.resize(wScale.size(), 0);
0073         for (unsigned int i = 0, n = wScale.size(); i < n; ++i)
0074           sumScale[i] += (w0 * wScale[i]);
0075       }
0076       if (!wPDF.empty()) {
0077         if (sumPDF.empty())
0078           sumPDF.resize(wPDF.size(), 0);
0079         for (unsigned int i = 0, n = wPDF.size(); i < n; ++i)
0080           sumPDF[i] += (w0 * wPDF[i]);
0081       }
0082       if (!wRwgt.empty()) {
0083         if (sumRwgt.empty())
0084           sumRwgt.resize(wRwgt.size(), 0);
0085         for (unsigned int i = 0, n = wRwgt.size(); i < n; ++i)
0086           sumRwgt[i] += (w0 * wRwgt[i]);
0087       }
0088       if (!wNamed.empty()) {
0089         if (sumNamed.empty())
0090           sumNamed.resize(wNamed.size(), 0);
0091         for (unsigned int i = 0, n = wNamed.size(); i < n; ++i)
0092           sumNamed[i] += (w0 * wNamed[i]);
0093       }
0094       incPSOnly(w0, wPS);
0095     }
0096 
0097     void merge(const Counter& other) {
0098       num += other.num;
0099       sumw += other.sumw;
0100       sumw2 += other.sumw2;
0101       if (sumScale.empty() && !other.sumScale.empty())
0102         sumScale.resize(other.sumScale.size(), 0);
0103       if (sumPDF.empty() && !other.sumPDF.empty())
0104         sumPDF.resize(other.sumPDF.size(), 0);
0105       if (sumRwgt.empty() && !other.sumRwgt.empty())
0106         sumRwgt.resize(other.sumRwgt.size(), 0);
0107       if (sumNamed.empty() && !other.sumNamed.empty())
0108         sumNamed.resize(other.sumNamed.size(), 0);
0109       if (sumPS.empty() && !other.sumPS.empty())
0110         sumPS.resize(other.sumPS.size(), 0);
0111       if (!other.sumScale.empty())
0112         for (unsigned int i = 0, n = sumScale.size(); i < n; ++i)
0113           sumScale[i] += other.sumScale[i];
0114       if (!other.sumPDF.empty())
0115         for (unsigned int i = 0, n = sumPDF.size(); i < n; ++i)
0116           sumPDF[i] += other.sumPDF[i];
0117       if (!other.sumRwgt.empty())
0118         for (unsigned int i = 0, n = sumRwgt.size(); i < n; ++i)
0119           sumRwgt[i] += other.sumRwgt[i];
0120       if (!other.sumNamed.empty())
0121         for (unsigned int i = 0, n = sumNamed.size(); i < n; ++i)
0122           sumNamed[i] += other.sumNamed[i];
0123       if (!other.sumPS.empty())
0124         for (unsigned int i = 0, n = sumPS.size(); i < n; ++i)
0125           sumPS[i] += other.sumPS[i];
0126     }
0127   };
0128 
0129   struct CounterMap {
0130     std::map<std::string, Counter> countermap;
0131     Counter* active_el = nullptr;
0132     std::string active_label = "";
0133     void merge(const CounterMap& other) {
0134       for (const auto& y : other.countermap)
0135         countermap[y.first].merge(y.second);
0136       active_el = nullptr;
0137     }
0138     void clear() {
0139       for (auto x : countermap)
0140         x.second.clear();
0141       active_el = nullptr;
0142       active_label = "";
0143     }
0144     void setLabel(std::string label) {
0145       active_el = &(countermap[label]);
0146       active_label = label;
0147     }
0148     void checkLabelSet() {
0149       if (!active_el)
0150         throw cms::Exception("LogicError", "Called CounterMap::get() before setting the active label\n");
0151     }
0152     Counter* get() {
0153       checkLabelSet();
0154       return active_el;
0155     }
0156     std::string& getLabel() {
0157       checkLabelSet();
0158       return active_label;
0159     }
0160   };
0161 
0162   ///  ---- RunCache object for dynamic choice of LHE IDs ----
0163   struct DynamicWeightChoice {
0164     // choice of LHE weights
0165     // ---- scale ----
0166     std::vector<std::string> scaleWeightIDs;
0167     std::string scaleWeightsDoc;
0168     // ---- pdf ----
0169     std::vector<std::string> pdfWeightIDs;
0170     std::string pdfWeightsDoc;
0171     // ---- rwgt ----
0172     std::vector<std::string> rwgtIDs;
0173     std::string rwgtWeightDoc;
0174   };
0175 
0176   struct DynamicWeightChoiceGenInfo {
0177     // choice of LHE weights
0178     // ---- scale ----
0179     std::vector<unsigned int> scaleWeightIDs;
0180     std::string scaleWeightsDoc;
0181     // ---- pdf ----
0182     std::vector<unsigned int> pdfWeightIDs;
0183     std::string pdfWeightsDoc;
0184     // ---- ps ----
0185     std::vector<unsigned int> defPSWeightIDs = {6, 7, 8, 9};
0186     std::vector<unsigned int> defPSWeightIDs_alt = {27, 5, 26, 4};
0187     bool matchPS_alt = false;
0188     std::vector<unsigned int> psWeightIDs;
0189     unsigned int psBaselineID = 1;
0190     std::string psWeightsDoc;
0191 
0192     void setMissingWeight(int idx) { psWeightIDs[idx] = (matchPS_alt) ? defPSWeightIDs_alt[idx] : defPSWeightIDs[idx]; }
0193 
0194     bool empty() const { return scaleWeightIDs.empty() && pdfWeightIDs.empty() && psWeightIDs.empty(); }
0195   };
0196 
0197   struct LumiCacheInfoHolder {
0198     CounterMap countermap;
0199     DynamicWeightChoiceGenInfo weightChoice;
0200     void clear() {
0201       countermap.clear();
0202       weightChoice = DynamicWeightChoiceGenInfo();
0203     }
0204   };
0205 
0206   float stof_fortrancomp(const std::string& str) {
0207     std::string::size_type match = str.find('d');
0208     if (match != std::string::npos) {
0209       std::string pre = str.substr(0, match);
0210       std::string post = str.substr(match + 1);
0211       return std::stof(pre) * std::pow(10.0f, std::stof(post));
0212     } else {
0213       return std::stof(str);
0214     }
0215   }
0216   ///  -------------- temporary objects --------------
0217   struct ScaleVarWeight {
0218     std::string wid, label;
0219     std::pair<float, float> scales;
0220     ScaleVarWeight(const std::string& id, const std::string& text, const std::string& muR, const std::string& muF)
0221         : wid(id), label(text), scales(stof_fortrancomp(muR), stof_fortrancomp(muF)) {}
0222     bool operator<(const ScaleVarWeight& other) {
0223       return (scales == other.scales ? wid < other.wid : scales < other.scales);
0224     }
0225   };
0226   struct PDFSetWeights {
0227     std::vector<std::string> wids;
0228     std::pair<unsigned int, unsigned int> lhaIDs;
0229     PDFSetWeights(const std::string& wid, unsigned int lhaID) : wids(1, wid), lhaIDs(lhaID, lhaID) {}
0230     bool operator<(const PDFSetWeights& other) const { return lhaIDs < other.lhaIDs; }
0231     void add(const std::string& wid, unsigned int lhaID) {
0232       wids.push_back(wid);
0233       lhaIDs.second = lhaID;
0234     }
0235     bool maybe_add(const std::string& wid, unsigned int lhaID) {
0236       if (lhaID == lhaIDs.second + 1) {
0237         lhaIDs.second++;
0238         wids.push_back(wid);
0239         return true;
0240       } else {
0241         return false;
0242       }
0243     }
0244   };
0245 }  // namespace
0246 
0247 class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<LumiCacheInfoHolder>,
0248                                                                edm::RunCache<DynamicWeightChoice>,
0249                                                                edm::RunSummaryCache<CounterMap>,
0250                                                                edm::EndRunProducer> {
0251 public:
0252   GenWeightsTableProducer(edm::ParameterSet const& params)
0253       : genTag_(consumes<GenEventInfoProduct>(params.getParameter<edm::InputTag>("genEvent"))),
0254         lheLabel_(params.getParameter<std::vector<edm::InputTag>>("lheInfo")),
0255         lheTag_(edm::vector_transform(lheLabel_,
0256                                       [this](const edm::InputTag& tag) { return mayConsume<LHEEventProduct>(tag); })),
0257         lheRunTag_(edm::vector_transform(
0258             lheLabel_, [this](const edm::InputTag& tag) { return mayConsume<LHERunInfoProduct, edm::InRun>(tag); })),
0259         genLumiInfoHeadTag_(
0260             mayConsume<GenLumiInfoHeader, edm::InLumi>(params.getParameter<edm::InputTag>("genLumiInfoHeader"))),
0261         namedWeightIDs_(params.getParameter<std::vector<std::string>>("namedWeightIDs")),
0262         namedWeightLabels_(params.getParameter<std::vector<std::string>>("namedWeightLabels")),
0263         lheWeightPrecision_(params.getParameter<int32_t>("lheWeightPrecision")),
0264         maxPdfWeights_(params.getParameter<uint32_t>("maxPdfWeights")),
0265         keepAllPSWeights_(params.getParameter<bool>("keepAllPSWeights")),
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 = &(streamCache(id)->weightChoice);
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 = "All PS weights (w_var / w_nominal)";
0523         } else {
0524           if (!psWeightWarning_.exchange(true))
0525             edm::LogWarning("LHETablesProducer")
0526                 << "GenLumiInfoHeader not found: Central PartonShower weights will fill with the 6-10th entries \n"
0527                 << "    This may incorrect for some mcs (madgraph 2.6.1 with its `isr:murfact=0.5` have a differnt "
0528                    "order )";
0529           for (std::size_t i = 6; i < 10; i++) {
0530             wPS.push_back(genWeights.at(i) / nominal);
0531           }
0532           psWeightDocStr =
0533               "PS weights (w_var / w_nominal);   [0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2"
0534               "[2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;";
0535         }
0536       } else {
0537         wPS.push_back(1.0);
0538         psWeightDocStr = "dummy PS weight (1.0) ";
0539       }
0540     }
0541   }
0542 
0543   // create an empty counter
0544   std::shared_ptr<DynamicWeightChoice> globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override {
0545     edm::Handle<LHERunInfoProduct> lheInfo;
0546 
0547     bool lheDebug = debugRun_.exchange(
0548         false);  // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
0549     auto weightChoice = std::make_shared<DynamicWeightChoice>();
0550 
0551     // getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
0552     //if (iRun.getByToken(lheRunTag_, lheInfo)) {
0553     for (const auto& lheLabel : lheLabel_) {
0554       iRun.getByLabel(lheLabel, lheInfo);
0555       if (lheInfo.isValid()) {
0556         break;
0557       }
0558     }
0559     if (lheInfo.isValid()) {
0560       std::vector<ScaleVarWeight> scaleVariationIDs;
0561       std::vector<PDFSetWeights> pdfSetWeightIDs;
0562       std::vector<std::string> lheReweighingIDs;
0563       bool isFirstGroup = true;
0564 
0565       std::regex weightgroupmg26x("<weightgroup\\s+(?:name|type)=\"(.*)\"\\s+combine=\"(.*)\"\\s*>");
0566       std::regex weightgroup("<weightgroup\\s+combine=\"(.*)\"\\s+(?:name|type)=\"(.*)\"\\s*>");
0567       std::regex weightgroupRwgt("<weightgroup\\s+(?:name|type)=\"(.*)\"\\s*>");
0568       std::regex endweightgroup("</weightgroup>");
0569       std::regex scalewmg26x(
0570           "<weight\\s+(?:.*\\s+)?id=\"(\\d+)\"\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:[mM][uU][rR]|renscfact)=\"("
0571           "\\S+)\"\\s+(?:[mM][uU][Ff]|facscfact)=\"(\\S+)\")(\\s+.*)?</weight>");
0572       std::regex scalewmg26xNew(
0573           "<weight\\s*((?:[mM][uU][fF]|facscfact)=\"(\\S+)\"\\s+(?:[mM][uU][Rr]|renscfact)=\"(\\S+)\").+id=\"(\\d+)\"(."
0574           "*)?</weight>");
0575 
0576       //<weight MUF="1.0" MUR="2.0" PDF="306000" id="1006"> MUR=2.0  </weight>
0577       std::regex scalew(
0578           "<weight\\s+(?:.*\\s+)?id=\"(\\d+|\\d+-NNLOPS)\">\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:mu[rR]|renscfact)"
0579           "=(\\S+)\\s+(?:mu[Ff]|facscfact)=(\\S+)(\\s+.*)?)</weight>");
0580       std::regex pdfw(
0581           "<weight\\s+id=\"(\\d+)\">\\s*(?:PDF set|lhapdf|PDF|pdfset)\\s*=\\s*(\\d+)\\s*(?:\\s.*)?</weight>");
0582       std::regex pdfwOld("<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*Member \\s*(\\d+)\\s*(?:.*)</weight>");
0583       std::regex pdfwmg26x(
0584           "<weight\\s+id=\"(\\d+)\"\\s*MUR=\"(?:\\S+)\"\\s*MUF=\"(?:\\S+)\"\\s*(?:PDF "
0585           "set|lhapdf|PDF|pdfset)\\s*=\\s*\"(\\d+)\"\\s*>\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</"
0586           "weight>");
0587       //<weightgroup combine="symmhessian+as" name="NNPDF31_nnlo_as_0118_mc_hessian_pdfas">
0588 
0589       //<weight MUF="1.0" MUR="1.0" PDF="325300" id="1048"> PDF=325300 MemberID=0 </weight>
0590       std::regex pdfwmg26xNew(
0591           "<weight\\s+MUF=\"(?:\\S+)\"\\s*MUR=\"(?:\\S+)\"\\s*PDF=\"(?:\\S+)\"\\s*id=\"(\\S+)\"\\s*>"
0592           "\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</"
0593           "weight>");
0594 
0595       std::regex rwgt("<weight\\s+id=\"(.+)\">(.+)?(</weight>)?");
0596       std::smatch groups;
0597       for (auto iter = lheInfo->headers_begin(), end = lheInfo->headers_end(); iter != end; ++iter) {
0598         if (iter->tag() != "initrwgt") {
0599           if (lheDebug)
0600             std::cout << "Skipping LHE header with tag" << iter->tag() << std::endl;
0601           continue;
0602         }
0603         if (lheDebug)
0604           std::cout << "Found LHE header with tag" << iter->tag() << std::endl;
0605         std::vector<std::string> lines = iter->lines();
0606         bool missed_weightgroup =
0607             false;  //Needed because in some of the samples ( produced with MG26X ) a small part of the header info is ordered incorrectly
0608         bool ismg26x = false;
0609         bool ismg26xNew = false;
0610         for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines;
0611              ++iLine) {  //First start looping through the lines to see which weightgroup pattern is matched
0612           boost::replace_all(lines[iLine], "&lt;", "<");
0613           boost::replace_all(lines[iLine], "&gt;", ">");
0614           if (std::regex_search(lines[iLine], groups, weightgroupmg26x)) {
0615             ismg26x = true;
0616           } else if (std::regex_search(lines[iLine], groups, scalewmg26xNew) ||
0617                      std::regex_search(lines[iLine], groups, pdfwmg26xNew)) {
0618             ismg26xNew = true;
0619           }
0620         }
0621         for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) {
0622           if (lheDebug)
0623             std::cout << lines[iLine];
0624           if (std::regex_search(lines[iLine], groups, ismg26x ? weightgroupmg26x : weightgroup)) {
0625             std::string groupname = groups.str(2);
0626             if (ismg26x)
0627               groupname = groups.str(1);
0628             if (lheDebug)
0629               std::cout << ">>> Looks like the beginning of a weight group for '" << groupname << "'" << std::endl;
0630             if (groupname.find("scale_variation") == 0 || groupname == "Central scale variation" || isFirstGroup) {
0631               if (lheDebug && groupname.find("scale_variation") != 0 && groupname != "Central scale variation")
0632                 std::cout << ">>> First weight is not scale variation, but assuming is the Central Weight" << std::endl;
0633               else if (lheDebug)
0634                 std::cout << ">>> Looks like scale variation for theory uncertainties" << std::endl;
0635               isFirstGroup = false;
0636               for (++iLine; iLine < nLines; ++iLine) {
0637                 if (lheDebug) {
0638                   std::cout << "    " << lines[iLine];
0639                 }
0640                 if (std::regex_search(
0641                         lines[iLine], groups, ismg26x ? scalewmg26x : (ismg26xNew ? scalewmg26xNew : scalew))) {
0642                   if (lheDebug)
0643                     std::cout << "    >>> Scale weight " << groups[1].str() << " for " << groups[3].str() << " , "
0644                               << groups[4].str() << " , " << groups[5].str() << std::endl;
0645                   if (ismg26xNew) {
0646                     scaleVariationIDs.emplace_back(groups.str(4), groups.str(1), groups.str(3), groups.str(2));
0647                   } else {
0648                     scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
0649                   }
0650                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0651                   if (lheDebug)
0652                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0653                   if (!missed_weightgroup) {
0654                     break;
0655                   } else
0656                     missed_weightgroup = false;
0657                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0658                   if (lheDebug)
0659                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0660                                  "of the group."
0661                               << std::endl;
0662                   if (ismg26x || ismg26xNew)
0663                     missed_weightgroup = true;
0664                   --iLine;  // rewind by one, and go back to the outer loop
0665                   break;
0666                 }
0667               }
0668             } else if (groupname == "PDF_variation" || groupname.find("PDF_variation ") == 0) {
0669               if (lheDebug)
0670                 std::cout << ">>> Looks like a new-style block of PDF weights for one or more pdfs" << std::endl;
0671               for (++iLine; iLine < nLines; ++iLine) {
0672                 if (lheDebug)
0673                   std::cout << "    " << lines[iLine];
0674                 if (std::regex_search(lines[iLine], groups, pdfw)) {
0675                   unsigned int lhaID = std::stoi(groups.str(2));
0676                   if (lheDebug)
0677                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
0678                               << std::endl;
0679                   if (pdfSetWeightIDs.empty() || !pdfSetWeightIDs.back().maybe_add(groups.str(1), lhaID)) {
0680                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0681                   }
0682                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0683                   if (lheDebug)
0684                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0685                   if (!missed_weightgroup) {
0686                     break;
0687                   } else
0688                     missed_weightgroup = false;
0689                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0690                   if (lheDebug)
0691                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0692                                  "of the group."
0693                               << std::endl;
0694                   if (ismg26x || ismg26xNew)
0695                     missed_weightgroup = true;
0696                   --iLine;  // rewind by one, and go back to the outer loop
0697                   break;
0698                 }
0699               }
0700             } else if (groupname == "PDF_variation1" || groupname == "PDF_variation2") {
0701               if (lheDebug)
0702                 std::cout << ">>> Looks like a new-style block of PDF weights for multiple pdfs" << std::endl;
0703               unsigned int lastid = 0;
0704               for (++iLine; iLine < nLines; ++iLine) {
0705                 if (lheDebug)
0706                   std::cout << "    " << lines[iLine];
0707                 if (std::regex_search(lines[iLine], groups, pdfw)) {
0708                   unsigned int id = std::stoi(groups.str(1));
0709                   unsigned int lhaID = std::stoi(groups.str(2));
0710                   if (lheDebug)
0711                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
0712                               << std::endl;
0713                   if (id != (lastid + 1) || pdfSetWeightIDs.empty()) {
0714                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0715                   } else {
0716                     pdfSetWeightIDs.back().add(groups.str(1), lhaID);
0717                   }
0718                   lastid = id;
0719                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0720                   if (lheDebug)
0721                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0722                   if (!missed_weightgroup) {
0723                     break;
0724                   } else
0725                     missed_weightgroup = false;
0726                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0727                   if (lheDebug)
0728                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0729                                  "of the group."
0730                               << std::endl;
0731                   if (ismg26x || ismg26xNew)
0732                     missed_weightgroup = true;
0733                   --iLine;  // rewind by one, and go back to the outer loop
0734                   break;
0735                 }
0736               }
0737             } else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) {
0738               if (lheDebug)
0739                 std::cout << ">>> Looks like an old-style PDF weight for an individual pdf" << std::endl;
0740               unsigned int firstLhaID = lhaNameToID_.find(groupname)->second;
0741               bool first = true;
0742               for (++iLine; iLine < nLines; ++iLine) {
0743                 if (lheDebug)
0744                   std::cout << "    " << lines[iLine];
0745                 if (std::regex_search(
0746                         lines[iLine], groups, ismg26x ? pdfwmg26x : (ismg26xNew ? pdfwmg26xNew : pdfwOld))) {
0747                   unsigned int member = 0;
0748                   if (!ismg26x && !ismg26xNew) {
0749                     member = std::stoi(groups.str(2));
0750                   } else if (ismg26xNew) {
0751                     if (!groups.str(3).empty()) {
0752                       member = std::stoi(groups.str(3));
0753                     }
0754                   } else {
0755                     if (!groups.str(4).empty()) {
0756                       member = std::stoi(groups.str(4));
0757                     }
0758                   }
0759                   unsigned int lhaID = member + firstLhaID;
0760                   if (lheDebug)
0761                     std::cout << "    >>> PDF weight " << groups.str(1) << " for " << member << " = " << lhaID
0762                               << std::endl;
0763                   //if (member == 0) continue; // let's keep also the central value for now
0764                   if (first) {
0765                     pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
0766                     first = false;
0767                   } else {
0768                     pdfSetWeightIDs.back().add(groups.str(1), lhaID);
0769                   }
0770                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0771                   if (lheDebug)
0772                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0773                   if (!missed_weightgroup) {
0774                     break;
0775                   } else
0776                     missed_weightgroup = false;
0777                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0778                   if (lheDebug)
0779                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0780                                  "of the group."
0781                               << std::endl;
0782                   if (ismg26x || ismg26xNew)
0783                     missed_weightgroup = true;
0784                   --iLine;  // rewind by one, and go back to the outer loop
0785                   break;
0786                 }
0787               }
0788             } else if (groupname == "mass_variation" || groupname == "sthw2_variation" ||
0789                        groupname == "width_variation") {
0790               if (lheDebug)
0791                 std::cout << ">>> Looks like an EW parameter weight" << std::endl;
0792               for (++iLine; iLine < nLines; ++iLine) {
0793                 if (lheDebug)
0794                   std::cout << "    " << lines[iLine];
0795                 if (std::regex_search(lines[iLine], groups, rwgt)) {
0796                   std::string rwgtID = groups.str(1);
0797                   if (lheDebug)
0798                     std::cout << "    >>> LHE reweighting weight: " << rwgtID << std::endl;
0799                   if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
0800                     // we're only interested in the beggining of the block
0801                     lheReweighingIDs.emplace_back(rwgtID);
0802                   }
0803                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0804                   if (lheDebug)
0805                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0806                 }
0807               }
0808             } else {
0809               for (++iLine; iLine < nLines; ++iLine) {
0810                 if (lheDebug)
0811                   std::cout << "    " << lines[iLine];
0812                 if (std::regex_search(lines[iLine], groups, endweightgroup)) {
0813                   if (lheDebug)
0814                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0815                   if (!missed_weightgroup) {
0816                     break;
0817                   } else
0818                     missed_weightgroup = false;
0819                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0820                   if (lheDebug)
0821                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0822                                  "of the group."
0823                               << std::endl;
0824                   if (ismg26x || ismg26xNew)
0825                     missed_weightgroup = true;
0826                   --iLine;  // rewind by one, and go back to the outer loop
0827                   break;
0828                 }
0829               }
0830             }
0831           } else if (std::regex_search(lines[iLine], groups, weightgroupRwgt)) {
0832             std::string groupname = groups.str(1);
0833             if (groupname.find("mg_reweighting") != std::string::npos) {
0834               if (lheDebug)
0835                 std::cout << ">>> Looks like a LHE weights for reweighting" << std::endl;
0836               for (++iLine; iLine < nLines; ++iLine) {
0837                 if (lheDebug)
0838                   std::cout << "    " << lines[iLine];
0839                 if (std::regex_search(lines[iLine], groups, rwgt)) {
0840                   std::string rwgtID = groups.str(1);
0841                   if (lheDebug)
0842                     std::cout << "    >>> LHE reweighting weight: " << rwgtID << std::endl;
0843                   if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
0844                     // we're only interested in the beggining of the block
0845                     lheReweighingIDs.emplace_back(rwgtID);
0846                   }
0847                 } else if (std::regex_search(lines[iLine], endweightgroup)) {
0848                   if (lheDebug)
0849                     std::cout << ">>> Looks like the end of a weight group" << std::endl;
0850                   if (!missed_weightgroup) {
0851                     break;
0852                   } else
0853                     missed_weightgroup = false;
0854                 } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
0855                   if (lheDebug)
0856                     std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
0857                                  "of the group."
0858                               << std::endl;
0859                   if (ismg26x)
0860                     missed_weightgroup = true;
0861                   --iLine;  // rewind by one, and go back to the outer loop
0862                   break;
0863                 }
0864               }
0865             }
0866           }
0867         }
0868         //std::cout << "============= END [ " << iter->tag() << " ] ============ \n\n" << std::endl;
0869 
0870         // ----- SCALE VARIATIONS -----
0871         std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
0872         if (lheDebug)
0873           std::cout << "Found " << scaleVariationIDs.size() << " scale variations: " << std::endl;
0874         std::stringstream scaleDoc;
0875         scaleDoc << "LHE scale variation weights (w_var / w_nominal); ";
0876         for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
0877           const auto& sw = scaleVariationIDs[isw];
0878           if (isw)
0879             scaleDoc << "; ";
0880           scaleDoc << "[" << isw << "] is " << sw.label;
0881           weightChoice->scaleWeightIDs.push_back(sw.wid);
0882           if (lheDebug)
0883             printf("    id %s: scales ren = % .2f  fact = % .2f  text = %s\n",
0884                    sw.wid.c_str(),
0885                    sw.scales.first,
0886                    sw.scales.second,
0887                    sw.label.c_str());
0888         }
0889         if (!scaleVariationIDs.empty())
0890           weightChoice->scaleWeightsDoc = scaleDoc.str();
0891 
0892         // ------ PDF VARIATIONS (take the preferred one) -----
0893         if (lheDebug) {
0894           std::cout << "Found " << pdfSetWeightIDs.size() << " PDF set errors: " << std::endl;
0895           for (const auto& pw : pdfSetWeightIDs)
0896             printf("lhaIDs %6d - %6d (%3lu weights: %s, ... )\n",
0897                    pw.lhaIDs.first,
0898                    pw.lhaIDs.second,
0899                    pw.wids.size(),
0900                    pw.wids.front().c_str());
0901         }
0902 
0903         // ------ LHE REWEIGHTING -------
0904         if (lheDebug) {
0905           std::cout << "Found " << lheReweighingIDs.size() << " reweighting weights" << std::endl;
0906         }
0907         std::copy(lheReweighingIDs.begin(), lheReweighingIDs.end(), std::back_inserter(weightChoice->rwgtIDs));
0908 
0909         std::stringstream pdfDoc;
0910         pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA IDs ";
0911         bool found = false;
0912         for (const auto& pw : pdfSetWeightIDs) {
0913           for (uint32_t lhaid : preferredPDFLHAIDs_) {
0914             if (pw.lhaIDs.first != lhaid && pw.lhaIDs.first != (lhaid + 1))
0915               continue;  // sometimes the first weight is not saved if that PDF is the nominal one for the sample
0916             if (pw.wids.size() == 1)
0917               continue;  // only consider error sets
0918             pdfDoc << pw.lhaIDs.first << " - " << pw.lhaIDs.second;
0919             weightChoice->pdfWeightIDs = pw.wids;
0920             if (maxPdfWeights_ < pw.wids.size()) {
0921               weightChoice->pdfWeightIDs.resize(maxPdfWeights_);  // drop some replicas
0922               pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
0923             }
0924             weightChoice->pdfWeightsDoc = pdfDoc.str();
0925             found = true;
0926             break;
0927           }
0928           if (found)
0929             break;
0930         }
0931       }
0932     }
0933     return weightChoice;
0934   }
0935 
0936   // create an empty counter
0937   std::unique_ptr<LumiCacheInfoHolder> beginStream(edm::StreamID) const override {
0938     return std::make_unique<LumiCacheInfoHolder>();
0939   }
0940   // inizialize to zero at begin run
0941   void streamBeginRun(edm::StreamID id, edm::Run const&, edm::EventSetup const&) const override {
0942     streamCache(id)->clear();
0943   }
0944   void streamBeginLuminosityBlock(edm::StreamID id,
0945                                   edm::LuminosityBlock const& lumiBlock,
0946                                   edm::EventSetup const& eventSetup) const override {
0947     auto counterMap = &(streamCache(id)->countermap);
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     std::string label;
0955     if (genLumiInfoHead.isValid()) {
0956       label = genLumiInfoHead->configDescription();
0957       boost::replace_all(label, "-", "_");
0958       boost::replace_all(label, "/", "_");
0959     }
0960     counterMap->setLabel(label);
0961 
0962     if (genLumiInfoHead.isValid()) {
0963       auto weightChoice = &(streamCache(id)->weightChoice);
0964 
0965       std::vector<ScaleVarWeight> scaleVariationIDs;
0966       std::vector<PDFSetWeights> pdfSetWeightIDs;
0967       weightChoice->psWeightIDs.clear();
0968 
0969       std::regex scalew("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+)\\,\\s+mur=(\\S+)\\smuf=(\\S+)");
0970       std::regex pdfw("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+),\\s+Member\\s+(\\d+)\\s+of\\ssets\\s+(\\w+\\b)");
0971       std::regex mainPSw("sr(Def|:murfac=)(Hi|Lo|_dn|_up|0.5|2.0)");
0972       std::smatch groups;
0973       auto weightNames = genLumiInfoHead->weightNames();
0974       std::unordered_map<std::string, uint32_t> knownPDFSetsFromGenInfo_;
0975       unsigned int weightIter = 0;
0976       for (const auto& line : weightNames) {
0977         if (std::regex_search(line, groups, scalew)) {  // scale variation
0978           auto id = groups.str(1);
0979           auto group = groups.str(2);
0980           auto mur = groups.str(3);
0981           auto muf = groups.str(4);
0982           if (group.find("Central scale variation") != std::string::npos)
0983             scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
0984         } else if (std::regex_search(line, groups, pdfw)) {  // PDF variation
0985           auto id = groups.str(1);
0986           auto group = groups.str(2);
0987           auto memberid = groups.str(3);
0988           auto pdfset = groups.str(4);
0989           if (group.find(pdfset) != std::string::npos) {
0990             if (knownPDFSetsFromGenInfo_.find(pdfset) == knownPDFSetsFromGenInfo_.end()) {
0991               knownPDFSetsFromGenInfo_[pdfset] = std::atoi(id.c_str());
0992               pdfSetWeightIDs.emplace_back(id, std::atoi(id.c_str()));
0993             } else
0994               pdfSetWeightIDs.back().add(id, std::atoi(id.c_str()));
0995           }
0996         } else if (line == "Baseline") {
0997           weightChoice->psBaselineID = weightIter;
0998         } else if (line.find("isr") != std::string::npos || line.find("fsr") != std::string::npos) {
0999           weightChoice->matchPS_alt = line.find("sr:") != std::string::npos;  // (f/i)sr: for new weights
1000           if (keepAllPSWeights_) {
1001             weightChoice->psWeightIDs.push_back(weightIter);  // PS variations
1002           } else if (std::regex_search(line, groups, mainPSw)) {
1003             if (weightChoice->psWeightIDs.empty())
1004               weightChoice->psWeightIDs = std::vector<unsigned int>(4, -1);
1005             int psIdx = (line.find("fsr") != std::string::npos) ? 1 : 0;
1006             psIdx += (groups.str(2) == "Hi" || groups.str(2) == "_up" || groups.str(2) == "2.0") ? 0 : 2;
1007             weightChoice->psWeightIDs[psIdx] = weightIter;
1008           }
1009         }
1010         weightIter++;
1011       }
1012       if (keepAllPSWeights_) {
1013         weightChoice->psWeightsDoc = "All PS weights (w_var / w_nominal)";
1014       } else if (weightChoice->psWeightIDs.size() == 4) {
1015         weightChoice->psWeightsDoc =
1016             "PS weights (w_var / w_nominal);   [0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2"
1017             "[2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;";
1018         for (int i = 0; i < 4; i++) {
1019           if (static_cast<int>(weightChoice->psWeightIDs[i]) == -1)
1020             weightChoice->setMissingWeight(i);
1021         }
1022       } else {
1023         weightChoice->psWeightsDoc = "dummy PS weight (1.0) ";
1024       }
1025 
1026       weightChoice->scaleWeightIDs.clear();
1027       weightChoice->pdfWeightIDs.clear();
1028 
1029       std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
1030       std::stringstream scaleDoc;
1031       scaleDoc << "LHE scale variation weights (w_var / w_nominal); ";
1032       for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
1033         const auto& sw = scaleVariationIDs[isw];
1034         if (isw)
1035           scaleDoc << "; ";
1036         scaleDoc << "[" << isw << "] is " << sw.label;
1037         weightChoice->scaleWeightIDs.push_back(std::atoi(sw.wid.c_str()));
1038       }
1039       if (!scaleVariationIDs.empty())
1040         weightChoice->scaleWeightsDoc = scaleDoc.str();
1041       std::stringstream pdfDoc;
1042       pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA names ";
1043       bool found = false;
1044       for (const auto& pw : pdfSetWeightIDs) {
1045         if (pw.wids.size() == 1)
1046           continue;  // only consider error sets
1047         for (const auto& wantedpdf : lhaNameToID_) {
1048           auto pdfname = wantedpdf.first;
1049           if (knownPDFSetsFromGenInfo_.find(pdfname) == knownPDFSetsFromGenInfo_.end())
1050             continue;
1051           uint32_t lhaid = knownPDFSetsFromGenInfo_.at(pdfname);
1052           if (pw.lhaIDs.first != lhaid)
1053             continue;
1054           pdfDoc << pdfname;
1055           for (const auto& x : pw.wids)
1056             weightChoice->pdfWeightIDs.push_back(std::atoi(x.c_str()));
1057           if (maxPdfWeights_ < pw.wids.size()) {
1058             weightChoice->pdfWeightIDs.resize(maxPdfWeights_);  // drop some replicas
1059             pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
1060           }
1061           weightChoice->pdfWeightsDoc = pdfDoc.str();
1062           found = true;
1063           break;
1064         }
1065         if (found)
1066           break;
1067       }
1068     }
1069   }
1070   // create an empty counter
1071   std::shared_ptr<CounterMap> globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override {
1072     return std::make_shared<CounterMap>();
1073   }
1074   // add this stream to the summary
1075   void streamEndRunSummary(edm::StreamID id,
1076                            edm::Run const&,
1077                            edm::EventSetup const&,
1078                            CounterMap* runCounterMap) const override {
1079     runCounterMap->merge(streamCache(id)->countermap);
1080   }
1081   // nothing to do per se
1082   void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {}
1083   // write the total to the run
1084   void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, CounterMap const* runCounterMap) const override {
1085     auto out = std::make_unique<nanoaod::MergeableCounterTable>();
1086 
1087     for (const auto& x : runCounterMap->countermap) {
1088       auto runCounter = &(x.second);
1089       std::string label = (!x.first.empty()) ? (std::string("_") + x.first) : "";
1090       std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : "";
1091 
1092       out->addInt("genEventCount" + label, "event count" + doclabel, runCounter->num);
1093       out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter->sumw);
1094       out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter->sumw2);
1095 
1096       double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1;
1097       auto sumScales = runCounter->sumScale;
1098       for (auto& val : sumScales)
1099         val *= norm;
1100       out->addVFloatWithNorm("LHEScaleSumw" + label,
1101                              "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel,
1102                              sumScales,
1103                              runCounter->sumw);
1104       auto sumPDFs = runCounter->sumPDF;
1105       for (auto& val : sumPDFs)
1106         val *= norm;
1107       out->addVFloatWithNorm("LHEPdfSumw" + label,
1108                              "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel,
1109                              sumPDFs,
1110                              runCounter->sumw);
1111       auto sumPS = runCounter->sumPS;
1112       for (auto& val : sumPS)
1113         val *= norm;
1114       out->addVFloatWithNorm("PSSumw" + label,
1115                              "Sum of genEventWeight * PSWeight[i], divided by genEventSumw" + doclabel,
1116                              sumPS,
1117                              runCounter->sumw);
1118       if (!runCounter->sumRwgt.empty()) {
1119         auto sumRwgts = runCounter->sumRwgt;
1120         for (auto& val : sumRwgts)
1121           val *= norm;
1122         out->addVFloatWithNorm("LHEReweightingSumw" + label,
1123                                "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel,
1124                                sumRwgts,
1125                                runCounter->sumw);
1126       }
1127       if (!runCounter->sumNamed.empty()) {  // it could be empty if there's no LHE info in the sample
1128         for (unsigned int i = 0, n = namedWeightLabels_.size(); i < n; ++i) {
1129           out->addFloatWithNorm(
1130               "LHESumw_" + namedWeightLabels_[i] + label,
1131               "Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[i] + ", divided by genEventSumw" + doclabel,
1132               runCounter->sumNamed[i] * norm,
1133               runCounter->sumw);
1134         }
1135       }
1136     }
1137     iRun.put(std::move(out));
1138   }
1139   // nothing to do here
1140   void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {}
1141 
1142   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1143     edm::ParameterSetDescription desc;
1144     desc.add<edm::InputTag>("genEvent", edm::InputTag("generator"))
1145         ->setComment("tag for the GenEventInfoProduct, to get the main weight");
1146     desc.add<edm::InputTag>("genLumiInfoHeader", edm::InputTag("generator"))
1147         ->setComment("tag for the GenLumiInfoProduct, to get the model string");
1148     desc.add<std::vector<edm::InputTag>>("lheInfo", std::vector<edm::InputTag>{{"externalLHEProducer"}, {"source"}})
1149         ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
1150 
1151     edm::ParameterSetDescription prefpdf;
1152     prefpdf.add<std::string>("name");
1153     prefpdf.add<uint32_t>("lhaid");
1154     desc.addVPSet("preferredPDFs", prefpdf, std::vector<edm::ParameterSet>())
1155         ->setComment(
1156             "LHA PDF Ids of the preferred PDF sets, in order of preference (the first matching one will be used)");
1157     desc.add<std::vector<std::string>>("namedWeightIDs")->setComment("set of LHA weight IDs for named LHE weights");
1158     desc.add<std::vector<std::string>>("namedWeightLabels")
1159         ->setComment("output names for the namedWeightIDs (in the same order)");
1160     desc.add<int32_t>("lheWeightPrecision")->setComment("Number of bits in the mantissa for LHE weights");
1161     desc.add<uint32_t>("maxPdfWeights")->setComment("Maximum number of PDF weights to save (to crop NN replicas)");
1162     desc.add<bool>("keepAllPSWeights")->setComment("Store all PS weights found");
1163     desc.addOptionalUntracked<bool>("debug")->setComment("dump out all LHE information for one event");
1164     descriptions.add("genWeightsTable", desc);
1165   }
1166 
1167 protected:
1168   const edm::EDGetTokenT<GenEventInfoProduct> genTag_;
1169   const std::vector<edm::InputTag> lheLabel_;
1170   const std::vector<edm::EDGetTokenT<LHEEventProduct>> lheTag_;
1171   const std::vector<edm::EDGetTokenT<LHERunInfoProduct>> lheRunTag_;
1172   const edm::EDGetTokenT<GenLumiInfoHeader> genLumiInfoHeadTag_;
1173 
1174   std::vector<uint32_t> preferredPDFLHAIDs_;
1175   std::unordered_map<std::string, uint32_t> lhaNameToID_;
1176   std::vector<std::string> namedWeightIDs_;
1177   std::vector<std::string> namedWeightLabels_;
1178   int lheWeightPrecision_;
1179   unsigned int maxPdfWeights_;
1180   bool keepAllPSWeights_;
1181 
1182   mutable std::atomic<bool> debug_, debugRun_, hasIssuedWarning_, psWeightWarning_;
1183 };
1184 
1185 #include "FWCore/Framework/interface/MakerMacros.h"
1186 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!