CMS 3D CMS Logo

GenWeightsTableProducer.cc
Go to the documentation of this file.
15 #include "boost/algorithm/string.hpp"
16 
17 #include <vector>
18 #include <unordered_map>
19 #include <iostream>
20 #include <regex>
21 
22 namespace {
24  struct Counter {
25  Counter() : num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumRwgt(), sumNamed(), sumPS() {}
26 
27  // the counters
28  long long num;
29  long double sumw;
30  long double sumw2;
31  std::vector<long double> sumPDF, sumScale, sumRwgt, sumNamed, sumPS;
32 
33  void clear() {
34  num = 0;
35  sumw = 0;
36  sumw2 = 0;
37  sumPDF.clear();
38  sumScale.clear();
39  sumRwgt.clear();
40  sumNamed.clear(), sumPS.clear();
41  }
42 
43  // inc the counters
44  void incGenOnly(double w) {
45  num++;
46  sumw += w;
47  sumw2 += (w * w);
48  }
49 
50  void incPSOnly(double w0, const std::vector<double>& wPS) {
51  if (!wPS.empty()) {
52  if (sumPS.empty())
53  sumPS.resize(wPS.size(), 0);
54  for (unsigned int i = 0, n = wPS.size(); i < n; ++i)
55  sumPS[i] += (w0 * wPS[i]);
56  }
57  }
58 
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) {
65  // add up weights
66  incGenOnly(w0);
67  // then add up variations
68  if (!wScale.empty()) {
69  if (sumScale.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]);
73  }
74  if (!wPDF.empty()) {
75  if (sumPDF.empty())
76  sumPDF.resize(wPDF.size(), 0);
77  for (unsigned int i = 0, n = wPDF.size(); i < n; ++i)
78  sumPDF[i] += (w0 * wPDF[i]);
79  }
80  if (!wRwgt.empty()) {
81  if (sumRwgt.empty())
82  sumRwgt.resize(wRwgt.size(), 0);
83  for (unsigned int i = 0, n = wRwgt.size(); i < n; ++i)
84  sumRwgt[i] += (w0 * wRwgt[i]);
85  }
86  if (!wNamed.empty()) {
87  if (sumNamed.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]);
91  }
92  incPSOnly(w0, wPS);
93  }
94 
95  void merge(const Counter& other) {
96  num += other.num;
97  sumw += other.sumw;
98  sumw2 += other.sumw2;
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];
124  }
125  };
126 
127  struct CounterMap {
128  std::map<std::string, Counter> countermap;
129  Counter* active_el = nullptr;
130  std::string active_label = "";
131  void merge(const CounterMap& other) {
132  for (const auto& y : other.countermap)
133  countermap[y.first].merge(y.second);
134  active_el = nullptr;
135  }
136  void clear() {
137  for (auto x : countermap)
138  x.second.clear();
139  active_el = nullptr;
140  active_label = "";
141  }
142  void setLabel(std::string label) {
143  active_el = &(countermap[label]);
144  active_label = label;
145  }
146  void checkLabelSet() {
147  if (!active_el)
148  throw cms::Exception("LogicError", "Called CounterMap::get() before setting the active label\n");
149  }
150  Counter* get() {
151  checkLabelSet();
152  return active_el;
153  }
154  std::string& getLabel() {
155  checkLabelSet();
156  return active_label;
157  }
158  };
159 
161  struct DynamicWeightChoice {
162  // choice of LHE weights
163  // ---- scale ----
164  std::vector<std::string> scaleWeightIDs;
165  std::string scaleWeightsDoc;
166  // ---- pdf ----
167  std::vector<std::string> pdfWeightIDs;
168  std::string pdfWeightsDoc;
169  // ---- rwgt ----
170  std::vector<std::string> rwgtIDs;
171  std::string rwgtWeightDoc;
172  };
173 
174  struct DynamicWeightChoiceGenInfo {
175  // choice of LHE weights
176  // ---- scale ----
177  std::vector<unsigned int> scaleWeightIDs;
178  std::string scaleWeightsDoc;
179  // ---- pdf ----
180  std::vector<unsigned int> pdfWeightIDs;
181  std::string pdfWeightsDoc;
182  // ---- ps ----
183  std::vector<unsigned int> psWeightIDs;
184  };
185 
186  struct LumiCacheInfoHolder {
187  CounterMap countermap;
188  DynamicWeightChoiceGenInfo weightChoice;
189  void clear() {
190  countermap.clear();
191  weightChoice = DynamicWeightChoiceGenInfo();
192  }
193  };
194 
195  float stof_fortrancomp(const std::string& str) {
196  std::string::size_type match = str.find("d");
197  if (match != std::string::npos) {
198  std::string pre = str.substr(0, match);
199  std::string post = str.substr(match + 1);
200  return std::stof(pre) * std::pow(10.0f, std::stof(post));
201  } else {
202  return std::stof(str);
203  }
204  }
206  struct ScaleVarWeight {
207  std::string wid, label;
208  std::pair<float, float> scales;
209  ScaleVarWeight(const std::string& id, const std::string& text, const std::string& muR, const std::string& muF)
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);
213  }
214  };
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; }
220  void add(const std::string& wid, unsigned int lhaID) {
221  wids.push_back(wid);
222  lhaIDs.second = lhaID;
223  }
224  bool maybe_add(const std::string& wid, unsigned int lhaID) {
225  if (lhaID == lhaIDs.second + 1) {
226  lhaIDs.second++;
227  wids.push_back(wid);
228  return true;
229  } else {
230  return false;
231  }
232  }
233  };
234 } // namespace
235 
236 class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<LumiCacheInfoHolder>,
237  edm::RunCache<DynamicWeightChoice>,
238  edm::RunSummaryCache<CounterMap>,
239  edm::EndRunProducer> {
240 public:
242  : genTag_(consumes<GenEventInfoProduct>(params.getParameter<edm::InputTag>("genEvent"))),
243  lheLabel_(params.getParameter<std::vector<edm::InputTag>>("lheInfo")),
244  lheTag_(edm::vector_transform(lheLabel_,
245  [this](const edm::InputTag& tag) { return mayConsume<LHEEventProduct>(tag); })),
246  lheRunTag_(edm::vector_transform(
247  lheLabel_, [this](const edm::InputTag& tag) { return mayConsume<LHERunInfoProduct, edm::InRun>(tag); })),
248  genLumiInfoHeadTag_(
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");
267  }
268  for (const edm::ParameterSet& pdfps : params.getParameter<std::vector<edm::ParameterSet>>("preferredPDFs")) {
269  const std::string& name = pdfps.getParameter<std::string>("name");
270  uint32_t lhaid = pdfps.getParameter<uint32_t>("lhaid");
271  preferredPDFLHAIDs_.push_back(lhaid);
272  lhaNameToID_[name] = lhaid;
273  lhaNameToID_[name + ".LHgrid"] = lhaid;
274  }
275  }
276 
278 
279  void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
280  // get my counter for weights
281  Counter* counter = streamCache(id)->countermap.get();
282 
283  // generator information (always available)
285  iEvent.getByToken(genTag_, genInfo);
286  double weight = genInfo->weight();
287 
288  // table for gen info, always available
289  auto out = std::make_unique<nanoaod::FlatTable>(1, "genWeight", true);
290  out->setDoc("generator weight");
291  out->addColumnValue<float>("", weight, "generator weight", nanoaod::FlatTable::FloatColumn);
292  iEvent.put(std::move(out));
293 
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 : "");
296  iEvent.put(std::move(outM), "genModel");
297  bool getLHEweightsFromGenInfo = !model_label.empty();
298 
299  // tables for LHE weights, may not be filled
300  std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab;
301  std::unique_ptr<nanoaod::FlatTable> genPSTab;
302 
304  for (const auto& lheTag : lheTag_) {
305  iEvent.getByToken(lheTag, lheInfo);
306  if (lheInfo.isValid()) {
307  break;
308  }
309  }
310  if (lheInfo.isValid()) {
311  if (getLHEweightsFromGenInfo)
312  edm::LogWarning("LHETablesProducer")
313  << "Found both a LHEEventProduct and a GenLumiInfoHeader: will only save weights from LHEEventProduct.\n";
314  // get the dynamic choice of weights
315  const DynamicWeightChoice* weightChoice = runCache(iEvent.getRun().index());
316  // go fill tables
317  fillLHEWeightTables(
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);
323  lheRwgtTab.reset(new nanoaod::FlatTable(1, "LHEReweightingWeights", true));
324  //lheNamedTab.reset(new nanoaod::FlatTable(1, "LHENamedWeights", true));
325  //genPSTab.reset(new nanoaod::FlatTable(1, "PSWeight", true));
326  } else {
327  // Still try to add the PS weights
328  fillOnlyPSWeightTable(counter, weight, *genInfo, genPSTab);
329  // make dummy values
330  lheScaleTab.reset(new nanoaod::FlatTable(1, "LHEScaleWeights", true));
331  lhePdfTab.reset(new nanoaod::FlatTable(1, "LHEPdfWeights", true));
332  lheRwgtTab.reset(new nanoaod::FlatTable(1, "LHEReweightingWeights", true));
333  lheNamedTab.reset(new nanoaod::FlatTable(1, "LHENamedWeights", true));
334  if (!hasIssuedWarning_.exchange(true)) {
335  edm::LogWarning("LHETablesProducer") << "No LHEEventProduct, so there will be no LHE Tables\n";
336  }
337  }
338 
339  iEvent.put(std::move(lheScaleTab), "LHEScale");
340  iEvent.put(std::move(lhePdfTab), "LHEPdf");
341  iEvent.put(std::move(lheRwgtTab), "LHEReweighting");
342  iEvent.put(std::move(lheNamedTab), "LHENamed");
343  iEvent.put(std::move(genPSTab), "PS");
344  }
345 
347  const DynamicWeightChoice* weightChoice,
348  double genWeight,
349  const LHEEventProduct& lheProd,
350  const GenEventInfoProduct& genProd,
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(
357  false); // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
358 
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;
362 
363  double w0 = lheProd.originalXWGTUP();
364 
365  std::vector<double> wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wRwgt(rwgtWeightIDs.size(), 1),
366  wNamed(namedWeightIDs_.size(), 1);
367  for (auto& weight : lheProd.weights()) {
368  if (lheDebug)
369  printf("Weight %+9.5f rel %+9.5f for id %s\n", weight.wgt, weight.wgt / w0, weight.id.c_str());
370  // now we do it slowly, can be optimized
371  auto mScale = std::find(scaleWeightIDs.begin(), scaleWeightIDs.end(), weight.id);
372  if (mScale != scaleWeightIDs.end())
373  wScale[mScale - scaleWeightIDs.begin()] = weight.wgt / w0;
374 
375  auto mPDF = std::find(pdfWeightIDs.begin(), pdfWeightIDs.end(), weight.id);
376  if (mPDF != pdfWeightIDs.end())
377  wPDF[mPDF - pdfWeightIDs.begin()] = weight.wgt / w0;
378 
379  auto mRwgt = std::find(rwgtWeightIDs.begin(), rwgtWeightIDs.end(), weight.id);
380  if (mRwgt != rwgtWeightIDs.end())
381  wRwgt[mRwgt - rwgtWeightIDs.begin()] = weight.wgt / w0;
382 
383  auto mNamed = std::find(namedWeightIDs_.begin(), namedWeightIDs_.end(), weight.id);
384  if (mNamed != namedWeightIDs_.end())
385  wNamed[mNamed - namedWeightIDs_.begin()] = weight.wgt / w0;
386  }
387 
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]; // Called 'Baseline' in GenLumiInfoHeader
392  for (unsigned int i = 6; i < 10; i++) {
393  wPS[i - 6] = (genProd.weights()[i]) / nominal;
394  }
395  }
396  outPS.reset(new nanoaod::FlatTable(wPS.size(), "PSWeight", false));
397  outPS->addColumn<float>("",
398  wPS,
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_);
404 
405  outScale.reset(new nanoaod::FlatTable(wScale.size(), "LHEScaleWeight", false));
406  outScale->addColumn<float>(
407  "", wScale, weightChoice->scaleWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
408 
409  outPdf.reset(new nanoaod::FlatTable(wPDF.size(), "LHEPdfWeight", false));
410  outPdf->addColumn<float>(
411  "", wPDF, weightChoice->pdfWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
412 
413  outRwgt.reset(new nanoaod::FlatTable(wRwgt.size(), "LHEReweightingWeight", false));
414  outRwgt->addColumn<float>(
415  "", wRwgt, weightChoice->rwgtWeightDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
416 
417  outNamed.reset(new nanoaod::FlatTable(1, "LHEWeight", true));
418  outNamed->addColumnValue<float>("originalXWGTUP",
419  lheProd.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],
424  wNamed[i],
425  "LHE weight for id " + namedWeightIDs_[i] + ", relative to nominal",
427  lheWeightPrecision_);
428  }
429 
430  counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS);
431  }
432 
434  const DynamicWeightChoiceGenInfo* weightChoice,
435  double genWeight,
436  const GenEventInfoProduct& genProd,
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;
444 
445  auto weights = genProd.weights();
446  double w0 = weights.at(1);
447  double originalXWGTUP = weights.at(1);
448  ;
449 
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);
455  }
456  if (!psWeightIDs.empty()) {
457  for (auto id : psWeightIDs)
458  wPS.push_back((weights.at(id)) / w0);
459  } else
460  wPS.push_back(1.0);
461 
462  outScale.reset(new nanoaod::FlatTable(wScale.size(), "LHEScaleWeight", false));
463  outScale->addColumn<float>(
464  "", wScale, weightChoice->scaleWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
465 
466  outPdf.reset(new nanoaod::FlatTable(wPDF.size(), "LHEPdfWeight", false));
467  outPdf->addColumn<float>(
468  "", wPDF, weightChoice->pdfWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
469 
470  outPS.reset(new nanoaod::FlatTable(wPS.size(), "PSWeight", false));
471  outPS->addColumn<float>("",
472  wPS,
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_);
478 
479  outNamed.reset(new nanoaod::FlatTable(1, "LHEWeight", true));
480  outNamed->addColumnValue<float>(
481  "originalXWGTUP", originalXWGTUP, "Nominal event weight in the LHE file", nanoaod::FlatTable::FloatColumn);
482  /*for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) {
483  outNamed->addColumnValue<float>(namedWeightLabels_[i], wNamed[i], "LHE weight for id "+namedWeightIDs_[i]+", relative to nominal", nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
484  }*/
485 
486  counter->incLHE(genWeight, wScale, wPDF, std::vector<double>(), std::vector<double>(), wPS);
487  }
488 
490  double genWeight,
491  const GenEventInfoProduct& genProd,
492  std::unique_ptr<nanoaod::FlatTable>& outPS) const {
493  int vectorSize = (genProd.weights().size() == 14 || genProd.weights().size() == 46) ? 4 : 1;
494 
495  std::vector<double> wPS(vectorSize, 1);
496  if (vectorSize > 1) {
497  double nominal = genProd.weights()[1]; // Called 'Baseline' in GenLumiInfoHeader
498  for (unsigned int i = 6; i < 10; i++) {
499  wPS[i - 6] = (genProd.weights()[i]) / nominal;
500  }
501  }
502 
503  outPS.reset(new nanoaod::FlatTable(wPS.size(), "PSWeight", false));
504  outPS->addColumn<float>("",
505  wPS,
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_);
511 
512  counter->incGenOnly(genWeight);
513  counter->incPSOnly(genWeight, wPS);
514  }
515 
516  // create an empty counter
517  std::shared_ptr<DynamicWeightChoice> globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override {
519 
520  bool lheDebug = debugRun_.exchange(
521  false); // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
522  auto weightChoice = std::make_shared<DynamicWeightChoice>();
523 
524  // getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
525  //if (iRun.getByToken(lheRunTag_, lheInfo)) {
526  for (const auto& lheLabel : lheLabel_) {
527  iRun.getByLabel(lheLabel, lheInfo);
528  if (lheInfo.isValid()) {
529  break;
530  }
531  }
532  if (lheInfo.isValid()) {
533  std::vector<ScaleVarWeight> scaleVariationIDs;
534  std::vector<PDFSetWeights> pdfSetWeightIDs;
535  std::vector<std::string> lheReweighingIDs;
536 
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>");
544  std::regex scalew(
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>");
547  std::regex pdfw(
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.*)?</"
553  "weight>");
554  std::regex rwgt("<weight\\s+id=\"(.+)\">(.+)?(</weight>)?");
555  std::smatch groups;
556  for (auto iter = lheInfo->headers_begin(), end = lheInfo->headers_end(); iter != end; ++iter) {
557  if (iter->tag() != "initrwgt") {
558  if (lheDebug)
559  std::cout << "Skipping LHE header with tag" << iter->tag() << std::endl;
560  continue;
561  }
562  if (lheDebug)
563  std::cout << "Found LHE header with tag" << iter->tag() << std::endl;
564  std::vector<std::string> lines = iter->lines();
565  bool missed_weightgroup =
566  false; //Needed because in some of the samples ( produced with MG26X ) a small part of the header info is ordered incorrectly
567  bool ismg26x = false;
568  for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines;
569  ++iLine) { //First start looping through the lines to see which weightgroup pattern is matched
570  boost::replace_all(lines[iLine], "&lt;", "<");
571  boost::replace_all(lines[iLine], "&gt;", ">");
572  if (std::regex_search(lines[iLine], groups, weightgroupmg26x)) {
573  ismg26x = true;
574  }
575  }
576  for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) {
577  if (lheDebug)
578  std::cout << lines[iLine];
579  if (std::regex_search(lines[iLine], groups, ismg26x ? weightgroupmg26x : weightgroup)) {
580  std::string groupname = groups.str(2);
581  if (ismg26x)
582  groupname = groups.str(1);
583  if (lheDebug)
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") {
586  if (lheDebug)
587  std::cout << ">>> Looks like scale variation for theory uncertainties" << std::endl;
588  for (++iLine; iLine < nLines; ++iLine) {
589  if (lheDebug)
590  std::cout << " " << lines[iLine];
591  if (std::regex_search(lines[iLine], groups, ismg26x ? scalewmg26x : scalew)) {
592  if (lheDebug)
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)) {
597  if (lheDebug)
598  std::cout << ">>> Looks like the end of a weight group" << std::endl;
599  if (!missed_weightgroup) {
600  break;
601  } else
602  missed_weightgroup = false;
603  } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
604  if (lheDebug)
605  std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
606  "of the group."
607  << std::endl;
608  if (ismg26x)
609  missed_weightgroup = true;
610  --iLine; // rewind by one, and go back to the outer loop
611  break;
612  }
613  }
614  } else if (groupname == "PDF_variation" || groupname.find("PDF_variation ") == 0) {
615  if (lheDebug)
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) {
618  if (lheDebug)
619  std::cout << " " << lines[iLine];
620  if (std::regex_search(lines[iLine], groups, pdfw)) {
621  unsigned int lhaID = std::stoi(groups.str(2));
622  if (lheDebug)
623  std::cout << " >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
624  << std::endl;
625  if (pdfSetWeightIDs.empty() || !pdfSetWeightIDs.back().maybe_add(groups.str(1), lhaID)) {
626  pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
627  }
628  } else if (std::regex_search(lines[iLine], endweightgroup)) {
629  if (lheDebug)
630  std::cout << ">>> Looks like the end of a weight group" << std::endl;
631  if (!missed_weightgroup) {
632  break;
633  } else
634  missed_weightgroup = false;
635  } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
636  if (lheDebug)
637  std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
638  "of the group."
639  << std::endl;
640  if (ismg26x)
641  missed_weightgroup = true;
642  --iLine; // rewind by one, and go back to the outer loop
643  break;
644  }
645  }
646  } else if (groupname == "PDF_variation1" || groupname == "PDF_variation2") {
647  if (lheDebug)
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) {
651  if (lheDebug)
652  std::cout << " " << lines[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));
656  if (lheDebug)
657  std::cout << " >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID
658  << std::endl;
659  if (id != (lastid + 1) || pdfSetWeightIDs.empty()) {
660  pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
661  } else {
662  pdfSetWeightIDs.back().add(groups.str(1), lhaID);
663  }
664  lastid = id;
665  } else if (std::regex_search(lines[iLine], endweightgroup)) {
666  if (lheDebug)
667  std::cout << ">>> Looks like the end of a weight group" << std::endl;
668  if (!missed_weightgroup) {
669  break;
670  } else
671  missed_weightgroup = false;
672  } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
673  if (lheDebug)
674  std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
675  "of the group."
676  << std::endl;
677  if (ismg26x)
678  missed_weightgroup = true;
679  --iLine; // rewind by one, and go back to the outer loop
680  break;
681  }
682  }
683  } else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) {
684  if (lheDebug)
685  std::cout << ">>> Looks like an old-style PDF weight for an individual pdf" << std::endl;
686  unsigned int firstLhaID = lhaNameToID_.find(groupname)->second;
687  bool first = true;
688  for (++iLine; iLine < nLines; ++iLine) {
689  if (lheDebug)
690  std::cout << " " << lines[iLine];
691  if (std::regex_search(lines[iLine], groups, ismg26x ? pdfwmg26x : pdfwOld)) {
692  unsigned int member = 0;
693  if (ismg26x == 0) {
694  member = std::stoi(groups.str(2));
695  } else {
696  if (!groups.str(4).empty()) {
697  member = std::stoi(groups.str(4));
698  }
699  }
700  unsigned int lhaID = member + firstLhaID;
701  if (lheDebug)
702  std::cout << " >>> PDF weight " << groups.str(1) << " for " << member << " = " << lhaID
703  << std::endl;
704  //if (member == 0) continue; // let's keep also the central value for now
705  if (first) {
706  pdfSetWeightIDs.emplace_back(groups.str(1), lhaID);
707  first = false;
708  } else {
709  pdfSetWeightIDs.back().add(groups.str(1), lhaID);
710  }
711  } else if (std::regex_search(lines[iLine], endweightgroup)) {
712  if (lheDebug)
713  std::cout << ">>> Looks like the end of a weight group" << std::endl;
714  if (!missed_weightgroup) {
715  break;
716  } else
717  missed_weightgroup = false;
718  } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
719  if (lheDebug)
720  std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
721  "of the group."
722  << std::endl;
723  if (ismg26x)
724  missed_weightgroup = true;
725  --iLine; // rewind by one, and go back to the outer loop
726  break;
727  }
728  }
729  } else {
730  for (++iLine; iLine < nLines; ++iLine) {
731  if (lheDebug)
732  std::cout << " " << lines[iLine];
733  if (std::regex_search(lines[iLine], groups, endweightgroup)) {
734  if (lheDebug)
735  std::cout << ">>> Looks like the end of a weight group" << std::endl;
736  if (!missed_weightgroup) {
737  break;
738  } else
739  missed_weightgroup = false;
740  } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
741  if (lheDebug)
742  std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
743  "of the group."
744  << std::endl;
745  if (ismg26x)
746  missed_weightgroup = true;
747  --iLine; // rewind by one, and go back to the outer loop
748  break;
749  }
750  }
751  }
752  } else if (std::regex_search(lines[iLine], groups, weightgroupRwgt)) {
753  std::string groupname = groups.str(1);
754  if (groupname == "mg_reweighting") {
755  if (lheDebug)
756  std::cout << ">>> Looks like a LHE weights for reweighting" << std::endl;
757  for (++iLine; iLine < nLines; ++iLine) {
758  if (lheDebug)
759  std::cout << " " << lines[iLine];
760  if (std::regex_search(lines[iLine], groups, rwgt)) {
761  std::string rwgtID = groups.str(1);
762  if (lheDebug)
763  std::cout << " >>> LHE reweighting weight: " << rwgtID << std::endl;
764  if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) {
765  // we're only interested in the beggining of the block
766  lheReweighingIDs.emplace_back(rwgtID);
767  }
768  } else if (std::regex_search(lines[iLine], endweightgroup)) {
769  if (lheDebug)
770  std::cout << ">>> Looks like the end of a weight group" << std::endl;
771  if (!missed_weightgroup) {
772  break;
773  } else
774  missed_weightgroup = false;
775  } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) {
776  if (lheDebug)
777  std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end "
778  "of the group."
779  << std::endl;
780  if (ismg26x)
781  missed_weightgroup = true;
782  --iLine; // rewind by one, and go back to the outer loop
783  break;
784  }
785  }
786  }
787  }
788  }
789  //std::cout << "============= END [ " << iter->tag() << " ] ============ \n\n" << std::endl;
790 
791  // ----- SCALE VARIATIONS -----
792  std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end());
793  if (lheDebug)
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];
799  if (isw)
800  scaleDoc << "; ";
801  scaleDoc << "[" << isw << "] is " << sw.label;
802  weightChoice->scaleWeightIDs.push_back(sw.wid);
803  if (lheDebug)
804  printf(" id %s: scales ren = % .2f fact = % .2f text = %s\n",
805  sw.wid.c_str(),
806  sw.scales.first,
807  sw.scales.second,
808  sw.label.c_str());
809  }
810  if (!scaleVariationIDs.empty())
811  weightChoice->scaleWeightsDoc = scaleDoc.str();
812 
813  // ------ PDF VARIATIONS (take the preferred one) -----
814  if (lheDebug) {
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",
818  pw.lhaIDs.first,
819  pw.lhaIDs.second,
820  pw.wids.size(),
821  pw.wids.front().c_str());
822  }
823 
824  // ------ LHE REWEIGHTING -------
825  if (lheDebug) {
826  std::cout << "Found " << lheReweighingIDs.size() << " reweighting weights" << std::endl;
827  }
828  std::copy(lheReweighingIDs.begin(), lheReweighingIDs.end(), std::back_inserter(weightChoice->rwgtIDs));
829 
830  std::stringstream pdfDoc;
831  pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA IDs ";
832  bool found = false;
833  for (uint32_t lhaid : preferredPDFLHAIDs_) {
834  for (const auto& pw : pdfSetWeightIDs) {
835  if (pw.lhaIDs.first != lhaid && pw.lhaIDs.first != (lhaid + 1))
836  continue; // sometimes the first weight is not saved if that PDF is the nominal one for the sample
837  if (pw.wids.size() == 1)
838  continue; // only consider error sets
839  pdfDoc << pw.lhaIDs.first << " - " << pw.lhaIDs.second;
840  weightChoice->pdfWeightIDs = pw.wids;
841  if (maxPdfWeights_ < pw.wids.size()) {
842  weightChoice->pdfWeightIDs.resize(maxPdfWeights_); // drop some replicas
843  pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
844  }
845  weightChoice->pdfWeightsDoc = pdfDoc.str();
846  found = true;
847  break;
848  }
849  if (found)
850  break;
851  }
852  }
853  }
854  return weightChoice;
855  }
856 
857  // create an empty counter
858  std::unique_ptr<LumiCacheInfoHolder> beginStream(edm::StreamID) const override {
859  return std::make_unique<LumiCacheInfoHolder>();
860  }
861  // inizialize to zero at begin run
862  void streamBeginRun(edm::StreamID id, edm::Run const&, edm::EventSetup const&) const override {
863  streamCache(id)->clear();
864  }
866  edm::LuminosityBlock const& lumiBlock,
867  edm::EventSetup const& eventSetup) const override {
868  auto counterMap = &(streamCache(id)->countermap);
869  edm::Handle<GenLumiInfoHeader> genLumiInfoHead;
870  lumiBlock.getByToken(genLumiInfoHeadTag_, genLumiInfoHead);
871  if (!genLumiInfoHead.isValid())
872  edm::LogWarning("LHETablesProducer")
873  << "No GenLumiInfoHeader product found, will not fill generator model string.\n";
875  if (genLumiInfoHead.isValid()) {
876  label = genLumiInfoHead->configDescription();
877  boost::replace_all(label, "-", "_");
878  boost::replace_all(label, "/", "_");
879  }
880  counterMap->setLabel(label);
881 
882  if (genLumiInfoHead.isValid()) {
883  auto weightChoice = &(streamCache(id)->weightChoice);
884 
885  std::vector<ScaleVarWeight> scaleVariationIDs;
886  std::vector<PDFSetWeights> pdfSetWeightIDs;
887  weightChoice->psWeightIDs.clear();
888 
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)");
891  std::smatch groups;
892  auto weightNames = genLumiInfoHead->weightNames();
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)) { // scale variation
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)) { // PDF variation
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()));
912  } else
913  pdfSetWeightIDs.back().add(id, std::atoi(id.c_str()));
914  }
915  } else if (line.find("isrDef") != std::string::npos ||
916  line.find("fsrDef") != std::string::npos) { // PS variation
917  weightChoice->psWeightIDs.push_back(weightIter);
918  }
919  weightIter++;
920  }
921 
922  weightChoice->scaleWeightIDs.clear();
923  weightChoice->pdfWeightIDs.clear();
924 
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];
930  if (isw)
931  scaleDoc << "; ";
932  scaleDoc << "[" << isw << "] is " << sw.label;
933  weightChoice->scaleWeightIDs.push_back(std::atoi(sw.wid.c_str()));
934  }
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 ";
939  bool found = false;
940  for (const auto& pw : pdfSetWeightIDs) {
941  if (pw.wids.size() == 1)
942  continue; // only consider error sets
943  for (auto wantedpdf : lhaNameToID_) {
944  auto pdfname = wantedpdf.first;
945  if (knownPDFSetsFromGenInfo_.find(pdfname) == knownPDFSetsFromGenInfo_.end())
946  continue;
947  uint32_t lhaid = knownPDFSetsFromGenInfo_.at(pdfname);
948  if (pw.lhaIDs.first != lhaid)
949  continue;
950  pdfDoc << pdfname;
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_); // drop some replicas
955  pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas";
956  }
957  weightChoice->pdfWeightsDoc = pdfDoc.str();
958  found = true;
959  break;
960  }
961  if (found)
962  break;
963  }
964  }
965  }
966  // create an empty counter
967  std::shared_ptr<CounterMap> globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override {
968  return std::make_shared<CounterMap>();
969  }
970  // add this stream to the summary
972  edm::Run const&,
973  edm::EventSetup const&,
974  CounterMap* runCounterMap) const override {
975  runCounterMap->merge(streamCache(id)->countermap);
976  }
977  // nothing to do per se
978  void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {}
979  // write the total to the run
980  void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, CounterMap const* runCounterMap) const override {
981  auto out = std::make_unique<nanoaod::MergeableCounterTable>();
982 
983  for (auto x : runCounterMap->countermap) {
984  auto runCounter = &(x.second);
985  std::string label = (!x.first.empty()) ? (std::string("_") + x.first) : "";
986  std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : "";
987 
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);
991 
992  double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1;
993  auto sumScales = runCounter->sumScale;
994  for (auto& val : sumScales)
995  val *= norm;
996  out->addVFloatWithNorm("LHEScaleSumw" + label,
997  "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel,
998  sumScales,
999  runCounter->sumw);
1000  auto sumPDFs = runCounter->sumPDF;
1001  for (auto& val : sumPDFs)
1002  val *= norm;
1003  out->addVFloatWithNorm("LHEPdfSumw" + label,
1004  "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel,
1005  sumPDFs,
1006  runCounter->sumw);
1007  if (!runCounter->sumRwgt.empty()) {
1008  auto sumRwgts = runCounter->sumRwgt;
1009  for (auto& val : sumRwgts)
1010  val *= norm;
1011  out->addVFloatWithNorm("LHEReweightingSumw" + label,
1012  "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel,
1013  sumRwgts,
1014  runCounter->sumw);
1015  }
1016  if (!runCounter->sumNamed.empty()) { // it could be empty if there's no LHE info in the sample
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,
1022  runCounter->sumw);
1023  }
1024  }
1025  }
1026  iRun.put(std::move(out));
1027  }
1028  // nothing to do here
1029  void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {}
1030 
1033  desc.add<edm::InputTag>("genEvent", edm::InputTag("generator"))
1034  ->setComment("tag for the GenEventInfoProduct, to get the main weight");
1035  desc.add<edm::InputTag>("genLumiInfoHeader", edm::InputTag("generator"))
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)");
1039 
1041  prefpdf.add<std::string>("name");
1042  prefpdf.add<uint32_t>("lhaid");
1043  desc.addVPSet("preferredPDFs", prefpdf, std::vector<edm::ParameterSet>())
1044  ->setComment(
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);
1053  }
1054 
1055 protected:
1057  const std::vector<edm::InputTag> lheLabel_;
1058  const std::vector<edm::EDGetTokenT<LHEEventProduct>> lheTag_;
1059  const std::vector<edm::EDGetTokenT<LHERunInfoProduct>> lheRunTag_;
1061 
1062  std::vector<uint32_t> preferredPDFLHAIDs_;
1063  std::unordered_map<std::string, uint32_t> lhaNameToID_;
1064  std::vector<std::string> namedWeightIDs_;
1065  std::vector<std::string> namedWeightLabels_;
1067  unsigned int maxPdfWeights_;
1068 
1069  mutable std::atomic<bool> debug_, debugRun_, hasIssuedWarning_;
1070 };
1071 
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
Definition: Run.h:280
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.
Definition: Event.h:125
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
const double w
Definition: UKUtility.cc:23
def copy(args, dbName)
void globalEndRun(edm::Run const &, edm::EventSetup const &) const override
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
bool getByToken(EDGetToken token, Handle< PROD > &result) const
const edm::EDGetTokenT< GenEventInfoProduct > genTag_
const std::vector< std::string > & weightNames() const
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
Definition: weight.py:1
headers_const_iterator headers_end() const
Run const & getRun() const
Definition: Event.cc:99
double weight() 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)
Definition: FindCaloHit.cc:20
uint16_t size_type
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
char const * label
std::function< unsigned int(align::ID)> Counter
std::shared_ptr< CounterMap > globalBeginRunSummary(edm::Run const &, edm::EventSetup const &) const override
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
headers_const_iterator headers_begin() const
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
Definition: matutil.cc:167
RunIndex index() const
Definition: Run.cc:21
double f[11][100]
#define end
Definition: vmac.h:39
void fillOnlyPSWeightTable(Counter *counter, double genWeight, const GenEventInfoProduct &genProd, std::unique_ptr< nanoaod::FlatTable > &outPS) const
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:74
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 >
Definition: transform.h:11
const std::string & configDescription() const
std::vector< std::string > namedWeightLabels_
GenWeightsTableProducer(edm::ParameterSet const &params)
void put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Run.h:108
void add(std::string const &label, ParameterSetDescription const &psetDescription)
const std::vector< edm::EDGetTokenT< LHEEventProduct > > lheTag_
bool operator<(DTCELinkId const &lhs, DTCELinkId const &rhs)
Definition: DTCELinkId.h:73
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_
HLT enums.
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.
Definition: Utils.h:10
std::vector< std::string > namedWeightIDs_
#define str(s)
const edm::EDGetTokenT< GenLumiInfoHeader > genLumiInfoHeadTag_
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
def move(src, dest)
Definition: eostools.py:511
Definition: Run.h:45
def merge(dictlist, TELL=False)
Definition: MatrixUtil.py:194