15 #include "boost/algorithm/string.hpp" 18 #include <unordered_map> 25 Counter() :
num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumRwgt(), sumNamed(), sumPS() {}
31 std::vector<long double> sumPDF, sumScale, sumRwgt, sumNamed, sumPS;
40 sumNamed.clear(), sumPS.clear();
44 void incGenOnly(
double w) {
50 void incPSOnly(
double w0,
const std::vector<double>& wPS) {
53 sumPS.resize(wPS.size(), 0);
54 for (
unsigned int i = 0,
n = wPS.size();
i <
n; ++
i)
55 sumPS[
i] += (w0 * wPS[
i]);
59 void incLHE(
double w0,
60 const std::vector<double>& wScale,
61 const std::vector<double>& wPDF,
62 const std::vector<double>& wRwgt,
63 const std::vector<double>& wNamed,
64 const std::vector<double>& wPS) {
68 if (!wScale.empty()) {
70 sumScale.resize(wScale.size(), 0);
71 for (
unsigned int i = 0,
n = wScale.size();
i <
n; ++
i)
72 sumScale[
i] += (w0 * wScale[
i]);
76 sumPDF.resize(wPDF.size(), 0);
77 for (
unsigned int i = 0, n = wPDF.size();
i <
n; ++
i)
78 sumPDF[
i] += (w0 * wPDF[
i]);
82 sumRwgt.resize(wRwgt.size(), 0);
83 for (
unsigned int i = 0, n = wRwgt.size();
i <
n; ++
i)
84 sumRwgt[
i] += (w0 * wRwgt[
i]);
86 if (!wNamed.empty()) {
88 sumNamed.resize(wNamed.size(), 0);
89 for (
unsigned int i = 0, n = wNamed.size();
i <
n; ++
i)
90 sumNamed[
i] += (w0 * wNamed[
i]);
99 if (sumScale.empty() && !other.sumScale.empty())
100 sumScale.resize(other.sumScale.size(), 0);
101 if (sumPDF.empty() && !other.sumPDF.empty())
102 sumPDF.resize(other.sumPDF.size(), 0);
103 if (sumRwgt.empty() && !other.sumRwgt.empty())
104 sumRwgt.resize(other.sumRwgt.size(), 0);
105 if (sumNamed.empty() && !other.sumNamed.empty())
106 sumNamed.resize(other.sumNamed.size(), 0);
107 if (sumPS.empty() && !other.sumPS.empty())
108 sumPS.resize(other.sumPS.size(), 0);
109 if (!other.sumScale.empty())
110 for (
unsigned int i = 0, n = sumScale.size();
i <
n; ++
i)
111 sumScale[
i] += other.sumScale[
i];
112 if (!other.sumPDF.empty())
113 for (
unsigned int i = 0, n = sumPDF.size();
i <
n; ++
i)
114 sumPDF[
i] += other.sumPDF[
i];
115 if (!other.sumRwgt.empty())
116 for (
unsigned int i = 0, n = sumRwgt.size();
i <
n; ++
i)
117 sumRwgt[
i] += other.sumRwgt[
i];
118 if (!other.sumNamed.empty())
119 for (
unsigned int i = 0, n = sumNamed.size();
i <
n; ++
i)
120 sumNamed[
i] += other.sumNamed[
i];
121 if (!other.sumPS.empty())
122 for (
unsigned int i = 0, n = sumPS.size();
i <
n; ++
i)
123 sumPS[
i] += other.sumPS[
i];
128 std::map<std::string, Counter> countermap;
131 void merge(
const CounterMap& other) {
132 for (
const auto&
y : other.countermap)
133 countermap[
y.first].merge(
y.second);
137 for (
auto x : countermap)
143 active_el = &(countermap[
label]);
144 active_label =
label;
146 void checkLabelSet() {
148 throw cms::Exception(
"LogicError",
"Called CounterMap::get() before setting the active label\n");
161 struct DynamicWeightChoice {
164 std::vector<std::string> scaleWeightIDs;
167 std::vector<std::string> pdfWeightIDs;
170 std::vector<std::string> rwgtIDs;
174 struct DynamicWeightChoiceGenInfo {
177 std::vector<unsigned int> scaleWeightIDs;
180 std::vector<unsigned int> pdfWeightIDs;
183 std::vector<unsigned int> psWeightIDs;
186 struct LumiCacheInfoHolder {
187 CounterMap countermap;
188 DynamicWeightChoiceGenInfo weightChoice;
191 weightChoice = DynamicWeightChoiceGenInfo();
197 if (match != std::string::npos) {
200 return std::stof(pre) *
std::pow(10.0
f, std::stof(post));
202 return std::stof(str);
206 struct ScaleVarWeight {
208 std::pair<float, float> scales;
210 : wid(id), label(text), scales(stof_fortrancomp(muR), stof_fortrancomp(muF)) {}
211 bool operator<(
const ScaleVarWeight& other) {
212 return (scales == other.scales ? wid < other.wid : scales < other.scales);
215 struct PDFSetWeights {
216 std::vector<std::string> wids;
217 std::pair<unsigned int, unsigned int> lhaIDs;
218 PDFSetWeights(
const std::string& wid,
unsigned int lhaID) : wids(1, wid), lhaIDs(lhaID, lhaID) {}
219 bool operator<(
const PDFSetWeights& other)
const {
return lhaIDs < other.lhaIDs; }
222 lhaIDs.second = lhaID;
224 bool maybe_add(
const std::string& wid,
unsigned int lhaID) {
225 if (lhaID == lhaIDs.second + 1) {
237 edm::RunCache<DynamicWeightChoice>,
238 edm::RunSummaryCache<CounterMap>,
239 edm::EndRunProducer> {
243 lheLabel_(params.getParameter<
std::vector<
edm::InputTag>>(
"lheInfo")),
245 [this](const
edm::InputTag&
tag) {
return mayConsume<LHEEventProduct>(
tag); })),
247 lheLabel_, [
this](
const edm::InputTag&
tag) {
return mayConsume<LHERunInfoProduct, edm::InRun>(
tag); })),
249 mayConsume<GenLumiInfoHeader, edm::InLumi>(params.getParameter<
edm::InputTag>(
"genLumiInfoHeader"))),
250 namedWeightIDs_(params.getParameter<std::vector<std::string>>(
"namedWeightIDs")),
251 namedWeightLabels_(params.getParameter<std::vector<std::string>>(
"namedWeightLabels")),
252 lheWeightPrecision_(params.getParameter<int32_t>(
"lheWeightPrecision")),
253 maxPdfWeights_(params.getParameter<uint32_t>(
"maxPdfWeights")),
254 debug_(params.getUntrackedParameter<
bool>(
"debug",
false)),
255 debugRun_(debug_.load()),
256 hasIssuedWarning_(
false) {
257 produces<nanoaod::FlatTable>();
258 produces<std::string>(
"genModel");
259 produces<nanoaod::FlatTable>(
"LHEScale");
260 produces<nanoaod::FlatTable>(
"LHEPdf");
261 produces<nanoaod::FlatTable>(
"LHEReweighting");
262 produces<nanoaod::FlatTable>(
"LHENamed");
263 produces<nanoaod::FlatTable>(
"PS");
264 produces<nanoaod::MergeableCounterTable, edm::Transition::EndRun>();
265 if (namedWeightIDs_.size() != namedWeightLabels_.size()) {
266 throw cms::Exception(
"Configuration",
"Size mismatch between namedWeightIDs & namedWeightLabels");
270 uint32_t
lhaid = pdfps.getParameter<uint32_t>(
"lhaid");
271 preferredPDFLHAIDs_.push_back(lhaid);
273 lhaNameToID_[name +
".LHgrid"] =
lhaid;
289 auto out = std::make_unique<nanoaod::FlatTable>(1,
"genWeight",
true);
290 out->setDoc(
"generator weight");
294 std::string model_label = streamCache(
id)->countermap.getLabel();
295 auto outM = std::make_unique<std::string>((!model_label.empty()) ?
std::string(
"GenModel_") + model_label :
"");
297 bool getLHEweightsFromGenInfo = !model_label.empty();
300 std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab;
301 std::unique_ptr<nanoaod::FlatTable> genPSTab;
304 for (
const auto& lheTag : lheTag_) {
311 if (getLHEweightsFromGenInfo)
313 <<
"Found both a LHEEventProduct and a GenLumiInfoHeader: will only save weights from LHEEventProduct.\n";
315 const DynamicWeightChoice* weightChoice = runCache(iEvent.
getRun().
index());
318 counter, weightChoice, weight, *lheInfo, *genInfo, lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab, genPSTab);
319 }
else if (getLHEweightsFromGenInfo) {
320 const DynamicWeightChoiceGenInfo* weightChoice = &(streamCache(
id)->weightChoice);
321 fillLHEPdfWeightTablesFromGenInfo(
322 counter, weightChoice, weight, *genInfo, lheScaleTab, lhePdfTab, lheNamedTab, genPSTab);
328 fillOnlyPSWeightTable(counter, weight, *genInfo, genPSTab);
334 if (!hasIssuedWarning_.exchange(
true)) {
335 edm::LogWarning(
"LHETablesProducer") <<
"No LHEEventProduct, so there will be no LHE Tables\n";
347 const DynamicWeightChoice* weightChoice,
351 std::unique_ptr<nanoaod::FlatTable>& outScale,
352 std::unique_ptr<nanoaod::FlatTable>& outPdf,
353 std::unique_ptr<nanoaod::FlatTable>& outRwgt,
354 std::unique_ptr<nanoaod::FlatTable>& outNamed,
355 std::unique_ptr<nanoaod::FlatTable>& outPS)
const {
356 bool lheDebug = debug_.exchange(
359 const std::vector<std::string>& scaleWeightIDs = weightChoice->scaleWeightIDs;
360 const std::vector<std::string>& pdfWeightIDs = weightChoice->pdfWeightIDs;
361 const std::vector<std::string>& rwgtWeightIDs = weightChoice->rwgtIDs;
365 std::vector<double> wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wRwgt(rwgtWeightIDs.size(), 1),
366 wNamed(namedWeightIDs_.size(), 1);
369 printf(
"Weight %+9.5f rel %+9.5f for id %s\n",
weight.wgt,
weight.wgt / w0,
weight.id.c_str());
371 auto mScale =
std::find(scaleWeightIDs.begin(), scaleWeightIDs.end(),
weight.id);
372 if (mScale != scaleWeightIDs.end())
373 wScale[mScale - scaleWeightIDs.begin()] =
weight.wgt / w0;
375 auto mPDF =
std::find(pdfWeightIDs.begin(), pdfWeightIDs.end(),
weight.id);
376 if (mPDF != pdfWeightIDs.end())
377 wPDF[mPDF - pdfWeightIDs.begin()] =
weight.wgt / w0;
379 auto mRwgt =
std::find(rwgtWeightIDs.begin(), rwgtWeightIDs.end(),
weight.id);
380 if (mRwgt != rwgtWeightIDs.end())
381 wRwgt[mRwgt - rwgtWeightIDs.begin()] =
weight.wgt / w0;
383 auto mNamed =
std::find(namedWeightIDs_.begin(), namedWeightIDs_.end(),
weight.id);
384 if (mNamed != namedWeightIDs_.end())
385 wNamed[mNamed - namedWeightIDs_.begin()] =
weight.wgt / w0;
388 int vectorSize = (genProd.
weights().size() == 14 || genProd.
weights().size() == 46) ? 4 : 1;
389 std::vector<double> wPS(vectorSize, 1);
390 if (vectorSize > 1) {
391 double nominal = genProd.
weights()[1];
392 for (
unsigned int i = 6;
i < 10;
i++) {
393 wPS[
i - 6] = (genProd.
weights()[
i]) / nominal;
397 outPS->addColumn<
float>(
"",
399 vectorSize > 1 ?
"PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 " 400 "FSR=0.5; [2] is ISR=2 FSR=1; [3] is ISR=1 FSR=2 " 401 :
"dummy PS weight (1.0) ",
403 lheWeightPrecision_);
406 outScale->addColumn<
float>(
410 outPdf->addColumn<
float>(
414 outRwgt->addColumn<
float>(
418 outNamed->addColumnValue<
float>(
"originalXWGTUP",
420 "Nominal event weight in the LHE file",
422 for (
unsigned int i = 0, n = wNamed.size();
i <
n; ++
i) {
423 outNamed->addColumnValue<
float>(namedWeightLabels_[
i],
425 "LHE weight for id " + namedWeightIDs_[
i] +
", relative to nominal",
427 lheWeightPrecision_);
430 counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS);
434 const DynamicWeightChoiceGenInfo* weightChoice,
437 std::unique_ptr<nanoaod::FlatTable>& outScale,
438 std::unique_ptr<nanoaod::FlatTable>& outPdf,
439 std::unique_ptr<nanoaod::FlatTable>& outNamed,
440 std::unique_ptr<nanoaod::FlatTable>& outPS)
const {
441 const std::vector<unsigned int>& scaleWeightIDs = weightChoice->scaleWeightIDs;
442 const std::vector<unsigned int>& pdfWeightIDs = weightChoice->pdfWeightIDs;
443 const std::vector<unsigned int>& psWeightIDs = weightChoice->psWeightIDs;
447 double originalXWGTUP =
weights.at(1);
450 std::vector<double> wScale, wPDF, wPS;
451 for (
auto id : scaleWeightIDs)
452 wScale.push_back(
weights.at(
id) / w0);
453 for (
auto id : pdfWeightIDs) {
454 wPDF.push_back(
weights.at(
id) / w0);
456 if (!psWeightIDs.empty()) {
457 for (
auto id : psWeightIDs)
458 wPS.push_back((
weights.at(
id)) / w0);
463 outScale->addColumn<
float>(
467 outPdf->addColumn<
float>(
471 outPS->addColumn<
float>(
"",
473 wPS.size() > 1 ?
"PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 " 474 "FSR=0.5; [2] is ISR=2 FSR=1; [3] is ISR=1 FSR=2 " 475 :
"dummy PS weight (1.0) ",
477 lheWeightPrecision_);
480 outNamed->addColumnValue<
float>(
486 counter->incLHE(genWeight, wScale, wPDF, std::vector<double>(), std::vector<double>(), wPS);
492 std::unique_ptr<nanoaod::FlatTable>& outPS)
const {
493 int vectorSize = (genProd.
weights().size() == 14 || genProd.
weights().size() == 46) ? 4 : 1;
495 std::vector<double> wPS(vectorSize, 1);
496 if (vectorSize > 1) {
497 double nominal = genProd.
weights()[1];
498 for (
unsigned int i = 6;
i < 10;
i++) {
499 wPS[
i - 6] = (genProd.
weights()[
i]) / nominal;
504 outPS->addColumn<
float>(
"",
506 vectorSize > 1 ?
"PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 " 507 "FSR=0.5; [2] is ISR=2 FSR=1; [3] is ISR=1 FSR=2 " 508 :
"dummy PS weight (1.0) ",
510 lheWeightPrecision_);
512 counter->incGenOnly(genWeight);
513 counter->incPSOnly(genWeight, wPS);
520 bool lheDebug = debugRun_.exchange(
522 auto weightChoice = std::make_shared<DynamicWeightChoice>();
526 for (
const auto& lheLabel : lheLabel_) {
533 std::vector<ScaleVarWeight> scaleVariationIDs;
534 std::vector<PDFSetWeights> pdfSetWeightIDs;
535 std::vector<std::string> lheReweighingIDs;
537 std::regex weightgroupmg26x(
"<weightgroup\\s+(?:name|type)=\"(.*)\"\\s+combine=\"(.*)\"\\s*>");
538 std::regex weightgroup(
"<weightgroup\\s+combine=\"(.*)\"\\s+(?:name|type)=\"(.*)\"\\s*>");
539 std::regex weightgroupRwgt(
"<weightgroup\\s+(?:name|type)=\"(.*)\"\\s*>");
540 std::regex endweightgroup(
"</weightgroup>");
541 std::regex scalewmg26x(
542 "<weight\\s+(?:.*\\s+)?id=\"(\\d+)\"\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:[mM][uU][rR]|renscfact)=\"(" 543 "\\S+)\"\\s+(?:[mM][uU][Ff]|facscfact)=\"(\\S+)\")(\\s+.*)?</weight>");
545 "<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:mu[rR]|renscfact)=(\\S+)\\s+(" 546 "?:mu[Ff]|facscfact)=(\\S+)(\\s+.*)?)</weight>");
548 "<weight\\s+id=\"(\\d+)\">\\s*(?:PDF set|lhapdf|PDF|pdfset)\\s*=\\s*(\\d+)\\s*(?:\\s.*)?</weight>");
549 std::regex pdfwOld(
"<weight\\s+(?:.*\\s+)?id=\"(\\d+)\">\\s*Member \\s*(\\d+)\\s*(?:.*)</weight>");
550 std::regex pdfwmg26x(
551 "<weight\\s+id=\"(\\d+)\"\\s*MUR=\"(?:\\S+)\"\\s*MUF=\"(?:\\S+)\"\\s*(?:PDF " 552 "set|lhapdf|PDF|pdfset)\\s*=\\s*\"(\\d+)\"\\s*>\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?</" 554 std::regex rwgt(
"<weight\\s+id=\"(.+)\">(.+)?(</weight>)?");
557 if (iter->tag() !=
"initrwgt") {
559 std::cout <<
"Skipping LHE header with tag" << iter->tag() << std::endl;
563 std::cout <<
"Found LHE header with tag" << iter->tag() << std::endl;
564 std::vector<std::string>
lines = iter->lines();
565 bool missed_weightgroup =
567 bool ismg26x =
false;
568 for (
unsigned int iLine = 0, nLines = lines.size(); iLine < nLines;
570 boost::replace_all(lines[iLine],
"<",
"<");
571 boost::replace_all(lines[iLine],
">",
">");
572 if (std::regex_search(lines[iLine], groups, weightgroupmg26x)) {
576 for (
unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) {
579 if (std::regex_search(lines[iLine], groups, ismg26x ? weightgroupmg26x : weightgroup)) {
582 groupname = groups.str(1);
584 std::cout <<
">>> Looks like the beginning of a weight group for '" << groupname <<
"'" << std::endl;
585 if (groupname.find(
"scale_variation") == 0 || groupname ==
"Central scale variation") {
587 std::cout <<
">>> Looks like scale variation for theory uncertainties" << std::endl;
588 for (++iLine; iLine < nLines; ++iLine) {
591 if (std::regex_search(lines[iLine], groups, ismg26x ? scalewmg26x : scalew)) {
593 std::cout <<
" >>> Scale weight " << groups[1].str() <<
" for " << groups[3].str() <<
" , " 594 << groups[4].str() <<
" , " << groups[5].str() << std::endl;
595 scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
596 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
598 std::cout <<
">>> Looks like the end of a weight group" << std::endl;
599 if (!missed_weightgroup) {
602 missed_weightgroup =
false;
603 }
else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
605 std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end " 609 missed_weightgroup =
true;
614 }
else if (groupname ==
"PDF_variation" || groupname.find(
"PDF_variation ") == 0) {
616 std::cout <<
">>> Looks like a new-style block of PDF weights for one or more pdfs" << std::endl;
617 for (++iLine; iLine < nLines; ++iLine) {
620 if (std::regex_search(lines[iLine], groups, pdfw)) {
621 unsigned int lhaID = std::stoi(groups.str(2));
623 std::cout <<
" >>> PDF weight " << groups.str(1) <<
" for " << groups.str(2) <<
" = " << lhaID
625 if (pdfSetWeightIDs.empty() || !pdfSetWeightIDs.back().maybe_add(groups.str(1), lhaID)) {
626 pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
628 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
630 std::cout <<
">>> Looks like the end of a weight group" << std::endl;
631 if (!missed_weightgroup) {
634 missed_weightgroup =
false;
635 }
else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
637 std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end " 641 missed_weightgroup =
true;
646 }
else if (groupname ==
"PDF_variation1" || groupname ==
"PDF_variation2") {
648 std::cout <<
">>> Looks like a new-style block of PDF weights for multiple pdfs" << std::endl;
649 unsigned int lastid = 0;
650 for (++iLine; iLine < nLines; ++iLine) {
653 if (std::regex_search(lines[iLine], groups, pdfw)) {
654 unsigned int id = std::stoi(groups.str(1));
655 unsigned int lhaID = std::stoi(groups.str(2));
657 std::cout <<
" >>> PDF weight " << groups.str(1) <<
" for " << groups.str(2) <<
" = " << lhaID
659 if (
id != (lastid + 1) || pdfSetWeightIDs.empty()) {
660 pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
662 pdfSetWeightIDs.back().add(groups.str(1), lhaID);
665 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
667 std::cout <<
">>> Looks like the end of a weight group" << std::endl;
668 if (!missed_weightgroup) {
671 missed_weightgroup =
false;
672 }
else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
674 std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end " 678 missed_weightgroup =
true;
683 }
else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) {
685 std::cout <<
">>> Looks like an old-style PDF weight for an individual pdf" << std::endl;
686 unsigned int firstLhaID = lhaNameToID_.find(groupname)->second;
688 for (++iLine; iLine < nLines; ++iLine) {
691 if (std::regex_search(lines[iLine], groups, ismg26x ? pdfwmg26x : pdfwOld)) {
692 unsigned int member = 0;
694 member = std::stoi(groups.str(2));
696 if (!groups.str(4).empty()) {
697 member = std::stoi(groups.str(4));
700 unsigned int lhaID = member + firstLhaID;
702 std::cout <<
" >>> PDF weight " << groups.str(1) <<
" for " << member <<
" = " << lhaID
706 pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
709 pdfSetWeightIDs.back().add(groups.str(1), lhaID);
711 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
713 std::cout <<
">>> Looks like the end of a weight group" << std::endl;
714 if (!missed_weightgroup) {
717 missed_weightgroup =
false;
718 }
else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
720 std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end " 724 missed_weightgroup =
true;
730 for (++iLine; iLine < nLines; ++iLine) {
733 if (std::regex_search(lines[iLine], groups, endweightgroup)) {
735 std::cout <<
">>> Looks like the end of a weight group" << std::endl;
736 if (!missed_weightgroup) {
739 missed_weightgroup =
false;
740 }
else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
742 std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end " 746 missed_weightgroup =
true;
752 }
else if (std::regex_search(lines[iLine], groups, weightgroupRwgt)) {
754 if (groupname ==
"mg_reweighting") {
756 std::cout <<
">>> Looks like a LHE weights for reweighting" << std::endl;
757 for (++iLine; iLine < nLines; ++iLine) {
760 if (std::regex_search(lines[iLine], groups, rwgt)) {
763 std::cout <<
" >>> LHE reweighting weight: " << rwgtID << std::endl;
764 if (
std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
766 lheReweighingIDs.emplace_back(rwgtID);
768 }
else if (std::regex_search(lines[iLine], endweightgroup)) {
770 std::cout <<
">>> Looks like the end of a weight group" << std::endl;
771 if (!missed_weightgroup) {
774 missed_weightgroup =
false;
775 }
else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
777 std::cout <<
">>> Looks like the beginning of a new weight group, I will assume I missed the end " 781 missed_weightgroup =
true;
792 std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
794 std::cout <<
"Found " << scaleVariationIDs.size() <<
" scale variations: " << std::endl;
795 std::stringstream scaleDoc;
796 scaleDoc <<
"LHE scale variation weights (w_var / w_nominal); ";
797 for (
unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
798 const auto& sw = scaleVariationIDs[isw];
801 scaleDoc <<
"[" << isw <<
"] is " << sw.label;
802 weightChoice->scaleWeightIDs.push_back(sw.wid);
804 printf(
" id %s: scales ren = % .2f fact = % .2f text = %s\n",
810 if (!scaleVariationIDs.empty())
811 weightChoice->scaleWeightsDoc = scaleDoc.str();
815 std::cout <<
"Found " << pdfSetWeightIDs.size() <<
" PDF set errors: " << std::endl;
816 for (
const auto& pw : pdfSetWeightIDs)
817 printf(
"lhaIDs %6d - %6d (%3lu weights: %s, ... )\n",
821 pw.wids.front().c_str());
826 std::cout <<
"Found " << lheReweighingIDs.size() <<
" reweighting weights" << std::endl;
828 std::copy(lheReweighingIDs.begin(), lheReweighingIDs.end(), std::back_inserter(weightChoice->rwgtIDs));
830 std::stringstream pdfDoc;
831 pdfDoc <<
"LHE pdf variation weights (w_var / w_nominal) for LHA IDs ";
833 for (uint32_t
lhaid : preferredPDFLHAIDs_) {
834 for (
const auto& pw : pdfSetWeightIDs) {
835 if (pw.lhaIDs.first !=
lhaid && pw.lhaIDs.first != (
lhaid + 1))
837 if (pw.wids.size() == 1)
839 pdfDoc << pw.lhaIDs.first <<
" - " << pw.lhaIDs.second;
840 weightChoice->pdfWeightIDs = pw.wids;
841 if (maxPdfWeights_ < pw.wids.size()) {
842 weightChoice->pdfWeightIDs.resize(maxPdfWeights_);
843 pdfDoc <<
", truncated to the first " << maxPdfWeights_ <<
" replicas";
845 weightChoice->pdfWeightsDoc = pdfDoc.str();
859 return std::make_unique<LumiCacheInfoHolder>();
863 streamCache(
id)->clear();
868 auto counterMap = &(streamCache(
id)->countermap);
870 lumiBlock.
getByToken(genLumiInfoHeadTag_, genLumiInfoHead);
871 if (!genLumiInfoHead.
isValid())
873 <<
"No GenLumiInfoHeader product found, will not fill generator model string.\n";
875 if (genLumiInfoHead.
isValid()) {
877 boost::replace_all(label,
"-",
"_");
878 boost::replace_all(label,
"/",
"_");
880 counterMap->setLabel(label);
882 if (genLumiInfoHead.
isValid()) {
883 auto weightChoice = &(streamCache(
id)->weightChoice);
885 std::vector<ScaleVarWeight> scaleVariationIDs;
886 std::vector<PDFSetWeights> pdfSetWeightIDs;
887 weightChoice->psWeightIDs.clear();
889 std::regex scalew(
"LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+)\\,\\s+mur=(\\S+)\\smuf=(\\S+)");
890 std::regex pdfw(
"LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+),\\s+Member\\s+(\\d+)\\s+of\\ssets\\s+(\\w+\\b)");
893 std::unordered_map<std::string, uint32_t> knownPDFSetsFromGenInfo_;
894 unsigned int weightIter = 0;
895 for (
auto line : weightNames) {
896 if (std::regex_search(
line, groups, scalew)) {
897 auto id = groups.str(1);
898 auto group = groups.str(2);
899 auto mur = groups.str(3);
900 auto muf = groups.str(4);
901 if (
group.find(
"Central scale variation") != std::string::npos)
902 scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4));
903 }
else if (std::regex_search(
line, groups, pdfw)) {
904 auto id = groups.str(1);
905 auto group = groups.str(2);
906 auto memberid = groups.str(3);
907 auto pdfset = groups.str(4);
908 if (
group.find(pdfset) != std::string::npos) {
909 if (knownPDFSetsFromGenInfo_.find(pdfset) == knownPDFSetsFromGenInfo_.end()) {
910 knownPDFSetsFromGenInfo_[pdfset] = std::atoi(
id.c_str());
911 pdfSetWeightIDs.emplace_back(
id, std::atoi(
id.c_str()));
913 pdfSetWeightIDs.back().add(
id, std::atoi(
id.c_str()));
915 }
else if (
line.find(
"isrDef") != std::string::npos ||
916 line.find(
"fsrDef") != std::string::npos) {
917 weightChoice->psWeightIDs.push_back(weightIter);
922 weightChoice->scaleWeightIDs.clear();
923 weightChoice->pdfWeightIDs.clear();
925 std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
926 std::stringstream scaleDoc;
927 scaleDoc <<
"LHE scale variation weights (w_var / w_nominal); ";
928 for (
unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) {
929 const auto& sw = scaleVariationIDs[isw];
932 scaleDoc <<
"[" << isw <<
"] is " << sw.label;
933 weightChoice->scaleWeightIDs.push_back(std::atoi(sw.wid.c_str()));
935 if (!scaleVariationIDs.empty())
936 weightChoice->scaleWeightsDoc = scaleDoc.str();
937 std::stringstream pdfDoc;
938 pdfDoc <<
"LHE pdf variation weights (w_var / w_nominal) for LHA names ";
940 for (
const auto& pw : pdfSetWeightIDs) {
941 if (pw.wids.size() == 1)
943 for (
auto wantedpdf : lhaNameToID_) {
944 auto pdfname = wantedpdf.first;
945 if (knownPDFSetsFromGenInfo_.find(pdfname) == knownPDFSetsFromGenInfo_.end())
947 uint32_t
lhaid = knownPDFSetsFromGenInfo_.at(pdfname);
948 if (pw.lhaIDs.first != lhaid)
951 for (
auto x : pw.wids)
952 weightChoice->pdfWeightIDs.push_back(std::atoi(
x.c_str()));
953 if (maxPdfWeights_ < pw.wids.size()) {
954 weightChoice->pdfWeightIDs.resize(maxPdfWeights_);
955 pdfDoc <<
", truncated to the first " << maxPdfWeights_ <<
" replicas";
957 weightChoice->pdfWeightsDoc = pdfDoc.str();
968 return std::make_shared<CounterMap>();
974 CounterMap* runCounterMap)
const override {
975 runCounterMap->merge(streamCache(
id)->countermap);
981 auto out = std::make_unique<nanoaod::MergeableCounterTable>();
983 for (
auto x : runCounterMap->countermap) {
984 auto runCounter = &(
x.second);
988 out->addInt(
"genEventCount" + label,
"event count" + doclabel, runCounter->num);
989 out->addFloat(
"genEventSumw" + label,
"sum of gen weights" + doclabel, runCounter->sumw);
990 out->addFloat(
"genEventSumw2" + label,
"sum of gen (weight^2)" + doclabel, runCounter->sumw2);
992 double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1;
993 auto sumScales = runCounter->sumScale;
994 for (
auto&
val : sumScales)
996 out->addVFloatWithNorm(
"LHEScaleSumw" + label,
997 "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel,
1000 auto sumPDFs = runCounter->sumPDF;
1001 for (
auto&
val : sumPDFs)
1003 out->addVFloatWithNorm(
"LHEPdfSumw" + label,
1004 "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel,
1007 if (!runCounter->sumRwgt.empty()) {
1008 auto sumRwgts = runCounter->sumRwgt;
1009 for (
auto&
val : sumRwgts)
1011 out->addVFloatWithNorm(
"LHEReweightingSumw" + label,
1012 "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel,
1016 if (!runCounter->sumNamed.empty()) {
1017 for (
unsigned int i = 0, n = namedWeightLabels_.size();
i <
n; ++
i) {
1018 out->addFloatWithNorm(
1019 "LHESumw_" + namedWeightLabels_[
i] + label,
1020 "Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[
i] +
", divided by genEventSumw" + doclabel,
1021 runCounter->sumNamed[
i] * norm,
1034 ->setComment(
"tag for the GenEventInfoProduct, to get the main weight");
1036 ->setComment(
"tag for the GenLumiInfoProduct, to get the model string");
1037 desc.
add<std::vector<edm::InputTag>>(
"lheInfo", std::vector<edm::InputTag>{{
"externalLHEProducer"}, {
"source"}})
1038 ->setComment(
"tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
1042 prefpdf.
add<uint32_t>(
"lhaid");
1043 desc.
addVPSet(
"preferredPDFs", prefpdf, std::vector<edm::ParameterSet>())
1045 "LHA PDF Ids of the preferred PDF sets, in order of preference (the first matching one will be used)");
1046 desc.
add<std::vector<std::string>>(
"namedWeightIDs")->setComment(
"set of LHA weight IDs for named LHE weights");
1047 desc.
add<std::vector<std::string>>(
"namedWeightLabels")
1048 ->setComment(
"output names for the namedWeightIDs (in the same order)");
1049 desc.
add<int32_t>(
"lheWeightPrecision")->setComment(
"Number of bits in the mantissa for LHE weights");
1050 desc.
add<uint32_t>(
"maxPdfWeights")->setComment(
"Maximum number of PDF weights to save (to crop NN replicas)");
1051 desc.
addOptionalUntracked<
bool>(
"debug")->setComment(
"dump out all LHE information for one event");
1052 descriptions.
add(
"genWeightsTable", desc);
1058 const std::vector<edm::EDGetTokenT<LHEEventProduct>>
lheTag_;
1059 const std::vector<edm::EDGetTokenT<LHERunInfoProduct>>
lheRunTag_;
const std::vector< edm::EDGetTokenT< LHERunInfoProduct > > lheRunTag_
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)
std::unique_ptr< LumiCacheInfoHolder > beginStream(edm::StreamID) const override
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)
void streamBeginLuminosityBlock(edm::StreamID id, edm::LuminosityBlock const &lumiBlock, edm::EventSetup const &eventSetup) const override
void globalEndRun(edm::Run const &, edm::EventSetup const &) const override
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
const edm::EDGetTokenT< GenEventInfoProduct > genTag_
void fillLHEPdfWeightTablesFromGenInfo(Counter *counter, const DynamicWeightChoiceGenInfo *weightChoice, double genWeight, const GenEventInfoProduct &genProd, std::unique_ptr< nanoaod::FlatTable > &outScale, std::unique_ptr< nanoaod::FlatTable > &outPdf, std::unique_ptr< nanoaod::FlatTable > &outNamed, std::unique_ptr< nanoaod::FlatTable > &outPS) const
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_
std::unordered_map< std::string, uint32_t > lhaNameToID_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
void fillLHEWeightTables(Counter *counter, const DynamicWeightChoice *weightChoice, double genWeight, const LHEEventProduct &lheProd, const GenEventInfoProduct &genProd, std::unique_ptr< nanoaod::FlatTable > &outScale, std::unique_ptr< nanoaod::FlatTable > &outPdf, std::unique_ptr< nanoaod::FlatTable > &outRwgt, std::unique_ptr< nanoaod::FlatTable > &outNamed, std::unique_ptr< nanoaod::FlatTable > &outPS) const
void globalEndRunSummary(edm::Run const &, edm::EventSetup const &, CounterMap *runCounterMap) const override
const std::vector< edm::InputTag > lheLabel_
const std::vector< WGT > & weights() const
std::function< unsigned int(align::ID)> Counter
std::shared_ptr< CounterMap > globalBeginRunSummary(edm::Run const &, edm::EventSetup const &) const override
#define DEFINE_FWK_MODULE(type)
headers_const_iterator headers_begin() const
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
void fillOnlyPSWeightTable(Counter *counter, double genWeight, const GenEventInfoProduct &genProd, std::unique_ptr< nanoaod::FlatTable > &outPS) const
ParameterDescriptionBase * add(U const &iLabel, T const &value)
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
auto vector_transform(std::vector< InputType > const &input, Function predicate) -> std::vector< typename std::remove_cv< typename std::remove_reference< decltype(predicate(input.front()))>::type >::type >
std::vector< std::string > namedWeightLabels_
~GenWeightsTableProducer() override
GenWeightsTableProducer(edm::ParameterSet const ¶ms)
void put(std::unique_ptr< PROD > product)
Put a new product.
void add(std::string const &label, ParameterSetDescription const &psetDescription)
unsigned int maxPdfWeights_
const std::vector< edm::EDGetTokenT< LHEEventProduct > > lheTag_
bool operator<(DTCELinkId const &lhs, DTCELinkId const &rhs)
void streamEndRunSummary(edm::StreamID id, edm::Run const &, edm::EventSetup const &, CounterMap *runCounterMap) const override
void globalEndRunProduce(edm::Run &iRun, edm::EventSetup const &, CounterMap const *runCounterMap) const override
std::vector< uint32_t > preferredPDFLHAIDs_
std::vector< double > & weights()
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_
const edm::EDGetTokenT< GenLumiInfoHeader > genLumiInfoHeadTag_
Power< A, B >::type pow(const A &a, const B &b)
def merge(dictlist, TELL=False)