Skip to content

Commit

Permalink
Merge pull request cms-sw#171 from adewit/fix-LHEWeights
Browse files Browse the repository at this point in the history 
Add PS weights to GenWeightsTableProducer
  • Loading branch information
@gpetruc
gpetruc authored May 23, 2018
2 parents d719e38 + 3115f26 commit 39122f8
Showing 1 changed file with 60 additions and 10 deletions.
70 changes: 60 additions & 10 deletions PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,24 +21,32 @@ namespace {
/// ---- Cache object for running sums of weights ----
struct Counter {
Counter() :
num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumNamed() {}
num(0), sumw(0), sumw2(0), sumPDF(), sumScale(), sumNamed(), sumPS() {}

// the counters
long long num;
long double sumw;
long double sumw2;
std::vector<long double> sumPDF, sumScale, sumNamed;
std::vector<long double> sumPDF, sumScale, sumNamed, sumPS;

void clear() {
num = 0; sumw = 0; sumw2 = 0;
sumPDF.clear(); sumScale.clear(); sumNamed.clear();
sumPDF.clear(); sumScale.clear(); sumNamed.clear(), sumPS.clear();
}

// inc the counters
void incGenOnly(double w) {
num++; sumw += w; sumw2 += (w*w);
}
void incLHE(double w0, const std::vector<double> & wScale, const std::vector<double> & wPDF, const std::vector<double> & wNamed) {

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

void incLHE(double w0, const std::vector<double> & wScale, const std::vector<double> & wPDF, const std::vector<double> & wNamed, const std::vector<double> & wPS) {
// add up weights
incGenOnly(w0);
// then add up variations
Expand All @@ -54,16 +62,19 @@ namespace {
if (sumNamed.empty()) sumNamed.resize(wNamed.size(), 0);
for (unsigned int i = 0, n = wNamed.size(); i < n; ++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 (sumNamed.empty() && !other.sumNamed.empty()) sumNamed.resize(other.sumNamed.size(),0);
if (sumPS.empty() && !other.sumPS.empty()) sumPS.resize(other.sumPS.size(),0);
for (unsigned int i = 0, n = sumScale.size(); i < n; ++i) sumScale[i] += other.sumScale[i];
for (unsigned int i = 0, n = sumPDF.size(); i < n; ++i) sumPDF[i] += other.sumPDF[i];
for (unsigned int i = 0, n = sumNamed.size(); i < n; ++i) sumNamed[i] += other.sumNamed[i];
for (unsigned int i = 0, n = sumPS.size(); i < n; ++i) sumPS[i] += other.sumPS[i];
}
};

Expand Down Expand Up @@ -135,6 +146,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
produces<nanoaod::FlatTable>("LHEScale");
produces<nanoaod::FlatTable>("LHEPdf");
produces<nanoaod::FlatTable>("LHENamed");
produces<nanoaod::FlatTable>("PS");
produces<nanoaod::MergeableCounterTable,edm::Transition::EndRun>();
if (namedWeightIDs_.size() != namedWeightLabels_.size()) {
throw cms::Exception("Configuration", "Size mismatch between namedWeightIDs & namedWeightLabels");
Expand All @@ -158,7 +170,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
edm::Handle<GenEventInfoProduct> genInfo;
iEvent.getByToken(genTag_, genInfo);
double weight = genInfo->weight();

// table for gen info, always available
auto out = std::make_unique<nanoaod::FlatTable>(1, "genWeight", true);
out->setDoc("generator weight");
Expand All @@ -167,16 +179,17 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<

// tables for LHE weights, may not be filled
std::unique_ptr<nanoaod::FlatTable> lheScaleTab, lhePdfTab, lheNamedTab;
std::unique_ptr<nanoaod::FlatTable> genPSTab;

edm::Handle<LHEEventProduct> lheInfo;
if (iEvent.getByToken(lheTag_, lheInfo)) {
// get the dynamic choice of weights
const DynamicWeightChoice * weightChoice = runCache(iEvent.getRun().index());
// go fill tables
fillLHEWeightTables(counter, weightChoice, weight, *lheInfo, lheScaleTab, lhePdfTab, lheNamedTab);
fillLHEWeightTables(counter, weightChoice, weight, *lheInfo, *genInfo, lheScaleTab, lhePdfTab, lheNamedTab, genPSTab);
} else {
// minimal book-keeping of weights
counter->incGenOnly(weight);
// Still try to add the PS weights
fillOnlyPSWeightTable(counter, weight, *genInfo, genPSTab);
// make dummy values
lheScaleTab.reset(new nanoaod::FlatTable(1, "LHEScaleWeights", true));
lhePdfTab.reset(new nanoaod::FlatTable(1, "LHEPdfWeights", true));
Expand All @@ -189,16 +202,19 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
iEvent.put(std::move(lheScaleTab), "LHEScale");
iEvent.put(std::move(lhePdfTab), "LHEPdf");
iEvent.put(std::move(lheNamedTab), "LHENamed");
iEvent.put(std::move(genPSTab), "PS");
}

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> & outNamed ) const
std::unique_ptr<nanoaod::FlatTable> & outNamed,
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)

Expand All @@ -221,6 +237,16 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
if (mNamed != namedWeightIDs_.end()) wNamed[mNamed-namedWeightIDs_.begin()] = weight.wgt/w0;
}

int vectorSize = genProd.weights().size() == 14 ? 4 : 1;
std::vector<double> wPS(vectorSize, 1);
if (vectorSize > 1 ) {
for (unsigned int i=6; i<10; i++){
wPS[i-6] = (genProd.weights()[i])/w0;
}
}
outPS.reset(new 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 FSR=0.5; [2] is ISR=2 FSR=1; [3] is ISR=1 FSR=2 " : "dummy PS weight (1.0) ", nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);

outScale.reset(new nanoaod::FlatTable(wScale.size(), "LHEScaleWeight", false));
outScale->addColumn<float>("", wScale, weightChoice->scaleWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);

Expand All @@ -233,9 +259,33 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
outNamed->addColumnValue<float>(namedWeightLabels_[i], wNamed[i], "LHE weight for id "+namedWeightIDs_[i]+", relative to nominal", nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
}

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

void fillOnlyPSWeightTable(
Counter * counter,
double genWeight,
const GenEventInfoProduct & genProd,
std::unique_ptr<nanoaod::FlatTable> & outPS ) const
{
int vectorSize = genProd.weights().size() == 14 ? 4 : 1;

std::vector<double> wPS(vectorSize, 1);
if (vectorSize > 1 ){
for (unsigned int i=6; i<10; i++){
wPS[i-6] = (genProd.weights()[i])/genWeight;
}
}

outPS.reset(new 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 FSR=0.5; [2] is ISR=2 FSR=1; [3] is ISR=1 FSR=2 " : "dummy PS weight (1.0) " , nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);

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;
Expand Down

0 comments on commit 39122f8

Please sign in to comment.