Skip to content

Commit

Permalink
give a bit more structure to weight counter class
Browse files Browse the repository at this point in the history
  • Loading branch information
@guitargeek
guitargeek committed Dec 20, 2019
1 parent d663d2f commit 45be54a
Show file tree
Hide file tree
Showing 5 changed files with 120 additions and 126 deletions.
191 changes: 94 additions & 97 deletions PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,39 +20,42 @@
#include <regex>

namespace {
/// ---- Cache object for running sums of weights ----
struct Counter {
Counter() : num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumRwgt(), sumNamed(), sumPS() {}

// the counters
long long num;
long double sumw;
long double sumw2;
std::vector<long double> sumPDF, sumScale, sumRwgt, sumNamed, sumPS;
void mergeSumVectors(std::vector<long double>& v1, std::vector<long double> const& v2) {
if (v1.empty() && !v2.empty())
v1.resize(v2.size(), 0);
if (!v2.empty())
for (unsigned int i = 0, n = v1.size(); i < n; ++i)
v1[i] += v2[i];
}

/// ---- Cache object for running sums of weights ----
class Counter {
public:
void clear() {
num = 0;
sumw = 0;
sumw2 = 0;
sumPDF.clear();
sumScale.clear();
sumRwgt.clear();
sumNamed.clear(), sumPS.clear();
num_ = 0;
sumw_ = 0;
sumw2_ = 0;
sumPDF_.clear();
sumScale_.clear();
sumRwgt_.clear();
sumNamed_.clear();
sumPS_.clear();
}

// inc the counters
void incGenOnly(double w) {
num++;
sumw += w;
sumw2 += (w * w);
num_++;
sumw_ += w;
sumw2_ += (w * w);
}

void incPSOnly(double w0, const std::vector<double>& wPS) {
if (!wPS.empty()) {
if (sumPS.empty())
sumPS.resize(wPS.size(), 0);
if (sumPS_.empty())
sumPS_.resize(wPS.size(), 0);
for (unsigned int i = 0, n = wPS.size(); i < n; ++i)
sumPS[i] += (w0 * wPS[i]);
sumPS_[i] += (w0 * wPS[i]);
}
}

Expand All @@ -66,62 +69,55 @@ namespace {
incGenOnly(w0);
// then add up variations
if (!wScale.empty()) {
if (sumScale.empty())
sumScale.resize(wScale.size(), 0);
if (sumScale_.empty())
sumScale_.resize(wScale.size(), 0);
for (unsigned int i = 0, n = wScale.size(); i < n; ++i)
sumScale[i] += (w0 * wScale[i]);
sumScale_[i] += (w0 * wScale[i]);
}
if (!wPDF.empty()) {
if (sumPDF.empty())
sumPDF.resize(wPDF.size(), 0);
if (sumPDF_.empty())
sumPDF_.resize(wPDF.size(), 0);
for (unsigned int i = 0, n = wPDF.size(); i < n; ++i)
sumPDF[i] += (w0 * wPDF[i]);
sumPDF_[i] += (w0 * wPDF[i]);
}
if (!wRwgt.empty()) {
if (sumRwgt.empty())
sumRwgt.resize(wRwgt.size(), 0);
if (sumRwgt_.empty())
sumRwgt_.resize(wRwgt.size(), 0);
for (unsigned int i = 0, n = wRwgt.size(); i < n; ++i)
sumRwgt[i] += (w0 * wRwgt[i]);
sumRwgt_[i] += (w0 * wRwgt[i]);
}
if (!wNamed.empty()) {
if (sumNamed.empty())
sumNamed.resize(wNamed.size(), 0);
if (sumNamed_.empty())
sumNamed_.resize(wNamed.size(), 0);
for (unsigned int i = 0, n = wNamed.size(); i < n; ++i)
sumNamed[i] += (w0 * wNamed[i]);
sumNamed_[i] += (w0 * wNamed[i]);
}
incPSOnly(w0, wPS);
}

void merge(const Counter& other) {
num += other.num;
sumw += other.sumw;
sumw2 += other.sumw2;
if (sumScale.empty() && !other.sumScale.empty())
sumScale.resize(other.sumScale.size(), 0);
if (sumPDF.empty() && !other.sumPDF.empty())
sumPDF.resize(other.sumPDF.size(), 0);
if (sumRwgt.empty() && !other.sumRwgt.empty())
sumRwgt.resize(other.sumRwgt.size(), 0);
if (sumNamed.empty() && !other.sumNamed.empty())
sumNamed.resize(other.sumNamed.size(), 0);
if (sumPS.empty() && !other.sumPS.empty())
sumPS.resize(other.sumPS.size(), 0);
if (!other.sumScale.empty())
for (unsigned int i = 0, n = sumScale.size(); i < n; ++i)
sumScale[i] += other.sumScale[i];
if (!other.sumPDF.empty())
for (unsigned int i = 0, n = sumPDF.size(); i < n; ++i)
sumPDF[i] += other.sumPDF[i];
if (!other.sumRwgt.empty())
for (unsigned int i = 0, n = sumRwgt.size(); i < n; ++i)
sumRwgt[i] += other.sumRwgt[i];
if (!other.sumNamed.empty())
for (unsigned int i = 0, n = sumNamed.size(); i < n; ++i)
sumNamed[i] += other.sumNamed[i];
if (!other.sumPS.empty())
for (unsigned int i = 0, n = sumPS.size(); i < n; ++i)
sumPS[i] += other.sumPS[i];
num_ += other.num_;
sumw_ += other.sumw_;
sumw2_ += other.sumw2_;

mergeSumVectors(sumScale_, other.sumScale_);
mergeSumVectors(sumPDF_, other.sumPDF_);
mergeSumVectors(sumRwgt_, other.sumRwgt_);
mergeSumVectors(sumNamed_, other.sumNamed_);
mergeSumVectors(sumPS_, other.sumPS_);
}

//private:
// the counters
long long num_ = 0;
long double sumw_ = 0;
long double sumw2_ = 0;

std::vector<long double> sumPDF_;
std::vector<long double> sumScale_;
std::vector<long double> sumRwgt_;
std::vector<long double> sumNamed_;
std::vector<long double> sumPS_;
};

struct CounterMap {
Expand Down Expand Up @@ -251,17 +247,15 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<

void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
// get my counter for weights
Counter* counter = streamCache(id)->get();
Counter& counter = *streamCache(id)->get();

// generator information (always available)
edm::Handle<GenEventInfoProduct> genInfo;
iEvent.getByToken(genTag_, genInfo);
double weight = genInfo->weight();
auto const& genInfo = iEvent.get(genTag_);

// table for gen info, always available
auto out = std::make_unique<nanoaod::FlatTable>(1, "genWeight", true);
out->setDoc("generator weight");
out->addColumnValue<float>("", weight, "generator weight", nanoaod::FlatTable::FloatColumn);
out->addColumnValue<float>("", genInfo.weight(), "generator weight", nanoaod::FlatTable::FloatColumn);
iEvent.put(std::move(out));

std::string model_label = streamCache(id)->getLabel();
Expand All @@ -280,12 +274,12 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
}
if (lheInfo.isValid()) {
// get the dynamic choice of weights
const DynamicWeightChoice* weightChoice = runCache(iEvent.getRun().index());
const DynamicWeightChoice& weightChoice = *runCache(iEvent.getRun().index());
// go fill tables
fillLHEWeightTables(counter, weightChoice, weight, *lheInfo, *genInfo, genPSTab);
fillLHEWeightTables(counter, weightChoice, genInfo.weight(), *lheInfo, genInfo, genPSTab);
} else {
// Still try to add the PS weights
fillOnlyPSWeightTable(counter, weight, *genInfo, genPSTab);
fillOnlyPSWeightTable(counter, genInfo.weight(), genInfo, genPSTab);
// make dummy values
if (!hasIssuedWarning_.exchange(true)) {
edm::LogWarning("LHETablesProducer") << "No LHEEventProduct, so there will be no LHE Tables\n";
Expand All @@ -295,23 +289,26 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
iEvent.put(std::move(genPSTab), "PS");
}

void fillLHEWeightTables(Counter* counter,
const DynamicWeightChoice* weightChoice,
void fillLHEWeightTables(Counter& counter,
const DynamicWeightChoice& weightChoice,
double genWeight,
const LHEEventProduct& lheProd,
const GenEventInfoProduct& genProd,
std::unique_ptr<nanoaod::FlatTable>& outPS) const {
bool lheDebug = debug_.exchange(
false); // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
// make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
bool lheDebug = debug_.exchange(false);

const std::vector<std::string>& scaleWeightIDs = weightChoice->scaleWeightIDs;
const std::vector<std::string>& pdfWeightIDs = weightChoice->pdfWeightIDs;
const std::vector<std::string>& rwgtWeightIDs = weightChoice->rwgtIDs;
const std::vector<std::string>& scaleWeightIDs = weightChoice.scaleWeightIDs;
const std::vector<std::string>& pdfWeightIDs = weightChoice.pdfWeightIDs;
const std::vector<std::string>& rwgtWeightIDs = weightChoice.rwgtIDs;

double w0 = lheProd.originalXWGTUP();

std::vector<double> wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wRwgt(rwgtWeightIDs.size(), 1),
wNamed(namedWeightIDs_.size(), 1);
std::vector<double> wScale(scaleWeightIDs.size(), 1);
std::vector<double> wPDF(pdfWeightIDs.size(), 1);
std::vector<double> wRwgt(rwgtWeightIDs.size(), 1);
std::vector<double> wNamed(namedWeightIDs_.size(), 1);

for (auto& weight : lheProd.weights()) {
if (lheDebug)
printf("Weight %+9.5f rel %+9.5f for id %s\n", weight.wgt, weight.wgt / w0, weight.id.c_str());
Expand Down Expand Up @@ -340,7 +337,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
wPS[i - 6] = (genProd.weights()[i]) / w0;
}
}
outPS.reset(new nanoaod::FlatTable(wPS.size(), "PSWeight", false));
outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
outPS->addColumn<float>("",
wPS,
vectorSize > 1 ? "PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 "
Expand All @@ -349,10 +346,10 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
nanoaod::FlatTable::FloatColumn,
lheWeightPrecision_);

counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS);
counter.incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS);
}

void fillOnlyPSWeightTable(Counter* counter,
void fillOnlyPSWeightTable(Counter& counter,
double genWeight,
const GenEventInfoProduct& genProd,
std::unique_ptr<nanoaod::FlatTable>& outPS) const {
Expand All @@ -365,7 +362,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
}
}

outPS.reset(new nanoaod::FlatTable(wPS.size(), "PSWeight", false));
outPS = std::make_unique<nanoaod::FlatTable>(wPS.size(), "PSWeight", false);
outPS->addColumn<float>("",
wPS,
vectorSize > 1 ? "PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 "
Expand All @@ -374,16 +371,16 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
nanoaod::FlatTable::FloatColumn,
lheWeightPrecision_);

counter->incGenOnly(genWeight);
counter->incPSOnly(genWeight, wPS);
counter.incGenOnly(genWeight);
counter.incPSOnly(genWeight, wPS);
}

// create an empty counter
std::shared_ptr<DynamicWeightChoice> globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override {
edm::Handle<LHERunInfoProduct> lheInfo;

bool lheDebug = debugRun_.exchange(
false); // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
// make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind)
bool lheDebug = debugRun_.exchange(false);
auto weightChoice = std::make_shared<DynamicWeightChoice>();

// getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
Expand Down Expand Up @@ -754,40 +751,40 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
auto out = std::make_unique<nanoaod::MergeableCounterTable>();

for (auto x : runCounterMap->countermap) {
auto runCounter = &(x.second);
auto& runCounter = x.second;
std::string label = std::string("_") + x.first;
std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : "";

out->addInt("genEventCount" + label, "event count" + doclabel, runCounter->num);
out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter->sumw);
out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter->sumw2);
out->addInt("genEventCount" + label, "event count" + doclabel, runCounter.num_);
out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter.sumw_);
out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter.sumw2_);

double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1;
auto sumScales = runCounter->sumScale;
double norm = runCounter.sumw_ ? 1.0 / runCounter.sumw_ : 1;
auto sumScales = runCounter.sumScale_;
for (auto& val : sumScales)
val *= norm;
out->addVFloat("LHEScaleSumw" + label,
"Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel,
sumScales);
auto sumPDFs = runCounter->sumPDF;
auto sumPDFs = runCounter.sumPDF_;
for (auto& val : sumPDFs)
val *= norm;
out->addVFloat(
"LHEPdfSumw" + label, "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel, sumPDFs);
if (!runCounter->sumRwgt.empty()) {
auto sumRwgts = runCounter->sumRwgt;
if (!runCounter.sumRwgt_.empty()) {
auto sumRwgts = runCounter.sumRwgt_;
for (auto& val : sumRwgts)
val *= norm;
out->addVFloat("LHEReweightingSumw" + label,
"Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel,
sumRwgts);
}
if (!runCounter->sumNamed.empty()) { // it could be empty if there's no LHE info in the sample
if (!runCounter.sumNamed_.empty()) { // it could be empty if there's no LHE info in the sample
for (unsigned int i = 0, n = namedWeightLabels_.size(); i < n; ++i) {
out->addFloat(
"LHESumw_" + namedWeightLabels_[i] + label,
"Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[i] + ", divided by genEventSumw" + doclabel,
runCounter->sumNamed[i] * norm);
runCounter.sumNamed_[i] * norm);
}
}
}
Expand Down
8 changes: 4 additions & 4 deletions PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,16 +94,16 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::RunCache<std
std::unordered_map<std::string, std::pair<std::string, std::vector<float>>> groupsWithWeights;
for (auto const& weight : lheInfo.weights()) {
auto& val = weightInfos[i].group ? groupsWithWeights[*weightInfos[i].group] : groupsWithWeights["ungrouped"];
if(val.first.empty()) {
val.first += ";id,text";
if (val.first.empty()) {
val.first += ";id,text";
}
val.first += ";" + weightInfos[i].id + "," + weightInfos[i].text;
val.second.push_back(weight.wgt / w0);
++i;
}
for (auto const& group : groupsWithWeights) {
if(std::find(weightgroups_.begin(), weightgroups_.end(), group.first) == weightgroups_.end()) {
continue;
if (std::find(weightgroups_.begin(), weightgroups_.end(), group.first) == weightgroups_.end()) {
continue;
}
std::string name = std::string("LHEWeight_") + group.first;
std::transform(name.begin(), name.end(), name.begin(), [](char ch) { return ch == ' ' ? '_' : ch; });
Expand Down
20 changes: 20 additions & 0 deletions PhysicsTools/NanoAOD/python/genWeightsTable_cfi.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
import FWCore.ParameterSet.Config as cms

genWeightsTable = cms.EDProducer("GenWeightsTableProducer",
genEvent = cms.InputTag("generator"),
genLumiInfoHeader = cms.InputTag("generator"),
lheInfo = cms.VInputTag(cms.InputTag("externalLHEProducer"), cms.InputTag("source")),
preferredPDFs = cms.VPSet( # see https://lhapdf.hepforge.org/pdfsets.html
cms.PSet( name = cms.string("PDF4LHC15_nnlo_30_pdfas"), lhaid = cms.uint32(91400) ),
cms.PSet( name = cms.string("NNPDF31_nnlo_hessian_pdfas"), lhaid = cms.uint32(306000) ),
cms.PSet( name = cms.string("NNPDF30_nlo_as_0118"), lhaid = cms.uint32(260000) ), # for some 92X samples. Note that the nominal weight, 260000, is not included in the LHE ...
cms.PSet( name = cms.string("NNPDF30_lo_as_0130"), lhaid = cms.uint32(262000) ), # some MLM 80X samples have only this (e.g. /store/mc/RunIISummer16MiniAODv2/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-madgraphMLM-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6_ext1-v2/120000/02A210D6-F5C3-E611-B570-008CFA197BD4.root )
cms.PSet( name = cms.string("NNPDF30_nlo_nf_4_pdfas"), lhaid = cms.uint32(292000) ), # some FXFX 80X samples have only this (e.g. WWTo1L1Nu2Q, WWTo4Q)
cms.PSet( name = cms.string("NNPDF30_nlo_nf_5_pdfas"), lhaid = cms.uint32(292200) ), # some FXFX 80X samples have only this (e.g. DYJetsToLL_Pt, WJetsToLNu_Pt, DYJetsToNuNu_Pt)
),
namedWeightIDs = cms.vstring(),
namedWeightLabels = cms.vstring(),
lheWeightPrecision = cms.int32(14),
maxPdfWeights = cms.uint32(150),
debug = cms.untracked.bool(False),
)
1 change: 1 addition & 0 deletions PhysicsTools/NanoAOD/python/nano_cff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
from PhysicsTools.NanoAOD.triggerObjects_cff import *
from PhysicsTools.NanoAOD.isotracks_cff import *
from PhysicsTools.NanoAOD.NanoAODEDMEventContent_cff import *
from PhysicsTools.NanoAOD.genWeightsTable_cfi import *

from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy
from Configuration.Eras.Modifier_run2_nanoAOD_94X2016_cff import run2_nanoAOD_94X2016
Expand Down
Loading

0 comments on commit 45be54a

Please sign in to comment.