diff --git a/DataFormats/NanoAOD/src/classes_def.xml b/DataFormats/NanoAOD/src/classes_def.xml index fb6c23b3646e9..ce490c5b6a357 100644 --- a/DataFormats/NanoAOD/src/classes_def.xml +++ b/DataFormats/NanoAOD/src/classes_def.xml @@ -9,6 +9,9 @@ + + + diff --git a/GeneratorInterface/Core/BuildFile.xml b/GeneratorInterface/Core/BuildFile.xml index 2d3cd1bb76d7f..eb5953775844c 100644 --- a/GeneratorInterface/Core/BuildFile.xml +++ b/GeneratorInterface/Core/BuildFile.xml @@ -10,6 +10,7 @@ + diff --git a/GeneratorInterface/Core/interface/GenWeightHelper.h b/GeneratorInterface/Core/interface/GenWeightHelper.h new file mode 100644 index 0000000000000..1a1724b36792e --- /dev/null +++ b/GeneratorInterface/Core/interface/GenWeightHelper.h @@ -0,0 +1,26 @@ +#ifndef GeneratorInterface_Core_GenWeightHelper_h +#define GeneratorInterface_Core_GenWeightHelper_h + +#include +#include +#include +#include +#include + +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoProduct.h" +#include "GeneratorInterface/Core/interface/WeightHelper.h" + +#include + +namespace gen { + class GenWeightHelper : public WeightHelper { + public: + GenWeightHelper(); + void parseWeightGroupsFromNames(std::vector weightNames); + }; +} // namespace gen + +#endif diff --git a/GeneratorInterface/Core/interface/LHEWeightHelper.h b/GeneratorInterface/Core/interface/LHEWeightHelper.h new file mode 100644 index 0000000000000..0418823095f2e --- /dev/null +++ b/GeneratorInterface/Core/interface/LHEWeightHelper.h @@ -0,0 +1,45 @@ +#ifndef GeneratorInterface_Core_LHEWeightHelper_h +#define GeneratorInterface_Core_LHEWeightHelper_h + +#include +#include +#include +#include +#include + +#include "SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" +#include "GeneratorInterface/Core/interface/WeightHelper.h" + +#include + +namespace gen { + class LHEWeightHelper : public WeightHelper { + public: + LHEWeightHelper() : WeightHelper(){}; + + enum ErrorType { SWAPHEADER, HTMLSTYLE, TRAILINGSTR, UNKNOWN }; + + void setHeaderLines(std::vector headerLines); + void parseWeights(); + bool isConsistent(); + void swapHeaders(); + void setFailIfInvalidXML(bool value) { failIfInvalidXML_ = value; } + + private: + std::vector headerLines_; + std::string weightgroupKet_ = ""; + bool failIfInvalidXML_ = false; + std::string parseGroupName(tinyxml2::XMLElement* el); + void addGroup(tinyxml2::XMLElement* inner, std::string groupName, int groupIndex, int& weightIndex); + bool parseLHE(tinyxml2::XMLDocument& xmlDoc); + tinyxml2::XMLError tryReplaceHtmlStyle(tinyxml2::XMLDocument& xmlDoc, std::string& fullHeader); + tinyxml2::XMLError tryRemoveTrailings(tinyxml2::XMLDocument& xmlDoc, std::string& fullHeader); + ErrorType findErrorType(std::string& fullHeader); + }; +} // namespace gen + +#endif diff --git a/GeneratorInterface/Core/interface/WeightHelper.h b/GeneratorInterface/Core/interface/WeightHelper.h new file mode 100644 index 0000000000000..517e57181d650 --- /dev/null +++ b/GeneratorInterface/Core/interface/WeightHelper.h @@ -0,0 +1,117 @@ +#ifndef GeneratorInterface_LHEInterface_WeightHelper_h +#define GeneratorInterface_LHEInterface_WeightHelper_h + +#include "DataFormats/Common/interface/OwnVector.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/WeightsInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h" +#include "LHAPDF/LHAPDF.h" +#include +#include +#include + +namespace gen { + struct ParsedWeight { + std::string id; + int index; + std::string groupname; + std::string content; + std::unordered_map attributes; + int wgtGroup_idx; + }; + + class WeightHelper { + public: + WeightHelper(); + edm::OwnVector weightGroups() { return weightGroups_; } + + template + std::unique_ptr weightProduct(std::vector weights, float w0); + + void setModel(std::string model) { model_ = model; } + void setGuessPSWeightIdx(bool guessPSWeightIdx) { + PartonShowerWeightGroupInfo::setGuessPSWeightIdx(guessPSWeightIdx); + } + void addUnassociatedGroup() { + weightGroups_.push_back(std::make_unique("unassociated")); + weightGroups_.back().setDescription("Weights with missing or invalid header meta data"); + } + int addWeightToProduct( + std::unique_ptr& product, double weight, std::string name, int weightNum, int groupIndex); + int findContainingWeightGroup(std::string wgtId, int weightIndex, int previousGroupIndex); + void setDebug(bool value) { debug_ = value; } + + protected: + // TODO: Make this only print from one thread a la + // https://github.com/kdlong/cmssw/blob/master/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc#L1069 + bool debug_ = false; + const unsigned int FIRST_PSWEIGHT_ENTRY = 2; + const unsigned int DEFAULT_PSWEIGHT_LENGTH = 46; + std::string model_; + std::vector parsedWeights_; + std::map currWeightAttributeMap_; + std::map currGroupAttributeMap_; + edm::OwnVector weightGroups_; + bool isScaleWeightGroup(const ParsedWeight& weight); + bool isMEParamWeightGroup(const ParsedWeight& weight); + bool isPdfWeightGroup(const ParsedWeight& weight); + bool isPartonShowerWeightGroup(const ParsedWeight& weight); + bool isOrphanPdfWeightGroup(ParsedWeight& weight); + void updateScaleInfo(gen::ScaleWeightGroupInfo& scaleGroup, const ParsedWeight& weight); + void updateMEParamInfo(const ParsedWeight& weight, int index); + void updatePdfInfo(gen::PdfWeightGroupInfo& pdfGroup, const ParsedWeight& weight); + void updatePartonShowerInfo(gen::PartonShowerWeightGroupInfo& psGroup, const ParsedWeight& weight); + void cleanupOrphanCentralWeight(); + bool splitPdfWeight(ParsedWeight& weight); + + int lhapdfId(const ParsedWeight& weight, gen::PdfWeightGroupInfo& pdfGroup); + std::string searchAttributes(const std::string& label, const ParsedWeight& weight) const; + std::string searchAttributesByTag(const std::string& label, const ParsedWeight& weight) const; + std::string searchAttributesByRegex(const std::string& label, const ParsedWeight& weight) const; + + // Possible names for the same thing + const std::unordered_map> attributeNames_ = { + {"muf", {"muF", "MUF", "muf", "facscfact"}}, + {"mur", {"muR", "MUR", "mur", "renscfact"}}, + {"pdf", {"PDF", "PDF set", "lhapdf", "pdf", "pdf set", "pdfset"}}, + {"dyn", {"DYN_SCALE"}}, + {"dyn_name", {"dyn_scale_choice"}}, + {"up", {"_up", "Hi"}}, + {"down", {"_dn", "Lo"}}, + {"me_variation", {"mass", "sthw2", "width"}}, + }; + void printWeights(); + std::unique_ptr buildGroup(ParsedWeight& weight); + void buildGroups(); + std::string searchString(const std::string& label, const std::string& name); + }; + + // Templated function (needed here because of plugins) + template + 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; + // This happens if there are no PS weights, so the weights vector contains only the central GEN weight. + // Just add an empty product (need for all cases or...?) + if (weights.size() > 1) { + for (const auto& weight : weights) { + if constexpr (std::is_same::value) + weightGroupIndex = addWeightToProduct(weightProduct, weight.wgt, weight.id, i, weightGroupIndex); + else if (std::is_same::value) + weightGroupIndex = addWeightToProduct(weightProduct, weight, std::to_string(i), i, weightGroupIndex); + i++; + } + } + return weightProduct; + } +} // namespace gen + +#endif diff --git a/GeneratorInterface/Core/plugins/BuildFile.xml b/GeneratorInterface/Core/plugins/BuildFile.xml index 40c9678cb4bc6..a8c9a22737fc6 100644 --- a/GeneratorInterface/Core/plugins/BuildFile.xml +++ b/GeneratorInterface/Core/plugins/BuildFile.xml @@ -3,6 +3,7 @@ + diff --git a/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc b/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc new file mode 100644 index 0000000000000..198917def4506 --- /dev/null +++ b/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc @@ -0,0 +1,95 @@ +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDProducer.h" +#include "FWCore/Framework/interface/LuminosityBlock.h" + +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h" + +#include "GeneratorInterface/Core/interface/GenWeightHelper.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include + +class GenWeightProductProducer : public edm::one::EDProducer { +public: + explicit GenWeightProductProducer(const edm::ParameterSet& iConfig); + ~GenWeightProductProducer() override; + +private: + std::vector weightNames_; + gen::GenWeightHelper weightHelper_; + edm::EDGetTokenT genLumiInfoToken_; + edm::EDGetTokenT genEventToken_; + const edm::EDGetTokenT genLumiInfoHeadTag_; + + void produce(edm::Event&, const edm::EventSetup&) override; + void beginLuminosityBlockProduce(edm::LuminosityBlock& lb, edm::EventSetup const& c) override; +}; + +// +// constructors and destructor +// +GenWeightProductProducer::GenWeightProductProducer(const edm::ParameterSet& iConfig) + : genLumiInfoToken_(consumes(iConfig.getParameter("genInfo"))), + genEventToken_(consumes(iConfig.getParameter("genInfo"))), + genLumiInfoHeadTag_( + mayConsume(iConfig.getParameter("genLumiInfoHeader"))) { + weightHelper_.setDebug(iConfig.getUntrackedParameter("debug", false)); + produces(); + produces(); + weightHelper_.setGuessPSWeightIdx(iConfig.getUntrackedParameter("guessPSWeightIdx", false)); +} + +GenWeightProductProducer::~GenWeightProductProducer() {} + +// ------------ method called to produce the data ------------ +void GenWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + edm::Handle genEventInfo; + iEvent.getByToken(genEventToken_, genEventInfo); + + float centralWeight = !genEventInfo->weights().empty() ? genEventInfo->weights().at(0) : 1.; + auto weightProduct = weightHelper_.weightProduct(genEventInfo->weights(), centralWeight); + iEvent.put(std::move(weightProduct)); +} + +void GenWeightProductProducer::beginLuminosityBlockProduce(edm::LuminosityBlock& iLumi, edm::EventSetup const& iSetup) { + edm::Handle genLumiInfoHead; + iLumi.getByToken(genLumiInfoHeadTag_, genLumiInfoHead); + if (genLumiInfoHead.isValid()) { + std::string label = genLumiInfoHead->configDescription(); + boost::replace_all(label, "-", "_"); + weightHelper_.setModel(label); + } + + edm::Handle genLumiInfoHandle; + iLumi.getByToken(genLumiInfoToken_, genLumiInfoHandle); + + weightNames_ = genLumiInfoHandle->weightNames(); + weightHelper_.parseWeightGroupsFromNames(weightNames_); + + auto weightInfoProduct = std::make_unique(); + if (weightHelper_.weightGroups().empty()) + weightHelper_.addUnassociatedGroup(); + + for (auto& weightGroup : weightHelper_.weightGroups()) { + weightInfoProduct->addWeightGroupInfo(std::unique_ptr(weightGroup.clone())); + } + iLumi.put(std::move(weightInfoProduct)); +} + +DEFINE_FWK_MODULE(GenWeightProductProducer); diff --git a/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc b/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc new file mode 100644 index 0000000000000..daa7b26bb874d --- /dev/null +++ b/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc @@ -0,0 +1,145 @@ +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDProducer.h" +#include "FWCore/Framework/interface/LuminosityBlock.h" + +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" + +#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h" +#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h" +#include "GeneratorInterface/Core/interface/LHEWeightHelper.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/transform.h" + +class LHEWeightProductProducer : public edm::one::EDProducer { +public: + explicit LHEWeightProductProducer(const edm::ParameterSet& iConfig); + ~LHEWeightProductProducer() override; + +private: + gen::LHEWeightHelper weightHelper_; + std::vector lheLabels_; + std::vector> lheEventTokens_; + std::vector> lheRunInfoTokens_; + std::vector> lheWeightInfoTokens_; + bool foundWeightProduct_ = false; + bool hasLhe_ = false; + + void produce(edm::Event&, const edm::EventSetup&) override; + void beginLuminosityBlockProduce(edm::LuminosityBlock& lumi, edm::EventSetup const& es) override; + void beginRun(edm::Run const& run, edm::EventSetup const& es) override; + void endRun(edm::Run const& run, edm::EventSetup const& es) override; +}; + +// TODO: Accept a vector of strings (source, externalLHEProducer) exit if neither are found +LHEWeightProductProducer::LHEWeightProductProducer(const edm::ParameterSet& iConfig) + : lheLabels_(iConfig.getParameter>("lheSourceLabels")), + lheEventTokens_(edm::vector_transform( + lheLabels_, [this](const std::string& tag) { return mayConsume(tag); })), + lheRunInfoTokens_(edm::vector_transform( + lheLabels_, [this](const std::string& tag) { return mayConsume(tag); })), + lheWeightInfoTokens_(edm::vector_transform( + lheLabels_, [this](const std::string& tag) { return mayConsume(tag); })) { + produces(); + produces(); + weightHelper_.setFailIfInvalidXML(iConfig.getUntrackedParameter("failIfInvalidXML", false)); + weightHelper_.setDebug(iConfig.getUntrackedParameter("debug", false)); + weightHelper_.setGuessPSWeightIdx(iConfig.getUntrackedParameter("guessPSWeightIdx", false)); +} + +LHEWeightProductProducer::~LHEWeightProductProducer() {} + +// ------------ method called to produce the data ------------ +void LHEWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + if (foundWeightProduct_ || !hasLhe_) + return; + + edm::Handle lheEventInfo; + for (auto& token : lheEventTokens_) { + iEvent.getByToken(token, lheEventInfo); + if (lheEventInfo.isValid()) { + break; + } + } + + auto weightProduct = weightHelper_.weightProduct(lheEventInfo->weights(), lheEventInfo->originalXWGTUP()); + iEvent.put(std::move(weightProduct)); +} + +void LHEWeightProductProducer::beginRun(edm::Run const& run, edm::EventSetup const& es) { + edm::Handle lheRunInfoHandle; + for (auto& label : lheLabels_) { + run.getByLabel(label, lheRunInfoHandle); + if (lheRunInfoHandle.isValid()) { + hasLhe_ = true; + break; + } + } + if (!hasLhe_) + return; + + typedef std::vector::const_iterator header_cit; + LHERunInfoProduct::Header headerWeightInfo; + for (header_cit iter = lheRunInfoHandle->headers_begin(); iter != lheRunInfoHandle->headers_end(); iter++) { + if (iter->tag() == "initrwgt") { + headerWeightInfo = *iter; + break; + } + } + + weightHelper_.setHeaderLines(headerWeightInfo.lines()); +} + +void LHEWeightProductProducer::endRun(edm::Run const& run, edm::EventSetup const& es) {} + +void LHEWeightProductProducer::beginLuminosityBlockProduce(edm::LuminosityBlock& lumi, edm::EventSetup const& es) { + edm::Handle lheWeightInfoHandle; + + for (auto& token : lheWeightInfoTokens_) { + lumi.getByToken(token, lheWeightInfoHandle); + if (lheWeightInfoHandle.isValid()) { + foundWeightProduct_ = true; + return; + } + } + + if (!hasLhe_) + return; + + try { + weightHelper_.parseWeights(); + } catch (cms::Exception& e) { + std::string error = e.what(); + error += + "\n NOTE: if you want to attempt to process this sample anyway, set failIfInvalidXML = False " + "in the configuration file\n. If you set this flag and the error persists, the issue " + " is fatal and must be solved at the LHE/gridpack level."; + throw cms::Exception("LHEWeightProductProducer") << error; + } + + if (weightHelper_.weightGroups().empty()) + weightHelper_.addUnassociatedGroup(); + + auto weightInfoProduct = std::make_unique(); + for (auto& weightGroup : weightHelper_.weightGroups()) { + weightInfoProduct->addWeightGroupInfo(std::unique_ptr(weightGroup.clone())); + } + lumi.put(std::move(weightInfoProduct)); +} + +DEFINE_FWK_MODULE(LHEWeightProductProducer); diff --git a/GeneratorInterface/Core/src/GenWeightHelper.cc b/GeneratorInterface/Core/src/GenWeightHelper.cc new file mode 100644 index 0000000000000..b06475941efd6 --- /dev/null +++ b/GeneratorInterface/Core/src/GenWeightHelper.cc @@ -0,0 +1,56 @@ +#include "GeneratorInterface/Core/interface/GenWeightHelper.h" +#include + +using namespace tinyxml2; + +namespace gen { + GenWeightHelper::GenWeightHelper() {} + + void GenWeightHelper::parseWeightGroupsFromNames(std::vector weightNames) { + parsedWeights_.clear(); + int index = 0; + int groupIndex = -1; + int showerGroupIndex = -1; + std::string curGroup = ""; + // If size is 1, it's just the central weight + if (weightNames.size() <= 1) + return; + + for (std::string weightName : weightNames) { + if (weightName.find("LHE") != std::string::npos) { + // Parse as usual, this is the SUSY workflow + std::vector info; + boost::split(info, weightName, boost::is_any_of(",")); + std::unordered_map attributes; + std::string text = info.back(); + info.pop_back(); + for (auto i : info) { + std::vector subInfo; + boost::split(subInfo, i, boost::is_any_of("=")); + if (subInfo.size() == 2) { + attributes[boost::algorithm::trim_copy(subInfo[0])] = boost::algorithm::trim_copy(subInfo[1]); + } + } + if (attributes["group"] != curGroup) { + curGroup = attributes["group"]; + groupIndex++; + } + // Gen Weights can't have an ID, because they are just a std::vector in the event + attributes["id"] = ""; + parsedWeights_.push_back({attributes["id"], index, curGroup, text, attributes, groupIndex}); + } else { + parsedWeights_.push_back( + {"", index, weightName, weightName, std::unordered_map(), groupIndex}); + if (isPartonShowerWeightGroup(parsedWeights_.back())) { + if (showerGroupIndex < 0) { + showerGroupIndex = ++groupIndex; + } + parsedWeights_.back().wgtGroup_idx = showerGroupIndex; // all parton showers are grouped together + } + } + index++; + } + buildGroups(); + printWeights(); + } +} // namespace gen diff --git a/GeneratorInterface/Core/src/LHEWeightHelper.cc b/GeneratorInterface/Core/src/LHEWeightHelper.cc new file mode 100644 index 0000000000000..290b8fd77f268 --- /dev/null +++ b/GeneratorInterface/Core/src/LHEWeightHelper.cc @@ -0,0 +1,186 @@ +#include "GeneratorInterface/Core/interface/LHEWeightHelper.h" +#include +#include +#include +#include + +using namespace tinyxml2; + +namespace gen { + void LHEWeightHelper::setHeaderLines(std::vector headerLines) { headerLines_ = headerLines; } + + bool LHEWeightHelper::parseLHE(tinyxml2::XMLDocument& xmlDoc) { + parsedWeights_.clear(); + + std::string fullHeader = boost::algorithm::join(headerLines_, ""); + if (debug_) + std::cout << "Full header is \n" << fullHeader << std::endl; + int xmlError = xmlDoc.Parse(fullHeader.c_str()); + + while (!isConsistent() || xmlError != 0) { + if (failIfInvalidXML_) { + xmlDoc.PrintError(); + throw cms::Exception("LHEWeightHelper") + << "The LHE header is not valid XML! Weight information was not properly parsed."; + } + + switch (findErrorType(fullHeader)) { + case ErrorType::SWAPHEADER: + swapHeaders(); + fullHeader = boost::algorithm::join(headerLines_, ""); + xmlError = xmlDoc.Parse(fullHeader.c_str()); + break; + case ErrorType::HTMLSTYLE: + xmlError = tryReplaceHtmlStyle(xmlDoc, fullHeader); + break; + case ErrorType::TRAILINGSTR: + xmlError = tryRemoveTrailings(xmlDoc, fullHeader); + break; + case ErrorType::UNKNOWN: + std::string error = + "Fatal error when parsing the LHE header. The header is not valid XML! Parsing error was "; + error += xmlDoc.ErrorStr(); + throw cms::Exception("LHEWeightHelper") << error; + } + } + + return true; + } + + void LHEWeightHelper::addGroup(tinyxml2::XMLElement* inner, std::string groupName, int groupIndex, int& weightIndex) { + if (debug_) + std::cout << " >> Found a weight inside the group. " << std::endl; + std::string text = ""; + if (inner->GetText()) + text = inner->GetText(); + + std::unordered_map attributes; + for (auto* att = inner->FirstAttribute(); att != nullptr; att = att->Next()) + attributes[att->Name()] = att->Value(); + if (debug_) + std::cout << " " << weightIndex << ": \"" << text << "\"" << std::endl; + parsedWeights_.push_back({inner->Attribute("id"), weightIndex++, groupName, text, attributes, groupIndex}); + } + + void LHEWeightHelper::parseWeights() { + tinyxml2::XMLDocument xmlDoc; + if (!parseLHE(xmlDoc)) { + return; + } + + int weightIndex = 0; + int groupIndex = 0; + for (auto* e = xmlDoc.RootElement(); e != nullptr; e = e->NextSiblingElement()) { + if (debug_) + std::cout << "XML element is " << e->Name() << std::endl; + std::string groupName = ""; + if (strcmp(e->Name(), "weight") == 0) { + if (debug_) + std::cout << "Found weight unmatched to group\n"; + addGroup(e, groupName, groupIndex, weightIndex); + } else if (strcmp(e->Name(), "weightgroup") == 0) { + groupName = parseGroupName(e); + if (debug_) + std::cout << ">>>> Found a weight group: " << groupName << std::endl; + for (auto inner = e->FirstChildElement("weight"); inner != nullptr; inner = inner->NextSiblingElement("weight")) + addGroup(inner, groupName, groupIndex, weightIndex); + } else + std::cout << "Found an invalid entry\n"; + groupIndex++; + } + buildGroups(); + if (debug_) + printWeights(); + } + + std::string LHEWeightHelper::parseGroupName(tinyxml2::XMLElement* el) { + std::vector nameAlts_ = {"name", "type"}; + for (const auto& nameAtt : nameAlts_) { + if (el->Attribute(nameAtt.c_str())) { + std::string groupName = el->Attribute(nameAtt.c_str()); + if (groupName.find('.') != std::string::npos) + groupName.erase(groupName.find('.'), groupName.size()); + return groupName; + } + } + bool hardFail = true; + if (hardFail) { + throw std::runtime_error("couldn't find groupname"); + } + return ""; + } + + bool LHEWeightHelper::isConsistent() { + int curLevel = 0; + + for (const auto& line : headerLines_) { + if (line.find("") != std::string::npos) { + curLevel--; + if (curLevel != 0) { + return false; + } + } + } + return curLevel == 0; + } + + void LHEWeightHelper::swapHeaders() { + int curLevel = 0; + int open = -1; + int close = -1; + for (size_t idx = 0; idx < headerLines_.size(); idx++) { + std::string line = headerLines_[idx]; + if (line.find("") != std::string::npos) { + curLevel--; + if (curLevel != 0) { + close = idx; + } + } + if (open > -1 && close > -1) { + std::swap(headerLines_[open], headerLines_[close]); + open = -1; + close = -1; + } + } + } + + tinyxml2::XMLError LHEWeightHelper::tryReplaceHtmlStyle(tinyxml2::XMLDocument& xmlDoc, std::string& fullHeader) { + // in case of > instead of < + boost::replace_all(fullHeader, "<", "<"); + boost::replace_all(fullHeader, ">", ">"); + + return xmlDoc.Parse(fullHeader.c_str()); + } + + tinyxml2::XMLError LHEWeightHelper::tryRemoveTrailings(tinyxml2::XMLDocument& xmlDoc, std::string& fullHeader) { + // delete extra strings after the last (occasionally contain '<' or '>') + std::size_t theLastKet = fullHeader.rfind(weightgroupKet_) + weightgroupKet_.length(); + fullHeader = fullHeader.substr(0, theLastKet); + + return xmlDoc.Parse(fullHeader.c_str()); + } + + LHEWeightHelper::ErrorType LHEWeightHelper::findErrorType(std::string& fullHeader) { + if (!isConsistent()) + return LHEWeightHelper::ErrorType::SWAPHEADER; + if (fullHeader.find("<") != std::string::npos || fullHeader.find(">") != std::string::npos) + return LHEWeightHelper::ErrorType::HTMLSTYLE; + + std::string trailingCand = + fullHeader.substr(fullHeader.rfind(weightgroupKet_) + std::string(weightgroupKet_).length()); + if (trailingCand.find('<') != std::string::npos || trailingCand.find('>') != std::string::npos) + return LHEWeightHelper::ErrorType::TRAILINGSTR; + + return LHEWeightHelper::ErrorType::UNKNOWN; + } +} // namespace gen diff --git a/GeneratorInterface/Core/src/WeightHelper.cc b/GeneratorInterface/Core/src/WeightHelper.cc new file mode 100644 index 0000000000000..234f041ec0cfa --- /dev/null +++ b/GeneratorInterface/Core/src/WeightHelper.cc @@ -0,0 +1,354 @@ +#include "GeneratorInterface/Core/interface/WeightHelper.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/Exception.h" +#include + +namespace gen { + WeightHelper::WeightHelper() { model_ = ""; } + + bool WeightHelper::isScaleWeightGroup(const ParsedWeight& weight) { + return (weight.groupname.find("scale_variation") != std::string::npos || + weight.groupname.find("Central scale variation") != std::string::npos); + } + + bool WeightHelper::isPdfWeightGroup(const ParsedWeight& weight) { + const std::string& name = weight.groupname; + if (name.find("PDF_variation") != std::string::npos) + return true; + return LHAPDF::lookupLHAPDFID(name) != -1; + } + + bool WeightHelper::isPartonShowerWeightGroup(const ParsedWeight& weight) { + const std::string& groupname = boost::to_lower_copy(weight.groupname); + std::vector psNames = {"isr", "fsr", "nominal", "baseline", "emission"}; + for (const auto& name : psNames) { + if (groupname.find(name) != std::string::npos) + return true; + } + return false; + } + + bool WeightHelper::isOrphanPdfWeightGroup(ParsedWeight& weight) { + std::pair pairLHA; + try { + pairLHA = LHAPDF::lookupPDF(stoi(searchAttributes("pdf", weight))); + } catch (...) { + return false; + } + + if (!pairLHA.first.empty() && pairLHA.second == 0) { + weight.groupname = std::string(pairLHA.first); + return true; + } else { + return false; + } + } + + bool WeightHelper::isMEParamWeightGroup(const ParsedWeight& weight) { + return (weight.groupname.find("mg_reweighting") != std::string::npos || + weight.groupname.find("variation") != std::string::npos); + // variation used for blanket of all variations, might need to change + } + + std::string WeightHelper::searchAttributes(const std::string& label, const ParsedWeight& weight) const { + std::string attribute = searchAttributesByTag(label, weight); + return attribute.empty() ? searchAttributesByRegex(label, weight) : attribute; + } + + std::string WeightHelper::searchAttributesByTag(const std::string& label, const ParsedWeight& weight) const { + auto& attributes = weight.attributes; + for (const auto& lab : attributeNames_.at(label)) { + if (attributes.find(lab) != attributes.end()) { + return boost::algorithm::trim_copy_if(attributes.at(lab), boost::is_any_of("\"")); + } + } + return ""; + } + + std::string WeightHelper::searchAttributesByRegex(const std::string& label, const ParsedWeight& weight) const { + auto& content = weight.content; + std::smatch match; + for (const auto& lab : attributeNames_.at(label)) { + std::regex floatExpr(lab + "\\s*=\\s*([0-9.]+(?:[eE][+-]?[0-9]+)?)"); + std::regex strExpr(lab + "\\s*=\\s*([^=]+)"); + if (std::regex_search(content, match, floatExpr)) { + return boost::algorithm::trim_copy(match.str(1)); + } else if (std::regex_search(content, match, strExpr)) { + return boost::algorithm::trim_copy(match.str(1)); + } + } + return ""; + } + + void WeightHelper::updateScaleInfo(gen::ScaleWeightGroupInfo& scaleGroup, const ParsedWeight& weight) { + std::string muRText = searchAttributes("mur", weight); + std::string muFText = searchAttributes("muf", weight); + std::string dynNumText = searchAttributes("dyn", weight); + float muR, muF; + try { + muR = std::stof(muRText); + muF = std::stof(muFText); + } catch (std::invalid_argument& e) { + if (debug_) + std::cout << "Tried to convert (" << muR << ", " << muF << ") to a int" << std::endl; + scaleGroup.setWeightIsCorrupt(); + return; + /// do something + } + + if (dynNumText.empty()) { + scaleGroup.setMuRMuFIndex(weight.index, weight.id, muR, muF); + } else { + std::string dynType = searchAttributes("dyn_name", weight); + try { + int dynNum = std::stoi(dynNumText); + scaleGroup.setDyn(weight.index, weight.id, muR, muF, dynNum, dynType); + } catch (std::invalid_argument& e) { + std::cout << "Tried to convert (" << dynNumText << ") a int" << std::endl; + scaleGroup.setWeightIsCorrupt(); + /// do something here + } + } + + if (scaleGroup.lhaid() == -1) { + std::string lhaidText = searchAttributes("pdf", weight); + try { + scaleGroup.setLhaid(std::stoi(lhaidText)); + } catch (std::invalid_argument& e) { + scaleGroup.setLhaid(-1); + // do something here + } + } + } + + int WeightHelper::lhapdfId(const ParsedWeight& weight, gen::PdfWeightGroupInfo& pdfGroup) { + std::string lhaidText = searchAttributes("pdf", weight); + if (!lhaidText.empty()) { + try { + return std::stoi(lhaidText); + } catch (std::invalid_argument& e) { + pdfGroup.setIsWellFormed(false); + } + } else if (!pdfGroup.lhaIds().empty()) { + return pdfGroup.lhaIds().back() + 1; + } else { + return LHAPDF::lookupLHAPDFID(weight.groupname); + } + return -1; + } + + void WeightHelper::updatePdfInfo(gen::PdfWeightGroupInfo& pdfGroup, const ParsedWeight& weight) { + int lhaid = lhapdfId(weight, pdfGroup); + if (pdfGroup.parentLhapdfId() < 0) { + int parentId = lhaid - LHAPDF::lookupPDF(lhaid).second; + pdfGroup.setParentLhapdfInfo(parentId); + pdfGroup.setUncertaintyType(gen::kUnknownUnc); + + std::string description = ""; + if (pdfGroup.uncertaintyType() == gen::kHessianUnc) + description += "Hessian "; + else if (pdfGroup.uncertaintyType() == gen::kMonteCarloUnc) + description += "Monte Carlo "; + description += "Uncertainty sets for LHAPDF set " + LHAPDF::lookupPDF(parentId).first; + description += " with LHAID = " + std::to_string(parentId); + description += "; "; + + pdfGroup.appendDescription(description); + } + // after setup parent info, add lhaid + pdfGroup.addLhaid(lhaid); + } + + void WeightHelper::updatePartonShowerInfo(gen::PartonShowerWeightGroupInfo& psGroup, const ParsedWeight& weight) { + if (psGroup.containedIds().size() == DEFAULT_PSWEIGHT_LENGTH) + psGroup.setIsWellFormed(true); + if (weight.content.find(':') != std::string::npos && weight.content.find('=') != std::string::npos) + psGroup.setNameIsPythiaSyntax(true); + } + + bool WeightHelper::splitPdfWeight(ParsedWeight& weight) { + if (weightGroups_[weight.wgtGroup_idx].weightType() == gen::WeightType::kPdfWeights) { + auto& pdfGroup = dynamic_cast(weightGroups_[weight.wgtGroup_idx]); + int lhaid = lhapdfId(weight, pdfGroup); + if (lhaid > 0 && !pdfGroup.isIdInParentSet(lhaid) && pdfGroup.parentLhapdfId() > 0) { + weightGroups_.push_back(*buildGroup(weight)); + weight.wgtGroup_idx++; + return true; + } + } + return false; + } + + void WeightHelper::cleanupOrphanCentralWeight() { + std::vector removeList; + for (auto it = weightGroups_.begin(); it < weightGroups_.end(); it++) { + if (it->weightType() != WeightType::kScaleWeights) + continue; + auto& baseWeight = dynamic_cast(*it); + if (baseWeight.containsCentralWeight()) + continue; + for (auto subIt = weightGroups_.begin(); subIt < it; subIt++) { + if (subIt->weightType() != WeightType::kPdfWeights) + continue; + auto& subWeight = dynamic_cast(*subIt); + if (subWeight.nIdsContained() == 1 && subWeight.parentLhapdfId() == baseWeight.lhaid()) { + removeList.push_back(subIt - weightGroups_.begin()); + auto info = subWeight.idsContained().at(0); + baseWeight.addContainedId(info.globalIndex, info.id, info.label, 1, 1); + } + } + } + std::sort(removeList.begin(), removeList.end(), std::greater()); + for (auto idx : removeList) { + weightGroups_.erase(weightGroups_.begin() + idx); + } + } + + int WeightHelper::addWeightToProduct( + std::unique_ptr& product, double weight, std::string name, int weightNum, int groupIndex) { + bool isUnassociated = false; + try { + groupIndex = findContainingWeightGroup(name, weightNum, groupIndex); + } catch (const cms::Exception& e) { + std::cerr << "WARNING: " << e.what() << std::endl; + isUnassociated = true; + + // initialize index to 0, prevent underflow when accessing vector + groupIndex = 0; + bool foundUnassocGroup = false; + while (!foundUnassocGroup && groupIndex < static_cast(weightGroups_.size())) { + auto& g = weightGroups_[groupIndex]; + if (g.weightType() == gen::WeightType::kUnknownWeights && g.name() == "unassociated") + foundUnassocGroup = true; + else + groupIndex++; + } + if (!foundUnassocGroup) { + addUnassociatedGroup(); + } + } + auto& group = weightGroups_[groupIndex]; + if (isUnassociated) { + group.addContainedId(weightNum, name, name); + } + int entry = !isUnassociated ? group.weightVectorEntry(name, weightNum) : group.nIdsContained(); + if (debug_) + std::cout << "Adding weight " << entry << " to group " << groupIndex; + product->addWeight(weight, groupIndex, entry); + return groupIndex; + } + + int WeightHelper::findContainingWeightGroup(std::string wgtId, int weightIndex, int previousGroupIndex) { + // Start search at previous index, under expectation of ordered weights + previousGroupIndex = previousGroupIndex >= 0 ? previousGroupIndex : 0; + for (int index = previousGroupIndex; index < std::min(index + 1, static_cast(weightGroups_.size())); index++) { + const gen::WeightGroupInfo& weightGroup = weightGroups_[index]; + if (weightGroup.indexInRange(weightIndex) && weightGroup.containsWeight(wgtId, weightIndex)) { + return static_cast(index); + } + } + + // Fall back to unordered search + int counter = 0; + for (const auto& weightGroup : weightGroups_) { + if (weightGroup.containsWeight(wgtId, weightIndex)) + return counter; + counter++; + } + // Needs to be properly handled + throw cms::Exception("Unmatched Generator weight! ID was " + wgtId + " index was " + std::to_string(weightIndex) + + "\nNot found in any of " + std::to_string(weightGroups_.size()) + " weightGroups."); + return -1; + } + + void WeightHelper::printWeights() { + // checks + for (auto& wgt : weightGroups_) { + std::cout << std::boolalpha << wgt.name() << " (" << wgt.firstId() << "-" << wgt.lastId() + << "): " << wgt.isWellFormed() << std::endl; + if (wgt.weightType() == gen::WeightType::kScaleWeights) { + auto& wgtScale = dynamic_cast(wgt); + std::cout << wgtScale.centralIndex() << " "; + std::cout << wgtScale.muR1muF2Index() << " "; + std::cout << wgtScale.muR1muF05Index() << " "; + std::cout << wgtScale.muR2muF1Index() << " "; + std::cout << wgtScale.muR2muF2Index() << " "; + std::cout << wgtScale.muR2muF05Index() << " "; + std::cout << wgtScale.muR05muF1Index() << " "; + std::cout << wgtScale.muR05muF2Index() << " "; + std::cout << wgtScale.muR05muF05Index() << " \n"; + for (auto name : wgtScale.dynNames()) { + std::cout << name << ": "; + std::cout << wgtScale.scaleIndex(1.0, 1.0, name) << " "; + std::cout << wgtScale.scaleIndex(1.0, 2.0, name) << " "; + std::cout << wgtScale.scaleIndex(1.0, 0.5, name) << " "; + std::cout << wgtScale.scaleIndex(2.0, 1.0, name) << " "; + std::cout << wgtScale.scaleIndex(2.0, 2.0, name) << " "; + std::cout << wgtScale.scaleIndex(2.0, 0.5, name) << " "; + std::cout << wgtScale.scaleIndex(0.5, 1.0, name) << " "; + std::cout << wgtScale.scaleIndex(0.5, 2.0, name) << " "; + std::cout << wgtScale.scaleIndex(0.5, 0.5, name) << "\n"; + } + + } else if (wgt.weightType() == gen::WeightType::kPdfWeights) { + std::cout << wgt.description() << "\n"; + } else if (wgt.weightType() == gen::WeightType::kPartonShowerWeights) { + auto& wgtPS = dynamic_cast(wgt); + std::vector labels = wgtPS.weightLabels(); + wgtPS.cacheWeightIndicesByLabel(); + wgtPS.printVariables(); + } + } + } + + std::unique_ptr WeightHelper::buildGroup(ParsedWeight& weight) { + if (debug_) { + std::cout << "Building group for weight group " << weight.groupname << " weight content is " << weight.content + << std::endl; + } + if (isScaleWeightGroup(weight)) + return std::make_unique(weight.groupname); + else if (isPdfWeightGroup(weight)) + return std::make_unique(weight.groupname); + else if (isMEParamWeightGroup(weight)) + return std::make_unique(weight.groupname); + else if (isPartonShowerWeightGroup(weight)) + return std::make_unique("shower"); + else if (isOrphanPdfWeightGroup(weight)) + return std::make_unique(weight.groupname); + + return std::make_unique(weight.groupname); + } + + void WeightHelper::buildGroups() { + weightGroups_.clear(); + int groupOffset = 0; + for (auto& weight : parsedWeights_) { + weight.wgtGroup_idx += groupOffset; + if (debug_) + std::cout << "Building group for weight " << weight.content << " group " << weight.groupname << " group index " + << weight.wgtGroup_idx << std::endl; + + int numGroups = static_cast(weightGroups_.size()); + if (weight.wgtGroup_idx == numGroups) { + weightGroups_.push_back(*buildGroup(weight)); + } else if (weight.wgtGroup_idx >= numGroups) + cms::Exception("Invalid group index " + std::to_string(weight.wgtGroup_idx)); + + // split PDF groups + if (splitPdfWeight(weight)) + groupOffset++; + + WeightGroupInfo& group = weightGroups_[weight.wgtGroup_idx]; + group.addContainedId(weight.index, weight.id, weight.content); + if (group.weightType() == gen::WeightType::kScaleWeights) + updateScaleInfo(dynamic_cast(group), weight); + else if (group.weightType() == gen::WeightType::kPdfWeights) + updatePdfInfo(dynamic_cast(group), weight); + else if (group.weightType() == gen::WeightType::kPartonShowerWeights) + updatePartonShowerInfo(dynamic_cast(group), weight); + } + cleanupOrphanCentralWeight(); + } + +} // namespace gen diff --git a/GeneratorInterface/Core/test/dumpWeightInfo.py b/GeneratorInterface/Core/test/dumpWeightInfo.py new file mode 100644 index 0000000000000..b7231c6945a33 --- /dev/null +++ b/GeneratorInterface/Core/test/dumpWeightInfo.py @@ -0,0 +1,44 @@ +from __future__ import print_function +from DataFormats.FWLite import Events,Handle,Runs,Lumis +import ROOT +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("infile", type=str, help="Input EDM file") +parser.add_argument("--source", type=str, help="product ID of weight product", default="externalLHEProducer") +args = parser.parse_args() + +lumis = Lumis(args.infile) +lumi = lumis.__iter__().next() +weightInfoHandle = Handle("GenWeightInfoProduct") +print("Trying to get weightInfo from lumi") +lumi.getByLabel(args.source, weightInfoHandle) +weightInfoProd = weightInfoHandle.product() + +events = Events(args.infile) +event = events.__iter__().next() +weightHandle = Handle("GenWeightProduct") +print("Trying to get weightProduct from event") +event.getByLabel(args.source, weightHandle) +weightInfo = weightHandle.product() +print("Number of weight groups in weightInfo is", len(weightInfo.weights())) +for j, weights in enumerate(weightInfo.weights()): + print("-"*10, "Looking at entry", j, "length is", len(weights),"-"*10) + matching = weightInfoProd.orderedWeightGroupInfo(j) + print("Group is", matching, "name is", matching.name(), "well formed?", matching.isWellFormed()) + print("Group description", matching.description()) + print(type(matching.weightType()), matching.weightType()) + if matching.weightType() == 's': + for var in [(x, y) for x in ["05", "1", "2"] for y in ["05", "1", "2"]]: + name = "muR%smuF%sIndex" % (var[0], var[1]) if not (var[0] == "1" and var[1] == "1") else "centralIndex" + print(name, getattr(matching, name)()) + elif matching.weightType() == 'P': + print("uncertaintyType", "Hessian" if matching.uncertaintyType() == ROOT.gen.kHessianUnc else "MC") + print("Has alphas? ", matching.hasAlphasVariations()) + print("Weights length?", len(weights), "Contained ids lenths?", len(matching.containedIds())) + print("-"*80) + for i,weight in enumerate(weights): + print(i, weight) + info = matching.weightMetaInfo(i) + print(" ID, localIndex, globalIndex, label, set:", info.id, info.localIndex, info.globalIndex, info.label, matching.name()) + print("-"*80) diff --git a/GeneratorInterface/Core/test/testGenWeightProducer_cfg.py b/GeneratorInterface/Core/test/testGenWeightProducer_cfg.py new file mode 100644 index 0000000000000..cde648a77d08c --- /dev/null +++ b/GeneratorInterface/Core/test/testGenWeightProducer_cfg.py @@ -0,0 +1,42 @@ +import FWCore.ParameterSet.Config as cms +from FWCore.ParameterSet.VarParsing import VarParsing + +options = VarParsing('analysis') +options.register( + "lheSource", + "externalLHEProducer", + VarParsing.multiplicity.singleton, + VarParsing.varType.string, + "LHE source (default externalLHEProducer)" +) +options.parseArguments() + +process = cms.Process("test") + +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring(options.inputFiles) +) + +process.maxEvents = cms.untracked.PSet(input=cms.untracked.int32(options.maxEvents)) + +process.out = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string('test.root'), + outputCommands = cms.untracked.vstring(['drop *', + 'keep GenWeightProduct_test*Weights*_*_*', + 'keep GenWeightInfoProduct_test*Weights*_*_*',]) +) + +# process.testLHEWeights = cms.EDProducer("LHEWeightProductProducer", +# lheSourceLabel = cms.string("externalLHEProducer"), +# failIfValidXML = cms.untracked.bool(True)) + +process.testGenWeights = cms.EDProducer("GenWeightProductProducer", + genInfo = cms.InputTag("generator"), + genLumiInfoHeader = cms.InputTag("generator")) + +#process.p = cms.Path(process.testLHEWeights*process.testGenWeights) +process.p = cms.Path(process.testGenWeights) + +process.output = cms.EndPath(process.out) +process.schedule = cms.Schedule(process.p,process.output) + diff --git a/GeneratorInterface/LHEInterface/plugins/BuildFile.xml b/GeneratorInterface/LHEInterface/plugins/BuildFile.xml index b96e39874c6f2..f43689ef148df 100644 --- a/GeneratorInterface/LHEInterface/plugins/BuildFile.xml +++ b/GeneratorInterface/LHEInterface/plugins/BuildFile.xml @@ -8,6 +8,7 @@ + @@ -15,7 +16,6 @@ - diff --git a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc index 63269fb35b56a..d50e95666d562 100644 --- a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc +++ b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc @@ -52,11 +52,17 @@ Description: [one line class summary] #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEXMLStringProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h" #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h" #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h" #include "GeneratorInterface/LHEInterface/interface/LHEReader.h" +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h" +#include "GeneratorInterface/Core/interface/LHEWeightHelper.h" + #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/Utilities/interface/RandomNumberGenerator.h" @@ -66,7 +72,8 @@ Description: [one line class summary] // class declaration // -class ExternalLHEProducer : public edm::one::EDProducer { +class ExternalLHEProducer + : public edm::one::EDProducer { public: explicit ExternalLHEProducer(const edm::ParameterSet& iConfig); @@ -76,6 +83,7 @@ class ExternalLHEProducer : public edm::one::EDProducer makeArgs(uint32_t nEvents, unsigned int nThreads, std::uint32_t seed) const; @@ -99,6 +107,7 @@ class ExternalLHEProducer : public edm::one::EDProducer> nPartonMapping_{}; std::unique_ptr reader_; + gen::LHEWeightHelper weightHelper_; std::shared_ptr runInfoLast_; std::shared_ptr runInfo_; std::shared_ptr partonLevel_; @@ -158,9 +167,11 @@ ExternalLHEProducer::ExternalLHEProducer(const edm::ParameterSet& iConfig) produces("LHEScriptOutput"); + produces(); produces(); produces(); produces(); + produces(); } // @@ -188,6 +199,10 @@ void ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe std::for_each(partonLevel_->weights().begin(), partonLevel_->weights().end(), std::bind(&LHEEventProduct::addWeight, product.get(), std::placeholders::_1)); + + auto weightProduct = weightHelper_.weightProduct(partonLevel_->weights(), partonLevel_->originalXWGTUP()); + iEvent.put(std::move(weightProduct)); + product->setScales(partonLevel_->scales()); if (nPartonMapping_.empty()) { product->setNpLO(partonLevel_->npLO()); @@ -422,6 +437,14 @@ std::vector ExternalLHEProducer::makeArgs(uint32_t nEvents, return args; } +void ExternalLHEProducer::beginLuminosityBlockProduce(edm::LuminosityBlock& lumi, edm::EventSetup const& es) { + auto weightInfoProduct = std::make_unique(); + for (auto& weightGroup : weightHelper_.weightGroups()) { + weightInfoProduct->addWeightGroupInfo(std::unique_ptr(weightGroup.clone())); + } + lumi.put(std::move(weightInfoProduct)); +} + // ------------ Close all the open file descriptors ------------ int ExternalLHEProducer::closeDescriptors(int preserve) const { int maxfd = 1024; diff --git a/GeneratorInterface/LHEInterface/test/createLHEFormatFromROOTFile.py b/GeneratorInterface/LHEInterface/test/createLHEFormatFromROOTFile.py new file mode 100644 index 0000000000000..33c60dbff7ca4 --- /dev/null +++ b/GeneratorInterface/LHEInterface/test/createLHEFormatFromROOTFile.py @@ -0,0 +1,107 @@ +import ROOT +import sys + + +""" +Usage +python createLHEFormatFromROOTFile.py inputfile outputfile pdgId_particle_to_undo_decay1 pdgId_particle_to_undo_decay2 pdgId_particle_to_undo_decay3 ... +""" + +args = sys.argv[:] + +class HEPPart(object): + def __init__(self,event, idx): + """ + Class to organize the description of a particle in the LHE event + event : whole event information (usually a entry in the input TTree) + idx : the index of the particle inside the LHE file + """ + for att in ["pt","eta","phi", "mass","lifetime","pdgId","status","spin","color1", "color2","mother1","mother2","incomingpz"]: + setattr(self, att, getattr(event, "LHEPart_"+att)[idx]) + self.setP4() + + def setP4(self): + self.p4 = ROOT.TLorentzVector() + if self.status != -1: + self.p4.SetPtEtaPhiM(self.pt, self.eta, self.phi, self.mass) + else: + self.p4.SetPxPyPzE(0.,0.,self.incomingpz, abs(self.incomingpz)) + def printPart(self): + """ Just to print it pretty """ + return " {pdg:d} {status:d} {mother1:d} {mother2:d} {color1:d} {color2:d} {px:e} {py:e} {pz:e} {energy:e} {mass:e} {time:e} {spin:e}\n".format(pdg=self.pdgId,status=self.status, mother1=self.mother1, mother2=self.mother2, color1=self.color1, color2=self.color2, px=self.p4.Px(), py=self.p4.Py(), pz=self.p4.Pz(), energy=self.p4.E(), mass=self.mass, time=self.lifetime, spin=self.spin) + + +class LHEPrinter(object): + def __init__(self,theFile,theTree,outputLHE,undoDecays=[], chunkers=-1, prDict={}): + """ + theFile : path to input root file with the whole LHEinformation + theTree : number of ttree inside the file, usually "Events" + outputLHE: name of output .lhe file + undoDecays: pdgId of particles whose decays we want to undo (i.e. W that are decayed with madspin, as the reweighting is called before madspin) + chunkers: process by chunks in case we want to later multithread the jobs. Argument is max size (in events) of the chunks + prDict : a dictionary indicating possible mismatchings between the ProcessID in the new gridpack and the old (matching the numbers at generate p p > x y z @0 in the run card) + """ + + self.prDict = prDict + self.outputLHE = outputLHE + self.fil = ROOT.TFile(theFile,"OPEN") + self.tree = self.fil.Get(theTree) + self.undoDecays = undoDecays + self.baseheader = "\n" + self.baseline1 = " {nparts} {prid} {weight} {scale} {aqed} {aqcd}\n" + self.baseender = "\n\n\n" + self.chunkers = chunkers + + def insideLoop(self): + """ Loop over all events and process the root file into a plain text LHE file""" + totalEvents = self.tree.GetEntries() + if self.chunkers == -1 or self.chunkers > totalEvents: + self.output = open(self.outputLHE,"w") + self.chunkers = totalEvents + 1 + else: + self.output = open(self.outputLHE+"chunk0","w") + chunk = 0 + + print "Processing %i events, please wait..."%totalEvents + iEv = 0 + pEv = 0 + chunk = 0 + for ev in self.tree: + iEv += 1 + pEv += 1 + if pEv >= self.chunkers: + pEv = 0 + chunk += 1 + self.output.close() + self.output = open(self.outputLHE+"chunk%i"%chunk,"w") + print "...Event %i/%i"%(iEv, totalEvents) + self.process(ev) + + def process(self, ev): + """First produce the global line like """ + self.output.write(self.baseheader.format(nplo = ord(str(ev.LHE_NpLO)) if ord(str(ev.LHE_NpLO)) != 255 else -1, npnlo = ord(str(ev.LHE_NpNLO)) if ord(str(ev.LHE_NpNLO)) != 255 else -1)) + + """Then we need to treat the whole thing to undo the madspin decays, update statuses and rewrite particle order""" + lhepart = [] + deletedIndexes = [] + for i in range(getattr(ev, "nLHEPart")): + testPart = HEPPart(ev,i) + testPart.mother1 = testPart.mother1 - sum([1*(testPart.mother1 > d) for d in deletedIndexes]) + testPart.mother2 = testPart.mother2 - sum([1*(testPart.mother2 > d) for d in deletedIndexes]) + if testPart.mother1 != 0: + if abs(lhepart[testPart.mother1-1].pdgId) in self.undoDecays: #If from something that decays after weighting just skip it and update particle indexes + deletedIndexes.append(i) + continue + if abs(testPart.pdgId) in self.undoDecays: + testPart.status = 1 + lhepart.append(testPart) + + """ Now we can compute properly the number of particles at LHE """ + self.output.write(self.baseline1.format(nparts=len(lhepart), prid=self.prDict[str(ord(str(ev.LHE_ProcessID)))], weight=ev.LHEWeight_originalXWGTUP, scale=ev.LHE_Scale,aqed=ev.LHE_AlphaQED,aqcd=ev.LHE_AlphaS)) + """ And save each particle information """ + for part in lhepart: + self.output.write(part.printPart()) + self.output.write(self.baseender) + +theP = LHEPrinter(args[1],"Events",args[2],undoDecays=[int(i) for i in args[4:]],chunkers=int(args[3]), prDict={str(i):i for i in range(1000)}) #PRDict by default set to not change anything as it is rare to use it +theP.insideLoop() diff --git a/PhysicsTools/NanoAOD/interface/GenWeightCounters.h b/PhysicsTools/NanoAOD/interface/GenWeightCounters.h new file mode 100644 index 0000000000000..35446d7172e78 --- /dev/null +++ b/PhysicsTools/NanoAOD/interface/GenWeightCounters.h @@ -0,0 +1,112 @@ +#ifndef GENWEIGHTSCOUNTERS_h +#define GENWEIGHTSCOUNTERS_h +#include +#include +#include +namespace genCounter { + + void mergeSumVectors(std::vector& v1, std::vector 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; + weightSumMap_.clear(); + } + + // inc only the gen counters + void incGenOnly(double w) { + num_++; + sumw_ += w; + sumw2_ += (w * w); + } + + void incLHE(double w0, const std::vector& weightV, const std::string& wName) { + //if new type of weight, create a map element + if (weightSumMap_.find(wName) == weightSumMap_.end()) { + std::vector temp; + weightSumMap_.insert({wName, temp}); + } + if (!weightV.empty()) { + if (weightSumMap_[wName].empty()) + weightSumMap_[wName].resize(weightV.size(), 0); + for (unsigned int i = 0, n = weightV.size(); i < n; ++i) + weightSumMap_[wName][i] += (w0 * weightV[i]); + } + // incGenOnly(w0); + //incPSOnly(w0, wPS); + } + + void mergeSumMap(const Counter& other) { + num_ += other.num_; + sumw_ += other.sumw_; + sumw2_ += other.sumw2_; + //if weightMap_ for "this" is empty, create map elements with empty + //vectors before merging + if (weightSumMap_.empty() && !other.weightSumMap_.empty()) { + for (auto& wmap : other.weightSumMap_) { + std::vector temp; + weightSumMap_.insert({wmap.first, temp}); + } + } + + for (auto& wmap : weightSumMap_) { + if (other.weightSumMap_.find(wmap.first) != other.weightSumMap_.end()) + mergeSumVectors(wmap.second, other.weightSumMap_.at(wmap.first)); + } + } + + //private: + // the counters + long long num_ = 0; + long double sumw_ = 0; + long double sumw2_ = 0; + std::map> weightSumMap_; + }; + + struct CounterMap { + std::map countermap; + Counter* active_el = nullptr; + std::string active_label = ""; + + void mergeSumMap(const CounterMap& other) { + for (const auto& y : other.countermap) { + countermap[y.first].mergeSumMap(y.second); + } + active_el = nullptr; + } + + void clear() { + for (auto x : countermap) + x.second.clear(); + } + + void setLabel(std::string label) { + active_el = &(countermap[label]); + active_label = label; + } + void checkLabelSet() { + if (!active_el) + throw cms::Exception("LogicError", "Called CounterMap::get() before setting the active label\n"); + } + Counter* get() { + checkLabelSet(); + return active_el; + } + std::string& getLabel() { + checkLabelSet(); + return active_label; + } + }; + +} // namespace genCounter +#endif diff --git a/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc b/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc index b4bd5cf399317..e069c4fb12c6b 100644 --- a/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc @@ -6,1066 +6,112 @@ #include "DataFormats/NanoAOD/interface/MergeableCounterTable.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" -#include "FWCore/Utilities/interface/transform.h" -#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" -#include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "boost/algorithm/string.hpp" - -#include - -#include -#include +#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h" +#include "FWCore/Utilities/interface/transform.h" +#include "PhysicsTools/NanoAOD/interface/GenWeightCounters.h" +#include #include -#include +#include 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 sumPDF, sumScale, sumRwgt, sumNamed, sumPS; - - void 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); - } - - void incPSOnly(double w0, const std::vector& 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& wScale, - const std::vector& wPDF, - const std::vector& wRwgt, - const std::vector& wNamed, - const std::vector& wPS) { - // add up weights - incGenOnly(w0); - // then add up variations - if (!wScale.empty()) { - if (sumScale.empty()) - sumScale.resize(wScale.size(), 0); - for (unsigned int i = 0, n = wScale.size(); i < n; ++i) - sumScale[i] += (w0 * wScale[i]); - } - if (!wPDF.empty()) { - if (sumPDF.empty()) - sumPDF.resize(wPDF.size(), 0); - for (unsigned int i = 0, n = wPDF.size(); i < n; ++i) - sumPDF[i] += (w0 * wPDF[i]); - } - if (!wRwgt.empty()) { - if (sumRwgt.empty()) - sumRwgt.resize(wRwgt.size(), 0); - for (unsigned int i = 0, n = wRwgt.size(); i < n; ++i) - sumRwgt[i] += (w0 * wRwgt[i]); - } - if (!wNamed.empty()) { - 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 (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]; - } - }; - - struct CounterMap { - std::map countermap; - Counter* active_el = nullptr; - std::string active_label = ""; - void merge(const CounterMap& other) { - for (const auto& y : other.countermap) - countermap[y.first].merge(y.second); - active_el = nullptr; - } - void clear() { - for (auto x : countermap) - x.second.clear(); - active_el = nullptr; - active_label = ""; - } - void setLabel(std::string label) { - active_el = &(countermap[label]); - active_label = label; - } - void checkLabelSet() { - if (!active_el) - throw cms::Exception("LogicError", "Called CounterMap::get() before setting the active label\n"); - } - Counter* get() { - checkLabelSet(); - return active_el; - } - std::string& getLabel() { - checkLabelSet(); - return active_label; - } - }; - - /// ---- RunCache object for dynamic choice of LHE IDs ---- - struct DynamicWeightChoice { - // choice of LHE weights - // ---- scale ---- - std::vector scaleWeightIDs; - std::string scaleWeightsDoc; - // ---- pdf ---- - std::vector pdfWeightIDs; - std::string pdfWeightsDoc; - // ---- rwgt ---- - std::vector rwgtIDs; - std::string rwgtWeightDoc; - }; - - struct DynamicWeightChoiceGenInfo { - // choice of LHE weights - // ---- scale ---- - std::vector scaleWeightIDs; - std::string scaleWeightsDoc; - // ---- pdf ---- - std::vector pdfWeightIDs; - std::string pdfWeightsDoc; - // ---- ps ---- - std::vector defPSWeightIDs = {6, 7, 8, 9}; - std::vector defPSWeightIDs_alt = {27, 5, 26, 4}; - bool matchPS_alt = false; - std::vector psWeightIDs; - unsigned int psBaselineID = 1; - std::string psWeightsDoc; - - void setMissingWeight(int idx) { psWeightIDs[idx] = (matchPS_alt) ? defPSWeightIDs_alt[idx] : defPSWeightIDs[idx]; } - - bool empty() const { return scaleWeightIDs.empty() && pdfWeightIDs.empty() && psWeightIDs.empty(); } - }; - - struct LumiCacheInfoHolder { - CounterMap countermap; - DynamicWeightChoiceGenInfo weightChoice; - void clear() { - countermap.clear(); - weightChoice = DynamicWeightChoiceGenInfo(); - } - }; - - float stof_fortrancomp(const std::string& str) { - std::string::size_type match = str.find('d'); - if (match != std::string::npos) { - std::string pre = str.substr(0, match); - std::string post = str.substr(match + 1); - return std::stof(pre) * std::pow(10.0f, std::stof(post)); - } else { - return std::stof(str); - } - } - /// -------------- temporary objects -------------- - struct ScaleVarWeight { - std::string wid, label; - std::pair scales; - ScaleVarWeight(const std::string& id, const std::string& text, const std::string& muR, const std::string& muF) - : wid(id), label(text), scales(stof_fortrancomp(muR), stof_fortrancomp(muF)) {} - bool operator<(const ScaleVarWeight& other) { - return (scales == other.scales ? wid < other.wid : scales < other.scales); - } - }; - struct PDFSetWeights { - std::vector wids; - std::pair lhaIDs; - PDFSetWeights(const std::string& wid, unsigned int lhaID) : wids(1, wid), lhaIDs(lhaID, lhaID) {} - bool operator<(const PDFSetWeights& other) const { return lhaIDs < other.lhaIDs; } - void add(const std::string& wid, unsigned int lhaID) { - wids.push_back(wid); - lhaIDs.second = lhaID; - } - bool maybe_add(const std::string& wid, unsigned int lhaID) { - if (lhaID == lhaIDs.second + 1) { - lhaIDs.second++; - wids.push_back(wid); - return true; - } else { - return false; - } - } - }; + typedef std::vector WeightGroupDataContainer; + typedef std::array, 2> WeightGroupsToStore; } // namespace +using CounterMap = genCounter::CounterMap; +using Counter = genCounter::Counter; -class GenWeightsTableProducer : public edm::global::EDProducer, - edm::RunCache, +class LHEWeightsTableProducer : public edm::global::EDProducer, + edm::StreamCache, edm::RunSummaryCache, edm::EndRunProducer> { public: - GenWeightsTableProducer(edm::ParameterSet const& params) - : genTag_(consumes(params.getParameter("genEvent"))), - lheLabel_(params.getParameter>("lheInfo")), - lheTag_(edm::vector_transform(lheLabel_, - [this](const edm::InputTag& tag) { return mayConsume(tag); })), - lheRunTag_(edm::vector_transform( - lheLabel_, [this](const edm::InputTag& tag) { return mayConsume(tag); })), - genLumiInfoHeadTag_( - mayConsume(params.getParameter("genLumiInfoHeader"))), - namedWeightIDs_(params.getParameter>("namedWeightIDs")), - namedWeightLabels_(params.getParameter>("namedWeightLabels")), - lheWeightPrecision_(params.getParameter("lheWeightPrecision")), - maxPdfWeights_(params.getParameter("maxPdfWeights")), - keepAllPSWeights_(params.getParameter("keepAllPSWeights")), - debug_(params.getUntrackedParameter("debug", false)), - debugRun_(debug_.load()), - hasIssuedWarning_(false), - psWeightWarning_(false) { - produces(); - produces("genModel"); - produces("LHEScale"); - produces("LHEPdf"); - produces("LHEReweighting"); - produces("LHENamed"); - produces("PS"); - produces(); - if (namedWeightIDs_.size() != namedWeightLabels_.size()) { - throw cms::Exception("Configuration", "Size mismatch between namedWeightIDs & namedWeightLabels"); - } - for (const edm::ParameterSet& pdfps : params.getParameter>("preferredPDFs")) { - const std::string& name = pdfps.getParameter("name"); - uint32_t lhaid = pdfps.getParameter("lhaid"); - preferredPDFLHAIDs_.push_back(lhaid); - lhaNameToID_[name] = lhaid; - lhaNameToID_[name + ".LHgrid"] = lhaid; - } - } - - ~GenWeightsTableProducer() override {} - - void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override { - // get my counter for weights - Counter* counter = streamCache(id)->countermap.get(); - - // generator information (always available) - edm::Handle genInfo; - iEvent.getByToken(genTag_, genInfo); - double weight = genInfo->weight(); - - // table for gen info, always available - auto out = std::make_unique(1, "genWeight", true); - out->setDoc("generator weight"); - out->addColumnValue("", weight, "generator weight"); - iEvent.put(std::move(out)); - - std::string model_label = streamCache(id)->countermap.getLabel(); - auto outM = std::make_unique((!model_label.empty()) ? std::string("GenModel_") + model_label : ""); - iEvent.put(std::move(outM), "genModel"); - bool getLHEweightsFromGenInfo = !model_label.empty(); - - // tables for LHE weights, may not be filled - std::unique_ptr lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab; - std::unique_ptr genPSTab; - - edm::Handle lheInfo; - for (const auto& lheTag : lheTag_) { - iEvent.getByToken(lheTag, lheInfo); - if (lheInfo.isValid()) { + LHEWeightsTableProducer(edm::ParameterSet const& params); + + void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; + //func changed//sroychow + void addWeightGroupToTable(std::map>& lheWeightTables, + std::map>& weightVecsizes, + std::map& weightlabels, + const char* typeName, + const WeightGroupDataContainer& weightInfos, + WeightsContainer& allWeights, + Counter& counter, + double genWeight) const; + WeightGroupDataContainer weightDataPerType(edm::Handle& weightsHandle, + gen::WeightType weightType, + int& maxStore) const; + + std::vector orderedScaleWeights(const std::vector& scaleWeights, + const gen::ScaleWeightGroupInfo& scaleGroup) const; + + std::vector preferredPSweights(const std::vector& psWeights, + const gen::PartonShowerWeightGroupInfo& pswV) const; + + //Lumiblock + std::shared_ptr globalBeginLuminosityBlock(edm::LuminosityBlock const& iLumi, + edm::EventSetup const&) const override { + // Set equal to the max number of groups + // subtrack 1 for each weight group you find + bool foundLheWeights = false; + edm::Handle lheWeightInfoHandle; + for (auto& token : lheWeightInfoTokens_) { + iLumi.getByToken(token, lheWeightInfoHandle); + if (lheWeightInfoHandle.isValid()) { + foundLheWeights = true; break; } } - - const auto genWeightChoice = &(streamCache(id)->weightChoice); - if (lheInfo.isValid()) { - if (getLHEweightsFromGenInfo && !hasIssuedWarning_.exchange(true)) - edm::LogWarning("LHETablesProducer") - << "Found both a LHEEventProduct and a GenLumiInfoHeader: will only save weights from LHEEventProduct.\n"; - // get the dynamic choice of weights - const DynamicWeightChoice* weightChoice = runCache(iEvent.getRun().index()); - // go fill tables - fillLHEWeightTables(counter, - weightChoice, - genWeightChoice, - weight, - *lheInfo, - *genInfo, - lheScaleTab, - lhePdfTab, - lheRwgtTab, - lheNamedTab, - genPSTab); - } else if (getLHEweightsFromGenInfo) { - fillLHEPdfWeightTablesFromGenInfo( - counter, genWeightChoice, weight, *genInfo, lheScaleTab, lhePdfTab, lheNamedTab, genPSTab); - lheRwgtTab = std::make_unique(1, "LHEReweightingWeights", true); - //lheNamedTab.reset(new nanoaod::FlatTable(1, "LHENamedWeights", true)); - //genPSTab.reset(new nanoaod::FlatTable(1, "PSWeight", true)); - } else { - // Still try to add the PS weights - fillOnlyPSWeightTable(counter, genWeightChoice, weight, *genInfo, genPSTab); - // make dummy values - lheScaleTab = std::make_unique(1, "LHEScaleWeights", true); - lhePdfTab = std::make_unique(1, "LHEPdfWeights", true); - lheRwgtTab = std::make_unique(1, "LHEReweightingWeights", true); - lheNamedTab = std::make_unique(1, "LHENamedWeights", true); - if (!hasIssuedWarning_.exchange(true)) { - edm::LogWarning("LHETablesProducer") << "No LHEEventProduct, so there will be no LHE Tables\n"; + edm::Handle genWeightInfoHandle; + iLumi.getByToken(genWeightInfoToken_, genWeightInfoHandle); + + std::unordered_map storePerType; + for (size_t i = 0; i < weightgroups_.size(); i++) + storePerType[weightgroups_.at(i)] = maxGroupsPerType_.at(i); + + WeightGroupsToStore weightsToStore; + for (auto weightType : gen::allWeightTypes) { + if (foundLheWeights) { + auto lheWeights = weightDataPerType(lheWeightInfoHandle, weightType, storePerType[weightType]); + for (auto& w : lheWeights) + weightsToStore.at(inLHE).push_back({w.index, std::move(w.group)}); } + auto genWeights = weightDataPerType(genWeightInfoHandle, weightType, storePerType[weightType]); + for (auto& w : genWeights) + weightsToStore.at(inGen).push_back({w.index, std::move(w.group)}); } - - 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"); - } - - void fillLHEWeightTables(Counter* counter, - const DynamicWeightChoice* weightChoice, - const DynamicWeightChoiceGenInfo* genWeightChoice, - double genWeight, - const LHEEventProduct& lheProd, - const GenEventInfoProduct& genProd, - std::unique_ptr& outScale, - std::unique_ptr& outPdf, - std::unique_ptr& outRwgt, - std::unique_ptr& outNamed, - std::unique_ptr& 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) - - const std::vector& scaleWeightIDs = weightChoice->scaleWeightIDs; - const std::vector& pdfWeightIDs = weightChoice->pdfWeightIDs; - const std::vector& rwgtWeightIDs = weightChoice->rwgtIDs; - - double w0 = lheProd.originalXWGTUP(); - - std::vector wScale(scaleWeightIDs.size(), 1), wPDF(pdfWeightIDs.size(), 1), wRwgt(rwgtWeightIDs.size(), 1), - 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()); - // now we do it slowly, can be optimized - auto mScale = std::find(scaleWeightIDs.begin(), scaleWeightIDs.end(), weight.id); - if (mScale != scaleWeightIDs.end()) - wScale[mScale - scaleWeightIDs.begin()] = weight.wgt / w0; - - auto mPDF = std::find(pdfWeightIDs.begin(), pdfWeightIDs.end(), weight.id); - if (mPDF != pdfWeightIDs.end()) - wPDF[mPDF - pdfWeightIDs.begin()] = weight.wgt / w0; - - auto mRwgt = std::find(rwgtWeightIDs.begin(), rwgtWeightIDs.end(), weight.id); - if (mRwgt != rwgtWeightIDs.end()) - wRwgt[mRwgt - rwgtWeightIDs.begin()] = weight.wgt / w0; - - auto mNamed = std::find(namedWeightIDs_.begin(), namedWeightIDs_.end(), weight.id); - if (mNamed != namedWeightIDs_.end()) - wNamed[mNamed - namedWeightIDs_.begin()] = weight.wgt / w0; - } - - std::vector wPS; - std::string psWeightDocStr; - setPSWeightInfo(genProd.weights(), genWeightChoice, wPS, psWeightDocStr); - - outPS = std::make_unique(wPS.size(), "PSWeight", false); - outPS->addColumn("", wPS, psWeightDocStr, lheWeightPrecision_); - - outScale = std::make_unique(wScale.size(), "LHEScaleWeight", false); - outScale->addColumn("", wScale, weightChoice->scaleWeightsDoc, lheWeightPrecision_); - - outPdf = std::make_unique(wPDF.size(), "LHEPdfWeight", false); - outPdf->addColumn("", wPDF, weightChoice->pdfWeightsDoc, lheWeightPrecision_); - - outRwgt = std::make_unique(wRwgt.size(), "LHEReweightingWeight", false); - outRwgt->addColumn("", wRwgt, weightChoice->rwgtWeightDoc, lheWeightPrecision_); - - outNamed = std::make_unique(1, "LHEWeight", true); - outNamed->addColumnValue("originalXWGTUP", lheProd.originalXWGTUP(), "Nominal event weight in the LHE file"); - for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) { - outNamed->addColumnValue(namedWeightLabels_[i], - wNamed[i], - "LHE weight for id " + namedWeightIDs_[i] + ", relative to nominal", - lheWeightPrecision_); - } - - counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS); - } - - void fillLHEPdfWeightTablesFromGenInfo(Counter* counter, - const DynamicWeightChoiceGenInfo* weightChoice, - double genWeight, - const GenEventInfoProduct& genProd, - std::unique_ptr& outScale, - std::unique_ptr& outPdf, - std::unique_ptr& outNamed, - std::unique_ptr& outPS) const { - const std::vector& scaleWeightIDs = weightChoice->scaleWeightIDs; - const std::vector& pdfWeightIDs = weightChoice->pdfWeightIDs; - - auto weights = genProd.weights(); - double w0 = (weights.size() > 1) ? weights.at(1) : 1.; - double originalXWGTUP = (weights.size() > 1) ? weights.at(1) : 1.; - - std::vector wScale, wPDF, wPS; - for (auto id : scaleWeightIDs) - wScale.push_back(weights.at(id) / w0); - for (auto id : pdfWeightIDs) { - wPDF.push_back(weights.at(id) / w0); - } - - std::string psWeightsDocStr; - setPSWeightInfo(genProd.weights(), weightChoice, wPS, psWeightsDocStr); - - outScale = std::make_unique(wScale.size(), "LHEScaleWeight", false); - outScale->addColumn("", wScale, weightChoice->scaleWeightsDoc, lheWeightPrecision_); - - outPdf = std::make_unique(wPDF.size(), "LHEPdfWeight", false); - outPdf->addColumn("", wPDF, weightChoice->pdfWeightsDoc, lheWeightPrecision_); - - outPS = std::make_unique(wPS.size(), "PSWeight", false); - outPS->addColumn("", wPS, psWeightsDocStr, lheWeightPrecision_); - - outNamed = std::make_unique(1, "LHEWeight", true); - outNamed->addColumnValue("originalXWGTUP", originalXWGTUP, "Nominal event weight in the LHE file"); - /*for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) { - outNamed->addColumnValue(namedWeightLabels_[i], wNamed[i], "LHE weight for id "+namedWeightIDs_[i]+", relative to nominal", lheWeightPrecision_); - }*/ - - counter->incLHE(genWeight, wScale, wPDF, std::vector(), std::vector(), wPS); - } - - void fillOnlyPSWeightTable(Counter* counter, - const DynamicWeightChoiceGenInfo* genWeightChoice, - double genWeight, - const GenEventInfoProduct& genProd, - std::unique_ptr& outPS) const { - std::vector wPS; - std::string psWeightDocStr; - setPSWeightInfo(genProd.weights(), genWeightChoice, wPS, psWeightDocStr); - outPS = std::make_unique(wPS.size(), "PSWeight", false); - outPS->addColumn("", wPS, psWeightDocStr, lheWeightPrecision_); - - counter->incGenOnly(genWeight); - counter->incPSOnly(genWeight, wPS); - } - - void setPSWeightInfo(const std::vector& genWeights, - const DynamicWeightChoiceGenInfo* genWeightChoice, - std::vector& wPS, - std::string& psWeightDocStr) const { - wPS.clear(); - // isRegularPSSet = keeping all weights and the weights are a usual size, ie - // all weights are PS weights (don't use header incase missing names) - bool isRegularPSSet = keepAllPSWeights_ && (genWeights.size() == 14 || genWeights.size() == 46); - if (!genWeightChoice->psWeightIDs.empty() && !isRegularPSSet) { - psWeightDocStr = genWeightChoice->psWeightsDoc; - double psNom = genWeights.at(genWeightChoice->psBaselineID); - for (auto wgtidx : genWeightChoice->psWeightIDs) { - wPS.push_back(genWeights.at(wgtidx) / psNom); - } - } else { - int vectorSize = - keepAllPSWeights_ ? (genWeights.size() - 2) : ((genWeights.size() == 14 || genWeights.size() == 46) ? 4 : 1); - - if (vectorSize > 1) { - double nominal = genWeights.at(1); // Called 'Baseline' in GenLumiInfoHeader - if (keepAllPSWeights_) { - for (int i = 0; i < vectorSize; i++) { - wPS.push_back(genWeights.at(i + 2) / nominal); - } - psWeightDocStr = "All PS weights (w_var / w_nominal)"; - } else { - if (!psWeightWarning_.exchange(true)) - edm::LogWarning("LHETablesProducer") - << "GenLumiInfoHeader not found: Central PartonShower weights will fill with the 6-10th entries \n" - << " This may incorrect for some mcs (madgraph 2.6.1 with its `isr:murfact=0.5` have a differnt " - "order )"; - for (std::size_t i = 6; i < 10; i++) { - wPS.push_back(genWeights.at(i) / nominal); - } - psWeightDocStr = - "PS weights (w_var / w_nominal); [0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2" - "[2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;"; - } - } else { - wPS.push_back(1.0); - psWeightDocStr = "dummy PS weight (1.0) "; - } - } - } - - // create an empty counter - std::shared_ptr globalBeginRun(edm::Run const& iRun, edm::EventSetup const&) const override { - edm::Handle 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) - auto weightChoice = std::make_shared(); - - // getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499) - //if (iRun.getByToken(lheRunTag_, lheInfo)) { - for (const auto& lheLabel : lheLabel_) { - iRun.getByLabel(lheLabel, lheInfo); - if (lheInfo.isValid()) { - break; - } - } - if (lheInfo.isValid()) { - std::vector scaleVariationIDs; - std::vector pdfSetWeightIDs; - std::vector lheReweighingIDs; - bool isFirstGroup = true; - - std::regex weightgroupmg26x(""); - std::regex weightgroup(""); - std::regex weightgroupRwgt(""); - std::regex endweightgroup(""); - std::regex scalewmg26x( - ""); - std::regex scalewmg26xNew( - ""); - - // MUR=2.0 - std::regex scalew( - "\\s*(?:lhapdf=\\d+|dyn=\\s*-?\\d+)?\\s*((?:mu[rR]|renscfact)" - "=(\\S+)\\s+(?:mu[Ff]|facscfact)=(\\S+)(\\s+.*)?)"); - std::regex pdfw( - "\\s*(?:PDF set|lhapdf|PDF|pdfset)\\s*=\\s*(\\d+)\\s*(?:\\s.*)?"); - std::regex pdfwOld("\\s*Member \\s*(\\d+)\\s*(?:.*)"); - std::regex pdfwmg26x( - "\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?"); - // - - // PDF=325300 MemberID=0 - std::regex pdfwmg26xNew( - "" - "\\s*(?:PDF=(\\d+)\\s*MemberID=(\\d+))?\\s*(?:\\s.*)?"); - - std::regex rwgt("(.+)?()?"); - std::smatch groups; - for (auto iter = lheInfo->headers_begin(), end = lheInfo->headers_end(); iter != end; ++iter) { - if (iter->tag() != "initrwgt") { - if (lheDebug) - std::cout << "Skipping LHE header with tag" << iter->tag() << std::endl; - continue; - } - if (lheDebug) - std::cout << "Found LHE header with tag" << iter->tag() << std::endl; - std::vector lines = iter->lines(); - bool missed_weightgroup = - false; //Needed because in some of the samples ( produced with MG26X ) a small part of the header info is ordered incorrectly - bool ismg26x = false; - bool ismg26xNew = false; - for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; - ++iLine) { //First start looping through the lines to see which weightgroup pattern is matched - boost::replace_all(lines[iLine], "<", "<"); - boost::replace_all(lines[iLine], ">", ">"); - if (std::regex_search(lines[iLine], groups, weightgroupmg26x)) { - ismg26x = true; - } else if (std::regex_search(lines[iLine], groups, scalewmg26xNew) || - std::regex_search(lines[iLine], groups, pdfwmg26xNew)) { - ismg26xNew = true; - } - } - for (unsigned int iLine = 0, nLines = lines.size(); iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << lines[iLine]; - if (std::regex_search(lines[iLine], groups, ismg26x ? weightgroupmg26x : weightgroup)) { - std::string groupname = groups.str(2); - if (ismg26x) - groupname = groups.str(1); - if (lheDebug) - std::cout << ">>> Looks like the beginning of a weight group for '" << groupname << "'" << std::endl; - if (groupname.find("scale_variation") == 0 || groupname == "Central scale variation" || isFirstGroup) { - if (lheDebug && groupname.find("scale_variation") != 0 && groupname != "Central scale variation") - std::cout << ">>> First weight is not scale variation, but assuming is the Central Weight" << std::endl; - else if (lheDebug) - std::cout << ">>> Looks like scale variation for theory uncertainties" << std::endl; - isFirstGroup = false; - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) { - std::cout << " " << lines[iLine]; - } - if (std::regex_search( - lines[iLine], groups, ismg26x ? scalewmg26x : (ismg26xNew ? scalewmg26xNew : scalew))) { - if (lheDebug) - std::cout << " >>> Scale weight " << groups[1].str() << " for " << groups[3].str() << " , " - << groups[4].str() << " , " << groups[5].str() << std::endl; - if (ismg26xNew) { - scaleVariationIDs.emplace_back(groups.str(4), groups.str(1), groups.str(3), groups.str(2)); - } else { - scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4)); - } - } else if (std::regex_search(lines[iLine], endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - if (!missed_weightgroup) { - break; - } else - missed_weightgroup = false; - } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end " - "of the group." - << std::endl; - if (ismg26x || ismg26xNew) - missed_weightgroup = true; - --iLine; // rewind by one, and go back to the outer loop - break; - } - } - } else if (groupname == "PDF_variation" || groupname.find("PDF_variation ") == 0) { - if (lheDebug) - std::cout << ">>> Looks like a new-style block of PDF weights for one or more pdfs" << std::endl; - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << " " << lines[iLine]; - if (std::regex_search(lines[iLine], groups, pdfw)) { - unsigned int lhaID = std::stoi(groups.str(2)); - if (lheDebug) - std::cout << " >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID - << std::endl; - if (pdfSetWeightIDs.empty() || !pdfSetWeightIDs.back().maybe_add(groups.str(1), lhaID)) { - pdfSetWeightIDs.emplace_back(groups.str(1), lhaID); - } - } else if (std::regex_search(lines[iLine], endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - if (!missed_weightgroup) { - break; - } else - missed_weightgroup = false; - } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end " - "of the group." - << std::endl; - if (ismg26x || ismg26xNew) - missed_weightgroup = true; - --iLine; // rewind by one, and go back to the outer loop - break; - } - } - } else if (groupname == "PDF_variation1" || groupname == "PDF_variation2") { - if (lheDebug) - std::cout << ">>> Looks like a new-style block of PDF weights for multiple pdfs" << std::endl; - unsigned int lastid = 0; - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << " " << lines[iLine]; - if (std::regex_search(lines[iLine], groups, pdfw)) { - unsigned int id = std::stoi(groups.str(1)); - unsigned int lhaID = std::stoi(groups.str(2)); - if (lheDebug) - std::cout << " >>> PDF weight " << groups.str(1) << " for " << groups.str(2) << " = " << lhaID - << std::endl; - if (id != (lastid + 1) || pdfSetWeightIDs.empty()) { - pdfSetWeightIDs.emplace_back(groups.str(1), lhaID); - } else { - pdfSetWeightIDs.back().add(groups.str(1), lhaID); - } - lastid = id; - } else if (std::regex_search(lines[iLine], endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - if (!missed_weightgroup) { - break; - } else - missed_weightgroup = false; - } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end " - "of the group." - << std::endl; - if (ismg26x || ismg26xNew) - missed_weightgroup = true; - --iLine; // rewind by one, and go back to the outer loop - break; - } - } - } else if (lhaNameToID_.find(groupname) != lhaNameToID_.end()) { - if (lheDebug) - std::cout << ">>> Looks like an old-style PDF weight for an individual pdf" << std::endl; - unsigned int firstLhaID = lhaNameToID_.find(groupname)->second; - bool first = true; - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << " " << lines[iLine]; - if (std::regex_search( - lines[iLine], groups, ismg26x ? pdfwmg26x : (ismg26xNew ? pdfwmg26xNew : pdfwOld))) { - unsigned int member = 0; - if (!ismg26x && !ismg26xNew) { - member = std::stoi(groups.str(2)); - } else if (ismg26xNew) { - if (!groups.str(3).empty()) { - member = std::stoi(groups.str(3)); - } - } else { - if (!groups.str(4).empty()) { - member = std::stoi(groups.str(4)); - } - } - unsigned int lhaID = member + firstLhaID; - if (lheDebug) - std::cout << " >>> PDF weight " << groups.str(1) << " for " << member << " = " << lhaID - << std::endl; - //if (member == 0) continue; // let's keep also the central value for now - if (first) { - pdfSetWeightIDs.emplace_back(groups.str(1), lhaID); - first = false; - } else { - pdfSetWeightIDs.back().add(groups.str(1), lhaID); - } - } else if (std::regex_search(lines[iLine], endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - if (!missed_weightgroup) { - break; - } else - missed_weightgroup = false; - } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end " - "of the group." - << std::endl; - if (ismg26x || ismg26xNew) - missed_weightgroup = true; - --iLine; // rewind by one, and go back to the outer loop - break; - } - } - } else if (groupname == "mass_variation" || groupname == "sthw2_variation" || - groupname == "width_variation") { - if (lheDebug) - std::cout << ">>> Looks like an EW parameter weight" << std::endl; - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << " " << lines[iLine]; - if (std::regex_search(lines[iLine], groups, rwgt)) { - std::string rwgtID = groups.str(1); - if (lheDebug) - std::cout << " >>> LHE reweighting weight: " << rwgtID << std::endl; - if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) { - // we're only interested in the beggining of the block - lheReweighingIDs.emplace_back(rwgtID); - } - } else if (std::regex_search(lines[iLine], endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - } - } - } else { - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << " " << lines[iLine]; - if (std::regex_search(lines[iLine], groups, endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - if (!missed_weightgroup) { - break; - } else - missed_weightgroup = false; - } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end " - "of the group." - << std::endl; - if (ismg26x || ismg26xNew) - missed_weightgroup = true; - --iLine; // rewind by one, and go back to the outer loop - break; - } - } - } - } else if (std::regex_search(lines[iLine], groups, weightgroupRwgt)) { - std::string groupname = groups.str(1); - if (groupname == "mg_reweighting") { - if (lheDebug) - std::cout << ">>> Looks like a LHE weights for reweighting" << std::endl; - for (++iLine; iLine < nLines; ++iLine) { - if (lheDebug) - std::cout << " " << lines[iLine]; - if (std::regex_search(lines[iLine], groups, rwgt)) { - std::string rwgtID = groups.str(1); - if (lheDebug) - std::cout << " >>> LHE reweighting weight: " << rwgtID << std::endl; - if (std::find(lheReweighingIDs.begin(), lheReweighingIDs.end(), rwgtID) == lheReweighingIDs.end()) { - // we're only interested in the beggining of the block - lheReweighingIDs.emplace_back(rwgtID); - } - } else if (std::regex_search(lines[iLine], endweightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the end of a weight group" << std::endl; - if (!missed_weightgroup) { - break; - } else - missed_weightgroup = false; - } else if (std::regex_search(lines[iLine], ismg26x ? weightgroupmg26x : weightgroup)) { - if (lheDebug) - std::cout << ">>> Looks like the beginning of a new weight group, I will assume I missed the end " - "of the group." - << std::endl; - if (ismg26x) - missed_weightgroup = true; - --iLine; // rewind by one, and go back to the outer loop - break; - } - } - } - } - } - //std::cout << "============= END [ " << iter->tag() << " ] ============ \n\n" << std::endl; - - // ----- SCALE VARIATIONS ----- - std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end()); - if (lheDebug) - std::cout << "Found " << scaleVariationIDs.size() << " scale variations: " << std::endl; - std::stringstream scaleDoc; - scaleDoc << "LHE scale variation weights (w_var / w_nominal); "; - for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) { - const auto& sw = scaleVariationIDs[isw]; - if (isw) - scaleDoc << "; "; - scaleDoc << "[" << isw << "] is " << sw.label; - weightChoice->scaleWeightIDs.push_back(sw.wid); - if (lheDebug) - printf(" id %s: scales ren = % .2f fact = % .2f text = %s\n", - sw.wid.c_str(), - sw.scales.first, - sw.scales.second, - sw.label.c_str()); - } - if (!scaleVariationIDs.empty()) - weightChoice->scaleWeightsDoc = scaleDoc.str(); - - // ------ PDF VARIATIONS (take the preferred one) ----- - if (lheDebug) { - std::cout << "Found " << pdfSetWeightIDs.size() << " PDF set errors: " << std::endl; - for (const auto& pw : pdfSetWeightIDs) - printf("lhaIDs %6d - %6d (%3lu weights: %s, ... )\n", - pw.lhaIDs.first, - pw.lhaIDs.second, - pw.wids.size(), - pw.wids.front().c_str()); - } - - // ------ LHE REWEIGHTING ------- - if (lheDebug) { - std::cout << "Found " << lheReweighingIDs.size() << " reweighting weights" << std::endl; - } - std::copy(lheReweighingIDs.begin(), lheReweighingIDs.end(), std::back_inserter(weightChoice->rwgtIDs)); - - std::stringstream pdfDoc; - pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA IDs "; - bool found = false; - for (const auto& pw : pdfSetWeightIDs) { - for (uint32_t lhaid : preferredPDFLHAIDs_) { - if (pw.lhaIDs.first != lhaid && pw.lhaIDs.first != (lhaid + 1)) - continue; // sometimes the first weight is not saved if that PDF is the nominal one for the sample - if (pw.wids.size() == 1) - continue; // only consider error sets - pdfDoc << pw.lhaIDs.first << " - " << pw.lhaIDs.second; - weightChoice->pdfWeightIDs = pw.wids; - if (maxPdfWeights_ < pw.wids.size()) { - weightChoice->pdfWeightIDs.resize(maxPdfWeights_); // drop some replicas - pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas"; - } - weightChoice->pdfWeightsDoc = pdfDoc.str(); - found = true; - break; - } - if (found) - break; - } - } - } - return weightChoice; + return std::make_shared(weightsToStore); } + // nothing to do here + void globalEndLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override {} // create an empty counter - std::unique_ptr beginStream(edm::StreamID) const override { - return std::make_unique(); - } + std::unique_ptr beginStream(edm::StreamID) const override { return std::make_unique(); } // inizialize to zero at begin run void streamBeginRun(edm::StreamID id, edm::Run const&, edm::EventSetup const&) const override { streamCache(id)->clear(); } + void streamBeginLuminosityBlock(edm::StreamID id, edm::LuminosityBlock const& lumiBlock, edm::EventSetup const& eventSetup) const override { - auto counterMap = &(streamCache(id)->countermap); + auto counterMap = streamCache(id); edm::Handle genLumiInfoHead; lumiBlock.getByToken(genLumiInfoHeadTag_, genLumiInfoHead); if (!genLumiInfoHead.isValid()) edm::LogWarning("LHETablesProducer") << "No GenLumiInfoHeader product found, will not fill generator model string.\n"; - - std::string label; - if (genLumiInfoHead.isValid()) { - label = genLumiInfoHead->configDescription(); - boost::replace_all(label, "-", "_"); - boost::replace_all(label, "/", "_"); - } - counterMap->setLabel(label); - - if (genLumiInfoHead.isValid()) { - auto weightChoice = &(streamCache(id)->weightChoice); - - std::vector scaleVariationIDs; - std::vector pdfSetWeightIDs; - weightChoice->psWeightIDs.clear(); - - std::regex scalew("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+)\\,\\s+mur=(\\S+)\\smuf=(\\S+)"); - std::regex pdfw("LHE,\\s+id\\s+=\\s+(\\d+),\\s+(.+),\\s+Member\\s+(\\d+)\\s+of\\ssets\\s+(\\w+\\b)"); - std::regex mainPSw("sr(Def|:murfac=)(Hi|Lo|_dn|_up|0.5|2.0)"); - std::smatch groups; - auto weightNames = genLumiInfoHead->weightNames(); - std::unordered_map knownPDFSetsFromGenInfo_; - unsigned int weightIter = 0; - for (const auto& line : weightNames) { - if (std::regex_search(line, groups, scalew)) { // scale variation - auto id = groups.str(1); - auto group = groups.str(2); - auto mur = groups.str(3); - auto muf = groups.str(4); - if (group.find("Central scale variation") != std::string::npos) - scaleVariationIDs.emplace_back(groups.str(1), groups.str(2), groups.str(3), groups.str(4)); - } else if (std::regex_search(line, groups, pdfw)) { // PDF variation - auto id = groups.str(1); - auto group = groups.str(2); - auto memberid = groups.str(3); - auto pdfset = groups.str(4); - if (group.find(pdfset) != std::string::npos) { - if (knownPDFSetsFromGenInfo_.find(pdfset) == knownPDFSetsFromGenInfo_.end()) { - knownPDFSetsFromGenInfo_[pdfset] = std::atoi(id.c_str()); - pdfSetWeightIDs.emplace_back(id, std::atoi(id.c_str())); - } else - pdfSetWeightIDs.back().add(id, std::atoi(id.c_str())); - } - } else if (line == "Baseline") { - weightChoice->psBaselineID = weightIter; - } else if (line.find("isr") != std::string::npos || line.find("fsr") != std::string::npos) { - weightChoice->matchPS_alt = line.find("sr:") != std::string::npos; // (f/i)sr: for new weights - if (keepAllPSWeights_) { - weightChoice->psWeightIDs.push_back(weightIter); // PS variations - } else if (std::regex_search(line, groups, mainPSw)) { - if (weightChoice->psWeightIDs.empty()) - weightChoice->psWeightIDs = std::vector(4, -1); - int psIdx = (line.find("fsr") != std::string::npos) ? 1 : 0; - psIdx += (groups.str(2) == "Hi" || groups.str(2) == "_up" || groups.str(2) == "2.0") ? 0 : 2; - weightChoice->psWeightIDs[psIdx] = weightIter; - } - } - weightIter++; - } - if (keepAllPSWeights_) { - weightChoice->psWeightsDoc = "All PS weights (w_var / w_nominal)"; - } else if (weightChoice->psWeightIDs.size() == 4) { - weightChoice->psWeightsDoc = - "PS weights (w_var / w_nominal); [0] is ISR=2 FSR=1; [1] is ISR=1 FSR=2" - "[2] is ISR=0.5 FSR=1; [3] is ISR=1 FSR=0.5;"; - for (int i = 0; i < 4; i++) { - if (static_cast(weightChoice->psWeightIDs[i]) == -1) - weightChoice->setMissingWeight(i); - } - } else { - weightChoice->psWeightsDoc = "dummy PS weight (1.0) "; - } - - weightChoice->scaleWeightIDs.clear(); - weightChoice->pdfWeightIDs.clear(); - - std::sort(scaleVariationIDs.begin(), scaleVariationIDs.end()); - std::stringstream scaleDoc; - scaleDoc << "LHE scale variation weights (w_var / w_nominal); "; - for (unsigned int isw = 0, nsw = scaleVariationIDs.size(); isw < nsw; ++isw) { - const auto& sw = scaleVariationIDs[isw]; - if (isw) - scaleDoc << "; "; - scaleDoc << "[" << isw << "] is " << sw.label; - weightChoice->scaleWeightIDs.push_back(std::atoi(sw.wid.c_str())); - } - if (!scaleVariationIDs.empty()) - weightChoice->scaleWeightsDoc = scaleDoc.str(); - std::stringstream pdfDoc; - pdfDoc << "LHE pdf variation weights (w_var / w_nominal) for LHA names "; - bool found = false; - for (const auto& pw : pdfSetWeightIDs) { - if (pw.wids.size() == 1) - continue; // only consider error sets - for (const auto& wantedpdf : lhaNameToID_) { - auto pdfname = wantedpdf.first; - if (knownPDFSetsFromGenInfo_.find(pdfname) == knownPDFSetsFromGenInfo_.end()) - continue; - uint32_t lhaid = knownPDFSetsFromGenInfo_.at(pdfname); - if (pw.lhaIDs.first != lhaid) - continue; - pdfDoc << pdfname; - for (const auto& x : pw.wids) - weightChoice->pdfWeightIDs.push_back(std::atoi(x.c_str())); - if (maxPdfWeights_ < pw.wids.size()) { - weightChoice->pdfWeightIDs.resize(maxPdfWeights_); // drop some replicas - pdfDoc << ", truncated to the first " << maxPdfWeights_ << " replicas"; - } - weightChoice->pdfWeightsDoc = pdfDoc.str(); - found = true; - break; - } - if (found) - break; - } - } + counterMap->setLabel(genLumiInfoHead.isValid() ? genLumiInfoHead->configDescription() : ""); + std::string label = genLumiInfoHead.isValid() ? counterMap->getLabel() : "NULL"; } // create an empty counter std::shared_ptr globalBeginRunSummary(edm::Run const&, edm::EventSetup const&) const override { @@ -1075,105 +121,314 @@ class GenWeightsTableProducer : public edm::global::EDProducermerge(streamCache(id)->countermap); - } + CounterMap* runCounterMap) const override; // nothing to do per se void globalEndRunSummary(edm::Run const&, edm::EventSetup const&, CounterMap* runCounterMap) const override {} // write the total to the run - void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const&, CounterMap const* runCounterMap) const override { - auto out = std::make_unique(); + void globalEndRunProduce(edm::Run& iRun, edm::EventSetup const& es, CounterMap const* runCounterMap) const override; + // nothing to do here + //void globalEndRun(edm::Run const& iRun, edm::EventSetup const& es, CounterMap* runCounterMap) const override {} - for (const auto& x : runCounterMap->countermap) { - auto runCounter = &(x.second); - std::string label = (!x.first.empty()) ? (std::string("_") + x.first) : ""; - std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : ""; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - 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); +protected: + const std::vector> lheWeightTokens_; + const std::vector> lheWeightInfoTokens_; + const edm::EDGetTokenT genWeightToken_; + const edm::EDGetTokenT genWeightInfoToken_; + const edm::EDGetTokenT genEventInfoToken_; + const edm::EDGetTokenT genLumiInfoHeadTag_; + const std::vector weightgroups_; + const std::vector maxGroupsPerType_; + const std::vector pdfIds_; + const std::unordered_map weightTypeNames_ = { + {gen::WeightType::kScaleWeights, "LHEScaleWeight"}, + {gen::WeightType::kPdfWeights, "LHEPdfWeight"}, + {gen::WeightType::kMEParamWeights, "MEParamWeight"}, + {gen::WeightType::kPartonShowerWeights, "PSWeight"}, + {gen::WeightType::kUnknownWeights, "UnknownWeight"}, + }; + //std::unordered_map weightGroupIndices_; + int lheWeightPrecision_; + bool keepAllPSWeights_; - double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1; - auto sumScales = runCounter->sumScale; - for (auto& val : sumScales) - val *= norm; - out->addVFloatWithNorm("LHEScaleSumw" + label, - "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel, - sumScales, - runCounter->sumw); - auto sumPDFs = runCounter->sumPDF; - for (auto& val : sumPDFs) - val *= norm; - out->addVFloatWithNorm("LHEPdfSumw" + label, - "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel, - sumPDFs, - runCounter->sumw); - if (!runCounter->sumRwgt.empty()) { - auto sumRwgts = runCounter->sumRwgt; - for (auto& val : sumRwgts) - val *= norm; - out->addVFloatWithNorm("LHEReweightingSumw" + label, - "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel, - sumRwgts, - runCounter->sumw); - } - 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->addFloatWithNorm( - "LHESumw_" + namedWeightLabels_[i] + label, - "Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[i] + ", divided by genEventSumw" + doclabel, - runCounter->sumNamed[i] * norm, - runCounter->sumw); - } - } + enum { inLHE, inGen }; +}; +//put back if needed; till now not used +LHEWeightsTableProducer::LHEWeightsTableProducer(edm::ParameterSet const& params) + : lheWeightTokens_( + edm::vector_transform(params.getParameter>("lheWeights"), + [this](const edm::InputTag& tag) { return mayConsume(tag); })), + lheWeightInfoTokens_(edm::vector_transform( + params.getParameter>("lheWeights"), + [this](const edm::InputTag& tag) { return mayConsume(tag); })), + genWeightToken_(consumes(params.getParameter("genWeights"))), + genWeightInfoToken_( + consumes(params.getParameter("genWeights"))), + genEventInfoToken_(consumes(params.getParameter("genEvent"))), + genLumiInfoHeadTag_( + mayConsume(params.getParameter("genLumiInfoHeader"))), + weightgroups_(edm::vector_transform(params.getParameter>("weightgroups"), + [](auto& c) { return gen::WeightType(c.at(0)); })), + maxGroupsPerType_(params.getParameter>("maxGroupsPerType")), + pdfIds_(params.getUntrackedParameter>("pdfIds", {})), + lheWeightPrecision_(params.getParameter("lheWeightPrecision")), + keepAllPSWeights_(params.getParameter("keepAllPSWeights")) { + if (weightgroups_.size() != maxGroupsPerType_.size()) + throw std::invalid_argument("Inputs 'weightgroups' and 'weightgroupNums' must have equal size"); + for (auto& wg : weightTypeNames_) { + produces(wg.second); + produces(wg.second + "sizes"); + } + produces("GENWeight"); + produces(); + produces("genModel"); +} + +void LHEWeightsTableProducer::produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + //access counter for weight sums + Counter& counter = *streamCache(id)->get(); + edm::Handle lheWeightHandle; + bool foundLheWeights = false; + for (auto& token : lheWeightTokens_) { + iEvent.getByToken(token, lheWeightHandle); + if (lheWeightHandle.isValid()) { + foundLheWeights = true; + break; } - iRun.put(std::move(out)); } - // nothing to do here - void globalEndRun(edm::Run const&, edm::EventSetup const&) const override {} + //Taken from genweight producer //Added sroychow + // generator information (always available) + auto const& genInfo = iEvent.get(genEventInfoToken_); + const double genWeight = genInfo.weight(); + // table for gen info, always available + auto outGeninfo = std::make_unique(1, "genWeight", true); + outGeninfo->setDoc("generator weight"); + outGeninfo->addColumnValue("", genInfo.weight(), "generator weight"); + iEvent.put(std::move(outGeninfo), "GENWeight"); + //this will take care of sum of genWeights + counter.incGenOnly(genWeight); + + std::string& model_label = streamCache(id)->getLabel(); + auto outM = std::make_unique((!model_label.empty()) ? std::string("GenModel_") + model_label : ""); + iEvent.put(std::move(outM), "genModel"); + + WeightsContainer lheWeights; + if (foundLheWeights) { + const GenWeightProduct* lheWeightProduct = lheWeightHandle.product(); + lheWeights = lheWeightProduct->weights(); + } - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - edm::ParameterSetDescription desc; - desc.add("genEvent", edm::InputTag("generator")) - ->setComment("tag for the GenEventInfoProduct, to get the main weight"); - desc.add("genLumiInfoHeader", edm::InputTag("generator")) - ->setComment("tag for the GenLumiInfoProduct, to get the model string"); - desc.add>("lheInfo", std::vector{{"externalLHEProducer"}, {"source"}}) - ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)"); + edm::Handle genWeightHandle; + iEvent.getByToken(genWeightToken_, genWeightHandle); + const GenWeightProduct* genWeightProduct = genWeightHandle.product(); + WeightsContainer genWeights = genWeightProduct->weights(); + + auto const& weightInfos = *luminosityBlockCache(iEvent.getLuminosityBlock().index()); + + //create a container with dummy weight vector + std::map> lheWeightTables; + std::map> weightVecsizes; + std::map weightlabels; + for (auto& wg : weightTypeNames_) { + lheWeightTables.insert(std::make_pair(wg.first, std::vector())); + weightVecsizes.insert(std::make_pair(wg.first, std::vector())); + weightlabels.insert(std::make_pair(wg.first, "")); + } + if (foundLheWeights) { + addWeightGroupToTable( + lheWeightTables, weightVecsizes, weightlabels, "LHE", weightInfos.at(inLHE), lheWeights, counter, genWeight); + } - edm::ParameterSetDescription prefpdf; - prefpdf.add("name"); - prefpdf.add("lhaid"); - desc.addVPSet("preferredPDFs", prefpdf, std::vector()) - ->setComment( - "LHA PDF Ids of the preferred PDF sets, in order of preference (the first matching one will be used)"); - desc.add>("namedWeightIDs")->setComment("set of LHA weight IDs for named LHE weights"); - desc.add>("namedWeightLabels") - ->setComment("output names for the namedWeightIDs (in the same order)"); - desc.add("lheWeightPrecision")->setComment("Number of bits in the mantissa for LHE weights"); - desc.add("maxPdfWeights")->setComment("Maximum number of PDF weights to save (to crop NN replicas)"); - desc.add("keepAllPSWeights")->setComment("Store all PS weights found"); - desc.addOptionalUntracked("debug")->setComment("dump out all LHE information for one event"); - descriptions.add("genWeightsTable", desc); + addWeightGroupToTable( + lheWeightTables, weightVecsizes, weightlabels, "Gen", weightInfos.at(inGen), genWeights, counter, genWeight); + + for (auto& wg : weightTypeNames_) { + std::string wname = wg.second; + auto& weightVec = lheWeightTables[wg.first]; + counter.incLHE(genWeight, weightVec, wname); + auto outTable = std::make_unique(weightVec.size(), wname, false); + outTable->addColumn("", weightVec, weightlabels[wg.first], lheWeightPrecision_); + + //now add the vector containing the sizes of alt sets + auto outTableSizes = + std::make_unique(weightVecsizes[wg.first].size(), wname + "_AltSetSizes", false); + outTableSizes->addColumn( + "", weightVecsizes[wg.first], "Sizes of weight arrays for weight type:" + wname, lheWeightPrecision_); + iEvent.put(std::move(outTable), wname); + iEvent.put(std::move(outTableSizes), wname + "sizes"); } +} + +/* + +*/ +void LHEWeightsTableProducer::addWeightGroupToTable(std::map>& lheWeightTables, + std::map>& weightVecsizes, + std::map& weightlabels, + const char* typeName, + const WeightGroupDataContainer& weightInfos, + WeightsContainer& allWeights, + Counter& counter, + double genWeight) const { + std::unordered_map typeCount = {}; + for (auto& type : gen::allWeightTypes) + typeCount[type] = 0; + + for (const auto& groupInfo : weightInfos) { + //std::string entryName = typeName; + gen::WeightType weightType = groupInfo.group->weightType(); + std::string name = weightTypeNames_.at(weightType); + std::string label = "[" + std::to_string(typeCount[weightType]) + "] " + groupInfo.group->description(); + label.append("["); + label.append(std::to_string(lheWeightTables[weightType].size())); //to append the start index of this set + label.append("]; "); + auto& weights = allWeights.at(groupInfo.index); + //std::cout << "Group name is " << groupInfo.group->name() << " is it wellFormed? " << groupInfo.group->isWellFormed() << std::endl; + if (weightType == gen::WeightType::kScaleWeights) { + if (groupInfo.group->isWellFormed()) { + const auto scaleGroup = *static_cast(groupInfo.group.get()); + std::cout << "They're well formed, will be ordered as expected\n"; + weights = orderedScaleWeights(weights, scaleGroup); + 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 { + std::cout << "NOT WELL FORMED!\n"; + size_t nstore = std::min(gen::ScaleWeightGroupInfo::MIN_SCALE_VARIATIONS, weights.size()); + weights = std::vector(weights.begin(), weights.begin() + nstore); + label.append("WARNING: Unexpected format found. Contains first " + std::to_string(nstore) + + " elements of weights vector, unordered"); + } + // TODO: Handle storeAllWeights and !isWellFormed + } else if (!keepAllPSWeights_ && weightType == gen::WeightType::kPartonShowerWeights && + groupInfo.group->isWellFormed()) { + const auto psGroup = *static_cast(groupInfo.group.get()); + weights = preferredPSweights(weights, psGroup); + label.append( + "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"); + } + //else + // label.append(groupInfo.group->description()); + lheWeightTables[weightType].insert(lheWeightTables[weightType].end(), weights.begin(), weights.end()); + weightVecsizes[weightType].emplace_back(weights.size()); -protected: - const edm::EDGetTokenT genTag_; - const std::vector lheLabel_; - const std::vector> lheTag_; - const std::vector> lheRunTag_; - const edm::EDGetTokenT genLumiInfoHeadTag_; + if (weightlabels[weightType].empty()) + weightlabels[weightType].append("[idx in AltSetSizes array] Name [start idx in weight array];\n"); - std::vector preferredPDFLHAIDs_; - std::unordered_map lhaNameToID_; - std::vector namedWeightIDs_; - std::vector namedWeightLabels_; - int lheWeightPrecision_; - unsigned int maxPdfWeights_; - bool keepAllPSWeights_; + weightlabels[weightType].append(label); + typeCount[weightType]++; + } +} + +WeightGroupDataContainer LHEWeightsTableProducer::weightDataPerType(edm::Handle& weightsHandle, + gen::WeightType weightType, + int& maxStore) const { + std::vector allgroups; + if (weightType == gen::WeightType::kPdfWeights && !pdfIds_.empty()) { + allgroups = weightsHandle->pdfGroupsWithIndicesByLHAIDs(pdfIds_); + } else + allgroups = weightsHandle->weightGroupsAndIndicesByType(weightType); + + int toStore = maxStore; + if (maxStore < 0 || static_cast(allgroups.size()) <= maxStore) { + // Modify size in case one type of weight is present in multiple products + maxStore -= allgroups.size(); + toStore = allgroups.size(); + } - mutable std::atomic debug_, debugRun_, hasIssuedWarning_, psWeightWarning_; -}; + WeightGroupDataContainer out; + for (int i = 0; i < toStore; i++) { + auto& group = allgroups.at(i); + gen::SharedWeightGroupData temp = {group.index, std::move(group.group)}; + out.push_back(temp); + } + return out; +} + +std::vector LHEWeightsTableProducer::orderedScaleWeights(const std::vector& scaleWeights, + const gen::ScaleWeightGroupInfo& scaleGroup) const { + std::vector weights; + weights.emplace_back(scaleWeights.at(scaleGroup.muR05muF05Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR05muF1Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR05muF2Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR1muF05Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.centralIndex())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR1muF2Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR2muF05Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR2muF1Index())); + weights.emplace_back(scaleWeights.at(scaleGroup.muR2muF2Index())); + + return weights; +} + +std::vector LHEWeightsTableProducer::preferredPSweights(const std::vector& psWeights, + const gen::PartonShowerWeightGroupInfo& pswV) const { + std::vector psTosave; + + double baseline = psWeights.at(pswV.weightIndexFromLabel("Baseline")); + psTosave.emplace_back(psWeights.at(pswV.variationIndex(true, true, gen::PSVarType::def)) / baseline); + psTosave.emplace_back(psWeights.at(pswV.variationIndex(false, true, gen::PSVarType::def)) / baseline); + psTosave.emplace_back(psWeights.at(pswV.variationIndex(true, false, gen::PSVarType::def)) / baseline); + psTosave.emplace_back(psWeights.at(pswV.variationIndex(false, false, gen::PSVarType::def)) / baseline); + return psTosave; +} + +void LHEWeightsTableProducer::streamEndRunSummary(edm::StreamID id, + edm::Run const&, + edm::EventSetup const&, + CounterMap* runCounterMap) const { + //this takes care for mergeing all the weight sums + runCounterMap->mergeSumMap(*streamCache(id)); +} + +void LHEWeightsTableProducer::globalEndRunProduce(edm::Run& iRun, + edm::EventSetup const&, + CounterMap const* runCounterMap) const { + auto out = std::make_unique(); + + for (auto x : runCounterMap->countermap) { + 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_); + + double norm = runCounter.sumw_ ? 1.0 / runCounter.sumw_ : 1; + //Sum from map + for (auto& sumw : runCounter.weightSumMap_) { + //Normalize with genEventSumw + for (auto& val : sumw.second) + val *= norm; + out->addVFloat(sumw.first + "Sumw" + label, + "Sum of genEventWeight *" + sumw.first + "[i]/genEventSumw" + doclabel, + sumw.second); + } + } + iRun.put(std::move(out)); +} +void LHEWeightsTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add>("lheWeights"); + desc.add>("lheInfo", std::vector{{"externalLHEProducer"}, {"source"}}) + ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)"); + //desc.add>("genWeights"); + desc.add("genWeights"); + desc.add("genEvent", edm::InputTag("generator")) + ->setComment("tag for the GenEventInfoProduct, to get the main weight"); + desc.add("genLumiInfoHeader", edm::InputTag("generator")) + ->setComment("tag for the GenLumiInfoProduct, to get the model string"); + desc.add>("weightgroups"); + desc.add>("maxGroupsPerType"); + desc.addOptionalUntracked>("pdfIds"); + desc.add("lheWeightPrecision", -1)->setComment("Number of bits in the mantissa for LHE weights"); + desc.add("keepAllPSWeights", false)->setComment("True:stores all 45 PS weights; False:saves preferred 4"); + descriptions.addDefault(desc); +} #include "FWCore/Framework/interface/MakerMacros.h" -DEFINE_FWK_MODULE(GenWeightsTableProducer); +DEFINE_FWK_MODULE(LHEWeightsTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/LHETablesProducer.cc b/PhysicsTools/NanoAOD/plugins/LHETablesProducer.cc index 735b78858e14b..877c5e8f7eb6f 100644 --- a/PhysicsTools/NanoAOD/plugins/LHETablesProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/LHETablesProducer.cc @@ -17,7 +17,8 @@ class LHETablesProducer : public edm::global::EDProducer<> { : lheTag_(edm::vector_transform(params.getParameter>("lheInfo"), [this](const edm::InputTag& tag) { return mayConsume(tag); })), precision_(params.getParameter("precision")), - storeLHEParticles_(params.getParameter("storeLHEParticles")) { + storeLHEParticles_(params.getParameter("storeLHEParticles")), + storeAllLHEInfo_(params.getParameter("storeAllLHEInfo")) { produces("LHE"); if (storeLHEParticles_) produces("LHEPart"); @@ -54,7 +55,9 @@ class LHETablesProducer : public edm::global::EDProducer<> { unsigned int lheNj = 0, lheNb = 0, lheNc = 0, lheNuds = 0, lheNglu = 0; double lheVpt = 0; double alphaS = 0; - + double alphaQED = 0; + double scale = 0; + int idproc = 0; const auto& hepeup = lheProd.hepeup(); const auto& pup = hepeup.PUP; int lep = -1, lepBar = -1, nu = -1, nuBar = -1; @@ -63,18 +66,33 @@ class LHETablesProducer : public edm::global::EDProducer<> { std::vector vals_phi; std::vector vals_mass; std::vector vals_pz; + std::vector vals_time; std::vector vals_pid; std::vector vals_status; std::vector vals_spin; + std::vector vals_col1; + std::vector vals_col2; + std::vector vals_mother1; + std::vector vals_mother2; alphaS = hepeup.AQCDUP; + alphaQED = hepeup.AQEDUP; + scale = hepeup.SCALUP; + idproc = hepeup.IDPRUP; for (unsigned int i = 0, n = pup.size(); i < n; ++i) { int status = hepeup.ISTUP[i]; int idabs = std::abs(hepeup.IDUP[i]); - if (status == 1 || status == -1) { + if (status == 1 || status == -1 || storeAllLHEInfo_) { TLorentzVector p4(pup[i][0], pup[i][1], pup[i][2], pup[i][3]); // x,y,z,t vals_pid.push_back(hepeup.IDUP[i]); vals_spin.push_back(hepeup.SPINUP[i]); vals_status.push_back(status); + if (storeAllLHEInfo_) { + vals_col1.push_back(hepeup.ICOLUP[i].first); + vals_col2.push_back(hepeup.ICOLUP[i].second); + vals_mother1.push_back(hepeup.MOTHUP[i].first); + vals_mother2.push_back(hepeup.MOTHUP[i].second); + vals_time.push_back(hepeup.VTIMUP[i]); + } if (status == -1) { vals_pt.push_back(0); vals_eta.push_back(0); @@ -143,7 +161,11 @@ class LHETablesProducer : public edm::global::EDProducer<> { out.addColumnValue("NpNLO", lheProd.npNLO(), "number of partons at NLO"); out.addColumnValue("NpLO", lheProd.npLO(), "number of partons at LO"); out.addColumnValue("AlphaS", alphaS, "Per-event alphaS"); - + if (storeAllLHEInfo_) { + out.addColumnValue("AlphaQED", alphaQED, "Per-event alphaQED"); + out.addColumnValue("Scale", scale, "Per-event scale"); + out.addColumnValue("ProcessID", idproc, "Process id (as in the card ordering)"); + } auto outPart = std::make_unique(vals_pt.size(), "LHEPart", false); outPart->addColumn("pt", vals_pt, "Pt of LHE particles", this->precision_); outPart->addColumn("eta", vals_eta, "Pseodorapidity of LHE particles", this->precision_); @@ -153,7 +175,13 @@ class LHETablesProducer : public edm::global::EDProducer<> { outPart->addColumn("pdgId", vals_pid, "PDG ID of LHE particles"); outPart->addColumn("status", vals_status, "LHE particle status; -1:incoming, 1:outgoing"); outPart->addColumn("spin", vals_spin, "Spin of LHE particles"); - + if (storeAllLHEInfo_) { + outPart->addColumn("color1", vals_col1, "First color index of LHE particles"); + outPart->addColumn("color2", vals_col2, "Second color index of LHE particles"); + outPart->addColumn("mother1", vals_mother1, "First mother index of LHE particles"); + outPart->addColumn("mother2", vals_mother2, "Second mother index of LHE particles"); + outPart->addColumn("lifetime", vals_time, "Own lifetime of LHE particles", this->precision_); + } return outPart; } @@ -164,6 +192,8 @@ class LHETablesProducer : public edm::global::EDProducer<> { desc.add("precision", -1)->setComment("precision on the 4-momenta of the LHE particles"); desc.add("storeLHEParticles", false) ->setComment("Whether we want to store the 4-momenta of the status 1 particles at LHE level"); + desc.add("storeAllLHEInfo", false) + ->setComment("Whether to store the whole set of intermediate LHE particles, not only the status +/-1 ones"); descriptions.add("lheInfoTable", desc); } @@ -171,6 +201,7 @@ class LHETablesProducer : public edm::global::EDProducer<> { const std::vector> lheTag_; const unsigned int precision_; const bool storeLHEParticles_; + const bool storeAllLHEInfo_; }; #include "FWCore/Framework/interface/MakerMacros.h" diff --git a/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc b/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc index 02198c149c17c..4a6999ed42733 100644 --- a/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc +++ b/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc @@ -48,7 +48,6 @@ class NanoAODOutputModule : public edm::one::OutputModule<> { public: NanoAODOutputModule(edm::ParameterSet const& pset); - ~NanoAODOutputModule() override; static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); @@ -120,6 +119,8 @@ class NanoAODOutputModule : public edm::one::OutputModule<> { } m_commonRunBranches; std::vector m_tables; + std::vector m_tableTokens; + std::vector m_tableVectorTokens; std::vector m_triggers; bool m_triggers_areSorted = false; std::vector m_evstrings; @@ -152,8 +153,6 @@ NanoAODOutputModule::NanoAODOutputModule(edm::ParameterSet const& pset) m_autoFlush(pset.getUntrackedParameter("autoFlush", -10000000)), m_processHistoryRegistry() {} -NanoAODOutputModule::~NanoAODOutputModule() {} - void NanoAODOutputModule::write(edm::EventForOutput const& iEvent) { //Get data from 'e' and write it to the file edm::Service jr; @@ -194,9 +193,29 @@ void NanoAODOutputModule::write(edm::EventForOutput const& iEvent) { m_commonBranches.fill(iEvent.id()); // fill all tables, starting from main tables and then doing extension tables + std::vector tables; + for (auto& tableToken : m_tableTokens) { + edm::Handle handle; + iEvent.getByToken(tableToken, handle); + tables.push_back(&(*handle)); + } + for (auto& tableToken : m_tableVectorTokens) { + edm::Handle> handle; + iEvent.getByToken(tableToken, handle); + for (auto const& table : *handle) { + tables.push_back(&table); + } + } + + if (m_tables.empty()) { + m_tables.resize(tables.size()); + } for (unsigned int extensions = 0; extensions <= 1; ++extensions) { - for (auto& t : m_tables) - t.fill(iEvent, *m_tree, extensions); + size_t iTable = 0; + for (auto& table : tables) { + m_tables[iTable].fill(*table, *m_tree, extensions); + ++iTable; + } } if (!m_triggers_areSorted) { // sort triggers/flags in inverse processHistory order, to save without any special label the most recent ones std::vector pnames; @@ -283,15 +302,18 @@ void NanoAODOutputModule::openFile(edm::FileBlock const&) { } /* Setup file structure here */ m_tables.clear(); + m_tableTokens.clear(); m_triggers.clear(); m_triggers_areSorted = false; m_evstrings.clear(); m_runTables.clear(); const auto& keeps = keptProducts(); for (const auto& keep : keeps[edm::InEvent]) { - if (keep.first->className() == "nanoaod::FlatTable") - m_tables.emplace_back(keep.first, keep.second); - else if (keep.first->className() == "edm::TriggerResults") { + if (keep.first->className() == "nanoaod::FlatTable") { + m_tableTokens.emplace_back(keep.second); + } else if (keep.first->className() == "std::vector") { + m_tableVectorTokens.emplace_back(keep.second); + } else if (keep.first->className() == "edm::TriggerResults") { m_triggers.emplace_back(keep.first, keep.second); } else if (keep.first->className() == "std::basic_string >" && keep.first->productInstanceName() == "genModel") { // friendlyClassName == "String" diff --git a/PhysicsTools/NanoAOD/plugins/TableOutputBranches.cc b/PhysicsTools/NanoAOD/plugins/TableOutputBranches.cc index 5744b329f7bbf..abf25ffaf831b 100644 --- a/PhysicsTools/NanoAOD/plugins/TableOutputBranches.cc +++ b/PhysicsTools/NanoAOD/plugins/TableOutputBranches.cc @@ -59,15 +59,12 @@ void TableOutputBranches::branch(TTree &tree) { } } -void TableOutputBranches::fill(const edm::EventForOutput &iEvent, TTree &tree, bool extensions) { +void TableOutputBranches::fill(const nanoaod::FlatTable &tab, TTree &tree, bool extensions) { if (m_extension != DontKnowYetIfMainOrExtension) { if (extensions != m_extension) return; // do nothing, wait to be called with the proper flag } - edm::Handle handle; - iEvent.getByToken(m_token, handle); - const nanoaod::FlatTable &tab = *handle; m_counter = tab.size(); m_singleton = tab.singleton(); if (!m_branchesBooked) { diff --git a/PhysicsTools/NanoAOD/plugins/TableOutputBranches.h b/PhysicsTools/NanoAOD/plugins/TableOutputBranches.h index df843bb6f4216..c8e4366e48539 100644 --- a/PhysicsTools/NanoAOD/plugins/TableOutputBranches.h +++ b/PhysicsTools/NanoAOD/plugins/TableOutputBranches.h @@ -4,31 +4,23 @@ #include #include #include -#include "FWCore/Framework/interface/EventForOutput.h" #include "DataFormats/NanoAOD/interface/FlatTable.h" #include "DataFormats/Provenance/interface/BranchDescription.h" #include "FWCore/Utilities/interface/EDGetToken.h" class TableOutputBranches { public: - TableOutputBranches(const edm::BranchDescription *desc, const edm::EDGetToken &token) - : m_token(token), m_extension(DontKnowYetIfMainOrExtension), m_branchesBooked(false) { - if (desc->className() != "nanoaod::FlatTable") - throw cms::Exception("Configuration", "NanoAODOutputModule can only write out nanoaod::FlatTable objects"); - } - void defineBranchesFromFirstEvent(const nanoaod::FlatTable &tab); void branch(TTree &tree); /// Fill the current table, if extensions == table.extension(). /// This parameter is used so that the fill is called first for non-extensions and then for extensions - void fill(const edm::EventForOutput &iEvent, TTree &tree, bool extensions); + void fill(const nanoaod::FlatTable &tab, TTree &tree, bool extensions); private: - edm::EDGetToken m_token; std::string m_baseName; bool m_singleton; - enum { IsMain = 0, IsExtension = 1, DontKnowYetIfMainOrExtension = 2 } m_extension; + enum { IsMain = 0, IsExtension = 1, DontKnowYetIfMainOrExtension = 2 } m_extension = DontKnowYetIfMainOrExtension; std::string m_doc; UInt_t m_counter; struct NamedBranchPtr { @@ -44,7 +36,7 @@ class TableOutputBranches { std::vector m_floatBranches; std::vector m_intBranches; std::vector m_uint8Branches; - bool m_branchesBooked; + bool m_branchesBooked = false; template void fillColumn(NamedBranchPtr &pair, const nanoaod::FlatTable &tab) { diff --git a/PhysicsTools/NanoAOD/python/NanoAODEDMEventContent_cff.py b/PhysicsTools/NanoAOD/python/NanoAODEDMEventContent_cff.py index d2e1f939d795e..6a23c316eb624 100644 --- a/PhysicsTools/NanoAOD/python/NanoAODEDMEventContent_cff.py +++ b/PhysicsTools/NanoAOD/python/NanoAODEDMEventContent_cff.py @@ -3,6 +3,7 @@ NanoAODEDMEventContent = cms.PSet( outputCommands = cms.untracked.vstring( 'drop *', + "keep *_lheWeightsTable_*_*", # event data "keep nanoaodFlatTable_*Table_*_*", # event data "keep edmTriggerResults_*_*_*", # event data "keep String_*_genModel_*", # generator model data @@ -19,12 +20,15 @@ compressionLevel = cms.untracked.int32(9), compressionAlgorithm = cms.untracked.string("LZMA"), ) - -NanoGenOutput = NanoAODEDMEventContent.outputCommands[:] -NanoGenOutput.remove("keep edmTriggerResults_*_*_*") - NANOAODGENEventContent = cms.PSet( + outputCommands = cms.untracked.vstring( + 'drop *', + "keep *_lheWeightsTable_*_*", # event data + "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 + ), compressionLevel = cms.untracked.int32(9), compressionAlgorithm = cms.untracked.string("LZMA"), - outputCommands = cms.untracked.vstring(NanoGenOutput) ) diff --git a/PhysicsTools/NanoAOD/python/genWeightsTable_cfi.py b/PhysicsTools/NanoAOD/python/genWeightsTable_cfi.py deleted file mode 100644 index e226f955e3f90..0000000000000 --- a/PhysicsTools/NanoAOD/python/genWeightsTable_cfi.py +++ /dev/null @@ -1,28 +0,0 @@ -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("NNPDF31_nnlo_hessian_pdfas"), lhaid = cms.uint32(306000) ), - cms.PSet( name = cms.string("NNPDF31_nnlo_as_0118_hessian"), lhaid = cms.uint32(304400) ), - cms.PSet( name = cms.string("NNPDF31_nnlo_as_0118_mc_hessian_pdfas"), lhaid = cms.uint32(325300) ), - cms.PSet( name = cms.string("NNPDF31_nnlo_as_0118_mc"), lhaid = cms.uint32(316200) ), - cms.PSet( name = cms.string("NNPDF31_nnlo_as_0118_nf_4_mc_hessian"), lhaid = cms.uint32(325500) ), - cms.PSet( name = cms.string("NNPDF31_nnlo_as_0118_nf_4"), lhaid = cms.uint32(320900) ), - 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) - cms.PSet( name = cms.string("PDF4LHC15_nnlo_30_pdfas"), lhaid = cms.uint32(91400) ), - cms.PSet( name = cms.string("PDF4LHC15_nlo_30_pdfas"), lhaid = cms.uint32(90400) ), - cms.PSet( name = cms.string("PDF4LHC15_nlo_30"), lhaid = cms.uint32(90900) ), - ), - namedWeightIDs = cms.vstring(), - namedWeightLabels = cms.vstring(), - lheWeightPrecision = cms.int32(14), - maxPdfWeights = cms.uint32(150), - keepAllPSWeights = cms.bool(False), - debug = cms.untracked.bool(False), -) diff --git a/PhysicsTools/NanoAOD/python/genWeights_cff.py b/PhysicsTools/NanoAOD/python/genWeights_cff.py new file mode 100644 index 0000000000000..ba5cd8c7bc4f1 --- /dev/null +++ b/PhysicsTools/NanoAOD/python/genWeights_cff.py @@ -0,0 +1,30 @@ +import FWCore.ParameterSet.Config as cms + +genWeights = cms.EDProducer("GenWeightProductProducer", + genInfo = cms.InputTag("generator"), + genLumiInfoHeader = cms.InputTag("generator")) + +lheWeights = cms.EDProducer("LHEWeightProductProducer", + lheSourceLabels = cms.vstring(["externalLHEProducer", "source"]), + failIfInvalidXML = cms.untracked.bool(True) +) + +genWeightsTable = cms.EDProducer( + "LHEWeightsTableProducer", + lheWeights = cms.VInputTag(["externalLHEProducer", "source", "lheWeights"]), + lheWeightPrecision = cms.int32(14), + genWeights = cms.InputTag("genWeights"), + # Warning: you can use a full string, but only the first character is read. + # Note also that the capitalization is important! For example, 'parton shower' + # must be lower case and 'PDF' must be capital + weightgroups = cms.vstring(['scale', 'PDF', 'matrix element', 'unknown', 'parton shower']), + # Max number of groups to store for each type above, -1 ==> store all found + maxGroupsPerType = cms.vint32([-1, -1, -1, -1, 1]), + # If empty or not specified, no critieria is applied to filter on LHAPDF IDs + #pdfIds = cms.untracked.vint32([91400, 306000, 260000]), + #unknownOnlyIfEmpty = cms.untracked.vstring(['scale', 'PDF']), + #unknownOnlyIfAllEmpty = cms.untracked.bool(False), + keepAllPSWeights = cms.bool(False) +) + +genWeightsTables = cms.Sequence(lheWeights*genWeights*genWeightsTable) diff --git a/PhysicsTools/NanoAOD/python/nano_cff.py b/PhysicsTools/NanoAOD/python/nano_cff.py index 89fbe3b2ca22f..b98faffe69fb2 100644 --- a/PhysicsTools/NanoAOD/python/nano_cff.py +++ b/PhysicsTools/NanoAOD/python/nano_cff.py @@ -11,8 +11,8 @@ from PhysicsTools.NanoAOD.ttbarCategorization_cff import * from PhysicsTools.NanoAOD.genparticles_cff import * from PhysicsTools.NanoAOD.particlelevel_cff import * -from PhysicsTools.NanoAOD.genWeightsTable_cfi import * from PhysicsTools.NanoAOD.genVertex_cff import * +from PhysicsTools.NanoAOD.genWeights_cff import * from PhysicsTools.NanoAOD.vertices_cff import * from PhysicsTools.NanoAOD.met_cff import * from PhysicsTools.NanoAOD.triggerObjects_cff import * @@ -102,7 +102,8 @@ lheInfoTable = cms.EDProducer("LHETablesProducer", lheInfo = cms.VInputTag(cms.InputTag("externalLHEProducer"), cms.InputTag("source")), precision = cms.int32(14), - storeLHEParticles = cms.bool(True) + storeLHEParticles = cms.bool(True), + storeAllLHEInfo = cms.bool(False), ) l1bits=cms.EDProducer("L1TriggerResultsConverter", src=cms.InputTag("gtStage2Digis"), legacyL1=cms.bool(False), diff --git a/PhysicsTools/NanoAOD/python/nanogen_cff.py b/PhysicsTools/NanoAOD/python/nanogen_cff.py index 8bb4f31e9322a..53978e512e11d 100644 --- a/PhysicsTools/NanoAOD/python/nanogen_cff.py +++ b/PhysicsTools/NanoAOD/python/nanogen_cff.py @@ -5,9 +5,9 @@ from PhysicsTools.NanoAOD.genparticles_cff import * from PhysicsTools.NanoAOD.particlelevel_cff import * from PhysicsTools.NanoAOD.lheInfoTable_cfi import * -from PhysicsTools.NanoAOD.genWeightsTable_cfi import * from PhysicsTools.NanoAOD.genVertex_cff import * from PhysicsTools.NanoAOD.common_cff import Var,CandVars +from PhysicsTools.NanoAOD.genWeights_cff import * nanoMetadata = cms.EDProducer("UniqueStringProducer", strings = cms.PSet( @@ -36,13 +36,14 @@ rivetProducerHTXS+ particleLevelTables+ metMCTable+ - genWeightsTable+ + genWeightsTables+ lheInfoTable ) def nanoGenCommonCustomize(process): process.rivetMetTable.extension = False process.lheInfoTable.storeLHEParticles = True + process.lheInfoTable.storeAllLHEInfo = True process.lheInfoTable.precision = 14 process.genWeightsTable.keepAllPSWeights = True process.genJetFlavourAssociation.jets = process.genJetTable.src diff --git a/PhysicsTools/NanoAOD/test/testNanoWeights.py b/PhysicsTools/NanoAOD/test/testNanoWeights.py new file mode 100644 index 0000000000000..2f288b8629c91 --- /dev/null +++ b/PhysicsTools/NanoAOD/test/testNanoWeights.py @@ -0,0 +1,29 @@ +import ROOT +import argparse + +def variableAndNumber(varName, tree): + countVar = "n"+varName + if not hasattr(tree, varName): + print("No variable %s found in file" % varName) + else: + count = getattr(tree, countVar) + var = getattr(tree, varName) + print("Found %i entries of %s in file" % (count, varName)) + branch = tree.GetBranch(varName) + print(" --> Desciption:%s" % branch.GetTitle()) + +parser = argparse.ArgumentParser() +parser.add_argument('inputFile', type=str, help='NanoAOD file to process') +parser.add_argument('--scan', action='store_true', help='Scan the weight values') +args = parser.parse_args() + +rtfile = ROOT.TFile(args.inputFile) +tree = rtfile.Get("Events") +tree.GetEntry(0) +variables = ["LHEScaleWeight", "LHEPdfWeight", "MEParamWeight", "UnknownWeight", "PSWeight", ] + +for varName in variables: + variableAndNumber(varName, tree) + +if args.scan: + tree.Scan(":".join(variables)) diff --git a/SimDataFormats/GeneratorProducts/BuildFile.xml b/SimDataFormats/GeneratorProducts/BuildFile.xml index b35ce8e95f1ff..b35e098a0995a 100644 --- a/SimDataFormats/GeneratorProducts/BuildFile.xml +++ b/SimDataFormats/GeneratorProducts/BuildFile.xml @@ -3,6 +3,7 @@ + diff --git a/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h b/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h new file mode 100644 index 0000000000000..cdccd169797c9 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h @@ -0,0 +1,53 @@ +#ifndef SimDataFormats_GeneratorProducts_GenWeightInfoProduct_h +#define SimDataFormats_GeneratorProducts_GenWeightInfoProduct_h + +#include +#include +#include +#include +#include + +//#include + +#include "DataFormats/Common/interface/OwnVector.h" +#include "SimDataFormats/GeneratorProducts/interface/LesHouches.h" +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" + +namespace gen { + struct WeightGroupData { + size_t index; + std::unique_ptr group; + }; + + struct SharedWeightGroupData { + size_t index; + std::shared_ptr group; + }; +} // namespace gen + +class GenWeightInfoProduct { +public: + GenWeightInfoProduct() {} + GenWeightInfoProduct(edm::OwnVector& weightGroups); + GenWeightInfoProduct(const GenWeightInfoProduct& other); + GenWeightInfoProduct(GenWeightInfoProduct&& other); + ~GenWeightInfoProduct() {} + GenWeightInfoProduct& operator=(const GenWeightInfoProduct& other); + GenWeightInfoProduct& operator=(GenWeightInfoProduct&& other); + + const edm::OwnVector& allWeightGroupsInfo() const; + std::unique_ptr containingWeightGroupInfo(int index) const; + std::unique_ptr orderedWeightGroupInfo(int index) const; + std::vector weightGroupsByType(gen::WeightType type) const; + std::vector weightGroupIndicesByType(gen::WeightType type) const; + std::vector weightGroupsAndIndicesByType(gen::WeightType type) const; + std::optional pdfGroupWithIndexByLHAID(int lhaid) const; + std::vector pdfGroupsWithIndicesByLHAIDs(const std::vector& lhaids) const; + void addWeightGroupInfo(std::unique_ptr info); + const int numberOfGroups() const { return weightGroupsInfo_.size(); } + +private: + edm::OwnVector weightGroupsInfo_; +}; + +#endif // GeneratorWeightInfo_LHEInterface_GenWeightInfoProduct_h diff --git a/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h b/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h new file mode 100644 index 0000000000000..601b37ca46bc9 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h @@ -0,0 +1,53 @@ +#ifndef SimDataFormats_GeneratorProducts_GenWeightProduct_h +#define SimDataFormats_GeneratorProducts_GenWeightProduct_h + +#include +#include +#include +#include + +#include "SimDataFormats/GeneratorProducts/interface/LesHouches.h" +#include "SimDataFormats/GeneratorProducts/interface/WeightsInfo.h" +#include "FWCore/Utilities/interface/Exception.h" + +typedef std::vector> WeightsContainer; + +class GenWeightProduct { +public: + GenWeightProduct() { + weightsVector_ = {}; + centralWeight_ = 1.; + } + GenWeightProduct(double w0) { + weightsVector_ = {}; + centralWeight_ = w0; + } + GenWeightProduct& operator=(GenWeightProduct&& other) { + weightsVector_ = std::move(other.weightsVector_); + centralWeight_ = other.centralWeight_; + return *this; + } + ~GenWeightProduct() {} + + void setNumWeightSets(int num) { weightsVector_.resize(num); } + void addWeightSet() { weightsVector_.push_back({}); } + void addWeight(double weight, int setEntry, int weightNum) { + if (weightsVector_.empty() && setEntry == 0) + addWeightSet(); + if (static_cast(weightsVector_.size()) <= setEntry) + throw cms::Exception("GenWeightProduct") << "Trying to add weight index outside the range of weights expected"; + auto& weights = weightsVector_.at(setEntry); + if (static_cast(weights.size()) <= weightNum) { + weights.resize(weightNum + 1); + } + 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/LHEEventProduct.h b/SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h index c2af82bf3ab90..b2b8e1358b08f 100644 --- a/SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h +++ b/SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h @@ -27,7 +27,7 @@ class LHEEventProduct { ~LHEEventProduct() = default; - void setPDF(const PDF &pdf) { pdf_.reset(new PDF(pdf)); } + void setPDF(const PDF &pdf) { pdf_ = std::make_unique(pdf); } void addWeight(const WGT &wgt) { weights_.push_back(wgt); } void addComment(const std::string &line) { comments_.push_back(line); } diff --git a/SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h new file mode 100644 index 0000000000000..2b6714b240542 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h @@ -0,0 +1,33 @@ +#ifndef SimDataFormats_GeneratorProducts_MEParamWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_MEParamWeightGroupInfo_h + +#include +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" + +namespace gen { + class MEParamWeightGroupInfo : public WeightGroupInfo { + public: + MEParamWeightGroupInfo() : WeightGroupInfo() { weightType_ = WeightType::kMEParamWeights; } + MEParamWeightGroupInfo(std::string header, std::string name) : WeightGroupInfo(header, name) { + weightType_ = WeightType::kMEParamWeights; + } + MEParamWeightGroupInfo(std::string header) : MEParamWeightGroupInfo(header, header) {} + ~MEParamWeightGroupInfo() override {} + void copy(const MEParamWeightGroupInfo& other); + MEParamWeightGroupInfo* clone() const override; + int getCentralIndex() { return centralIdx; } + int getVariationIndex(int sig) { return massValue.at(numSigma + sig).second; } + double getVariationValue(int sig) { return massValue.at(numSigma + sig).first; } + void updateWeight(int globalIndex, std::string id, double weight); + void updateWeight(int globalIndex, std::string id, std::vector& content); + + private: + std::unordered_map> splitContent; + std::vector> massValue; + double central; + int centralIdx; + int numSigma; + }; +} // namespace gen + +#endif // SimDataFormats_GeneratorProducts_MEParamWeightGroupInfo_h diff --git a/SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h new file mode 100644 index 0000000000000..feb8face8658b --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h @@ -0,0 +1,59 @@ +#ifndef SimDataFormats_GeneratorProducts_PartonShowerWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_PartonShowerWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" + +namespace gen { + enum class PSVarType { muR, cNS, con, def, red, alphaS, LAST }; + enum class PSSplittingType { combined, g2gg, x2xg, g2qq, q2qg }; + + class PartonShowerWeightGroupInfo : public WeightGroupInfo { + public: + PartonShowerWeightGroupInfo() : PartonShowerWeightGroupInfo("") {} + PartonShowerWeightGroupInfo(std::string header, std::string name); + PartonShowerWeightGroupInfo(std::string header) : PartonShowerWeightGroupInfo(header, header) {} + PartonShowerWeightGroupInfo(const PartonShowerWeightGroupInfo &other) { copy(other); } + ~PartonShowerWeightGroupInfo() override {} + void copy(const PartonShowerWeightGroupInfo &other); + PartonShowerWeightGroupInfo *clone() const override; + void setNameIsPythiaSyntax(bool isPythiaSyntax) { nameIsPythiaSyntax_ = isPythiaSyntax; } + bool nameIsPythiaSyntax() const { return nameIsPythiaSyntax_; } + int variationIndex(bool isISR, bool isUp, PSVarType variationType, PSSplittingType splittingType) const; + std::string variationName(bool isISR, bool isUp, PSVarType variationType, PSSplittingType splittingType) const; + int variationIndex(bool isISR, bool isUp, PSVarType variationType) const; + static void setGuessPSWeightIdx(bool guessPSWeightIdx) { guessPSWeightIdx_ = guessPSWeightIdx; } + int psWeightIdxGuess(const std::string &varName) const; + void printVariables() const; + + private: + bool nameIsPythiaSyntax_ = false; + static inline bool guessPSWeightIdx_ = false; + + const std::vector expectedOrderPythiaSyntax_ = { + "fsr:murfac=0.707", "fsr:murfac=1.414", "fsr:murfac=0.5", "fsr:murfac=2.0", + "fsr:murfac=0.25", "fsr:murfac=4.0", "fsr:g2gg:murfac=0.5", "fsr:g2gg:murfac=2.0", + "fsr:g2qq:murfac=0.5", "fsr:g2qq:murfac=2.0", "fsr:q2qg:murfac=0.5", "fsr:q2qg:murfac=2.0", + "fsr:x2xg:murfac=0.5", "fsr:x2xg:murfac=2.0", "fsr:g2gg:cns=-2.0", "fsr:g2gg:cns=2.0", + "fsr:g2qq:cns=-2.0", "fsr:g2qq:cns=2.0", "fsr:q2qg:cns=-2.0", "fsr:q2qg:cns=2.0", + "fsr:x2xg:cns=-2.0", "fsr:x2xg:cns=2.0", "isr:murfac=0.707", "isr:murfac=1.414", + "isr:murfac=0.5", "isr:murfac=2.0", "isr:murfac=0.25", "isr:murfac=4.0", + "isr:g2gg:murfac=0.5", "isr:g2gg:murfac=2.0", "isr:g2qq:murfac=0.5", "isr:g2qq:murfac=2.0", + "isr:q2qg:murfac=0.5", "isr:q2qg:murfac=2.0", "isr:x2xg:murfac=0.5", "isr:x2xg:murfac=2.0", + "isr:g2gg:cns=-2.0", "isr:g2gg:cns=2.0", "isr:g2qq:cns=-2.0", "isr:g2qq:cns=2.0", + "isr:q2qg:cns=-2.0", "isr:q2qg:cns=2.0", "isr:x2xg:cns=-2.0", "isr:x2xg:cns=2.0", + }; + const std::vector expectedOrder_ = { + "isrRedHi", "fsrRedHi", "isrRedLo", "fsrRedLo", "isrDefHi", + "fsrDefHi", "isrDefLo", "fsrDefLo", "isrConHi", "fsrConHi", + "isrConLo", "fsrConLo", "fsr_G2GG_muR_dn", "fsr_G2GG_muR_up", "fsr_G2QQ_muR_dn", + "fsr_G2QQ_muR_up", "fsr_Q2QG_muR_dn", "fsr_Q2QG_muR_up", "fsr_X2XG_muR_dn", "fsr_X2XG_muR_up", + "fsr_G2GG_cNS_dn", "fsr_G2GG_cNS_up", "fsr_G2QQ_cNS_dn", "fsr_G2QQ_cNS_up", "fsr_Q2QG_cNS_dn", + "fsr_Q2QG_cNS_up", "fsr_X2XG_cNS_dn", "fsr_X2XG_cNS_up", "isr_G2GG_muR_dn", "isr_G2GG_muR_up", + "isr_G2QQ_muR_dn", "isr_G2QQ_muR_up", "isr_Q2QG_muR_dn", "isr_Q2QG_muR_up", "isr_X2XG_muR_dn", + "isr_X2XG_muR_up", "isr_G2GG_cNS_dn", "isr_G2GG_cNS_up", "isr_G2QQ_cNS_dn", "isr_G2QQ_cNS_up", + "isr_Q2QG_cNS_dn", "isr_Q2QG_cNS_up", "isr_X2XG_cNS_dn", "isr_X2XG_cNS_up", + }; + }; +} // namespace gen + +#endif diff --git a/SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h new file mode 100644 index 0000000000000..5c8214fa9eb35 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h @@ -0,0 +1,62 @@ +#ifndef SimDataFormats_GeneratorProducts_PdfWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_PdfWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include "LHAPDF/LHAPDF.h" +#include +#include +#include + +namespace gen { + enum PdfUncertaintyType { + kHessianUnc, + kMonteCarloUnc, + kVariationSet, + kUnknownUnc, + }; + + class PdfWeightGroupInfo : public WeightGroupInfo { + private: + PdfUncertaintyType uncertaintyType_; + bool hasAlphasVars_; + int alphasUpIndex_; + int alphasDownIndex_; + int parentLhapdfId_ = -1; + size_t parentLhapdfSize_ = -1; + std::string parentLhapdfError_; + std::vector lhaids_; + int parentLhapdfId(int lhaid) const { return lhaid - LHAPDF::lookupPDF(lhaid).second; } + + public: + PdfWeightGroupInfo() : WeightGroupInfo() { weightType_ = WeightType::kPdfWeights; } + PdfWeightGroupInfo(std::string header, std::string name) : WeightGroupInfo(header, name) { + weightType_ = WeightType::kPdfWeights; + } + PdfWeightGroupInfo(std::string header) : WeightGroupInfo(header) { weightType_ = WeightType::kPdfWeights; } + PdfWeightGroupInfo(const PdfWeightGroupInfo& other) { copy(other); } + ~PdfWeightGroupInfo() override {} + void copy(const PdfWeightGroupInfo& other); + PdfWeightGroupInfo* clone() const override; + + void setUncertaintyType(PdfUncertaintyType uncertaintyType) { uncertaintyType_ = uncertaintyType; } + void setHasAlphasVariations(bool hasAlphasVars) { hasAlphasVars_ = hasAlphasVars; } + void setAlphasUpIndex(int alphasUpIndex) { alphasUpIndex_ = alphasUpIndex; } + void setAlphasDownIndex(int alphasDownIndex) { alphasDownIndex_ = alphasDownIndex; } + PdfUncertaintyType uncertaintyType() const { return uncertaintyType_; } + bool hasAlphasVariations() const { return hasAlphasVars_; } + void addLhaid(int lhaid); + std::vector& lhaIds() { return lhaids_; } + + bool isIdInParentSet(int lhaid) const { return parentLhapdfId_ == parentLhapdfId(lhaid); } + int parentLhapdfId() const { return parentLhapdfId_; } + void setParentLhapdfInfo(int lhaid); + + // need to remove + bool containsLhapdfId(int lhaid) const { return isIdInParentSet(lhaid); } + + int alphasUpIndex() const { return alphasUpIndex_; } + int alphasDownIndex() const { return alphasDownIndex_; } + }; +} // namespace gen + +#endif // SimDataFormats_GeneratorProducts_PdfWeightGroupInfo_h diff --git a/SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h new file mode 100644 index 0000000000000..7bb460f57d210 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h @@ -0,0 +1,93 @@ +#ifndef SimDataFormats_GeneratorProducts_ScaleWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_ScaleWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include +#include + +namespace gen { + class ScaleWeightGroupInfo : public WeightGroupInfo { + private: + bool isFunctionalFormVar_; + std::vector muIndices_; + bool containsCentral_ = false; + int lhaid_ = -1; + bool weightIsCorrupt_ = false; + // Dyn_scale + std::vector dynNames_; + std::vector> dynVec_; + + inline int indexFromMus(int muR, int muF) const { return 3 * muR + muF; } + inline bool isValidValue(float mu) const { return mu == 0.5 || mu == 1.0 || mu == 2.0; } + + enum scaleIndices { + muR0p5_muF0p5_idx = 0, + muR0p5_muF1_idx = 1, + muR0p5_muF2_idx = 2, + muR1_muF0p5_idx = 3, + Central_idx = 4, + muR1_muF2_idx = 5, + muR2_muF0p5_idx = 6, + muR2_muF1_idx = 7, + muR2_muF2_idx = 8 + }; + + public: + static const unsigned int MIN_SCALE_VARIATIONS = 9; + ScaleWeightGroupInfo() : ScaleWeightGroupInfo("") {} + ScaleWeightGroupInfo(std::string header, std::string name) + : WeightGroupInfo(header, name), muIndices_(MIN_SCALE_VARIATIONS, -1), dynVec_(MIN_SCALE_VARIATIONS) { + weightType_ = WeightType::kScaleWeights; + isFunctionalFormVar_ = false; + } + ScaleWeightGroupInfo(std::string header) : ScaleWeightGroupInfo(header, header) {} + ScaleWeightGroupInfo(const ScaleWeightGroupInfo& other) { copy(other); } + ~ScaleWeightGroupInfo() override {} + void copy(const ScaleWeightGroupInfo& other); + ScaleWeightGroupInfo* clone() const override; + bool containsCentralWeight() const { return containsCentral_; } + void addContainedId(int globalIndex, std::string id, std::string label, float muR, float muF); + void setWeightIsCorrupt() { + isWellFormed_ = false; + weightIsCorrupt_ = true; + } + + void setMuRMuFIndex(int globalIndex, std::string id, float muR, float muF); + void setDyn(int globalIndex, std::string id, float muR, float muF, size_t dynNum, std::string dynName); + int lhaid() { return lhaid_; } + void setLhaid(int lhaid) { lhaid_ = lhaid; } + // Is a variation of the functional form of the dynamic scale + bool isFunctionalFormVariation(); + void setIsFunctionalFormVariation(bool functionalVar) { isFunctionalFormVar_ = functionalVar; } + int centralIndex() const { return muIndices_.at(Central_idx); } + int muR1muF2Index() const { return muIndices_.at(muR1_muF2_idx); } + int muR1muF05Index() const { return muIndices_.at(muR1_muF0p5_idx); } + int muR2muF05Index() const { return muIndices_.at(muR2_muF0p5_idx); } + int muR2muF1Index() const { return muIndices_.at(muR2_muF1_idx); } + int muR2muF2Index() const { return muIndices_.at(muR2_muF2_idx); } + int muR05muF05Index() const { return muIndices_.at(muR0p5_muF0p5_idx); } + int muR05muF1Index() const { return muIndices_.at(muR0p5_muF1_idx); } + int muR05muF2Index() const { return muIndices_.at(muR0p5_muF2_idx); } + + // dynweight version + size_t centralIndex(std::string& dynName) const { return scaleIndex(Central_idx, dynName); } + size_t muR1muF2Index(std::string& dynName) const { return scaleIndex(muR1_muF2_idx, dynName); } + size_t muR1muF05Index(std::string& dynName) const { return scaleIndex(muR1_muF0p5_idx, dynName); } + size_t muR2muF05Index(std::string& dynName) const { return scaleIndex(muR2_muF0p5_idx, dynName); } + size_t muR2muF1Index(std::string& dynName) const { return scaleIndex(muR2_muF1_idx, dynName); } + size_t muR2muF2Index(std::string& dynName) const { return scaleIndex(muR2_muF2_idx, dynName); } + size_t muR05muF05Index(std::string& dynName) const { return scaleIndex(muR0p5_muF0p5_idx, dynName); } + size_t muR05muF1Index(std::string& dynName) const { return scaleIndex(muR0p5_muF1_idx, dynName); } + size_t muR05muF2Index(std::string& dynName) const { return scaleIndex(muR0p5_muF2_idx, dynName); } + + size_t scaleIndex(float muR, float muF, size_t dynNum) const; + size_t scaleIndex(float muR, float muF, std::string& dynName) const; + size_t scaleIndex(int index, std::string& dynName) const; + size_t scaleIndex(float muR, float muF) const; + + size_t scaleIndex(int index, size_t dynNum) const { return dynVec_.at(index).at(dynNum); } + std::vector dynNames() const; + }; +} // namespace gen + +#endif diff --git a/SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h new file mode 100644 index 0000000000000..3513cc70841dc --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h @@ -0,0 +1,24 @@ +#ifndef SimDataFormats_GeneratorProducts_UnknownWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_UnknownWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" + +namespace gen { + class UnknownWeightGroupInfo : public WeightGroupInfo { + public: + UnknownWeightGroupInfo() : WeightGroupInfo() { weightType_ = WeightType::kUnknownWeights; } + UnknownWeightGroupInfo(std::string header, std::string name) : WeightGroupInfo(header, name) { + weightType_ = WeightType::kUnknownWeights; + isWellFormed_ = false; + } + UnknownWeightGroupInfo(std::string header) : WeightGroupInfo(header) { + weightType_ = WeightType::kUnknownWeights; + isWellFormed_ = false; + } + ~UnknownWeightGroupInfo() override {} + void copy(const UnknownWeightGroupInfo& other); + UnknownWeightGroupInfo* clone() const override; + }; +} // namespace gen + +#endif // SimDataFormats_GeneratorProducts_UnknownWeightGroupInfo_h diff --git a/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h new file mode 100644 index 0000000000000..e9058546d584c --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h @@ -0,0 +1,103 @@ +#ifndef SimDataFormats_GeneratorProducts_WeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_WeightGroupInfo_h + +/** \class PdfInfo + * + */ +#include +#include +#include +#include +#include + +namespace gen { + struct WeightMetaInfo { + size_t globalIndex; + size_t localIndex; + std::string id; + std::string label; + + bool operator==(const WeightMetaInfo& other) { + return (other.globalIndex == globalIndex && other.localIndex == localIndex && other.id == id && + other.label == label); + } + }; + + enum class WeightType : char { + kPdfWeights = 'P', + kScaleWeights = 's', + kMEParamWeights = 'm', + kPartonShowerWeights = 'p', + kUnknownWeights = 'u', + }; + + const std::array allWeightTypes = {{ + WeightType::kPdfWeights, + WeightType::kScaleWeights, + WeightType::kMEParamWeights, + WeightType::kPartonShowerWeights, + WeightType::kUnknownWeights, + }}; + + class WeightGroupInfo { + public: + WeightGroupInfo() : WeightGroupInfo("") {} + WeightGroupInfo(std::string header, std::string name) + : isWellFormed_(false), headerEntry_(header), name_(name), firstId_(-1), lastId_(-1) {} + WeightGroupInfo(std::string header) + : isWellFormed_(false), headerEntry_(header), name_(header), firstId_(-1), lastId_(-1) {} + WeightGroupInfo(const WeightGroupInfo& other) { copy(other); } + WeightGroupInfo& operator=(const WeightGroupInfo& other) { + copy(other); + return *this; + } + virtual ~WeightGroupInfo(){}; + void copy(const WeightGroupInfo& other); + virtual WeightGroupInfo* clone() const; + WeightMetaInfo weightMetaInfo(int weightEntry) const; + WeightMetaInfo weightMetaInfoByGlobalIndex(std::string wgtId, int weightEntry) const; + int weightVectorEntry(std::string& wgtId) const; + bool containsWeight(std::string& wgtId, int weightEntry) const; + int weightVectorEntry(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 appendDescription(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; } + // Call before doing lots of searches by label + void cacheWeightIndicesByLabel(); + + 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 + void setIsWellFormed(bool wellFormed) { isWellFormed_ = wellFormed; } + bool isWellFormed() const { return isWellFormed_; } + int weightIndexFromLabel(std::string weightLabel) const; + std::vector weightLabels() const; + + protected: + bool isWellFormed_ = false; + std::string headerEntry_; + std::string name_; + std::string description_; + WeightType weightType_; + std::vector idsContained_; + int firstId_; + int lastId_; + std::unordered_map weightLabelsToIndices_; + }; +} // namespace gen + +#endif // SimDataFormats_GeneratorProducts_WeightGroupInfo_h diff --git a/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc b/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc new file mode 100644 index 0000000000000..f54ced13fc6d7 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc @@ -0,0 +1,105 @@ +#include +#include + +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "FWCore/Utilities/interface/Exception.h" + +GenWeightInfoProduct::GenWeightInfoProduct(edm::OwnVector& weightGroups) { + weightGroupsInfo_ = weightGroups; +} + +GenWeightInfoProduct& GenWeightInfoProduct::operator=(const GenWeightInfoProduct& other) { + weightGroupsInfo_ = other.weightGroupsInfo_; + return *this; +} + +GenWeightInfoProduct& GenWeightInfoProduct::operator=(GenWeightInfoProduct&& other) { + weightGroupsInfo_ = std::move(other.weightGroupsInfo_); + return *this; +} + +const edm::OwnVector& GenWeightInfoProduct::allWeightGroupsInfo() const { + return weightGroupsInfo_; +} + +std::unique_ptr GenWeightInfoProduct::containingWeightGroupInfo(int index) const { + for (const auto& weightGroup : weightGroupsInfo_) { + if (weightGroup.indexInRange(index)) + return std::unique_ptr(weightGroup.clone()); + } + throw cms::Exception("GenWeightInfoProduct") << "No weight group found containing the weight index requested"; +} + +std::unique_ptr GenWeightInfoProduct::orderedWeightGroupInfo(int weightGroupIndex) const { + if (weightGroupIndex >= static_cast(weightGroupsInfo_.size())) + throw cms::Exception("GenWeightInfoProduct") + << "Weight index requested is outside the range of weights in the product"; + return std::unique_ptr(weightGroupsInfo_[weightGroupIndex].clone()); +} + +std::vector GenWeightInfoProduct::weightGroupsAndIndicesByType(gen::WeightType type) const { + std::vector matchingGroups; + for (size_t i = 0; i < weightGroupsInfo_.size(); i++) { + const gen::WeightGroupInfo& group = weightGroupsInfo_[i]; + if (weightGroupsInfo_[i].weightType() == type) + matchingGroups.push_back({i, std::unique_ptr(group.clone())}); + } + return matchingGroups; +} + +std::vector GenWeightInfoProduct::weightGroupsByType(gen::WeightType type) const { + std::vector matchingGroups; + for (size_t i = 0; i < weightGroupsInfo_.size(); i++) { + const gen::WeightGroupInfo& group = weightGroupsInfo_[i]; + if (weightGroupsInfo_[i].weightType() == type) + matchingGroups.push_back({i, std::unique_ptr(group.clone())}); + } + return matchingGroups; +} + +std::optional GenWeightInfoProduct::pdfGroupWithIndexByLHAID(int lhaid) const { + std::vector pdfGroups = weightGroupsAndIndicesByType(gen::WeightType::kPdfWeights); + + auto matchingPdfSet = std::find_if(pdfGroups.begin(), pdfGroups.end(), [lhaid](gen::WeightGroupData& data) { + auto pdfGroup = std::unique_ptr( + static_cast(data.group.release())); + return pdfGroup->containsLhapdfId(lhaid); + }); + + return matchingPdfSet == pdfGroups.end() + ? std::nullopt + : std::optional({matchingPdfSet->index, std::move(matchingPdfSet->group)}); +} + +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 = std::unique_ptr( + static_cast(data.group.release())); + return pdfGroup->containsLhapdfId(lhaid); + }); + if (matchingPdfSet != pdfGroups.end()) { + pdfGroups.push_back({matchingPdfSet->index, std::move(matchingPdfSet->group)}); + } + } + + return pdfGroups; +} + +std::vector GenWeightInfoProduct::weightGroupIndicesByType(gen::WeightType type) const { + std::vector matchingGroupIndices; + for (size_t i = 0; i < weightGroupsInfo_.size(); i++) { + if (weightGroupsInfo_[i].weightType() == type) + matchingGroupIndices.push_back(i); + } + return matchingGroupIndices; +} + +void GenWeightInfoProduct::addWeightGroupInfo(std::unique_ptr info) { + weightGroupsInfo_.push_back(std::move(info)); +} diff --git a/SimDataFormats/GeneratorProducts/src/MEParamWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/MEParamWeightGroupInfo.cc new file mode 100644 index 0000000000000..f3053ef5a1624 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/MEParamWeightGroupInfo.cc @@ -0,0 +1,27 @@ +#include "SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h" +#include + +namespace gen { + void MEParamWeightGroupInfo::copy(const MEParamWeightGroupInfo& other) { WeightGroupInfo::copy(other); } + + MEParamWeightGroupInfo* MEParamWeightGroupInfo::clone() const { return new MEParamWeightGroupInfo(*this); } + + void MEParamWeightGroupInfo::updateWeight(int globalIndex, std::string id, double weight) { + size_t localIndex = weightMetaInfoByGlobalIndex(id, globalIndex).localIndex; + auto lower = + std::lower_bound(massValue.begin(), massValue.end(), std::make_pair(weight, std::numeric_limits::min())); + massValue.insert(lower, std::make_pair(weight, localIndex)); + isWellFormed_ = massValue.size() % 2 == 1; + if (isWellFormed_) { + numSigma = massValue.size() / 2; + central = massValue.at(centralIdx).first; + centralIdx = massValue.at(centralIdx).second; + } + } + + void MEParamWeightGroupInfo::updateWeight(int globalIndex, std::string id, std::vector& content) { + size_t localIndex = weightMetaInfoByGlobalIndex(id, globalIndex).localIndex; + splitContent[localIndex] = content; + } + +} // namespace gen diff --git a/SimDataFormats/GeneratorProducts/src/PartonShowerWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/PartonShowerWeightGroupInfo.cc new file mode 100644 index 0000000000000..d54d00212e7b1 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/PartonShowerWeightGroupInfo.cc @@ -0,0 +1,121 @@ +#include "SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h" +#include +#include + +namespace gen { + PartonShowerWeightGroupInfo::PartonShowerWeightGroupInfo(std::string header, std::string name) + : WeightGroupInfo(header, name) { + weightType_ = WeightType::kPartonShowerWeights; + } + + void PartonShowerWeightGroupInfo::copy(const PartonShowerWeightGroupInfo& other) { + WeightGroupInfo::copy(other); + nameIsPythiaSyntax_ = other.nameIsPythiaSyntax_; + } + + PartonShowerWeightGroupInfo* PartonShowerWeightGroupInfo::clone() const { + return new PartonShowerWeightGroupInfo(*this); + } + + int PartonShowerWeightGroupInfo::variationIndex(bool isISR, bool isUp, PSVarType variationType) const { + return variationIndex(isISR, isUp, variationType, PSSplittingType::combined); + } + + int PartonShowerWeightGroupInfo::variationIndex(bool isISR, + bool isUp, + PSVarType variationType, + PSSplittingType splittingType) const { + std::string varName = variationName(isISR, isUp, variationType, splittingType); + int wgtIdx = weightIndexFromLabel(varName); + // Guess PS idx if not in label list + if (wgtIdx == -1 && guessPSWeightIdx_) + wgtIdx = psWeightIdxGuess(varName); + return wgtIdx; + } + + std::string PartonShowerWeightGroupInfo::variationName(bool isISR, + bool isUp, + PSVarType variationType, + PSSplittingType splittingType) const { + std::string label = isISR ? "isr" : "fsr"; + + // if ((variationType == PSVarType::con || variationType == PSVarType::def || variationType == PSVarType::red) && + // splittingType != PSSplittingType::combined) + // throw std::invalid_argument("VariationType must be muR or CNS if subprocess is specified"); + + if (nameIsPythiaSyntax_) { + // Splitting + if (splittingType == PSSplittingType::g2gg) + label += ":g2gg"; + else if (splittingType == PSSplittingType::g2qq) + label += ":g2qq"; + else if (splittingType == PSSplittingType::x2xg) + label += ":x2xg"; + else if (splittingType == PSSplittingType::q2qg) + label += ":q2qg"; + // type + if (variationType == PSVarType::con) + label += isUp ? ":murfac=4.0" : ":murfac=0.25"; + else if (variationType == PSVarType::def || variationType == PSVarType::muR) + label += isUp ? ":murfac=2.0" : ":murfac=0.5"; + else if (variationType == PSVarType::red) + label += isUp ? ":murfac=1.414" : ":murfac=0.707"; + else if (variationType == PSVarType::cNS) + label += isUp ? ":cns=2.0" : ":cns=-2.0"; + } else { + // Splitting + if (splittingType == PSSplittingType::g2gg) + label += "_G2GG_"; + else if (splittingType == PSSplittingType::g2qq) + label += "_G2QQ_"; + else if (splittingType == PSSplittingType::x2xg) + label += "_X2XG_"; + else if (splittingType == PSSplittingType::q2qg) + label += "_Q2QG_"; + // type + if (variationType == PSVarType::con) + label += "Con"; + else if (variationType == PSVarType::def) + label += "Def"; + else if (variationType == PSVarType::muR) + label += "muR"; + else if (variationType == PSVarType::red) + label += "Red"; + else if (variationType == PSVarType::cNS) + label += "cNS"; + // Up/Down + if (splittingType != PSSplittingType::combined) { + label += isUp ? "_up" : "_dn"; + } else + label += isUp ? "Hi" : "Lo"; + } + return label; + } + + int PartonShowerWeightGroupInfo::psWeightIdxGuess(const std::string& varName) const { + int wgtIdx; + if (nameIsPythiaSyntax_) { + auto wgtIter = std::find(expectedOrderPythiaSyntax_.begin(), expectedOrderPythiaSyntax_.end(), varName); + wgtIdx = wgtIter - expectedOrderPythiaSyntax_.begin() + 2; + } else { + auto wgtIter = std::find(expectedOrder_.begin(), expectedOrder_.end(), varName); + wgtIdx = wgtIter - expectedOrder_.begin() + 2; + } + if (wgtIdx >= (int)containedIds().size()) + wgtIdx = -1; + return wgtIdx; + } + + void PartonShowerWeightGroupInfo::printVariables() const { + const auto& variations = (nameIsPythiaSyntax_) ? expectedOrderPythiaSyntax_ : expectedOrder_; + for (const auto& varName : variations) { + int wgtIdx = weightIndexFromLabel(varName); + // Guess PS idx if not in label list + if (wgtIdx == -1 && guessPSWeightIdx_) + wgtIdx = psWeightIdxGuess(varName); + if (wgtIdx != -1) + std::cout << varName << " : " << wgtIdx << std::endl; + } + } + +} // namespace gen diff --git a/SimDataFormats/GeneratorProducts/src/PdfWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/PdfWeightGroupInfo.cc new file mode 100644 index 0000000000000..818ad16a78ef2 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/PdfWeightGroupInfo.cc @@ -0,0 +1,29 @@ +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" + +namespace gen { + void PdfWeightGroupInfo::copy(const PdfWeightGroupInfo& other) { + uncertaintyType_ = other.uncertaintyType(); + hasAlphasVars_ = other.hasAlphasVariations(); + alphasUpIndex_ = other.alphasDownIndex(); + alphasDownIndex_ = other.alphasDownIndex(); + parentLhapdfId_ = other.parentLhapdfId(); + WeightGroupInfo::copy(other); + } + + PdfWeightGroupInfo* PdfWeightGroupInfo::clone() const { return new PdfWeightGroupInfo(*this); } + void PdfWeightGroupInfo::addLhaid(int lhaid) { + lhaids_.push_back(lhaid); + if (lhaids_.size() == parentLhapdfSize_) + setIsWellFormed(true); + else + setIsWellFormed(false); + } + + void PdfWeightGroupInfo::setParentLhapdfInfo(int lhaid) { + parentLhapdfId_ = lhaid; + LHAPDF::PDFSet pdfSet(LHAPDF::lookupPDF(lhaid).first); + parentLhapdfSize_ = pdfSet.size(); + parentLhapdfError_ = pdfSet.errorType(); + } + +} // namespace gen diff --git a/SimDataFormats/GeneratorProducts/src/ScaleWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/ScaleWeightGroupInfo.cc new file mode 100644 index 0000000000000..955b4f5c8e260 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/ScaleWeightGroupInfo.cc @@ -0,0 +1,104 @@ +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include +#include + +namespace gen { + void ScaleWeightGroupInfo::copy(const ScaleWeightGroupInfo& other) { + muIndices_ = other.muIndices_; + dynVec_ = other.dynVec_; + dynNames_ = other.dynNames_; + WeightGroupInfo::copy(other); + } + + ScaleWeightGroupInfo* ScaleWeightGroupInfo::clone() const { return new ScaleWeightGroupInfo(*this); } + + void ScaleWeightGroupInfo::addContainedId(int globalIndex, std::string id, std::string label, float muR, float muF) { + int idxDiff = firstId_ - globalIndex; + if (idxDiff > 0) { + for (auto& entry : muIndices_) { + entry += idxDiff; + } + for (auto& subVec : dynVec_) { + for (auto& entry : subVec) { + entry += idxDiff; + } + } + } + WeightGroupInfo::addContainedId(globalIndex, id, label); + setMuRMuFIndex(globalIndex, id, muR, muF); + } + + void ScaleWeightGroupInfo::setMuRMuFIndex(int globalIndex, std::string id, float muR, float muF) { + auto info = weightMetaInfoByGlobalIndex(id, globalIndex); + int index = indexFromMus(muR, muF); + if (!(isValidValue(muR) && isValidValue(muF))) { + setWeightIsCorrupt(); + return; + } + if (index == Central_idx) + containsCentral_ = true; + muIndices_[index] = info.localIndex; + + for (int muidx : muIndices_) { + if (muidx == -1) + return; + } + isWellFormed_ = !weightIsCorrupt_; + } + + void ScaleWeightGroupInfo::setDyn( + int globalIndex, std::string id, float muR, float muF, size_t dynNum, std::string dynName) { + auto info = weightMetaInfoByGlobalIndex(id, globalIndex); + int index = indexFromMus(muR, muF); + // resize if too small + if (dynVec_.at(index).size() < dynNum + 1) { + for (auto& dynIt : dynVec_) + dynIt.resize(dynNum + 1); + dynNames_.resize(dynNum + 1); + } + + if (dynNames_.at(dynNum).empty()) + dynNames_[dynNum] = dynName; + dynVec_[index][dynNum] = info.localIndex; + } + + size_t ScaleWeightGroupInfo::scaleIndex(float muR, float muF, std::string& dynName) const { + auto it = std::find(dynNames_.begin(), dynNames_.end(), dynName); + if (it == dynNames_.end()) + return -1; + else + return scaleIndex(muR, muF, it - dynNames_.begin()); + } + + size_t ScaleWeightGroupInfo::scaleIndex(int index, std::string& dynName) const { + auto it = std::find(dynNames_.begin(), dynNames_.end(), dynName); + if (it == dynNames_.end()) + return -1; + else + return scaleIndex(index, it - dynNames_.begin()); + } + + size_t ScaleWeightGroupInfo::scaleIndex(float muR, float muF, size_t dynNum) const { + ; + if (!(isValidValue(muR) && isValidValue(muF)) || dynNum + 1 > dynNames_.size()) + return -1; + else + return scaleIndex(indexFromMus(muR, muF), dynNum); + } + size_t ScaleWeightGroupInfo::scaleIndex(float muR, float muF) const { + if (!(isValidValue(muR) && isValidValue(muF))) + return -1; + else + return muIndices_.at(indexFromMus(muR, muF)); + } + + std::vector ScaleWeightGroupInfo::dynNames() const { + std::vector returnVec; + for (const auto& item : dynNames_) { + if (!item.empty()) + returnVec.push_back(item); + } + return returnVec; + } + +} // namespace gen diff --git a/SimDataFormats/GeneratorProducts/src/UnknownWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/UnknownWeightGroupInfo.cc new file mode 100644 index 0000000000000..8c33d59eff35b --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/UnknownWeightGroupInfo.cc @@ -0,0 +1,5 @@ +#include "SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h" + +namespace gen { + UnknownWeightGroupInfo* UnknownWeightGroupInfo::clone() const { return new UnknownWeightGroupInfo(*this); } +} // namespace gen diff --git a/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc new file mode 100644 index 0000000000000..90ad7c284aca5 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc @@ -0,0 +1,116 @@ +#include +#include +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include "FWCore/Utilities/interface/Exception.h" + +#include + +namespace gen { + void WeightGroupInfo::copy(const WeightGroupInfo& other) { + isWellFormed_ = other.isWellFormed_; + headerEntry_ = other.headerEntry_; + name_ = other.name_; + description_ = other.description_; + weightType_ = other.weightType_; + idsContained_ = other.idsContained_; + firstId_ = other.firstId_; + lastId_ = other.lastId_; + } + + WeightGroupInfo* WeightGroupInfo::clone() const { + throw cms::Exception("WeightGroupInfo") + << "WeightGroupInfo is abstract, so it's clone() method can't be implemented."; + } + + WeightMetaInfo WeightGroupInfo::weightMetaInfo(int weightEntry) const { return idsContained_.at(weightEntry); } + + WeightMetaInfo WeightGroupInfo::weightMetaInfoByGlobalIndex(std::string wgtId, int weightEntry) const { + if (wgtId.empty()) + wgtId = std::to_string(weightEntry); + int entry = weightVectorEntry(wgtId, weightEntry); + if (entry < 0 || entry >= static_cast(idsContained_.size())) + throw cms::Exception("WeightGroupInfo") + << "Weight entry " << std::to_string(weightEntry) << " is not a member of group " << name_ + << ". \n firstID = " << std::to_string(firstId_) << " lastId = " << std::to_string(lastId_); + return idsContained_.at(entry); + } + + int WeightGroupInfo::weightVectorEntry(std::string& wgtId) const { return weightVectorEntry(wgtId, 0); } + + bool WeightGroupInfo::containsWeight(std::string& wgtId, int weightEntry) const { + return weightVectorEntry(wgtId, weightEntry) != -1; + } + + int WeightGroupInfo::weightVectorEntry(std::string& wgtId, int weightEntry) const { + if (wgtId.empty()) + wgtId = std::to_string(weightEntry); + // First try ordered search + size_t orderedEntry = weightEntry - firstId_; + if (indexInRange(weightEntry) && orderedEntry < idsContained_.size()) { + if (wgtId.empty() || idsContained_.at(orderedEntry).id == wgtId) { + return orderedEntry; + } + } + // Fall back to search on ID + else if (!wgtId.empty()) { + auto it = std::find_if( + idsContained_.begin(), idsContained_.end(), [wgtId](const WeightMetaInfo& w) { return w.id == wgtId; }); + if (it != idsContained_.end()) + return std::distance(idsContained_.begin(), it); + } + return -1; + } + + void WeightGroupInfo::addContainedId(int weightEntry, std::string id, std::string label = "") { + if (id.empty()) + id = std::to_string(weightEntry); + + if (firstId_ == -1 || weightEntry < firstId_) { + firstId_ = weightEntry; + for (auto& entry : idsContained_) // Reset if indices need to be shifted + entry.localIndex++; + } + if (weightEntry > lastId_) + lastId_ = weightEntry; + + size_t localIndex = std::min(weightEntry - firstId_, static_cast(idsContained_.size())); + WeightMetaInfo info = {static_cast(weightEntry), localIndex, id, label}; + // logic to insert for all cases e.g. inserting in the middle of the vector + if (localIndex == idsContained_.size()) + idsContained_.emplace_back(info); + else + idsContained_.insert(idsContained_.begin() + localIndex, info); + } + + std::vector WeightGroupInfo::containedIds() const { return idsContained_; } + + bool WeightGroupInfo::indexInRange(int index) const { return (index <= lastId_ && index >= firstId_); } + + void WeightGroupInfo::cacheWeightIndicesByLabel() { + for (const auto& weight : idsContained_) + weightLabelsToIndices_[weight.label] = weight.localIndex; + } + + std::vector WeightGroupInfo::weightLabels() const { + std::vector labels; + labels.reserve(idsContained_.size()); + for (const auto& weight : idsContained_) + labels.push_back(weight.label); + return labels; + } + + int WeightGroupInfo::weightIndexFromLabel(std::string weightLabel) const { + if (!weightLabelsToIndices_.empty()) { + if (weightLabelsToIndices_.find(weightLabel) != weightLabelsToIndices_.end()) + return static_cast(weightLabelsToIndices_.at(weightLabel)); + return -1; + } + + auto it = std::find_if( + idsContained_.begin(), idsContained_.end(), [weightLabel](const auto& w) { return weightLabel == w.label; }); + if (it == idsContained_.end()) + return -1; + return std::distance(idsContained_.begin(), it); + } + +} // namespace gen diff --git a/SimDataFormats/GeneratorProducts/src/classes.h b/SimDataFormats/GeneratorProducts/src/classes.h index f4ef658588213..4562fa5097487 100644 --- a/SimDataFormats/GeneratorProducts/src/classes.h +++ b/SimDataFormats/GeneratorProducts/src/classes.h @@ -8,6 +8,14 @@ #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEXMLStringProduct.h" #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" diff --git a/SimDataFormats/GeneratorProducts/src/classes_def.xml b/SimDataFormats/GeneratorProducts/src/classes_def.xml index 665a75ab79e8f..5cf7b347c409d 100644 --- a/SimDataFormats/GeneratorProducts/src/classes_def.xml +++ b/SimDataFormats/GeneratorProducts/src/classes_def.xml @@ -199,6 +199,8 @@ + + @@ -227,9 +229,33 @@ + + + + + + + + + + + + + + + + + + + + + + + +