15 #include <unordered_map> 23 num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumNamed() {}
29 std::vector<long double> sumPDF, sumScale, sumNamed;
32 num = 0; sumw = 0; sumw2 = 0;
33 sumPDF.clear(); sumScale.clear(); sumNamed.clear();
37 void incGenOnly(
double w) {
38 num++; sumw +=
w; sumw2 += (w*
w);
40 void incLHE(
double w0,
const std::vector<double> & wScale,
const std::vector<double> & wPDF,
const std::vector<double> & wNamed) {
44 if (!wScale.empty()) {
45 if (sumScale.empty()) sumScale.resize(wScale.size(), 0);
46 for (
unsigned int i = 0,
n = wScale.size();
i <
n; ++
i) sumScale[
i] += (w0 * wScale[
i]);
49 if (sumPDF.empty()) sumPDF.resize(wPDF.size(), 0);
50 for (
unsigned int i = 0, n = wPDF.size();
i <
n; ++
i) sumPDF[
i] += (w0 * wPDF[
i]);
52 if (!wNamed.empty()) {
53 if (sumNamed.empty()) sumNamed.resize(wNamed.size(), 0);
54 for (
unsigned int i = 0, n = wNamed.size();
i <
n; ++
i) sumNamed[
i] += (w0 * wNamed[
i]);
59 num += other.num; sumw += other.sumw; sumw2 += other.sumw2;
60 if (sumScale.empty() && !other.sumScale.empty()) sumScale.resize(other.sumScale.size(),0);
61 if (sumPDF.empty() && !other.sumPDF.empty()) sumPDF.resize(other.sumPDF.size(),0);
62 if (sumNamed.empty() && !other.sumNamed.empty()) sumNamed.resize(other.sumNamed.size(),0);
63 for (
unsigned int i = 0, n = sumScale.size();
i <
n; ++
i) sumScale[
i] += other.sumScale[
i];
64 for (
unsigned int i = 0, n = sumPDF.size();
i <
n; ++
i) sumPDF[
i] += other.sumPDF[
i];
65 for (
unsigned int i = 0, n = sumNamed.size();
i <
n; ++
i) sumNamed[
i] += other.sumNamed[
i];
70 struct DynamicWeightChoice {
73 std::vector<std::string> scaleWeightIDs;
76 std::vector<std::string> pdfWeightIDs;
82 if (match != std::string::npos) {
85 return std::stof(pre) *
std::pow(10.0
f, std::stof(post));
87 return std::stof(str);
91 struct ScaleVarWeight {
93 std::pair<float,float> scales;
95 wid(id), label(text), scales(stof_fortrancomp(muR), stof_fortrancomp(muF)) {}
96 bool operator<(
const ScaleVarWeight & other) {
return (scales == other.scales ? wid < other.wid : scales < other.scales); }
98 struct PDFSetWeights {
99 std::vector<std::string> wids;
100 std::pair<unsigned int,unsigned int> lhaIDs;
101 PDFSetWeights(
const std::string & wid,
unsigned int lhaID) : wids(1,wid), lhaIDs(lhaID,lhaID) {}
102 bool operator<(
const PDFSetWeights & other)
const {
return lhaIDs < other.lhaIDs; }
105 lhaIDs.second = lhaID;
107 bool maybe_add(
const std::string & wid,
unsigned int lhaID) {
108 if (lhaID == lhaIDs.second+1) {
123 lheLabel_(params.getParameter<
edm::InputTag>(
"lheInfo")),
126 namedWeightIDs_(params.getParameter<
std::vector<
std::
string>>(
"namedWeightIDs")),
127 namedWeightLabels_(params.getParameter<
std::vector<
std::
string>>(
"namedWeightLabels")),
128 lheWeightPrecision_(params.getParameter<int32_t>(
"lheWeightPrecision")),
129 maxPdfWeights_(params.getParameter<uint32_t>(
"maxPdfWeights")),
130 debug_(params.getUntrackedParameter<
bool>(
"debug",
false)), debugRun_(debug_.
load()),
131 hasIssuedWarning_(
false)
133 produces<nanoaod::FlatTable>();
134 produces<nanoaod::FlatTable>(
"LHEScale");
135 produces<nanoaod::FlatTable>(
"LHEPdf");
136 produces<nanoaod::FlatTable>(
"LHENamed");
137 produces<nanoaod::MergeableCounterTable,edm::Transition::EndRun>();
138 if (namedWeightIDs_.size() != namedWeightLabels_.size()) {
139 throw cms::Exception(
"Configuration",
"Size mismatch between namedWeightIDs & namedWeightLabels");
143 uint32_t
lhaid = pdfps.getParameter<uint32_t>(
"lhaid");
144 preferredPDFLHAIDs_.push_back(lhaid);
146 lhaNameToID_[name+
".LHgrid"] =
lhaid;
162 auto out = std::make_unique<nanoaod::FlatTable>(1,
"genWeight",
true);
163 out->setDoc(
"generator weight");
168 std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheNamedTab;
173 const DynamicWeightChoice * weightChoice = runCache(iEvent.
getRun().
index());
175 fillLHEWeightTables(counter, weightChoice, weight, *lheInfo, lheScaleTab, lhePdfTab, lheNamedTab);
178 counter->incGenOnly(weight);
183 if (!hasIssuedWarning_.exchange(
true)) {
184 edm::LogWarning(
"LHETablesProducer") <<
"No LHEEventProduct, so there will be no LHE Tables\n";
195 const DynamicWeightChoice * weightChoice,
198 std::unique_ptr<nanoaod::FlatTable> & outScale,
199 std::unique_ptr<nanoaod::FlatTable> & outPdf,
200 std::unique_ptr<nanoaod::FlatTable> & outNamed )
const 202 bool lheDebug = debug_.exchange(
false);
204 const std::vector<std::string> & scaleWeightIDs = weightChoice->scaleWeightIDs;
205 const std::vector<std::string> & pdfWeightIDs = weightChoice->pdfWeightIDs;
209 std::vector<double> wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wNamed(namedWeightIDs_.size(), 1);
211 if (lheDebug) printf(
"Weight %+9.5f rel %+9.5f for id %s\n",
weight.wgt,
weight.wgt/w0,
weight.id.c_str());
213 auto mScale =
std::find(scaleWeightIDs.begin(), scaleWeightIDs.end(),
weight.id);
214 if (mScale != scaleWeightIDs.end()) wScale[mScale-scaleWeightIDs.begin()] =
weight.wgt/w0;
216 auto mPDF =
std::find(pdfWeightIDs.begin(), pdfWeightIDs.end(),
weight.id);
217 if (mPDF != pdfWeightIDs.end()) wPDF[mPDF-pdfWeightIDs.begin()] =
weight.wgt/w0;
219 auto mNamed =
std::find(namedWeightIDs_.begin(), namedWeightIDs_.end(),
weight.id);
220 if (mNamed != namedWeightIDs_.end()) wNamed[mNamed-namedWeightIDs_.begin()] =
weight.wgt/w0;
231 for (
unsigned int i = 0, n = wNamed.size();
i <
n; ++
i) {
232 outNamed->addColumnValue<
float>(namedWeightLabels_[
i], wNamed[
i],
"LHE weight for id "+namedWeightIDs_[
i]+
", relative to nominal",
nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
235 counter->incLHE(genWeight, wScale, wPDF, wNamed);
242 bool lheDebug = debugRun_.exchange(
false);
243 auto weightChoice = std::make_shared<DynamicWeightChoice>();
248 std::vector<ScaleVarWeight> scaleVariationIDs;
249 std::vector<PDFSetWeights> pdfSetWeightIDs;
251 std::regex weightgroup(
"<weightgroup\\s+combine=\"(.*)\"\\s+(?:name|type)=\"(.*)\"\\s*>");
252 std::regex endweightgroup(
"</weightgroup>");
253 std::regex scalew(
"<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:mu[rR]|renscfact)=(\\S+)\\s+(?:mu[Ff]|facscfact)=(\\S+)(\\s+.*)?)</weight>");
254 std::regex pdfw(
"<weight\\s+id=\"(\\d+)\">\\s*(?:PDF set|lhapdf|PDF|pdfset)\\s*=\\s*(\\d+)\\s*(?:\\s.*)?</weight>");
255 std::regex pdfwOld(
"<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*Member \\s*(\\d+)\\s*(?:.*)</weight>");
258 if (iter->tag() !=
"initrwgt") {
259 if (lheDebug)
std::cout <<
"Skipping LHE header with tag" << iter->tag() << std::endl;
262 if (lheDebug)
std::cout <<
"Found LHE header with tag" << iter->tag() << std::endl;
263 const std::vector<std::string> &
lines = iter->lines();
264 for (
unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) {
266 if (std::regex_search(lines[iLine], groups, weightgroup)) {
268 if (lheDebug)
std::cout <<
">>> Looks like the beginning of a weight group for '" << groupname <<
"'" << std::endl;
269 if (groupname.find(
"scale_variation") == 0 || groupname ==
"Central scale variation") {
270 if (lheDebug)
std::cout <<
">>> Looks like scale variation for theory uncertainties" << std::endl;
271 for ( ++iLine; iLine < nLines; ++iLine) {
272 if (lheDebug)
std::cout <<
" " << lines[iLine];
273 if (std::regex_search(lines[iLine], groups, scalew)) {
274 if (lheDebug)
std::cout <<
" >>> Scale weight " << groups[1].str() <<
" for " << groups[3].str() <<
" , " << groups[4].str() <<
" , " << groups[5].str() << std::endl;
275 scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
276 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
277 if (lheDebug)
std::cout <<
">>> Looks like the end of a weight group" << std::endl;
279 }
else if (std::regex_search(lines[iLine], weightgroup)) {
280 if (lheDebug)
std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end of the group." << std::endl;
285 }
else if (groupname ==
"PDF_variation" || groupname.find(
"PDF_variation ") == 0) {
286 if (lheDebug)
std::cout <<
">>> Looks like a new-style block of PDF weights for one or more pdfs" << std::endl;
287 for ( ++iLine; iLine < nLines; ++iLine) {
288 if (lheDebug)
std::cout <<
" " << lines[iLine];
289 if (std::regex_search(lines[iLine], groups, pdfw)) {
290 unsigned int lhaID = std::stoi(groups.str(2));
291 if (lheDebug)
std::cout <<
" >>> PDF weight " << groups.str(1) <<
" for " << groups.str(2) <<
" = " << lhaID << std::endl;
292 if (pdfSetWeightIDs.empty() || ! pdfSetWeightIDs.back().maybe_add(groups.str(1),lhaID)) {
293 pdfSetWeightIDs.emplace_back(groups.str(1),lhaID);
295 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
296 if (lheDebug)
std::cout <<
">>> Looks like the end of a weight group" << std::endl;
298 }
else if (std::regex_search(lines[iLine], weightgroup)) {
299 if (lheDebug)
std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end of the group." << std::endl;
304 }
else if (groupname ==
"PDF_variation1" || groupname ==
"PDF_variation2") {
305 if (lheDebug)
std::cout <<
">>> Looks like a new-style block of PDF weights for multiple pdfs" << std::endl;
306 unsigned int lastid = 0;
307 for ( ++iLine; iLine < nLines; ++iLine) {
308 if (lheDebug)
std::cout <<
" " << lines[iLine];
309 if (std::regex_search(lines[iLine], groups, pdfw)) {
310 unsigned int id = std::stoi(groups.str(1));
311 unsigned int lhaID = std::stoi(groups.str(2));
312 if (lheDebug)
std::cout <<
" >>> PDF weight " << groups.str(1) <<
" for " << groups.str(2) <<
" = " << lhaID << std::endl;
313 if (
id != (lastid+1) || pdfSetWeightIDs.empty()) {
314 pdfSetWeightIDs.emplace_back(groups.str(1),lhaID);
316 pdfSetWeightIDs.back().add(groups.str(1),lhaID);
319 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
320 if (lheDebug)
std::cout <<
">>> Looks like the end of a weight group" << std::endl;
322 }
else if (std::regex_search(lines[iLine], weightgroup)) {
323 if (lheDebug)
std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end of the group." << std::endl;
328 }
else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) {
329 if (lheDebug)
std::cout <<
">>> Looks like an old-style PDF weight for an individual pdf" << std::endl;
330 unsigned int firstLhaID = lhaNameToID_.find(groupname)->second;
332 for ( ++iLine; iLine < nLines; ++iLine) {
333 if (lheDebug)
std::cout <<
" " << lines[iLine];
334 if (std::regex_search(lines[iLine], groups, pdfwOld)) {
335 unsigned int member = std::stoi(groups.str(2));
336 unsigned int lhaID = member+firstLhaID;
337 if (lheDebug)
std::cout <<
" >>> PDF weight " << groups.str(1) <<
" for " << groups.str(2) <<
" = " << lhaID << std::endl;
340 pdfSetWeightIDs.emplace_back(groups.str(1),lhaID);
343 pdfSetWeightIDs.back().add(groups.str(1),lhaID);
345 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
346 if (lheDebug)
std::cout <<
">>> Looks like the end of a weight group" << std::endl;
348 }
else if (std::regex_search(lines[iLine], weightgroup)) {
349 if (lheDebug)
std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end of the group." << std::endl;
355 for ( ++iLine; iLine < nLines; ++iLine) {
356 if (lheDebug)
std::cout <<
" " << lines[iLine];
357 if (std::regex_search(lines[iLine], groups, endweightgroup)) {
358 if (lheDebug)
std::cout <<
">>> Looks like the end of a weight group" << std::endl;
360 }
else if (std::regex_search(lines[iLine], weightgroup)) {
361 if (lheDebug)
std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end of the group." << std::endl;
372 std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
373 if (lheDebug)
std::cout <<
"Found " << scaleVariationIDs.size() <<
" scale variations: " << std::endl;
374 std::stringstream scaleDoc; scaleDoc <<
"LHE scale variation weights (w_var / w_nominal); ";
375 for (
unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
376 const auto & sw = scaleVariationIDs[isw];
377 if (isw) scaleDoc <<
"; ";
378 scaleDoc <<
"[" << isw <<
"] is " << sw.label;
379 weightChoice->scaleWeightIDs.push_back(sw.wid);
380 if (lheDebug) printf(
" id %s: scales ren = % .2f fact = % .2f text = %s\n", sw.wid.c_str(), sw.scales.first, sw.scales.second, sw.label.c_str());
382 if (!scaleVariationIDs.empty()) weightChoice->scaleWeightsDoc = scaleDoc.str();
386 std::cout <<
"Found " << pdfSetWeightIDs.size() <<
" PDF set errors: " << std::endl;
387 for (
const auto & pw : pdfSetWeightIDs) printf(
"lhaIDs %6d - %6d (%3lu weights: %s, ... )\n", pw.lhaIDs.first, pw.lhaIDs.second, pw.wids.size(), pw.wids.front().c_str());
390 std::stringstream pdfDoc; pdfDoc <<
"LHE pdf variation weights (w_var / w_nominal) for LHA IDs ";
392 for (uint32_t
lhaid : preferredPDFLHAIDs_) {
393 for (
const auto & pw : pdfSetWeightIDs) {
394 if (pw.lhaIDs.first !=
lhaid && pw.lhaIDs.first != (
lhaid+1))
continue;
395 if (pw.wids.size() == 1)
continue;
396 pdfDoc << pw.lhaIDs.first <<
" - " << pw.lhaIDs.second;
397 weightChoice->pdfWeightIDs = pw.wids;
398 if (maxPdfWeights_ < pw.wids.size()) {
399 weightChoice->pdfWeightIDs.resize(maxPdfWeights_);
400 pdfDoc <<
", truncated to the first " << maxPdfWeights_ <<
" replicas";
402 weightChoice->pdfWeightsDoc = pdfDoc.str();
415 return std::make_unique<Counter>();
419 streamCache(
id)->clear();
423 return std::make_shared<Counter>();
427 runCounter->merge(*streamCache(
id));
434 auto out = std::make_unique<nanoaod::MergeableCounterTable>();
435 out->addInt(
"genEventCount",
"event count", runCounter->num);
436 out->addFloat(
"genEventSumw",
"sum of gen weights", runCounter->sumw);
437 out->addFloat(
"genEventSumw2",
"sum of gen (weight^2)", runCounter->sumw2);
439 double norm = runCounter->sumw ? 1.0/runCounter->sumw : 1;
440 auto sumScales = runCounter->sumScale;
for (
auto &
val : sumScales)
val *= norm;
441 out->addVFloat(
"LHEScaleSumw",
"Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw", sumScales);
442 auto sumPDFs = runCounter->sumPDF;
for (
auto &
val : sumPDFs)
val *= norm;
443 out->addVFloat(
"LHEPdfSumw",
"Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw", sumPDFs);
444 if (!runCounter->sumNamed.empty()) {
445 for (
unsigned int i = 0, n = namedWeightLabels_.size();
i <
n; ++
i) {
446 out->addFloat(
"LHESumw_"+namedWeightLabels_[
i],
"Sum of genEventWeight * LHEWeight_"+namedWeightLabels_[i]+
", divided by genEventSumw", runCounter->sumNamed[i] * norm);
457 desc.
add<
edm::InputTag>(
"lheInfo",
edm::InputTag(
"externalLHEProducer"))->setComment(
"tag for the LHE information (LHEEventProduct and LHERunInfoProduct)");
461 prefpdf.
add<uint32_t>(
"lhaid");
462 desc.
addVPSet(
"preferredPDFs", prefpdf, std::vector<edm::ParameterSet>())->
setComment(
"LHA PDF Ids of the preferred PDF sets, in order of preference (the first matching one will be used)");
463 desc.
add<std::vector<std::string>>(
"namedWeightIDs")->setComment(
"set of LHA weight IDs for named LHE weights");
464 desc.
add<std::vector<std::string>>(
"namedWeightLabels")->setComment(
"output names for the namedWeightIDs (in the same order)");
465 desc.
add<int32_t>(
"lheWeightPrecision")->setComment(
"Number of bits in the mantissa for LHE weights");
466 desc.
add<uint32_t>(
"maxPdfWeights")->setComment(
"Maximum number of PDF weights to save (to crop NN replicas)");
467 desc.
addOptionalUntracked<
bool>(
"debug")->setComment(
"dump out all LHE information for one event");
468 descriptions.
add(
"genWeightsTable", desc);
double originalXWGTUP() const
T getParameter(std::string const &) const
bool getByLabel(std::string const &label, Handle< PROD > &result) const
void setComment(std::string const &value)
void streamBeginRun(edm::StreamID id, edm::Run const &, edm::EventSetup const &) const override
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
ParameterDescriptionBase * addVPSet(U const &iLabel, ParameterSetDescription const &validator, std::vector< ParameterSet > const &defaults)
std::unordered_map< std::string, uint32_t > lhaNameToID_
void globalEndRun(edm::Run const &, edm::EventSetup const &) const override
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
void globalEndRunProduce(edm::Run &iRun, edm::EventSetup const &, Counter const *runCounter) const override
#define DEFINE_FWK_MODULE(type)
void streamEndRunSummary(edm::StreamID id, edm::Run const &, edm::EventSetup const &, Counter *runCounter) const override
const edm::EDGetTokenT< GenEventInfoProduct > genTag_
std::shared_ptr< DynamicWeightChoice > globalBeginRun(edm::Run const &iRun, edm::EventSetup const &) const override
headers_const_iterator headers_end() const
Run const & getRun() const
std::atomic< bool > hasIssuedWarning_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
std::unique_ptr< Counter > beginStream(edm::StreamID) const override
const std::vector< WGT > & weights() const
std::function< unsigned int(align::ID)> Counter
bool operator<(const FedChannelConnection &, const FedChannelConnection &)
headers_const_iterator headers_begin() const
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
const edm::EDGetTokenT< LHEEventProduct > lheTag_
const edm::InputTag lheLabel_
ParameterDescriptionBase * add(U const &iLabel, T const &value)
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
std::vector< std::string > namedWeightLabels_
~GenWeightsTableProducer() override
GenWeightsTableProducer(edm::ParameterSet const ¶ms)
std::shared_ptr< Counter > globalBeginRunSummary(edm::Run const &, edm::EventSetup const &) const override
void put(std::unique_ptr< PROD > product)
Put a new product.
void add(std::string const &label, ParameterSetDescription const &psetDescription)
unsigned int maxPdfWeights_
const edm::EDGetTokenT< LHERunInfoProduct > lheRunTag_
std::vector< uint32_t > preferredPDFLHAIDs_
void produce(edm::StreamID id, edm::Event &iEvent, const edm::EventSetup &iSetup) const override
ParameterDescriptionBase * addOptionalUntracked(U const &iLabel, T const &value)
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
std::vector< std::string > namedWeightIDs_
void globalEndRunSummary(edm::Run const &, edm::EventSetup const &, Counter *runCounter) const override
Power< A, B >::type pow(const A &a, const B &b)
void fillLHEWeightTables(Counter *counter, const DynamicWeightChoice *weightChoice, double genWeight, const LHEEventProduct &lheProd, std::unique_ptr< nanoaod::FlatTable > &outScale, std::unique_ptr< nanoaod::FlatTable > &outPdf, std::unique_ptr< nanoaod::FlatTable > &outNamed) const
def merge(dictlist, TELL=False)