diff --git a/GeneratorInterface/Core/interface/WeightHelper.h b/GeneratorInterface/Core/interface/WeightHelper.h index 64c284fb4d7b5..5e3102aa60ea4 100644 --- a/GeneratorInterface/Core/interface/WeightHelper.h +++ b/GeneratorInterface/Core/interface/WeightHelper.h @@ -21,8 +21,8 @@ namespace gen { edm::OwnVector weightGroups() { return weightGroups_; } - std::unique_ptr weightProduct(std::vector); - std::unique_ptr weightProduct(std::vector); + std::unique_ptr weightProduct(std::vector, float w0); + std::unique_ptr weightProduct(std::vector, float w0); void setGroupInfo(); bool currentGroupIsScale(); bool currentGroupIsMEParam(); diff --git a/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc b/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc index 2cbb6748a08de..104083173e6be 100644 --- a/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc +++ b/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc @@ -43,10 +43,8 @@ class GenWeightProductProducer : public edm::one::EDProducer(edm::InputTag("generator"))), - //iConfig.getUntrackedParameter("lheSource", edm::InputTag("externalLHEProducer")))), - genEventToken_(consumes(edm::InputTag("generator"))) - //iConfig.getUntrackedParameter("lheSource", edm::InputTag("externalLHEProducer")))) + genLumiInfoToken_(consumes(iConfig.getParameter("genInfo"))), + genEventToken_(consumes(iConfig.getParameter("genInfo"))) { produces(); produces(); @@ -63,14 +61,12 @@ void GenWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { edm::Handle genEventInfo; iEvent.getByToken(genEventToken_, genEventInfo); - // Read weights from LHEEventProduct - auto weightProduct = weightHelper_.weightProduct(genEventInfo->weights()); + + float centralWeight = genEventInfo->weights().size() > 0 ? genEventInfo->weights().at(0) : 1.; + auto weightProduct = weightHelper_.weightProduct(genEventInfo->weights(), centralWeight); iEvent.put(std::move(weightProduct)); } -//void -//GenWeightProductProducer::endLuminosityBlockProduce(edm::LuminosityBlock& iLumi, edm::EventSetup const& iSetup) {} - void GenWeightProductProducer::beginLuminosityBlockProduce(edm::LuminosityBlock& iLumi, edm::EventSetup const& iSetup) { if (weightNames_.size() == 0) { diff --git a/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc b/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc index 0c98c2832d1a4..a4f0836314460 100644 --- a/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc +++ b/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc @@ -70,7 +70,7 @@ LHEWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe edm::Handle lheEventInfo; iEvent.getByToken(lheEventToken_, lheEventInfo); // Read weights from LHEEventProduct - auto weightProduct = weightHelper_.weightProduct(lheEventInfo->weights()); + auto weightProduct = weightHelper_.weightProduct(lheEventInfo->weights(), lheEventInfo->originalXWGTUP()); iEvent.put(std::move(weightProduct)); } diff --git a/GeneratorInterface/Core/src/LHEWeightHelper.cc b/GeneratorInterface/Core/src/LHEWeightHelper.cc index 55952e28210ed..1b3a7f4ee1039 100644 --- a/GeneratorInterface/Core/src/LHEWeightHelper.cc +++ b/GeneratorInterface/Core/src/LHEWeightHelper.cc @@ -136,6 +136,8 @@ namespace gen { else { weightGroups_.push_back(std::make_unique(name)); } + auto& group = weightGroups_.back(); + group.setDescription("Test description"); } else if (std::regex_match(headerLine, std::regex(".*.*\n*"))) { std::string fullTag = headerLine; diff --git a/GeneratorInterface/Core/src/WeightHelper.cc b/GeneratorInterface/Core/src/WeightHelper.cc index 091feab7d6265..fe96bd7298735 100644 --- a/GeneratorInterface/Core/src/WeightHelper.cc +++ b/GeneratorInterface/Core/src/WeightHelper.cc @@ -85,8 +85,8 @@ namespace gen { } // TODO: Could probably recycle this code better - std::unique_ptr WeightHelper::weightProduct(std::vector weights) { - auto weightProduct = std::make_unique(); + std::unique_ptr WeightHelper::weightProduct(std::vector weights, float w0) { + auto weightProduct = std::make_unique(w0); weightProduct->setNumWeightSets(weightGroups_.size()); int weightGroupIndex = 0; for (unsigned int i = 0; i < weights.size(); i++) { @@ -95,8 +95,8 @@ namespace gen { return std::move(weightProduct); } - std::unique_ptr WeightHelper::weightProduct(std::vector weights) { - auto weightProduct = std::make_unique(); + std::unique_ptr WeightHelper::weightProduct(std::vector weights, float w0) { + auto weightProduct = std::make_unique(w0); weightProduct->setNumWeightSets(weightGroups_.size()); int weightGroupIndex = 0; int i = 0; diff --git a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc index 88d0d3cb9ed08..49aaa321f9702 100644 --- a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc +++ b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc @@ -197,7 +197,8 @@ void ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe product.get(), _1)); // Should also zero out the weights in the GenInfoProduct - auto weightProduct = weightHelper_.weightProduct(partonLevel->weights()); + auto weightProduct = weightHelper_.weightProduct(partonLevel->weights(), + partonLevel->originalXWGTUP()); iEvent.put(std::move(weightProduct)); product->setScales(partonLevel->scales()); diff --git a/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc b/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc index f200d31ded731..cc4dc25809135 100644 --- a/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc @@ -21,13 +21,13 @@ #include namespace { - typedef std::pair, std::vector> WeightGroupsToStore; + typedef std::vector WeightGroupDataContainer; + typedef std::array, 2> WeightGroupsToStore; } // namespace class LHEWeightsTableProducer : public edm::global::EDProducer> { public: LHEWeightsTableProducer(edm::ParameterSet const& params) : - lheToken_(consumes(params.getParameter("lheInfo"))), lheWeightToken_(consumes(params.getParameter("lheWeights"))), lheWeightInfoToken_(consumes(params.getParameter("lheWeights"))), genWeightToken_(consumes(params.getParameter("genWeights"))), @@ -44,8 +44,6 @@ class LHEWeightsTableProducer : public edm::global::EDProducer lheWeightHandle; iEvent.getByToken(lheWeightToken_, lheWeightHandle); const GenWeightProduct* lheWeightProduct = lheWeightHandle.product(); @@ -57,17 +55,16 @@ class LHEWeightsTableProducer : public edm::global::EDProducerweights(); auto lheWeightTables = std::make_unique>(); - double w0 = lheInfo.originalXWGTUP(); auto const& weightInfos = *luminosityBlockCache(iEvent.getLuminosityBlock().index()); - addWeightGroupToTable(lheWeightTables, "LHE", weightInfos.first, lheWeights, w0); - addWeightGroupToTable(lheWeightTables, "GEN", weightInfos.second, genWeights, w0); + addWeightGroupToTable(lheWeightTables, "LHE", weightInfos.at(inLHE), lheWeights); + addWeightGroupToTable(lheWeightTables, "Gen", weightInfos.at(inGen), genWeights); iEvent.put(std::move(lheWeightTables)); } void addWeightGroupToTable(std::unique_ptr>& lheWeightTables, const char* typeName, - const std::vector& weightInfos, WeightsContainer& lheWeights, double w0) const { + const WeightGroupDataContainer& weightInfos, WeightsContainer& allWeights) const { size_t typeCount = 0; gen::WeightType previousType = gen::WeightType::kUnknownWeights; @@ -76,9 +73,23 @@ class LHEWeightsTableProducer : public edm::global::EDProducerweightType(); if (previousType != weightType) typeCount = 0; + std::string name = weightTypeNames_.at(weightType); - auto weights = weightType != gen::WeightType::kScaleWeights ? normalizedWeights(lheWeights, groupInfo, w0) : - orderedScaleWeights(lheWeights, groupInfo, w0); + std::string label = groupInfo.group->name(); + + auto& weights = allWeights.at(groupInfo.index); + label.append("; "); + if (weightType == gen::WeightType::kScaleWeights && groupInfo.group->isWellFormed() + && groupInfo.group->nIdsContained() < 10) { + weights = orderedScaleWeights(weights, + dynamic_cast(groupInfo.group)); + label.append("[1] is mur=0.5 muf=1; [2] is mur=0.5 muf=2; [3] is mur=1 muf=0.5 ;" + " [4] is mur=1 muf=1; [5] is mur=1 muf=2; [6] is mur=2 muf=0.5;" + " [7] is mur=2 muf=1 ; [8] is mur=2 muf=2)"); + } + else + label.append(groupInfo.group->description()); + entryName.append(name); entryName.append("Weight"); if (typeCount > 0) { @@ -88,7 +99,7 @@ class LHEWeightsTableProducer : public edm::global::EDProduceremplace_back(weights.size(), entryName, false); lheWeightTables->back().addColumn( - "", weights, groupInfo.group->name(), nanoaod::FlatTable::FloatColumn, lheWeightPrecision_); + "", weights, label, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_); typeCount++; previousType = weightType; @@ -112,33 +123,25 @@ class LHEWeightsTableProducer : public edm::global::EDProducer(weightsToStore); } - std::vector weightDataPerType(edm::Handle& weightsHandle, + WeightGroupDataContainer weightDataPerType(edm::Handle& weightsHandle, gen::WeightType weightType, int& maxStore) const { - std::vector group; - if (weightType == gen::WeightType::kPdfWeights) { - for (auto lhaid : pdfIds_) { - if (auto pdfGroup = weightsHandle->pdfGroupWithIndexByLHAID(lhaid)) { - group.push_back(pdfGroup.value()); - maxStore--; - if (maxStore == 0) - break; - } - } - return group; + WeightGroupDataContainer group; + if (weightType == gen::WeightType::kPdfWeights && pdfIds_.size() > 0) { + group = weightsHandle->pdfGroupsWithIndicesByLHAIDs(pdfIds_); } - - group = weightsHandle->weightGroupsAndIndicesByType(weightType); + else + group = weightsHandle->weightGroupsAndIndicesByType(weightType); if (maxStore < 0 || static_cast(group.size()) <= maxStore) { // Modify size in case one type of weight is present in multiple products @@ -149,33 +152,20 @@ class LHEWeightsTableProducer : public edm::global::EDProducer normalizedWeights(WeightsContainer& lheWeights, const gen::WeightGroupData& meGroupInfo, double w0) const { - std::vector normalizedWeights; - for (const auto& weight : lheWeights.at(meGroupInfo.index)) - normalizedWeights.push_back(weight/w0); - return normalizedWeights; - } + std::vector orderedScaleWeights(const std::vector& scaleWeights, const gen::ScaleWeightGroupInfo* scaleGroup) const { + + std::vector weights; + weights.push_back(scaleWeights.at(scaleGroup->muR05muF05Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR05muF1Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR05muF2Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR1muF05Index())); + weights.push_back(scaleWeights.at(scaleGroup->centralIndex())); + weights.push_back(scaleWeights.at(scaleGroup->muR1muF2Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR2muF05Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR2muF1Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR2muF2Index())); - std::vector orderedScaleWeights(WeightsContainer& lheWeights, const gen::WeightGroupData& scaleGroupInfo, double w0) const { - const gen::ScaleWeightGroupInfo* scaleGroup = static_cast(scaleGroupInfo.group); - size_t scaleGroupIndex = scaleGroupInfo.index; - - std::vector normalizedWeights; - std::vector& scaleWeights = lheWeights.at(scaleGroupIndex); - - // nano ordering of mur=0.5 muf=0.5 ; [1] is mur=0.5 muf=1 ; [2] is mur=0.5 muf=2 ; [3] is mur=1 muf=0.5 ; - // [4] is mur=1 muf=1 ; [5] is mur=1 muf=2 ; [6] is mur=2 muf=0.5 ; [7] is mur=2 muf=1 ; [8] is mur=2 muf=2 * - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR05muF05Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR05muF1Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR05muF2Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR1muF05Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->centralIndex())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR1muF2Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR2muF05Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR2muF1Index())/w0); - normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR2muF2Index())/w0); - - return normalizedWeights; + return weights; } // nothing to do here @@ -183,8 +173,6 @@ class LHEWeightsTableProducer : public edm::global::EDProducer("lheInfo", {"externalLHEProducer"}) - ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)"); desc.add("lheWeights"); desc.add("genWeights"); desc.add>("weightgroups"); @@ -213,6 +201,7 @@ class LHEWeightsTableProducer : public edm::global::EDProducer weightGroupIndices_; int lheWeightPrecision_; + enum {inLHE, inGen}; }; #include "FWCore/Framework/interface/MakerMacros.h" diff --git a/PhysicsTools/NanoAOD/python/nanogen_cff.py b/PhysicsTools/NanoAOD/python/nanogen_cff.py index b86f63905d6c9..aed174dbf7a2b 100644 --- a/PhysicsTools/NanoAOD/python/nanogen_cff.py +++ b/PhysicsTools/NanoAOD/python/nanogen_cff.py @@ -6,11 +6,11 @@ from PhysicsTools.NanoAOD.lheInfoTable_cfi import * from PhysicsTools.NanoAOD.genWeightsTable_cfi import * -genWeights = cms.EDProducer("GenWeightProductProducer") +genWeights = cms.EDProducer("GenWeightProductProducer", + genInfo = cms.InputTag("generator")) lheWeightsTable = cms.EDProducer( "LHEWeightsTableProducer", - lheInfo = cms.InputTag("externalLHEProducer"), lheWeights = cms.InputTag("externalLHEProducer"), genWeights = cms.InputTag("genWeights"), # Warning: you can use a full string, but only the first character is read. @@ -18,6 +18,7 @@ # must be lower case and 'PDF' must be capital weightgroups = cms.vstring(['scale', 'PDF', 'matrix element', 'unknown', 'parton shower']), maxGroupsPerType = cms.vint32([1, -1, 1, 2, 1]), + # If empty or not specified, all pdf weights are stored pdfIds = cms.vint32([91400, 306000, 260000]), lheWeightPrecision = cms.int32(14), ) diff --git a/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h b/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h index 06a3ff7907287..f990ed1931ea5 100644 --- a/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h +++ b/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h @@ -36,7 +36,8 @@ class GenWeightInfoProduct { std::vector weightGroupsByType(gen::WeightType type) const; std::vector weightGroupIndicesByType(gen::WeightType type) const; std::vector weightGroupsAndIndicesByType(gen::WeightType type) const; - std::optional pdfGroupWithIndexByLHAID(size_t lhaid) const; + std::optional pdfGroupWithIndexByLHAID(int lhaid) const; + std::vector pdfGroupsWithIndicesByLHAIDs(const std::vector& lhaids) const; void addWeightGroupInfo(gen::WeightGroupInfo* info); const int numberOfGroups() const { return weightGroupsInfo_.size(); } diff --git a/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h b/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h index e09f2541f8601..4675e17c6c80d 100644 --- a/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h +++ b/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h @@ -13,10 +13,12 @@ typedef std::vector> WeightsContainer; class GenWeightProduct { public: - GenWeightProduct() { weightsVector_ = {}; } + GenWeightProduct() { weightsVector_ = {}; centralWeight_ = 1.; } + GenWeightProduct(double w0) { weightsVector_ = {}; centralWeight_ = w0; } GenWeightProduct& operator=(GenWeightProduct&& other) { - weightsVector_ = std::move(other.weightsVector_); - return *this; + weightsVector_ = std::move(other.weightsVector_); + centralWeight_ = other.centralWeight_; + return *this; } ~GenWeightProduct() {} @@ -31,12 +33,14 @@ class GenWeightProduct { if (static_cast(weights.size()) <= weightNum) { weights.resize(weightNum+1); } - weights[weightNum] = weight; + weights[weightNum] = weight/centralWeight_; } const WeightsContainer& weights() const { return weightsVector_; } + double centralWeight() const { return centralWeight_; } private: WeightsContainer weightsVector_; + double centralWeight_; }; #endif // GeneratorEvent_LHEInterface_GenWeightProduct_h diff --git a/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h index fc71436bd5ed6..aa097b399cf12 100644 --- a/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h +++ b/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h @@ -25,7 +25,7 @@ namespace gen { kUnknownWeights = 'u', }; - const std::array allGenWeightTypes = {{WeightType::kPdfWeights, + const std::array allWeightTypes = {{WeightType::kPdfWeights, WeightType::kScaleWeights, WeightType::kMEParamWeights, WeightType::kPartonShowerWeights, WeightType::kUnknownWeights, }}; @@ -50,22 +50,25 @@ namespace gen { WeightMetaInfo weightMetaInfo(int weightEntry) const; WeightMetaInfo weightMetaInfoByGlobalIndex(std::string wgtId, int weightEntry) const; int weightVectorEntry(const std::string& wgtId) const; - int containsWeight(const std::string& wgtId, int weightEntry) const; + bool containsWeight(const std::string& wgtId, int weightEntry) const; int weightVectorEntry(const std::string& wgtId, int weightEntry) const; void addContainedId(int weightEntry, std::string id, std::string label); std::vector containedIds() const; bool indexInRange(int index) const; void setName(std::string name) { name_ = name; } + void setDescription(std::string description) { description_ = description; } void setHeaderEntry(std::string header) { headerEntry_ = header; } void setWeightType(WeightType type) { weightType_ = type; } void setFirstId(int firstId) { firstId_ = firstId; } void setLastId(int lastId) { lastId_ = lastId; } std::string name() const { return name_; } + std::string description() const { return description_; } std::string headerEntry() const { return headerEntry_; } WeightType weightType() const { return weightType_; } std::vector idsContained() const { return idsContained_; } + size_t nIdsContained() const { return idsContained_.size(); } int firstId() const { return firstId_; } int lastId() const { return lastId_; } // Store whether the group was fully parsed succesfully @@ -76,6 +79,7 @@ namespace gen { bool isWellFormed_; std::string headerEntry_; std::string name_; + std::string description_; WeightType weightType_; std::vector idsContained_; int firstId_; diff --git a/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc b/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc index 2b08fce3ef09e..e16bbebcc27da 100644 --- a/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc +++ b/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc @@ -54,7 +54,7 @@ std::vector GenWeightInfoProduct::weightGroupsByType(gen: return matchingGroups; } -std::optional GenWeightInfoProduct::pdfGroupWithIndexByLHAID(size_t lhaid) const { +std::optional GenWeightInfoProduct::pdfGroupWithIndexByLHAID(int lhaid) const { auto pdfGroups = weightGroupsAndIndicesByType(gen::WeightType::kPdfWeights); auto matchingPdfSet = std::find_if(pdfGroups.begin(), pdfGroups.end(), @@ -65,6 +65,22 @@ std::optional GenWeightInfoProduct::pdfGroupWithIndexByLHA return matchingPdfSet != pdfGroups.end() ? std::optional(*matchingPdfSet) : std::nullopt; } +std::vector GenWeightInfoProduct::pdfGroupsWithIndicesByLHAIDs(const std::vector& lhaids) const { + auto pdfGroups = weightGroupsAndIndicesByType(gen::WeightType::kPdfWeights); + std::vector groups; + + for (auto lhaid : lhaids) { + auto matchingPdfSet = std::find_if(pdfGroups.begin(), pdfGroups.end(), + [lhaid] (gen::WeightGroupData& data) { + auto pdfGroup = dynamic_cast(data.group); + return pdfGroup->containsLhapdfId(lhaid); + }); + if (matchingPdfSet != pdfGroups.end()) + pdfGroups.push_back(*matchingPdfSet); + } + + return pdfGroups; +} std::vector GenWeightInfoProduct::weightGroupIndicesByType(gen::WeightType type) const { std::vector matchingGroupIndices; diff --git a/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc index aeb877bfb7b29..23c22b2342838 100644 --- a/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc +++ b/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc @@ -10,6 +10,7 @@ namespace gen { isWellFormed_ = true; headerEntry_ = other.headerEntry(); name_ = other.name(); + description_ = other.description(); weightType_ = other.weightType(); idsContained_ = other.idsContained(); firstId_ = other.firstId(); @@ -37,7 +38,7 @@ namespace gen { return weightVectorEntry(wgtId, 0); } - int WeightGroupInfo::containsWeight(const std::string& wgtId, int weightEntry) const { + bool WeightGroupInfo::containsWeight(const std::string& wgtId, int weightEntry) const { return weightVectorEntry(wgtId, weightEntry) != -1; }