Skip to content

Commit

Permalink
New LHEWeightsTableProducer
Browse files Browse the repository at this point in the history
  • Loading branch information
guitargeek committed Dec 19, 2019
1 parent 5b0a828 commit 70cdb97
Show file tree
Hide file tree
Showing 4 changed files with 254 additions and 44 deletions.
1 change: 1 addition & 0 deletions PhysicsTools/NanoAOD/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
<use name="DataFormats/Candidate"/>
<use name="DataFormats/NanoAOD"/>
<use name="boost"/>
<use name="tinyxml2"/>
<export>
<lib name="1"/>
</export>
45 changes: 1 addition & 44 deletions PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -233,10 +233,6 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
hasIssuedWarning_(false) {
produces<nanoaod::FlatTable>();
produces<std::string>("genModel");
produces<nanoaod::FlatTable>("LHEScale");
produces<nanoaod::FlatTable>("LHEPdf");
produces<nanoaod::FlatTable>("LHEReweighting");
produces<nanoaod::FlatTable>("LHENamed");
produces<nanoaod::FlatTable>("PS");
produces<nanoaod::MergeableCounterTable, edm::Transition::EndRun>();
if (namedWeightIDs_.size() != namedWeightLabels_.size()) {
Expand Down Expand Up @@ -273,7 +269,6 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
iEvent.put(std::move(outM), "genModel");

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

edm::Handle<LHEEventProduct> lheInfo;
Expand All @@ -287,25 +282,16 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
// get the dynamic choice of weights
const DynamicWeightChoice* weightChoice = runCache(iEvent.getRun().index());
// go fill tables
fillLHEWeightTables(
counter, weightChoice, weight, *lheInfo, *genInfo, lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab, genPSTab);
fillLHEWeightTables(counter, weightChoice, weight, *lheInfo, *genInfo, genPSTab);
} else {
// 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));
lheRwgtTab.reset(new nanoaod::FlatTable(1, "LHEReweightingWeights", true));
lheNamedTab.reset(new nanoaod::FlatTable(1, "LHENamedWeights", true));
if (!hasIssuedWarning_.exchange(true)) {
edm::LogWarning("LHETablesProducer") << "No LHEEventProduct, so there will be no LHE Tables\n";
}
}

iEvent.put(std::move(lheScaleTab), "LHEScale");
iEvent.put(std::move(lhePdfTab), "LHEPdf");
iEvent.put(std::move(lheRwgtTab), "LHEReweighting");
iEvent.put(std::move(lheNamedTab), "LHENamed");
iEvent.put(std::move(genPSTab), "PS");
}

Expand All @@ -314,10 +300,6 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
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 {
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 Down Expand Up @@ -367,31 +349,6 @@ class GenWeightsTableProducer : public edm::global::EDProducer<edm::StreamCache<
nanoaod::FlatTable::FloatColumn,
lheWeightPrecision_);

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

outPdf.reset(new nanoaod::FlatTable(wPDF.size(), "LHEPdfWeight", false));
outPdf->addColumn<float>(
"", wPDF, weightChoice->pdfWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);

outRwgt.reset(new nanoaod::FlatTable(wRwgt.size(), "LHEReweightingWeight", false));
outRwgt->addColumn<float>(
"", wRwgt, weightChoice->rwgtWeightDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);

outNamed.reset(new nanoaod::FlatTable(1, "LHEWeight", true));
outNamed->addColumnValue<float>("originalXWGTUP",
lheProd.originalXWGTUP(),
"Nominal event weight in the LHE file",
nanoaod::FlatTable::FloatColumn);
for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) {
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, wRwgt, wNamed, wPS);
}

Expand Down
138 changes: 138 additions & 0 deletions PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/NanoAOD/interface/FlatTable.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <optional>
#include <tinyxml2.h>

namespace {

// join vector of strings to one string
std::string join(const std::vector<std::string>& vec, const char* delim) {
std::stringstream res;
std::copy(vec.begin(), vec.end(), std::ostream_iterator<std::string>(res, delim));
return res.str();
}

struct LHEWeightInfo {
std::string id;
std::string text;
std::optional<std::string> group = std::nullopt;
};


std::vector<LHEWeightInfo> getLHEWeightInfos(LHERunInfoProduct const& lheInfo) {
std::vector<LHEWeightInfo> out;

for (auto iter = lheInfo.headers_begin(), end = lheInfo.headers_end(); iter != end; ++iter) {
if (iter->tag() != "initrwgt") {
continue;
}

tinyxml2::XMLDocument xmlDoc;
xmlDoc.Parse(("<root>" + join(iter->lines(), "") + "</root>").c_str());
tinyxml2::XMLElement* root = xmlDoc.FirstChildElement("root");

for (auto* e = root->FirstChildElement(); e != nullptr; e = e->NextSiblingElement()) {
if (strcmp(e->Name(), "weight") == 0) {
// we are here if there is a weight that does not belong to any group
std::string text = "";
if(e->GetText()) text = e->GetText();
out.push_back({e->Attribute("id"), text});
}
if (strcmp(e->Name(), "weightgroup") == 0) {
std::string groupName = e->Attribute("name");
for (auto* inner = e->FirstChildElement("weight"); inner != nullptr;
inner = inner->NextSiblingElement("weight")) {
// we are here if there is a weight in a weightgroup
std::string text = "";
if(inner->GetText()) text = inner->GetText();
out.push_back({inner->Attribute("id"), text, groupName});
}
}
}
}

return out;
}

}

class LHEWeightsTableProducer : public edm::global::EDProducer<edm::RunCache<std::vector<LHEWeightInfo>>> {
public:
LHEWeightsTableProducer( edm::ParameterSet const & params ) :
lheInputTag_(params.getParameter<edm::InputTag>("lheInfo")),
lheToken_(consumes<LHEEventProduct>(params.getParameter<edm::InputTag>("lheInfo"))),
lheWeightPrecision_(params.getParameter<int32_t>("lheWeightPrecision"))
{
consumes<LHERunInfoProduct, edm::InRun>(lheInputTag_);
produces<nanoaod::FlatTable>();
}

void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
// tables for LHE weights, may not be filled
auto const& lheInfo = iEvent.get(lheToken_);

auto lheWeightsTab = std::make_unique<nanoaod::FlatTable>(1, "LHEWeight", true);

auto const& weightInfos = *runCache(iEvent.getRun().index());

double w0 = lheInfo.originalXWGTUP();

int i = 0;
if(lheInfo.weights().size() != weightInfos.size()) {
throw cms::Exception("LogicError", "Weight labels don't have the same size as weights!\n");
}
for(auto const& weight : lheInfo.weights()) {
if(!weightInfos[i].group) {
// for now we ignore the ungrouped weights
continue;
}
auto const& label = (*weightInfos[i].group) + "__" + weightInfos[i].id;
lheWeightsTab->addColumnValue<float>(label, weight.wgt / w0, weightInfos[i].text + " (w_var / w_nominal)", nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
++i;
}

iEvent.put(std::move(lheWeightsTab));
}

std::shared_ptr<std::vector<LHEWeightInfo>> globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override {
edm::Handle<LHERunInfoProduct> lheInfo;

auto weightChoice = std::make_shared<std::vector<LHEWeightInfo>>();

// getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499)
iRun.getByLabel(lheInputTag_, lheInfo);
if (lheInfo.isValid()) {
getLHEWeightInfos(*lheInfo).swap(*weightChoice);
}

return weightChoice;
}

// nothing to do here
void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {}

static void fillDescriptions(edm::ConfigurationDescriptions & descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("lheInfo", {"externalLHEProducer"})->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
desc.add<int32_t>("lheWeightPrecision", -1)->setComment("Number of bits in the mantissa for LHE weights");
descriptions.addDefault(desc);
}


protected:
const edm::InputTag lheInputTag_;
const edm::EDGetTokenT<LHEEventProduct> lheToken_;
int lheWeightPrecision_;
};

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(LHEWeightsTableProducer);
114 changes: 114 additions & 0 deletions PhysicsTools/NanoAOD/python/nanogen_cff.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
from PhysicsTools.NanoAOD.common_cff import *
from PhysicsTools.NanoAOD.taus_cff import *
from PhysicsTools.NanoAOD.jets_cff import *
from PhysicsTools.NanoAOD.globals_cff import *
from PhysicsTools.NanoAOD.genparticles_cff import *
from PhysicsTools.NanoAOD.particlelevel_cff import *
from PhysicsTools.NanoAOD.lheInfoTable_cfi import *

# duplicate definition with nano_cff right now
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)),
# for some 92X samples. Note that the nominal weight, 260000, is not included in the LHE ...
cms.PSet(name = cms.string("NNPDF30_nlo_as_0118"), lhaid = cms.uint32(260000)),
# 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_lo_as_0130"), lhaid = cms.uint32(262000)),
# some FXFX 80X samples have only this (e.g. WWTo1L1Nu2Q, WWTo4Q)
cms.PSet(name = cms.string("NNPDF30_nlo_nf_4_pdfas"), lhaid = cms.uint32(292000)),
# some FXFX 80X samples have only this (e.g. DYJetsToLL_Pt, WJetsToLNu_Pt, DYJetsToNuNu_Pt)
cms.PSet(name = cms.string("NNPDF30_nlo_nf_5_pdfas"), lhaid = cms.uint32(292200)),
),
namedWeightIDs = cms.vstring(),
namedWeightLabels = cms.vstring(),
lheWeightPrecision = cms.int32(14),
maxPdfWeights = cms.uint32(150),
debug = cms.untracked.bool(False),
)

lheWeightsTable = cms.EDProducer(
"LHEWeightsTableProducer",
lheInfo = cms.InputTag("externalLHEProducer"),
lheWeightPrecision = cms.int32(14),
)

nanoMetadata = cms.EDProducer("UniqueStringProducer",
strings = cms.PSet(
tag = cms.string("untagged"),
)
)

metGenTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
src = cms.InputTag("genMetTrue"),
name = cms.string("GenMET"),
doc = cms.string("Gen MET"),
singleton = cms.bool(True),
extension = cms.bool(False),
variables = cms.PSet(
pt = Var("pt", float, doc="pt", precision=10),
phi = Var("phi", float, doc="phi", precision=10),
),
)

nanogenSequence = cms.Sequence(
nanoMetadata+
particleLevel+
genJetTable+
patJetPartons+
genJetFlavourAssociation+
genJetFlavourTable+
genJetAK8Table+
genJetAK8FlavourAssociation+
genJetAK8FlavourTable+
tauGenJets+
tauGenJetsSelectorAllHadrons+
genVisTaus+
genVisTauTable+
genTable+
genWeightsTable+
lheWeightsTable+
genParticleTables+
tautagger+
rivetProducerHTXS+
particleLevelTables+
metGenTable+
lheInfoTable
)

NANOAODGENoutput = cms.OutputModule("NanoAODOutputModule",
compressionAlgorithm = cms.untracked.string('LZMA'),
compressionLevel = cms.untracked.int32(9),
dataset = cms.untracked.PSet(
dataTier = cms.untracked.string('NANOAODSIM'),
filterName = cms.untracked.string('')
),
fileName = cms.untracked.string('nanogen.root'),
outputCommands = cms.untracked.vstring(
'drop *',
"keep nanoaodFlatTable_*Table_*_*", # event data
"keep String_*_genModel_*", # generator model data
"keep nanoaodMergeableCounterTable_*Table_*_*", # accumulated per/run or per/lumi data
"keep nanoaodUniqueString_nanoMetadata_*_*", # basic metadata
)
)

def customizeNanoGEN(process):
process.genParticleTable.src = "genParticles"
process.patJetPartons.particles = "genParticles"
process.particleLevel.src = "generatorSmeared"
process.rivetProducerHTXS.HepMCCollection = "generatorSmeared"

process.genJetTable.src = "ak4GenJetsNoNu"
process.genJetFlavourAssociation.jets = process.genJetTable.src
process.genJetFlavourTable.src = process.genJetTable.src
process.genJetFlavourTable.jetFlavourInfos = "genJetFlavourAssociation"
process.genJetAK8Table.src = "ak8GenJetsNoNu"
process.genJetAK8FlavourAssociation.jets = process.genJetAK8Table.src
process.genJetAK8FlavourTable.src = process.genJetAK8Table.src
process.tauGenJets.GenParticles = "genParticles"
process.genVisTaus.srcGenParticles = "genParticles"

0 comments on commit 70cdb97

Please sign in to comment.