diff --git a/Configuration/Applications/python/ConfigBuilder.py b/Configuration/Applications/python/ConfigBuilder.py index 608651f121385..375ac04f59e7d 100644 --- a/Configuration/Applications/python/ConfigBuilder.py +++ b/Configuration/Applications/python/ConfigBuilder.py @@ -938,6 +938,7 @@ def define_Configs(self): self.RECOSIMDefaultCFF="Configuration/StandardSequences/RecoSim_cff" self.PATDefaultCFF="Configuration/StandardSequences/PAT_cff" self.NANODefaultCFF="PhysicsTools/NanoAOD/nano_cff" + self.NANOGENDefaultCFF="PhysicsTools/NanoAOD/nanogen_cff" self.EIDefaultCFF=None self.SKIMDefaultCFF="Configuration/StandardSequences/Skims_cff" self.POSTRECODefaultCFF="Configuration/StandardSequences/PostRecoGenerator_cff" @@ -987,6 +988,7 @@ def define_Configs(self): self.REPACKDefaultSeq='DigiToRawRepack' self.PATDefaultSeq='miniAOD' self.PATGENDefaultSeq='miniGEN' + self.NANOGENDefaultSeq='nanogenSequence' self.NANODefaultSeq='nanoSequence' self.EVTCONTDefaultCFF="Configuration/EventContent/EventContent_cff" @@ -1690,6 +1692,15 @@ def prepare_NANO(self, sequence = "nanoAOD"): self._options.customise_commands = self._options.customise_commands + " \n" self._options.customise_commands = self._options.customise_commands + "process.unpackedPatTrigger.triggerResults= cms.InputTag( 'TriggerResults::"+self._options.hltProcess+"' )\n" + def prepare_NANOGEN(self, sequence = "nanoAOD"): + ''' Enrich the schedule with NANO ''' + self.loadDefaultOrSpecifiedCFF(sequence,self.NANOGENDefaultCFF) + self.scheduleSequence(sequence.split('.')[-1],'nanoAOD_step') + custom = "customizeNanoGEN" + if self._options.runUnscheduled: + self._options.customisation_file_unsch.insert(0, '.'.join([self.NANOGENDefaultCFF, custom])) + else: + self._options.customisation_file.insert(0, '.'.join([self.NANOGENDefaultCFF, custom])) def prepare_EI(self, sequence = None): ''' Enrich the schedule with event interpretation ''' diff --git a/DataFormats/NanoAOD/src/classes_def.xml b/DataFormats/NanoAOD/src/classes_def.xml index 32c077349ff95..24d6da6064c73 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 4257af0e29ebe..07d1b0adc6194 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..4836ac90b5332 --- /dev/null +++ b/GeneratorInterface/Core/interface/GenWeightHelper.h @@ -0,0 +1,30 @@ +#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); + private: + }; +} + +#endif + + diff --git a/GeneratorInterface/Core/interface/LHEWeightHelper.h b/GeneratorInterface/Core/interface/LHEWeightHelper.h new file mode 100644 index 0000000000000..12f973ca45fcc --- /dev/null +++ b/GeneratorInterface/Core/interface/LHEWeightHelper.h @@ -0,0 +1,33 @@ +#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() {}; + void setHeaderLines(std::vector headerLines); + void parseWeights(); + void buildGroups(); + std::unique_ptr buildGroup(const ParsedWeight& weight); + private: + std::vector headerLines_; + }; +} + +#endif + diff --git a/GeneratorInterface/Core/interface/WeightHelper.h b/GeneratorInterface/Core/interface/WeightHelper.h new file mode 100644 index 0000000000000..66cfbfd547155 --- /dev/null +++ b/GeneratorInterface/Core/interface/WeightHelper.h @@ -0,0 +1,65 @@ +#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/ScaleWeightGroupInfo.h" +#include "SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h" +#include + +namespace gen { + struct PdfSetInfo { + std::string name; + int lhapdfId; + PdfUncertaintyType uncertaintyType; + }; + + struct ParsedWeight { + std::string id; + size_t index; + std::string groupname; + std::string content; + std::unordered_map attributes; + }; + + class WeightHelper { + public: + WeightHelper(); + edm::OwnVector weightGroups() { + return weightGroups_; + } + std::unique_ptr weightProduct(std::vector, float w0); + std::unique_ptr weightProduct(std::vector, float w0); + void setModel(std::string model) { model_ = model; } + int addWeightToProduct(std::unique_ptr& product, double weight, std::string name, int weightNum, int groupIndex); + int findContainingWeightGroup(std::string wgtId, int weightIndex, int previousGroupIndex); + protected: + std::string model_; + std::vector parsedWeights_; + const std::vector pdfSetsInfo; + std::map currWeightAttributeMap_; + std::map currGroupAttributeMap_; + edm::OwnVector weightGroups_; + bool isScaleWeightGroup(const ParsedWeight& weight); + bool isMEParamWeightGroup(const ParsedWeight& weight); + bool isPdfWeightGroup(const ParsedWeight& weight); + void updateScaleInfo(const ParsedWeight& weight); + void updatePdfInfo(const ParsedWeight& weight); + std::string searchAttributes(const std::string& label, const ParsedWeight& weight); + + // Possible names for the same thing + const std::unordered_map> attributeNames_ = { + {"muf", {"muR", "MUR", "muf","facscfact"}}, + {"mur", {"muF", "MUF", "mur","renscfact"}}, + {"pdf", {"PDF", "PDF set", "lhapdf", "pdf", "pdf set", "pdfset"}} + }; + }; +} + +#endif + diff --git a/GeneratorInterface/Core/plugins/BuildFile.xml b/GeneratorInterface/Core/plugins/BuildFile.xml index af12a8119d35c..83ce8d6a7c13d 100644 --- a/GeneratorInterface/Core/plugins/BuildFile.xml +++ b/GeneratorInterface/Core/plugins/BuildFile.xml @@ -6,7 +6,12 @@ + + + + + + + - - diff --git a/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc b/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc new file mode 100644 index 0000000000000..98075bedd6a3b --- /dev/null +++ b/GeneratorInterface/Core/plugins/GenWeightProductProducer.cc @@ -0,0 +1,99 @@ +#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"))) +{ + produces(); + produces(); +} + + +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().size() > 0 ? 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); + } + + if (weightNames_.size() == 0) { + edm::Handle genLumiInfoHandle; + iLumi.getByToken(genLumiInfoToken_, genLumiInfoHandle); + + weightNames_ = genLumiInfoHandle->weightNames(); + weightHelper_.parseWeightGroupsFromNames(weightNames_); + } + auto weightInfoProduct = std::make_unique(); + for (auto& weightGroup : weightHelper_.weightGroups()) { + weightInfoProduct->addWeightGroupInfo(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..f259540d79b59 --- /dev/null +++ b/GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc @@ -0,0 +1,114 @@ +#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" + +class LHEWeightProductProducer : public edm::one::EDProducer { +public: + explicit LHEWeightProductProducer(const edm::ParameterSet& iConfig); + ~LHEWeightProductProducer() override; + +private: + gen::LHEWeightHelper weightHelper_; + std::string lheLabel_; + edm::EDGetTokenT lheRunInfoToken_; + edm::EDGetTokenT lheEventToken_; + const edm::EDGetTokenT lheWeightInfoToken_; + bool foundWeightProduct_; + + 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) : + lheLabel_(iConfig.getParameter("lheSourceLabel")), + lheRunInfoToken_(consumes(lheLabel_)), + lheEventToken_(consumes(lheLabel_)), + lheWeightInfoToken_(consumes(lheLabel_)) +{ + produces(); + produces(); +} + + +LHEWeightProductProducer::~LHEWeightProductProducer() +{ +} + + +// ------------ method called to produce the data ------------ +void +LHEWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + if (foundWeightProduct_) + return; + + edm::Handle lheEventInfo; + iEvent.getByToken(lheEventToken_, lheEventInfo); + // Read weights from LHEEventProduct + 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; + run.getByLabel(lheLabel_, lheRunInfoHandle); + + 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; + } + + //weightHelper_.parseWeightGroupsFromHeader(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; + lumi.getByToken(lheWeightInfoToken_, lheWeightInfoHandle); + if (lheWeightInfoHandle.isValid()) { + foundWeightProduct_ = true; + return; + } + + auto weightInfoProduct = std::make_unique(); + for (auto& weightGroup : weightHelper_.weightGroups()) { + weightInfoProduct->addWeightGroupInfo(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..74bad0d66781a --- /dev/null +++ b/GeneratorInterface/Core/src/GenWeightHelper.cc @@ -0,0 +1,35 @@ +#include "GeneratorInterface/Core/interface/GenWeightHelper.h" +#include + +using namespace tinyxml2; + +namespace gen { + GenWeightHelper::GenWeightHelper() { + } + + void + GenWeightHelper::parseWeightGroupsFromNames(std::vector weightNames) { + int index = 0; + + 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(",")); + weightGroups_.push_back(new gen::UnknownWeightGroupInfo("inGen")); + } + // Working on the not-so-nice assumption that all non-LHE gen weights are PS weights + else if (weightGroups_.size() == 0) { + weightGroups_.push_back(new gen::PartonShowerWeightGroupInfo("shower")); + } + auto& group = weightGroups_.back(); + // No IDs for Gen weights + group.addContainedId(index++, "", weightName); + } + } +} + + diff --git a/GeneratorInterface/Core/src/LHEWeightHelper.cc b/GeneratorInterface/Core/src/LHEWeightHelper.cc new file mode 100644 index 0000000000000..7a882e07ea84e --- /dev/null +++ b/GeneratorInterface/Core/src/LHEWeightHelper.cc @@ -0,0 +1,78 @@ +#include "GeneratorInterface/Core/interface/LHEWeightHelper.h" +#include +#include + +using namespace tinyxml2; + +namespace gen { + void LHEWeightHelper::setHeaderLines(std::vector headerLines) { + headerLines_ = headerLines; + } + + void LHEWeightHelper::parseWeights() { + tinyxml2::XMLDocument xmlDoc; + xmlDoc.Parse(("" + boost::algorithm::join(headerLines_, "") + "").c_str()); + tinyxml2::XMLElement* root = xmlDoc.FirstChildElement("root"); + + size_t weightIndex = 0; + for (auto* e = root->FirstChildElement(); e != nullptr; e = e->NextSiblingElement()) { + std::string groupName = ""; + if (strcmp(e->Name(), "weight") == 0) { + // we are here if there is a weight that does not belong to any group + std::string text = ""; + if (e->GetText()) + text = e->GetText(); + parsedWeights_.push_back({e->Attribute("id"), weightIndex++, groupName, text}); + } + if (strcmp(e->Name(), "weightgroup") == 0) { + groupName = e->Attribute("name"); + for (auto* inner = e->FirstChildElement("weight"); inner != nullptr; + inner = inner->NextSiblingElement("weight")) { + // we are here if there is a weight in a weightgroup + std::string text = ""; + if (inner->GetText()) + text = inner->GetText(); + std::unordered_map attributes; + for (auto* att = inner->FirstAttribute(); att != nullptr; att = att->Next()) + attributes[att->Name()] = att->Value(); + parsedWeights_.push_back({inner->Attribute("id"), weightIndex++, groupName, text, attributes}); + } + } + } + buildGroups(); + } + + //void LHEWeightHelper::setGroupPriority() { + // groupCheckOrder_ = {gen::WeightType::kScaleWeights, gen::WeightType::kPdfWeights, + // gen::WeightType::kMEParamWeights + // }; + //} + + void LHEWeightHelper::buildGroups() { + std::string currentGroupName; + for (const auto& weight : parsedWeights_) { + if (weight.groupname != currentGroupName) { + weightGroups_.push_back(*buildGroup(weight)); + } + currentGroupName = weight.groupname; + WeightGroupInfo& group = weightGroups_.back(); + group.addContainedId(weight.index, weight.id, weight.content); + + if (group.weightType() == gen::WeightType::kScaleWeights) + updateScaleInfo(weight); + else if (group.weightType() == gen::WeightType::kPdfWeights) + updatePdfInfo(weight); + } + } + + std::unique_ptr LHEWeightHelper::buildGroup(const ParsedWeight& weight) { + 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); + + return std::make_unique(weight.groupname); + } +} diff --git a/GeneratorInterface/Core/src/WeightHelper.cc b/GeneratorInterface/Core/src/WeightHelper.cc new file mode 100644 index 0000000000000..b178134a7b156 --- /dev/null +++ b/GeneratorInterface/Core/src/WeightHelper.cc @@ -0,0 +1,191 @@ +#include "GeneratorInterface/Core/interface/WeightHelper.h" + +namespace gen { + WeightHelper::WeightHelper() : + pdfSetsInfo({ + // In principle this can be parsed from $LHAPDF_DATA_PATH/pdfsets.index, + // but do we really want to do that? Can also just hardcode a subset... + // TODO: Actually we can just take this from LHAPDF + {"NNPDF31_nnlo_hessian_pdfas", 306000, kHessianUnc}, + {"NNPDF31_nnlo_as_0118", 303600, kMonteCarloUnc}, + {"NNPDF31_nlo_as_0118", 303400, kMonteCarloUnc}, + {"NNPDF31_nlo_hessian_pdfas", 305800, kHessianUnc}, + {"NNPDF31_nnlo_as_0108", 322500, kVariationSet}, + {"NNPDF31_nnlo_as_0110", 322700, kVariationSet}, + {"NNPDF31_nnlo_as_0112", 322900, kVariationSet}, + {"NNPDF31_nnlo_as_0114", 323100, kVariationSet}, + {"NNPDF31_nnlo_as_0117", 323300, kVariationSet}, + {"NNPDF31_nnlo_as_0119", 323500, kVariationSet}, + {"NNPDF31_nnlo_as_0122", 323700, kVariationSet}, + {"NNPDF31_nnlo_as_0124", 323900, kVariationSet}, + {"NNPDF31_nnlo_as_0118_nf_4_mc_hessian", 325500, kHessianUnc}, + {"NNPDF31_nlo_as_0118_nf_4", 320500, kMonteCarloUnc}, + {"NNPDF31_nnlo_as_0118_nf_4", 320900, kMonteCarloUnc}, + {"NNPDF30_nlo_nf_5_pdfas", 292200, kMonteCarloUnc}, + {"NNPDF30_nnlo_nf_5_pdfas", 292600, kMonteCarloUnc}, + {"NNPDF30_nnlo_nf_4_pdfas", 292400, kMonteCarloUnc}, + {"NNPDF30_nlo_nf_4_pdfas", 292000, kMonteCarloUnc}, + {"NNPDF30_lo_as_0130", 263000, kMonteCarloUnc}, + {"NNPDF30_lo_as_0118", 262000, kMonteCarloUnc}, + {"CT14nnlo", 13000, kHessianUnc}, + {"CT14nlo", 13100, kHessianUnc}, + {"CT14nnlo_as_0116", 13065, kVariationSet}, + {"CT14nnlo_as_0120", 13069, kVariationSet}, + {"CT14nlo_as_0116", 13163, kVariationSet}, + {"CT14nlo_as_0120", 13167, kVariationSet}, + {"CT14lo", 13200, kVariationSet}, + {"MMHT2014nlo68clas118", 25200, kHessianUnc}, + {"MMHT2014nnlo68cl", 25300, kHessianUnc}, + {"MMHT2014lo68cl", 25000, kHessianUnc}, + {"PDF4LHC15_nlo_100_pdfas", 90200, kMonteCarloUnc}, + {"PDF4LHC15_nnlo_100_pdfas", 91200, kMonteCarloUnc}, + {"PDF4LHC15_nlo_30_pdfas", 90400, kMonteCarloUnc}, + {"PDF4LHC15_nnlo_30_pdfas", 91400, kMonteCarloUnc}, + {"ABMP16als118_5_nnlo", 42780, kHessianUnc}, + {"HERAPDF20_NLO_EIG", 61130, kHessianUnc}, + {"HERAPDF20_NNLO_EIG", 61200, kHessianUnc}, + {"HERAPDF20_NLO_VAR", 61130, kHessianUnc}, + {"HERAPDF20_NNLO_VAR", 61230, kHessianUnc}, + {"CT14qed_inc_proton", 13400, kHessianUnc}, + {"LUXqed17_plus_PDF4LHC15_nnlo_100", 82200, kMonteCarloUnc}, + }) + { 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 std::find_if(pdfSetsInfo.begin(), pdfSetsInfo.end(), + [name] (const PdfSetInfo& setInfo) { return setInfo.name == name; }) != pdfSetsInfo.end(); + } + + bool WeightHelper::isMEParamWeightGroup(const ParsedWeight& weight) { + return (weight.groupname.find("mg_reweighting") != std::string::npos); + } + + std::string WeightHelper::searchAttributes(const std::string& label, const ParsedWeight& weight) { + for (const auto& lab : attributeNames_.at(label)) { + auto& attributes = weight.attributes; + if (attributes.find(lab) != attributes.end()) { + return boost::algorithm::trim_copy_if(attributes.at(lab), boost::is_any_of("\"")); + } + } + // Next regexes + return ""; + } + + void WeightHelper::updateScaleInfo(const ParsedWeight& weight) { + auto& group = weightGroups_.back(); + auto& scaleGroup = dynamic_cast(group); + std::string muRText = searchAttributes("mur", weight); + std::string muFText = searchAttributes("mur", weight); + if (muRText.empty() || muFText.empty()) { + scaleGroup.setIsWellFormed(false); + return; + } + + try { + float muR = std::stof(muRText); + float muF = std::stof(muFText); + scaleGroup.setMuRMuFIndex(weight.index, weight.id, muR, muF); + } + catch(std::invalid_argument& e) { + scaleGroup.setIsWellFormed(false); + } + } + + void WeightHelper::updatePdfInfo(const ParsedWeight& weight) { + auto& pdfGroup = dynamic_cast(weightGroups_.back()); + std::string lhaidText = searchAttributes("pdf", weight); + int lhaid = 0; + if (!lhaidText.empty()) { + try { + lhaid = std::stoi(lhaidText); + } + catch(std::invalid_argument& e) { + pdfGroup.setIsWellFormed(false); + } + + if (!pdfGroup.containsParentLhapdfId(lhaid, weight.index)) { + std::string description; + auto pdfInfo = std::find_if(pdfSetsInfo.begin(), pdfSetsInfo.end(), + [lhaid] (const PdfSetInfo& setInfo) { return setInfo.lhapdfId == lhaid; }); + if (pdfInfo != pdfSetsInfo.end()) { + pdfGroup.setUncertaintyType(pdfInfo->uncertaintyType); + if (pdfInfo->uncertaintyType == gen::kHessianUnc) + description = "Hessian "; + else if (pdfInfo->uncertaintyType == gen::kMonteCarloUnc) + description = "Monte Carlo "; + description += "Uncertainty sets for LHAPDF set " + pdfInfo->name; + description += " with LHAID = " + std::to_string(lhaid); + } + description = "Uncertainty sets for LHAPDF set with LHAID = " + std::to_string(lhaid); + pdfGroup.addLhapdfId(lhaid, weight.index); + pdfGroup.setDescription(description); + } + } + else + pdfGroup.setIsWellFormed(false); + } + + // TODO: Could probably recycle this code better + std::unique_ptr WeightHelper::weightProduct(std::vector weights, float w0) { + auto weightProduct = std::make_unique(w0); + weightProduct->setNumWeightSets(weightGroups_.size()); + int weightGroupIndex = 0; + for (unsigned int i = 0; i < weights.size(); i++) { + addWeightToProduct(weightProduct, weights.at(i), "", i, weightGroupIndex); + } + return std::move(weightProduct); + } + + 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; + for (const auto& weight : weights) { + weightGroupIndex = addWeightToProduct(weightProduct, weight.wgt, weight.id, i++, weightGroupIndex); + } + return std::move(weightProduct); + } + + int WeightHelper::addWeightToProduct(std::unique_ptr& product, + double weight, std::string name, int weightNum, int groupIndex) { + groupIndex = findContainingWeightGroup(name, weightNum, groupIndex); + auto group = weightGroups_[groupIndex]; + int entry = group.weightVectorEntry(name, weightNum); + 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 (auto weightGroup : weightGroups_) { + if (weightGroup.containsWeight(wgtId, weightIndex)) + return counter; + counter++; + } + // Needs to be properly handled + throw std::range_error("Unmatched Generator weight! ID was " + wgtId + " index was " + std::to_string(weightIndex) + + "\nNot found in any of " + std::to_string(weightGroups_.size()) + " weightGroups."); + } +} + diff --git a/GeneratorInterface/Core/test/dumpWeightInfo.py b/GeneratorInterface/Core/test/dumpWeightInfo.py new file mode 100644 index 0000000000000..2ee498c9fb383 --- /dev/null +++ b/GeneratorInterface/Core/test/dumpWeightInfo.py @@ -0,0 +1,45 @@ +from DataFormats.FWLite import Events,Handle,Runs,Lumis +import ROOT +sources = ["externalLHEProducer"] +#sources = ["testLHEWeights", "testGenWeights"] + +#for filename in ["HIG-RunIIFall18wmLHEGS-00509.root"," HIG-RunIIFall18wmLHEGS-00509_ordered.root","HIG-RunIIFall18wmLHEGS-00509_unordered.root"]: +for filename in ["HIG-RunIIFall18wmLHEGS-00509.root"]: +#for filename in ["test.root"]: + for source in sources: + lumis = Lumis(filename) + lumi = lumis.__iter__().next() + weightInfoHandle = Handle("GenWeightInfoProduct") + lumi.getByLabel(source, weightInfoHandle) + weightInfoProd = weightInfoHandle.product() + + events = Events(filename) + event = events.__iter__().next() + weightHandle = Handle("GenWeightProduct") + event.getByLabel(source, weightHandle) + weightInfo = weightHandle.product() + print weightInfoProd.allWeightGroupsInfo(), len(weightInfoProd.allWeightGroupsInfo()) + 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 "contains LHAPDFIds", len(matching.lhapdfIdsContained()), [i for i in matching.lhapdfIdsContained()], + 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..0c0925dc2e46e --- /dev/null +++ b/GeneratorInterface/Core/test/testGenWeightProducer_cfg.py @@ -0,0 +1,41 @@ +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")) + +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 e585d6c7a8f48..96e1834824f70 100644 --- a/GeneratorInterface/LHEInterface/plugins/BuildFile.xml +++ b/GeneratorInterface/LHEInterface/plugins/BuildFile.xml @@ -17,9 +17,3 @@ - - - - - - diff --git a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc index e1d6be38f8359..eb2ffa677c3fa 100644 --- a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc +++ b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc @@ -47,12 +47,15 @@ Description: [one line class summary] #include "SimDataFormats/GeneratorProducts/interface/LesHouches.h" #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEXMLStringProduct.h" #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h" #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h" #include "GeneratorInterface/LHEInterface/interface/LHEReader.h" +#include "GeneratorInterface/Core/interface/LHEWeightHelper.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/Utilities/interface/RandomNumberGenerator.h" @@ -63,7 +66,9 @@ 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); ~ExternalLHEProducer() override; @@ -74,10 +79,12 @@ class ExternalLHEProducer : public edm::one::EDProducer readOutput(); void nextEvent(); @@ -95,13 +102,14 @@ class ExternalLHEProducer : public edm::one::EDProducer> nPartonMapping_{}; - std::unique_ptr reader_; - std::shared_ptr runInfoLast; - std::shared_ptr runInfo; - std::shared_ptr partonLevel; - boost::ptr_deque runInfoProducts; - bool wasMerged; - + std::unique_ptr reader_; + gen::LHEWeightHelper weightHelper_; + std::shared_ptr runInfoLast; + std::shared_ptr runInfo; + std::shared_ptr partonLevel; + boost::ptr_deque runInfoProducts; + bool wasMerged; + class FileCloseSentry : private boost::noncopyable { public: explicit FileCloseSentry(int fd) : fd_(fd){}; @@ -151,8 +159,10 @@ ExternalLHEProducer::ExternalLHEProducer(const edm::ParameterSet& iConfig) produces("LHEScriptOutput"); produces(); + produces(); produces(); produces(); + produces(); } ExternalLHEProducer::~ExternalLHEProducer() {} @@ -181,7 +191,14 @@ void ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe } std::for_each(partonLevel->weights().begin(), partonLevel->weights().end(), - boost::bind(&LHEEventProduct::addWeight, product.get(), _1)); + boost::bind(&LHEEventProduct::addWeight, + product.get(), _1)); + + // Should also zero out the weights in the GenInfoProduct + auto weightProduct = weightHelper_.weightProduct(partonLevel->weights(), + partonLevel->originalXWGTUP()); + iEvent.put(std::move(weightProduct)); + product->setScales(partonLevel->scales()); if (nPartonMapping_.empty()) { product->setNpLO(partonLevel->npLO()); @@ -246,31 +263,31 @@ void ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe } // ------------ method called when starting to processes a run ------------ -void ExternalLHEProducer::beginRunProduce(edm::Run& run, edm::EventSetup const& es) { - // pass the number of events as previous to last argument +void +ExternalLHEProducer::beginRunProduce(edm::Run& run, edm::EventSetup const& es) +{ - std::ostringstream eventStream; - eventStream << nEvents_; - // args_.push_back(eventStream.str()); - args_.insert(args_.begin() + 1, eventStream.str()); + // pass the number of events as previous to last argument + std::ostringstream eventStream; + eventStream << nEvents_; + // args_.push_back(eventStream.str()); + args_.insert(args_.begin() + 1, eventStream.str()); - // pass the random number generator seed as last argument + // pass the random number generator seed as last argument - edm::Service rng; + edm::Service rng; - if (!rng.isAvailable()) { - throw cms::Exception("Configuration") - << "The ExternalLHEProducer module requires the RandomNumberGeneratorService\n" - "which is not present in the configuration file. You must add the service\n" - "in the configuration file if you want to run ExternalLHEProducer"; - } - std::ostringstream randomStream; - randomStream << rng->mySeed(); - // args_.push_back(randomStream.str()); - args_.insert(args_.begin() + 2, randomStream.str()); + if ( ! rng.isAvailable()) { + throw cms::Exception("Configuration") + << "The ExternalLHEProducer module requires the RandomNumberGeneratorService\n" + "which is not present in the configuration file. You must add the service\n" + "in the configuration file if you want to run ExternalLHEProducer"; + } + std::ostringstream randomStream; + randomStream << rng->mySeed(); + args_.insert(args_.begin() + 2, randomStream.str()); - // args_.emplace_back(std::to_string(nThreads_)); - args_.insert(args_.begin() + 3, std::to_string(nThreads_)); + args_.insert(args_.begin() + 3, std::to_string(nThreads_)); for (unsigned int iArg = 0; iArg < args_.size(); iArg++) { LogDebug("LHEInputArgs") << "arg [" << iArg << "] = " << args_[iArg]; @@ -293,14 +310,14 @@ void ExternalLHEProducer::beginRunProduce(edm::Run& run, edm::EventSetup const& p->fillCompressedContent(instream, 0.25 * insize); instream.close(); } - run.put(std::move(p), "LHEScriptOutput"); + run.put(std::move(p), "LHEScriptOutput"); - // LHE C++ classes translation - // (read back uncompressed file from disk in streaming mode again to save memory) + // LHE C++ classes translation + // (read back uncompressed file from disk in streaming mode again to save memory) - std::vector infiles(1, outputFile_); - unsigned int skip = 0; - reader_ = std::make_unique(infiles, skip); + std::vector infiles(1, outputFile_); + unsigned int skip = 0; + reader_ = std::make_unique(infiles, skip); nextEvent(); if (runInfoLast) { @@ -322,8 +339,36 @@ void ExternalLHEProducer::beginRunProduce(edm::Run& run, edm::EventSetup const& runInfo.reset(); } + + nextEvent(); + if (runInfoLast) { + runInfo = runInfoLast; + + std::unique_ptr product(new LHERunInfoProduct(*runInfo->getHEPRUP())); + std::for_each(runInfo->getHeaders().begin(), + runInfo->getHeaders().end(), + boost::bind(&LHERunInfoProduct::addHeader, + product.get(), _1)); + std::for_each(runInfo->getComments().begin(), + runInfo->getComments().end(), + boost::bind(&LHERunInfoProduct::addComment, + product.get(), _1)); + + // keep a copy around in case of merging + runInfoProducts.push_back(new LHERunInfoProduct(*product)); + wasMerged = false; + + run.put(std::move(product)); + + //weightHelper_.parseWeightGroupsFromHeader(runInfo->findHeader("initrwgt")); + weightHelper_.setHeaderLines(runInfo->findHeader("initrwgt")); + weightHelper_.parseWeights(); + + runInfo.reset(); + } } + // ------------ method called when ending the processing of a run ------------ void ExternalLHEProducer::endRunProduce(edm::Run& run, edm::EventSetup const& es) { if (!runInfoProducts.empty()) { @@ -333,18 +378,25 @@ void ExternalLHEProducer::endRunProduce(edm::Run& run, edm::EventSetup const& es nextEvent(); if (partonLevel) { - throw edm::Exception(edm::errors::EventGenerationFailure) - << "Error in ExternalLHEProducer::endRunProduce(). " - << "Event loop is over, but there are still lhe events to process." - << "This could happen if lhe file contains more events than requested. This is never expected to happen."; - } - - reader_.reset(); + throw edm::Exception(edm::errors::EventGenerationFailure) << "Error in ExternalLHEProducer::endRunProduce(). " + << "Event loop is over, but there are still lhe events to process." + << "This could happen if lhe file contains more events than requested. This is never expected to happen."; + } + + reader_.reset(); if (unlink(outputFile_.c_str())) { - throw cms::Exception("OutputDeleteError") << "Unable to delete original script output file " << outputFile_ - << " (errno=" << errno << ", " << strerror(errno) << ")."; + throw cms::Exception("OutputDeleteError") << "Unable to delete original script output file " << outputFile_ << " (errno=" << errno << ", " << strerror(errno) << ")."; + } +} + +void +ExternalLHEProducer::beginLuminosityBlockProduce(edm::LuminosityBlock& lumi, edm::EventSetup const& es) { + auto weightInfoProduct = std::make_unique(); + for (auto& weightGroup : weightHelper_.weightGroups()) { + weightInfoProduct->addWeightGroupInfo(weightGroup.clone()); } + lumi.put(std::move(weightInfoProduct)); } // ------------ Close all the open file descriptors ------------ @@ -473,6 +525,7 @@ void ExternalLHEProducer::executeScript() { } } + // ------------ Read the output script ------------ #define BUFSIZE 4096 std::unique_ptr ExternalLHEProducer::readOutput() { diff --git a/GeneratorInterface/LHEInterface/plugins/LHESource.cc b/GeneratorInterface/LHEInterface/plugins/LHESource.cc index bc8f77ce549ea..a29960d0a0a1b 100644 --- a/GeneratorInterface/LHEInterface/plugins/LHESource.cc +++ b/GeneratorInterface/LHEInterface/plugins/LHESource.cc @@ -23,6 +23,7 @@ #include "SimDataFormats/GeneratorProducts/interface/LesHouches.h" #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h" #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h" @@ -105,6 +106,7 @@ void LHESource::readRun_(edm::RunPrincipal& runPrincipal) { runPrincipal.fillRunPrincipal(processHistoryRegistryForUpdate()); putRunInfoProduct(runPrincipal); + putWeightInfoProduct(runPrincipal); } void LHESource::readLuminosityBlock_(edm::LuminosityBlockPrincipal& lumiPrincipal) { @@ -120,7 +122,28 @@ void LHESource::putRunInfoProduct(edm::RunPrincipal& iRunPrincipal) { } } -bool LHESource::setRunAndEventInfo(edm::EventID&, edm::TimeValue_t&, edm::EventAuxiliary::ExperimentType&) { +void LHESource::putWeightInfoProduct(edm::RunPrincipal& iRunPrincipal) { + if (runInfoProductLast_) { + auto product = std::make_unique(); + gen::WeightGroupInfo scaleInfo( + "" + ); + scaleInfo.setWeightType(gen::WeightType::kScaleWeights); + + gen::WeightGroupInfo cenPdfInfo( + "" + ); + cenPdfInfo.setWeightType(gen::WeightType::kPdfWeights); + + product->addWeightGroupInfo(&scaleInfo); + product->addWeightGroupInfo(&cenPdfInfo); + std::unique_ptr rdp(new edm::Wrapper(std::move(product))); + //iRunPrincipal.put(lheProvenanceHelper_.weightProductBranchDescription_, std::move(rdp)); + } +} + +bool LHESource::setRunAndEventInfo(edm::EventID&, edm::TimeValue_t&, edm::EventAuxiliary::ExperimentType&) +{ nextEvent(); if (!partonLevel_) { // We just finished an input file. See if there is another. diff --git a/GeneratorInterface/LHEInterface/plugins/LHESource.h b/GeneratorInterface/LHEInterface/plugins/LHESource.h index 05208fdf03a6c..46bad67ce4f4c 100644 --- a/GeneratorInterface/LHEInterface/plugins/LHESource.h +++ b/GeneratorInterface/LHEInterface/plugins/LHESource.h @@ -50,7 +50,8 @@ class LHESource : public edm::ProducerSourceFromFiles { void nextEvent(); void putRunInfoProduct(edm::RunPrincipal&); - void fillRunInfoProduct(lhef::LHERunInfo const&, LHERunInfoProduct&); + void putWeightInfoProduct(edm::RunPrincipal&); + void fillRunInfoProduct(lhef::LHERunInfo const&, LHERunInfoProduct& ); std::unique_ptr reader_; diff --git a/GeneratorInterface/LHEInterface/test/test_Weights_cfg.py b/GeneratorInterface/LHEInterface/test/test_Weights_cfg.py new file mode 100644 index 0000000000000..b70c90b3ca2d0 --- /dev/null +++ b/GeneratorInterface/LHEInterface/test/test_Weights_cfg.py @@ -0,0 +1,85 @@ +# Auto generated configuration file +# using: +# Revision: 1.19 +# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v +# with command line options: Configuration/GenProduction/python/HIG-RunIIFall18wmLHEGS-00509-fragment.py --fileout file:HIG-RunIIFall18wmLHEGS-00509.root --mc --eventcontent LHE --datatier LHE --conditions 102X_upgrade2018_realistic_v11 --step LHE --python_filename HIG-RunIIFall18wmLHEGS-00509_1_cfg.py --no_exec +import FWCore.ParameterSet.Config as cms + + + +process = cms.Process('LHE') + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10) +) + +# Input source +process.source = cms.Source("EmptySource") + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('Configuration/GenProduction/python/HIG-RunIIFall18wmLHEGS-00509-fragment.py nevts:10'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition +process.LHEoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('LHE'), + filterName = cms.untracked.string('') + ), + fileName = cms.untracked.string('file:HIG-RunIIFall18wmLHEGS-00509.root'), + #outputCommands = process.LHEEventContent.outputCommands, + outputCommands = cms.untracked.vstring('keep *', + 'drop ME*_MEtoEDM*_*_*'), + splitLevel = cms.untracked.int32(0) +) + +# Additional output definition + +process.externalLHEProducer = cms.EDProducer("ExternalLHEProducer", + #args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/WLLJJ_WToLNu_EWK_4F_MLL-60_slc6_amd64_gcc481_CMSSW_7_1_30_tarball_Dummy.tgz'), + args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/VVV_aQGCfs_dummy.tgz'), + nEvents = cms.untracked.uint32(10), + numberOfParameters = cms.uint32(1), + outputFile = cms.string('cmsgrid_final.lhe'), + scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh') +) + +# args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/DarkMatter_MonoZPrime_V_Mx50_Mv500_gDMgQ1_LO_slc6_amd64_gcc481_CMSSW_7_1_30_tarball_Dummy.tgz'), +# args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/WLLJJ_WToLNu_EWK_4F_MLL-60_slc6_amd64_gcc481_CMSSW_7_1_30_tarball_Dummy.tgz'), +# args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/ZZ_4L_NNPDF30_13TeV_tarballDummy.tar.gz'), +# args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/ZZ_slc6_amd64_gcc630_CMSSW_9_3_0_ZZTo4L2017_pdf306000Dummy.tgz'), + + + +# Path and EndPath definitions +process.lhe_step = cms.Path(process.externalLHEProducer) +process.endjob_step = cms.EndPath(process.endOfProcess) +process.LHEoutput_step = cms.EndPath(process.LHEoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.lhe_step,process.endjob_step,process.LHEoutput_step) +from PhysicsTools.PatAlgos.tools.helpers import associatePatAlgosToolsTask +associatePatAlgosToolsTask(process) + + + + +# Customisation from command line + +# Add early deletion of temporary data products to reduce peak memory need +from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete +process = customiseEarlyDelete(process) +# End adding early deletion diff --git a/PhysicsTools/NanoAOD/BuildFile.xml b/PhysicsTools/NanoAOD/BuildFile.xml index c36ed4184624a..1d6f1f6d5dccb 100644 --- a/PhysicsTools/NanoAOD/BuildFile.xml +++ b/PhysicsTools/NanoAOD/BuildFile.xml @@ -5,6 +5,7 @@ + diff --git a/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc b/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc index 77c8989defef1..6b678099e98b8 100644 --- a/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/GenWeightsTableProducer.cc @@ -20,39 +20,42 @@ #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 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; - sumPDF.clear(); - sumScale.clear(); - sumRwgt.clear(); - sumNamed.clear(), sumPS.clear(); + num_ = 0; + sumw_ = 0; + sumw2_ = 0; + sumPDF_.clear(); + sumScale_.clear(); + sumRwgt_.clear(); + sumNamed_.clear(); + sumPS_.clear(); } // inc the counters void incGenOnly(double w) { - num++; - sumw += w; - sumw2 += (w * w); + num_++; + sumw_ += w; + sumw2_ += (w * w); } void incPSOnly(double w0, const std::vector& wPS) { if (!wPS.empty()) { - if (sumPS.empty()) - sumPS.resize(wPS.size(), 0); + if (sumPS_.empty()) + sumPS_.resize(wPS.size(), 0); for (unsigned int i = 0, n = wPS.size(); i < n; ++i) - sumPS[i] += (w0 * wPS[i]); + sumPS_[i] += (w0 * wPS[i]); } } @@ -66,62 +69,55 @@ namespace { incGenOnly(w0); // then add up variations if (!wScale.empty()) { - if (sumScale.empty()) - sumScale.resize(wScale.size(), 0); + if (sumScale_.empty()) + sumScale_.resize(wScale.size(), 0); for (unsigned int i = 0, n = wScale.size(); i < n; ++i) - sumScale[i] += (w0 * wScale[i]); + sumScale_[i] += (w0 * wScale[i]); } if (!wPDF.empty()) { - if (sumPDF.empty()) - sumPDF.resize(wPDF.size(), 0); + if (sumPDF_.empty()) + sumPDF_.resize(wPDF.size(), 0); for (unsigned int i = 0, n = wPDF.size(); i < n; ++i) - sumPDF[i] += (w0 * wPDF[i]); + sumPDF_[i] += (w0 * wPDF[i]); } if (!wRwgt.empty()) { - if (sumRwgt.empty()) - sumRwgt.resize(wRwgt.size(), 0); + if (sumRwgt_.empty()) + sumRwgt_.resize(wRwgt.size(), 0); for (unsigned int i = 0, n = wRwgt.size(); i < n; ++i) - sumRwgt[i] += (w0 * wRwgt[i]); + sumRwgt_[i] += (w0 * wRwgt[i]); } if (!wNamed.empty()) { - if (sumNamed.empty()) - sumNamed.resize(wNamed.size(), 0); + if (sumNamed_.empty()) + sumNamed_.resize(wNamed.size(), 0); for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) - sumNamed[i] += (w0 * wNamed[i]); + sumNamed_[i] += (w0 * wNamed[i]); } incPSOnly(w0, wPS); } void merge(const Counter& other) { - num += other.num; - sumw += other.sumw; - sumw2 += other.sumw2; - if (sumScale.empty() && !other.sumScale.empty()) - sumScale.resize(other.sumScale.size(), 0); - if (sumPDF.empty() && !other.sumPDF.empty()) - sumPDF.resize(other.sumPDF.size(), 0); - if (sumRwgt.empty() && !other.sumRwgt.empty()) - sumRwgt.resize(other.sumRwgt.size(), 0); - if (sumNamed.empty() && !other.sumNamed.empty()) - sumNamed.resize(other.sumNamed.size(), 0); - if (sumPS.empty() && !other.sumPS.empty()) - sumPS.resize(other.sumPS.size(), 0); - if (!other.sumScale.empty()) - for (unsigned int i = 0, n = sumScale.size(); i < n; ++i) - sumScale[i] += other.sumScale[i]; - if (!other.sumPDF.empty()) - for (unsigned int i = 0, n = sumPDF.size(); i < n; ++i) - sumPDF[i] += other.sumPDF[i]; - if (!other.sumRwgt.empty()) - for (unsigned int i = 0, n = sumRwgt.size(); i < n; ++i) - sumRwgt[i] += other.sumRwgt[i]; - if (!other.sumNamed.empty()) - for (unsigned int i = 0, n = sumNamed.size(); i < n; ++i) - sumNamed[i] += other.sumNamed[i]; - if (!other.sumPS.empty()) - for (unsigned int i = 0, n = sumPS.size(); i < n; ++i) - sumPS[i] += other.sumPS[i]; + num_ += other.num_; + sumw_ += other.sumw_; + sumw2_ += other.sumw2_; + + mergeSumVectors(sumScale_, other.sumScale_); + mergeSumVectors(sumPDF_, other.sumPDF_); + mergeSumVectors(sumRwgt_, other.sumRwgt_); + mergeSumVectors(sumNamed_, other.sumNamed_); + mergeSumVectors(sumPS_, other.sumPS_); } + + //private: + // the counters + long long num_ = 0; + long double sumw_ = 0; + long double sumw2_ = 0; + + std::vector sumPDF_; + std::vector sumScale_; + std::vector sumRwgt_; + std::vector sumNamed_; + std::vector sumPS_; }; struct CounterMap { @@ -233,10 +229,6 @@ class GenWeightsTableProducer : public edm::global::EDProducer(); produces("genModel"); - produces("LHEScale"); - produces("LHEPdf"); - produces("LHEReweighting"); - produces("LHENamed"); produces("PS"); produces(); if (namedWeightIDs_.size() != namedWeightLabels_.size()) { @@ -255,17 +247,15 @@ class GenWeightsTableProducer : public edm::global::EDProducerget(); + Counter& counter = *streamCache(id)->get(); // generator information (always available) - edm::Handle genInfo; - iEvent.getByToken(genTag_, genInfo); - double weight = genInfo->weight(); + auto const& genInfo = iEvent.get(genTag_); // table for gen info, always available auto out = std::make_unique(1, "genWeight", true); out->setDoc("generator weight"); - out->addColumnValue("", weight, "generator weight", nanoaod::FlatTable::FloatColumn); + out->addColumnValue("", genInfo.weight(), "generator weight", nanoaod::FlatTable::FloatColumn); iEvent.put(std::move(out)); std::string model_label = streamCache(id)->getLabel(); @@ -273,7 +263,6 @@ class GenWeightsTableProducer : public edm::global::EDProducer lheScaleTab, lhePdfTab, lheRwgtTab, lheNamedTab; std::unique_ptr genPSTab; edm::Handle lheInfo; @@ -285,51 +274,41 @@ class GenWeightsTableProducer : public edm::global::EDProducer& 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) + // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind) + bool lheDebug = debug_.exchange(false); - const std::vector& scaleWeightIDs = weightChoice->scaleWeightIDs; - const std::vector& pdfWeightIDs = weightChoice->pdfWeightIDs; - const std::vector& rwgtWeightIDs = weightChoice->rwgtIDs; + 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); + std::vector wScale(scaleWeightIDs.size(), 1); + std::vector wPDF(pdfWeightIDs.size(), 1); + std::vector wRwgt(rwgtWeightIDs.size(), 1); + std::vector 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()); @@ -358,7 +337,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer(wPS.size(), "PSWeight", false); outPS->addColumn("", wPS, vectorSize > 1 ? "PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 " @@ -367,35 +346,10 @@ class GenWeightsTableProducer : public edm::global::EDProduceraddColumn( - "", wScale, weightChoice->scaleWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_); - - outPdf.reset(new nanoaod::FlatTable(wPDF.size(), "LHEPdfWeight", false)); - outPdf->addColumn( - "", wPDF, weightChoice->pdfWeightsDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_); - - outRwgt.reset(new nanoaod::FlatTable(wRwgt.size(), "LHEReweightingWeight", false)); - outRwgt->addColumn( - "", wRwgt, weightChoice->rwgtWeightDoc, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_); - - outNamed.reset(new nanoaod::FlatTable(1, "LHEWeight", true)); - outNamed->addColumnValue("originalXWGTUP", - lheProd.originalXWGTUP(), - "Nominal event weight in the LHE file", - nanoaod::FlatTable::FloatColumn); - for (unsigned int i = 0, n = wNamed.size(); i < n; ++i) { - outNamed->addColumnValue(namedWeightLabels_[i], - wNamed[i], - "LHE weight for id " + namedWeightIDs_[i] + ", relative to nominal", - nanoaod::FlatTable::FloatColumn, - lheWeightPrecision_); - } - - counter->incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS); + counter.incLHE(genWeight, wScale, wPDF, wRwgt, wNamed, wPS); } - void fillOnlyPSWeightTable(Counter* counter, + void fillOnlyPSWeightTable(Counter& counter, double genWeight, const GenEventInfoProduct& genProd, std::unique_ptr& outPS) const { @@ -408,7 +362,7 @@ class GenWeightsTableProducer : public edm::global::EDProducer(wPS.size(), "PSWeight", false); outPS->addColumn("", wPS, vectorSize > 1 ? "PS weights (w_var / w_nominal); [0] is ISR=0.5 FSR=1; [1] is ISR=1 " @@ -417,16 +371,16 @@ class GenWeightsTableProducer : public edm::global::EDProducerincGenOnly(genWeight); - counter->incPSOnly(genWeight, wPS); + counter.incGenOnly(genWeight); + counter.incPSOnly(genWeight, wPS); } // 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) + // make sure only the first thread dumps out this (even if may still be mixed up with other output, but nevermind) + bool lheDebug = debugRun_.exchange(false); auto weightChoice = std::make_shared(); // getByToken throws since we're not in the endRun (see https://github.com/cms-sw/cmssw/pull/18499) @@ -797,40 +751,40 @@ class GenWeightsTableProducer : public edm::global::EDProducer(); for (auto x : runCounterMap->countermap) { - auto runCounter = &(x.second); + auto& runCounter = x.second; std::string label = std::string("_") + x.first; std::string doclabel = (!x.first.empty()) ? (std::string(", for model label ") + x.first) : ""; - out->addInt("genEventCount" + label, "event count" + doclabel, runCounter->num); - out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter->sumw); - out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter->sumw2); + out->addInt("genEventCount" + label, "event count" + doclabel, runCounter.num_); + out->addFloat("genEventSumw" + label, "sum of gen weights" + doclabel, runCounter.sumw_); + out->addFloat("genEventSumw2" + label, "sum of gen (weight^2)" + doclabel, runCounter.sumw2_); - double norm = runCounter->sumw ? 1.0 / runCounter->sumw : 1; - auto sumScales = runCounter->sumScale; + double norm = runCounter.sumw_ ? 1.0 / runCounter.sumw_ : 1; + auto sumScales = runCounter.sumScale_; for (auto& val : sumScales) val *= norm; out->addVFloat("LHEScaleSumw" + label, "Sum of genEventWeight * LHEScaleWeight[i], divided by genEventSumw" + doclabel, sumScales); - auto sumPDFs = runCounter->sumPDF; + auto sumPDFs = runCounter.sumPDF_; for (auto& val : sumPDFs) val *= norm; out->addVFloat( "LHEPdfSumw" + label, "Sum of genEventWeight * LHEPdfWeight[i], divided by genEventSumw" + doclabel, sumPDFs); - if (!runCounter->sumRwgt.empty()) { - auto sumRwgts = runCounter->sumRwgt; + if (!runCounter.sumRwgt_.empty()) { + auto sumRwgts = runCounter.sumRwgt_; for (auto& val : sumRwgts) val *= norm; out->addVFloat("LHEReweightingSumw" + label, "Sum of genEventWeight * LHEReweightingWeight[i], divided by genEventSumw" + doclabel, sumRwgts); } - if (!runCounter->sumNamed.empty()) { // it could be empty if there's no LHE info in the sample + if (!runCounter.sumNamed_.empty()) { // it could be empty if there's no LHE info in the sample for (unsigned int i = 0, n = namedWeightLabels_.size(); i < n; ++i) { out->addFloat( "LHESumw_" + namedWeightLabels_[i] + label, "Sum of genEventWeight * LHEWeight_" + namedWeightLabels_[i] + ", divided by genEventSumw" + doclabel, - runCounter->sumNamed[i] * norm); + runCounter.sumNamed_[i] * norm); } } } diff --git a/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc b/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc new file mode 100644 index 0000000000000..5587cbd3d31e6 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc @@ -0,0 +1,220 @@ +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/NanoAOD/interface/FlatTable.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.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 "FWCore/Utilities/interface/transform.h" + + +#include +#include +#include + +namespace { + typedef std::vector WeightGroupDataContainer; + typedef std::array, 2> WeightGroupsToStore; +} // namespace + +class LHEWeightsTableProducer : public edm::global::EDProducer> { +public: + 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"))), + 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")) { + if (weightgroups_.size() != maxGroupsPerType_.size()) + throw std::invalid_argument("Inputs 'weightgroups' and 'weightgroupNums' must have equal size"); + produces>(); + } + + void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override { + edm::Handle lheWeightHandle; + for (auto& token : lheWeightTokens_) { + iEvent.getByToken(token, lheWeightHandle); + if (lheWeightHandle.isValid()) { + break; + } + } + + const GenWeightProduct* lheWeightProduct = lheWeightHandle.product(); + WeightsContainer lheWeights = lheWeightProduct->weights(); + + edm::Handle genWeightHandle; + iEvent.getByToken(genWeightToken_, genWeightHandle); + const GenWeightProduct* genWeightProduct = genWeightHandle.product(); + WeightsContainer genWeights = genWeightProduct->weights(); + + auto lheWeightTables = std::make_unique>(); + auto const& weightInfos = *luminosityBlockCache(iEvent.getLuminosityBlock().index()); + + addWeightGroupToTable(lheWeightTables, "LHE", weightInfos.at(inLHE), lheWeights); + addWeightGroupToTable(lheWeightTables, "Gen", weightInfos.at(inGen), genWeights); + + iEvent.put(std::move(lheWeightTables)); + } + + void addWeightGroupToTable(std::unique_ptr>& lheWeightTables, const char* typeName, + const WeightGroupDataContainer& weightInfos, WeightsContainer& allWeights) const { + size_t typeCount = 0; + gen::WeightType previousType = gen::WeightType::kUnknownWeights; + + for (const auto& groupInfo : weightInfos) { + std::string entryName = typeName; + gen::WeightType weightType = groupInfo.group->weightType(); + if (previousType != weightType) + typeCount = 0; + + std::string name = weightTypeNames_.at(weightType); + std::string label = groupInfo.group->name(); + + auto& weights = allWeights.at(groupInfo.index); + label.append("; "); + if (weightType == gen::WeightType::kScaleWeights && groupInfo.group->isWellFormed() + && groupInfo.group->nIdsContained() < 10) { + weights = orderedScaleWeights(weights, + dynamic_cast(groupInfo.group)); + label.append("[1] is mur=0.5 muf=1; [2] is mur=0.5 muf=2; [3] is mur=1 muf=0.5 ;" + " [4] is mur=1 muf=1; [5] is mur=1 muf=2; [6] is mur=2 muf=0.5;" + " [7] is mur=2 muf=1 ; [8] is mur=2 muf=2)"); + } + else + label.append(groupInfo.group->description()); + + entryName.append(name); + entryName.append("Weight"); + if (typeCount > 0) { + entryName.append("AltSet"); + entryName.append(std::to_string(typeCount)); + } + + lheWeightTables->emplace_back(weights.size(), entryName, false); + lheWeightTables->back().addColumn( + "", weights, label, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_); + + typeCount++; + previousType = weightType; + } + } + + + 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 + edm::Handle lheWeightInfoHandle; + for (auto& token : lheWeightInfoTokens_) { + iLumi.getByToken(token, lheWeightInfoHandle); + if (lheWeightInfoHandle.isValid()) { + break; + } + } + + 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) { + auto lheWeights = weightDataPerType(lheWeightInfoHandle, weightType, storePerType[weightType]); + weightsToStore.at(inLHE).insert(weightsToStore.at(inLHE).end(), lheWeights.begin(), lheWeights.end()); + + auto genWeights = weightDataPerType(genWeightInfoHandle, weightType, storePerType[weightType]); + weightsToStore.at(inGen).insert(weightsToStore.at(inGen).end(), genWeights.begin(), genWeights.end()); + } + + return std::make_shared(weightsToStore); + } + + WeightGroupDataContainer weightDataPerType(edm::Handle& weightsHandle, + gen::WeightType weightType, int& maxStore) const { + WeightGroupDataContainer group; + if (weightType == gen::WeightType::kPdfWeights && pdfIds_.size() > 0) { + group = weightsHandle->pdfGroupsWithIndicesByLHAIDs(pdfIds_); + } + else + group = weightsHandle->weightGroupsAndIndicesByType(weightType); + + if (maxStore < 0 || static_cast(group.size()) <= maxStore) { + // Modify size in case one type of weight is present in multiple products + maxStore -= group.size(); + return group; + } + return std::vector(group.begin(), group.begin()+maxStore); + } + + + std::vector orderedScaleWeights(const std::vector& scaleWeights, const gen::ScaleWeightGroupInfo* scaleGroup) const { + + std::vector weights; + weights.push_back(scaleWeights.at(scaleGroup->muR05muF05Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR05muF1Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR05muF2Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR1muF05Index())); + weights.push_back(scaleWeights.at(scaleGroup->centralIndex())); + weights.push_back(scaleWeights.at(scaleGroup->muR1muF2Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR2muF05Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR2muF1Index())); + weights.push_back(scaleWeights.at(scaleGroup->muR2muF2Index())); + + return weights; + } + + // nothing to do here + virtual void globalEndLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add>("lheWeights"); + //desc.add>("genWeights"); + desc.add("genWeights"); + desc.add>("weightgroups"); + desc.add>("maxGroupsPerType"); + desc.addOptionalUntracked>("pdfIds"); + desc.add("lheWeightPrecision", -1)->setComment("Number of bits in the mantissa for LHE weights"); + descriptions.addDefault(desc); + } + +protected: + const edm::EDGetTokenT lheToken_; + const std::vector> lheWeightTokens_; + const std::vector> lheWeightInfoTokens_; + const edm::EDGetTokenT genWeightToken_; + const edm::EDGetTokenT genWeightInfoToken_; + const std::vector weightgroups_; + const std::vector maxGroupsPerType_; + const std::vector pdfIds_; + const std::unordered_map weightTypeNames_ = { + {gen::WeightType::kScaleWeights, "Scale"}, + {gen::WeightType::kPdfWeights, "Pdf"}, + {gen::WeightType::kMEParamWeights, "MEParam"}, + {gen::WeightType::kPartonShowerWeights, "PartonShower"}, + {gen::WeightType::kUnknownWeights, "Unknown"}, + }; + + //std::unordered_map weightGroupIndices_; + int lheWeightPrecision_; + enum {inLHE, inGen}; +}; +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(LHEWeightsTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc b/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc index 4499a23ef8cfb..8f6464c095d6b 100644 --- a/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc +++ b/PhysicsTools/NanoAOD/plugins/NanoAODOutputModule.cc @@ -46,7 +46,6 @@ class NanoAODOutputModule : public edm::one::OutputModule<> { public: NanoAODOutputModule(edm::ParameterSet const& pset); - ~NanoAODOutputModule() override; static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); @@ -118,6 +117,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; std::vector m_evstrings; @@ -149,8 +150,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; @@ -191,9 +190,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; + } } // fill triggers for (auto& t : m_triggers) @@ -270,14 +289,17 @@ void NanoAODOutputModule::openFile(edm::FileBlock const&) { } /* Setup file structure here */ m_tables.clear(); + m_tableTokens.clear(); m_triggers.clear(); 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 1c72a63cb4fdb..4f411b12dfe55 100644 --- a/PhysicsTools/NanoAOD/plugins/TableOutputBranches.cc +++ b/PhysicsTools/NanoAOD/plugins/TableOutputBranches.cc @@ -57,15 +57,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/genWeightsTable_cfi.py b/PhysicsTools/NanoAOD/python/genWeightsTable_cfi.py new file mode 100644 index 0000000000000..15c2ea05d814b --- /dev/null +++ b/PhysicsTools/NanoAOD/python/genWeightsTable_cfi.py @@ -0,0 +1,20 @@ +import FWCore.ParameterSet.Config as cms + +genWeightsTable = cms.EDProducer("GenWeightsTableProducer", + genEvent = cms.InputTag("generator"), + genLumiInfoHeader = cms.InputTag("generator"), + lheInfo = cms.VInputTag(cms.InputTag("externalLHEProducer"), cms.InputTag("source")), + preferredPDFs = cms.VPSet( # see https://lhapdf.hepforge.org/pdfsets.html + cms.PSet( name = cms.string("PDF4LHC15_nnlo_30_pdfas"), lhaid = cms.uint32(91400) ), + cms.PSet( name = cms.string("NNPDF31_nnlo_hessian_pdfas"), lhaid = cms.uint32(306000) ), + cms.PSet( name = cms.string("NNPDF30_nlo_as_0118"), lhaid = cms.uint32(260000) ), # for some 92X samples. Note that the nominal weight, 260000, is not included in the LHE ... + cms.PSet( name = cms.string("NNPDF30_lo_as_0130"), lhaid = cms.uint32(262000) ), # some MLM 80X samples have only this (e.g. /store/mc/RunIISummer16MiniAODv2/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-madgraphMLM-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6_ext1-v2/120000/02A210D6-F5C3-E611-B570-008CFA197BD4.root ) + cms.PSet( name = cms.string("NNPDF30_nlo_nf_4_pdfas"), lhaid = cms.uint32(292000) ), # some FXFX 80X samples have only this (e.g. WWTo1L1Nu2Q, WWTo4Q) + cms.PSet( name = cms.string("NNPDF30_nlo_nf_5_pdfas"), lhaid = cms.uint32(292200) ), # some FXFX 80X samples have only this (e.g. DYJetsToLL_Pt, WJetsToLNu_Pt, DYJetsToNuNu_Pt) + ), + namedWeightIDs = cms.vstring(), + namedWeightLabels = cms.vstring(), + lheWeightPrecision = cms.int32(14), + maxPdfWeights = cms.uint32(150), + debug = cms.untracked.bool(False), +) diff --git a/PhysicsTools/NanoAOD/python/nano_cff.py b/PhysicsTools/NanoAOD/python/nano_cff.py index 466b4c4554040..bfc9951ebb39e 100644 --- a/PhysicsTools/NanoAOD/python/nano_cff.py +++ b/PhysicsTools/NanoAOD/python/nano_cff.py @@ -16,6 +16,7 @@ from PhysicsTools.NanoAOD.triggerObjects_cff import * from PhysicsTools.NanoAOD.isotracks_cff import * from PhysicsTools.NanoAOD.NanoAODEDMEventContent_cff import * +from PhysicsTools.NanoAOD.genWeightsTable_cfi import * from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy from Configuration.Eras.Modifier_run2_nanoAOD_94X2016_cff import run2_nanoAOD_94X2016 diff --git a/PhysicsTools/NanoAOD/python/nanogen_cff.py b/PhysicsTools/NanoAOD/python/nanogen_cff.py new file mode 100644 index 0000000000000..66c938ac408ae --- /dev/null +++ b/PhysicsTools/NanoAOD/python/nanogen_cff.py @@ -0,0 +1,111 @@ +from PhysicsTools.NanoAOD.taus_cff import * +from PhysicsTools.NanoAOD.jets_cff import * +from PhysicsTools.NanoAOD.globals_cff import * +from PhysicsTools.NanoAOD.genparticles_cff import * +from PhysicsTools.NanoAOD.particlelevel_cff import * +from PhysicsTools.NanoAOD.lheInfoTable_cfi import * +from PhysicsTools.NanoAOD.genWeightsTable_cfi import * + +genWeights = cms.EDProducer("GenWeightProductProducer", + genInfo = cms.InputTag("generator"), + genLumiInfoHeader = cms.InputTag("generator")) + +lheWeights = cms.EDProducer("LHEWeightProductProducer", + lheSourceLabel = cms.string("externalLHEProducer")) + +lheWeightsTable = cms.EDProducer( + "LHEWeightsTableProducer", + lheWeights = cms.VInputTag(["externalLHEProducer", "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, 2, 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), +) + +nanoMetadata = cms.EDProducer("UniqueStringProducer", + strings = cms.PSet( + tag = cms.string("untagged"), + ) +) + +metGenTable = cms.EDProducer("SimpleCandidateFlatTableProducer", + src = cms.InputTag("genMetTrue"), + name = cms.string("GenMET"), + doc = cms.string("Gen MET"), + singleton = cms.bool(True), + extension = cms.bool(False), + variables = cms.PSet( + pt = Var("pt", float, doc="pt", precision=10), + phi = Var("phi", float, doc="phi", precision=10), + ), +) + +nanogenSequence = cms.Sequence( + genWeights+ + lheWeights+ + nanoMetadata+ + particleLevel+ + genJetTable+ + patJetPartons+ + genJetFlavourAssociation+ + genJetFlavourTable+ + genJetAK8Table+ + genJetAK8FlavourAssociation+ + genJetAK8FlavourTable+ + tauGenJets+ + tauGenJetsSelectorAllHadrons+ + genVisTaus+ + genVisTauTable+ + genTable+ + genWeightsTable+ + lheWeightsTable+ + genParticleTables+ + tautagger+ + rivetProducerHTXS+ + particleLevelTables+ + metGenTable+ + lheInfoTable +) + +NANOAODGENoutput = cms.OutputModule("NanoAODOutputModule", + compressionAlgorithm = cms.untracked.string('LZMA'), + compressionLevel = cms.untracked.int32(9), + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('NANOAODSIM'), + filterName = cms.untracked.string('') + ), + fileName = cms.untracked.string('nanogen.root'), + outputCommands = cms.untracked.vstring( + 'drop *', + "keep *_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 + ) +) + +def customizeNanoGEN(process): + process.genParticleTable.src = "genParticles" + process.patJetPartons.particles = "genParticles" + process.particleLevel.src = "generatorSmeared" + process.rivetProducerHTXS.HepMCCollection = "generatorSmeared" + + process.genJetTable.src = "ak4GenJetsNoNu" + process.genJetFlavourAssociation.jets = process.genJetTable.src + process.genJetFlavourTable.src = process.genJetTable.src + process.genJetFlavourTable.jetFlavourInfos = "genJetFlavourAssociation" + process.genJetAK8Table.src = "ak8GenJetsNoNu" + process.genJetAK8FlavourAssociation.jets = process.genJetAK8Table.src + process.genJetAK8FlavourTable.src = process.genJetAK8Table.src + process.tauGenJets.GenParticles = "genParticles" + process.genVisTaus.srcGenParticles = "genParticles" + return process diff --git a/PhysicsTools/NanoAOD/test/nanoGen_test_cfg.py b/PhysicsTools/NanoAOD/test/nanoGen_test_cfg.py new file mode 100644 index 0000000000000..bf6124b75a66f --- /dev/null +++ b/PhysicsTools/NanoAOD/test/nanoGen_test_cfg.py @@ -0,0 +1,169 @@ +import FWCore.ParameterSet.Config as cms + + + +process = cms.Process('GEN') + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.Generator_cff') +process.load('IOMC.EventVertexGenerators.VtxSmearedRealistic50ns13TeVCollision_cfi') +process.load('GeneratorInterface.Core.genFilterSummary_cff') +process.load('PhysicsTools.NanoAOD.nanogen_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10), + output = cms.optional.untracked.allowed(cms.int32,cms.PSet) +) + +# Input source +process.source = cms.Source("EmptySource") + +process.options = cms.untracked.PSet( + FailPath = cms.untracked.vstring(), + IgnoreCompletely = cms.untracked.vstring(), + Rethrow = cms.untracked.vstring(), + SkipEvent = cms.untracked.vstring(), + allowUnscheduled = cms.obsolete.untracked.bool, + canDeleteEarly = cms.untracked.vstring(), + emptyRunLumiMode = cms.obsolete.untracked.string, + eventSetup = cms.untracked.PSet( + forceNumberOfConcurrentIOVs = cms.untracked.PSet( + + ), + numberOfConcurrentIOVs = cms.untracked.uint32(1) + ), + fileMode = cms.untracked.string('FULLMERGE'), + forceEventSetupCacheClearOnNewRun = cms.untracked.bool(False), + makeTriggerResults = cms.obsolete.untracked.bool, + numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1), + numberOfConcurrentRuns = cms.untracked.uint32(1), + numberOfStreams = cms.untracked.uint32(0), + numberOfThreads = cms.untracked.uint32(1), + printDependencies = cms.untracked.bool(False), + sizeOfStackForThreadsInKB = cms.optional.untracked.uint32, + throwIfIllegalParameter = cms.untracked.bool(True), + wantSummary = cms.untracked.bool(False) +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('Configuration/GenProduction/python/SMP-RunIISummer15wmLHEGS-00344-fragment.py nevts:10'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Other statements +process.genstepfilter.triggerConditions=cms.vstring("generation_step") + +process.generator = cms.EDFilter("Pythia8HadronizerFilter", + PythiaParameters = cms.PSet( + parameterSets = cms.vstring( + 'pythia8CommonSettings', + 'pythia8CUEP8M1Settings', + 'pythia8PowhegEmissionVetoSettings', + 'pythia8PSweightsSettings', + 'processParameters' + ), + processParameters = cms.vstring( + 'POWHEG:nFinal = 2', + 'ParticleDecays:allowPhotonRadiation = on', + 'TimeShower:QEDshowerByL = on' + ), + pythia8CUEP8M1Settings = cms.vstring( + 'Tune:pp 14', + 'Tune:ee 7', + 'MultipartonInteractions:pT0Ref=2.4024', + 'MultipartonInteractions:ecmPow=0.25208', + 'MultipartonInteractions:expPow=1.6' + ), + pythia8CommonSettings = cms.vstring( + 'Tune:preferLHAPDF = 2', + 'Main:timesAllowErrors = 10000', + 'Check:epTolErr = 0.01', + 'Beams:setProductionScalesFromLHEF = off', + 'SLHA:keepSM = on', + 'SLHA:minMassSM = 1000.', + 'ParticleDecays:limitTau0 = on', + 'ParticleDecays:tau0Max = 10', + 'ParticleDecays:allowPhotonRadiation = on' + ), + pythia8PSweightsSettings = cms.vstring( + 'UncertaintyBands:doVariations = on', + 'UncertaintyBands:List = {isrRedHi isr:muRfac=0.707,fsrRedHi fsr:muRfac=0.707,isrRedLo isr:muRfac=1.414,fsrRedLo fsr:muRfac=1.414,isrDefHi isr:muRfac=0.5,fsrDefHi fsr:muRfac=0.5,isrDefLo isr:muRfac=2.0,fsrDefLo fsr:muRfac=2.0,isrConHi isr:muRfac=0.25,fsrConHi fsr:muRfac=0.25,isrConLo isr:muRfac=4.0,fsrConLo fsr:muRfac=4.0,fsr_G2GG_muR_dn fsr:G2GG:muRfac=0.5,fsr_G2GG_muR_up fsr:G2GG:muRfac=2.0,fsr_G2QQ_muR_dn fsr:G2QQ:muRfac=0.5,fsr_G2QQ_muR_up fsr:G2QQ:muRfac=2.0,fsr_Q2QG_muR_dn fsr:Q2QG:muRfac=0.5,fsr_Q2QG_muR_up fsr:Q2QG:muRfac=2.0,fsr_X2XG_muR_dn fsr:X2XG:muRfac=0.5,fsr_X2XG_muR_up fsr:X2XG:muRfac=2.0,fsr_G2GG_cNS_dn fsr:G2GG:cNS=-2.0,fsr_G2GG_cNS_up fsr:G2GG:cNS=2.0,fsr_G2QQ_cNS_dn fsr:G2QQ:cNS=-2.0,fsr_G2QQ_cNS_up fsr:G2QQ:cNS=2.0,fsr_Q2QG_cNS_dn fsr:Q2QG:cNS=-2.0,fsr_Q2QG_cNS_up fsr:Q2QG:cNS=2.0,fsr_X2XG_cNS_dn fsr:X2XG:cNS=-2.0,fsr_X2XG_cNS_up fsr:X2XG:cNS=2.0,isr_G2GG_muR_dn isr:G2GG:muRfac=0.5,isr_G2GG_muR_up isr:G2GG:muRfac=2.0,isr_G2QQ_muR_dn isr:G2QQ:muRfac=0.5,isr_G2QQ_muR_up isr:G2QQ:muRfac=2.0,isr_Q2QG_muR_dn isr:Q2QG:muRfac=0.5,isr_Q2QG_muR_up isr:Q2QG:muRfac=2.0,isr_X2XG_muR_dn isr:X2XG:muRfac=0.5,isr_X2XG_muR_up isr:X2XG:muRfac=2.0,isr_G2GG_cNS_dn isr:G2GG:cNS=-2.0,isr_G2GG_cNS_up isr:G2GG:cNS=2.0,isr_G2QQ_cNS_dn isr:G2QQ:cNS=-2.0,isr_G2QQ_cNS_up isr:G2QQ:cNS=2.0,isr_Q2QG_cNS_dn isr:Q2QG:cNS=-2.0,isr_Q2QG_cNS_up isr:Q2QG:cNS=2.0,isr_X2XG_cNS_dn isr:X2XG:cNS=-2.0,isr_X2XG_cNS_up isr:X2XG:cNS=2.0}', + 'UncertaintyBands:nFlavQ = 4', + 'UncertaintyBands:MPIshowers = on', + 'UncertaintyBands:overSampleFSR = 10.0', + 'UncertaintyBands:overSampleISR = 10.0', + 'UncertaintyBands:FSRpTmin2Fac = 20', + 'UncertaintyBands:ISRpTmin2Fac = 1' + ), + pythia8PowhegEmissionVetoSettings = cms.vstring( + 'POWHEG:veto = 1', + 'POWHEG:pTdef = 1', + 'POWHEG:emitted = 0', + 'POWHEG:pTemt = 0', + 'POWHEG:pThard = 0', + 'POWHEG:vetoCount = 100', + 'SpaceShower:pTmaxMatch = 2', + 'TimeShower:pTmaxMatch = 2' + ) + ), + comEnergy = cms.double(13000.0), + filterEfficiency = cms.untracked.double(1.0), + maxEventsToPrint = cms.untracked.int32(1), + pythiaHepMCVerbosity = cms.untracked.bool(False), + pythiaPylistVerbosity = cms.untracked.int32(1) +) + +process.externalLHEProducer = cms.EDProducer("ExternalLHEProducer", + args = cms.vstring('/afs/cern.ch/user/k/kelong/work/public/DummyGridpacks/VVV_aQGCfs_dummy.tgz'), + nEvents = cms.untracked.uint32(10), + numberOfParameters = cms.uint32(1), + outputFile = cms.string('cmsgrid_final.lhe'), + scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh') +) + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc', '') + +process.ProductionFilterSequence = cms.Sequence(process.generator) + +# Path and EndPath definitions +process.lhe_step = cms.Path(process.externalLHEProducer) +process.generation_step = cms.Path(process.pgen) +process.genfiltersummary_step = cms.EndPath(process.genFilterSummary) + +# Path and EndPath definitions +process.nanoGEN_step = cms.Path(process.genWeights*process.nanogenSequence) +process.endjob_step = cms.EndPath(process.endOfProcess) +process.NANOAODGENoutput_step = cms.EndPath(process.NANOAODGENoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.lhe_step,process.generation_step,process.genfiltersummary_step,process.nanoGEN_step,process.endjob_step,process.NANOAODGENoutput_step) + +# filter all path with the production filter sequence +for path in process.paths: + if path in ['lhe_step']: continue + getattr(process,path).insert(0, process.ProductionFilterSequence) + + +process.RandomNumberGeneratorService.externalLHEProducer.initialSeed=int(96) +# Customisation from command line +from PhysicsTools.NanoAOD.nanogen_cff import customizeNanoGEN +process = customizeNanoGEN(process) + + +# Add early deletion of temporary data products to reduce peak memory need +from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete +process = customiseEarlyDelete(process) +# End adding early deletion + diff --git a/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h b/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h new file mode 100644 index 0000000000000..f990ed1931ea5 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h @@ -0,0 +1,51 @@ +#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; + const gen::WeightGroupInfo* group; + }; +} // namespace + +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; + const gen::WeightGroupInfo* containingWeightGroupInfo(int index) const; + const gen::WeightGroupInfo* 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(gen::WeightGroupInfo* 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..4675e17c6c80d --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h @@ -0,0 +1,47 @@ +#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" + +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_.size() == 0 && setEntry == 0) + addWeightSet(); + if (static_cast(weightsVector_.size()) <= setEntry) + throw std::domain_error("Out of range weight"); + 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/MEParamWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h new file mode 100644 index 0000000000000..3b2dd09d37609 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h @@ -0,0 +1,23 @@ +#ifndef SimDataFormats_GeneratorProducts_MEParamWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_MEParamWeightGroupInfo_h + +#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) {} + virtual ~MEParamWeightGroupInfo() override {} + void copy(const MEParamWeightGroupInfo &other); + virtual MEParamWeightGroupInfo* clone() const override; + }; +} + +#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..ea6c1b1fba683 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h @@ -0,0 +1,29 @@ +#ifndef SimDataFormats_GeneratorProducts_PartonShowerWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_PartonShowerWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" + +namespace gen { + class PartonShowerWeightGroupInfo : public WeightGroupInfo { + public: + PartonShowerWeightGroupInfo() : PartonShowerWeightGroupInfo("") {} + PartonShowerWeightGroupInfo(std::string header, std::string name) : + WeightGroupInfo(header, name) { + weightType_ = WeightType::kPartonShowerWeights; + } + PartonShowerWeightGroupInfo(std::string header) : + PartonShowerWeightGroupInfo(header, header) { } + PartonShowerWeightGroupInfo(const PartonShowerWeightGroupInfo &other) { + copy(other); + } + virtual ~PartonShowerWeightGroupInfo() override {} + void copy(const PartonShowerWeightGroupInfo &other); + virtual PartonShowerWeightGroupInfo* clone() const override; + + // Is a variation of the functional form of the dynamic scale + }; +} + +#endif + + diff --git a/SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h new file mode 100644 index 0000000000000..18dc88fdbd067 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.h @@ -0,0 +1,74 @@ +#ifndef SimDataFormats_GeneratorProducts_PdfWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_PdfWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" +#include +#include + +namespace gen { + enum PdfUncertaintyType { + kHessianUnc, + kMonteCarloUnc, + kVariationSet, + kUnknownUnc, + }; + + class PdfWeightGroupInfo : public WeightGroupInfo { + private: + PdfUncertaintyType uncertaintyType_; + bool hasAlphasVars_; + int alphasUpIndex_; + int alphasDownIndex_; + std::vector> lhapdfIdsContained_; + 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); + } + virtual ~PdfWeightGroupInfo() override {} + void copy(const PdfWeightGroupInfo &other); + virtual 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_; } + bool containsMultipleSets() const { return lhapdfIdsContained_.size() > 1; } + bool containsParentLhapdfId(int lhaid, int globalIndex) const { + if (indexOfLhapdfId(lhaid) != -1) + return true; + int parentid = lhaid - (globalIndex - firstId_); + return indexOfLhapdfId(parentid) != -1; + } + bool containsLhapdfId(int lhaid) const { return indexOfLhapdfId(lhaid) != -1; } + int indexOfLhapdfId(int lhaid) const { + for (const auto& id : lhapdfIdsContained_) { + if (id.first == lhaid) + return id.second; + } + return -1; + } + int alphasUpIndex() const { return alphasUpIndex_; } + int alphasDownIndex() const { return alphasDownIndex_; } + void addLhapdfId(int lhaid, int globalIndex) { + lhapdfIdsContained_.push_back(std::make_pair(lhaid, globalIndex-firstId_)); + } + std::vector lhapdfIdsContained() const { + std::vector lhaids; + for (const auto& id : lhapdfIdsContained_) { + lhaids.push_back(id.first); + } + return lhaids; + } + }; +} + + +#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..71b6a6b0a987a --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h @@ -0,0 +1,65 @@ +#ifndef SimDataFormats_GeneratorProducts_ScaleWeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_ScaleWeightGroupInfo_h + +#include "SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h" + +namespace gen { + class ScaleWeightGroupInfo : public WeightGroupInfo { + private: + bool isFunctionalFormVar_; + size_t centralIndex_; + size_t muR1muF2Index_; + size_t muR1muF05Index_; + size_t muR2muF05Index_; + size_t muR2muF1Index_; + size_t muR2muF2Index_; + size_t muR05muF05Index_; + size_t muR05muF1Index_; + size_t muR05muF2Index_; + public: + ScaleWeightGroupInfo() : ScaleWeightGroupInfo("") {} + ScaleWeightGroupInfo(std::string header, std::string name) : + WeightGroupInfo(header, name) { + weightType_ = WeightType::kScaleWeights; + isFunctionalFormVar_ = false; + centralIndex_ = 0; + muR1muF2Index_ = 0; + muR1muF05Index_ = 0; + muR2muF05Index_ = 0; + muR2muF1Index_ = 0; + muR2muF2Index_ = 0; + muR2muF05Index_ = 0; + muR05muF05Index_ = 0; + muR05muF1Index_ = 0; + muR05muF2Index_ = 0; + } + ScaleWeightGroupInfo(std::string header) : + ScaleWeightGroupInfo(header, header) { } + ScaleWeightGroupInfo(const ScaleWeightGroupInfo &other) { + copy(other); + } + virtual ~ScaleWeightGroupInfo() override {} + void copy(const ScaleWeightGroupInfo &other); + virtual ScaleWeightGroupInfo* clone() const override; + + void setMuRMuFIndex(WeightMetaInfo& info, float muR, float muF); + void setMuRMuFIndex(int globalIndex, std::string id, float muR, float muF); + void addContainedId(int weightEntry, std::string id, std::string label, float muR, float muF); + + // Is a variation of the functional form of the dynamic scale + bool isFunctionalFormVariation(); + void setIsFunctionalFormVariation(bool functionalVar) {isFunctionalFormVar_ = functionalVar; } + size_t centralIndex() const {return centralIndex_; } + size_t muR1muF2Index() const { return muR1muF2Index_; } + size_t muR1muF05Index() const { return muR1muF05Index_; } + size_t muR2muF05Index() const { return muR2muF05Index_; } + size_t muR2muF1Index() const { return muR2muF1Index_; } + size_t muR2muF2Index() const { return muR2muF2Index_; } + size_t muR05muF05Index() const { return muR05muF05Index_; } + size_t muR05muF1Index() const { return muR05muF1Index_; } + size_t muR05muF2Index() const { return muR05muF2Index_; } + }; +} + +#endif + diff --git a/SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h b/SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h new file mode 100644 index 0000000000000..66152d8840434 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h @@ -0,0 +1,22 @@ +#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;} + virtual ~UnknownWeightGroupInfo() override {} + void copy(const UnknownWeightGroupInfo &other); + virtual UnknownWeightGroupInfo* clone() const override; + }; +} + +#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..aa097b399cf12 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/interface/WeightGroupInfo.h @@ -0,0 +1,91 @@ +#ifndef SimDataFormats_GeneratorProducts_WeightGroupInfo_h +#define SimDataFormats_GeneratorProducts_WeightGroupInfo_h + +/** \class PdfInfo + * + */ +#include +#include +#include +#include + +namespace gen { + struct WeightMetaInfo { + size_t globalIndex; + size_t localIndex; + std::string id; + std::string 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): + headerEntry_(header), name_(name), firstId_(-1), lastId_(-1) {} + WeightGroupInfo(std::string header): + 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(const std::string& wgtId) const; + bool containsWeight(const std::string& wgtId, int weightEntry) const; + int weightVectorEntry(const std::string& wgtId, int weightEntry) const; + void addContainedId(int weightEntry, std::string id, std::string label); + std::vector containedIds() const; + bool indexInRange(int index) const; + + void setName(std::string name) { name_ = name; } + void setDescription(std::string description) { description_ = description; } + void setHeaderEntry(std::string header) { headerEntry_ = header; } + void setWeightType(WeightType type) { weightType_ = type; } + void setFirstId(int firstId) { firstId_ = firstId; } + void setLastId(int lastId) { lastId_ = lastId; } + + std::string name() const { return name_; } + std::string description() const { return description_; } + std::string headerEntry() const { return headerEntry_; } + WeightType weightType() const { return weightType_; } + std::vector idsContained() const { return idsContained_; } + size_t nIdsContained() const { return idsContained_.size(); } + int firstId() const { return firstId_; } + int lastId() const { return lastId_; } + // Store whether the group was fully parsed succesfully + void setIsWellFormed(bool wellFormed) { isWellFormed_ = wellFormed; } + bool isWellFormed() const { return isWellFormed_; } + + protected: + bool isWellFormed_; + std::string headerEntry_; + std::string name_; + std::string description_; + WeightType weightType_; + std::vector idsContained_; + int firstId_; + int lastId_; + }; +} + +#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..e16bbebcc27da --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/GenWeightInfoProduct.cc @@ -0,0 +1,96 @@ +#include +#include + +#include "SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/PdfWeightGroupInfo.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_; +} + +const gen::WeightGroupInfo* GenWeightInfoProduct::containingWeightGroupInfo(int index) const { + for (const auto& weightGroup : weightGroupsInfo_) { + if (weightGroup.indexInRange(index)) + return &weightGroup; + } + throw std::domain_error("Failed to find containing weight group"); +} + +const gen::WeightGroupInfo* GenWeightInfoProduct::orderedWeightGroupInfo(int weightGroupIndex) const { + if (weightGroupIndex >= static_cast(weightGroupsInfo_.size())) + throw std::range_error("Weight index out of range!"); + return &weightGroupsInfo_[weightGroupIndex]; +} + +std::vector GenWeightInfoProduct::weightGroupsAndIndicesByType(gen::WeightType type) const { + std::vector matchingGroups; + for (size_t i = 0; i < weightGroupsInfo_.size(); i++) { + if (weightGroupsInfo_[i].weightType() == type) + matchingGroups.push_back({i, weightGroupsInfo_[i].clone()}); + } + return matchingGroups; +} + +std::vector GenWeightInfoProduct::weightGroupsByType(gen::WeightType type) const { + std::vector matchingGroups; + for (size_t i = 0; i < weightGroupsInfo_.size(); i++) { + if (weightGroupsInfo_[i].weightType() == type) + matchingGroups.push_back(weightGroupsInfo_[i].clone()); + } + return matchingGroups; +} + +std::optional GenWeightInfoProduct::pdfGroupWithIndexByLHAID(int lhaid) const { + auto pdfGroups = weightGroupsAndIndicesByType(gen::WeightType::kPdfWeights); + + auto matchingPdfSet = std::find_if(pdfGroups.begin(), pdfGroups.end(), + [lhaid] (gen::WeightGroupData& data) { + auto pdfGroup = dynamic_cast(data.group); + return pdfGroup->containsLhapdfId(lhaid); + }); + return matchingPdfSet != pdfGroups.end() ? std::optional(*matchingPdfSet) : std::nullopt; +} + +std::vector GenWeightInfoProduct::pdfGroupsWithIndicesByLHAIDs(const std::vector& lhaids) const { + auto pdfGroups = weightGroupsAndIndicesByType(gen::WeightType::kPdfWeights); + std::vector groups; + + for (auto lhaid : lhaids) { + auto matchingPdfSet = std::find_if(pdfGroups.begin(), pdfGroups.end(), + [lhaid] (gen::WeightGroupData& data) { + auto pdfGroup = dynamic_cast(data.group); + return pdfGroup->containsLhapdfId(lhaid); + }); + if (matchingPdfSet != pdfGroups.end()) + pdfGroups.push_back(*matchingPdfSet); + } + + return pdfGroups; +} + +std::vector GenWeightInfoProduct::weightGroupIndicesByType(gen::WeightType type) const { + std::vector matchingGroupIndices; + for (size_t i = 0; i < weightGroupsInfo_.size(); i++) { + if (weightGroupsInfo_[i].weightType() == type) + matchingGroupIndices.push_back(i); + } + return matchingGroupIndices; +} + +void GenWeightInfoProduct::addWeightGroupInfo(gen::WeightGroupInfo* info) { + weightGroupsInfo_.push_back(info); +} diff --git a/SimDataFormats/GeneratorProducts/src/MEParamWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/MEParamWeightGroupInfo.cc new file mode 100644 index 0000000000000..8d75b537cf890 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/MEParamWeightGroupInfo.cc @@ -0,0 +1,12 @@ +#include "SimDataFormats/GeneratorProducts/interface/MEParamWeightGroupInfo.h" + +namespace gen { + void MEParamWeightGroupInfo::copy(const MEParamWeightGroupInfo &other) { + WeightGroupInfo::copy(other); + } + + MEParamWeightGroupInfo* MEParamWeightGroupInfo::clone() const { + return new MEParamWeightGroupInfo(*this); + } +} + diff --git a/SimDataFormats/GeneratorProducts/src/PartonShowerWeights.cc b/SimDataFormats/GeneratorProducts/src/PartonShowerWeights.cc new file mode 100644 index 0000000000000..893db2f20433e --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/PartonShowerWeights.cc @@ -0,0 +1,12 @@ +#include "SimDataFormats/GeneratorProducts/interface/PartonShowerWeightGroupInfo.h" + +namespace gen { + void PartonShowerWeightGroupInfo::copy(const PartonShowerWeightGroupInfo &other) { + WeightGroupInfo::copy(other); + } + + PartonShowerWeightGroupInfo* PartonShowerWeightGroupInfo::clone() const { + return new PartonShowerWeightGroupInfo(*this); + } +} + diff --git a/SimDataFormats/GeneratorProducts/src/PdfWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/PdfWeightGroupInfo.cc new file mode 100644 index 0000000000000..5dfbf814f323e --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/PdfWeightGroupInfo.cc @@ -0,0 +1,16 @@ +#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(); + lhapdfIdsContained_ = other.lhapdfIdsContained_; + WeightGroupInfo::copy(other); + } + + PdfWeightGroupInfo* PdfWeightGroupInfo::clone() const { + return new PdfWeightGroupInfo(*this); + } +} diff --git a/SimDataFormats/GeneratorProducts/src/ScaleWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/ScaleWeightGroupInfo.cc new file mode 100644 index 0000000000000..575871de6f101 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/ScaleWeightGroupInfo.cc @@ -0,0 +1,57 @@ +#include "SimDataFormats/GeneratorProducts/interface/ScaleWeightGroupInfo.h" +#include +#include + +namespace gen { + void ScaleWeightGroupInfo::copy(const ScaleWeightGroupInfo &other) { + centralIndex_ = other.centralIndex_; + muR1muF2Index_ = other.muR1muF2Index_; + muR1muF05Index_ = other.muR1muF05Index_; + muR2muF05Index_ = other.muR2muF05Index_; + muR2muF1Index_ = other.muR2muF1Index_; + muR2muF2Index_ = other.muR2muF2Index_; + muR05muF05Index_ = other.muR2muF05Index_; + muR05muF1Index_ = other.muR05muF1Index_; + muR05muF2Index_ = other.muR05muF2Index_; + 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) { + WeightGroupInfo::addContainedId(globalIndex, id, label); + setMuRMuFIndex(globalIndex, id, muR, muF); + } + + void ScaleWeightGroupInfo::setMuRMuFIndex(int globalIndex, std::string id, float muR, float muF) { + auto metaInfo = weightMetaInfoByGlobalIndex(id, globalIndex); + setMuRMuFIndex(metaInfo, muR, muF); + } + + void ScaleWeightGroupInfo::setMuRMuFIndex(WeightMetaInfo& info, float muR, float muF) { + if (muR == 0.5 && muF == 0.5) + muR05muF05Index_ = info.localIndex; + else if (muR == 0.5 && muF == 1.0) + muR05muF1Index_ = info.localIndex; + else if (muR == 0.5 && muF == 2.0) + muR05muF2Index_ = info.localIndex; + else if (muR == 1.0 && muF == 0.5) + muR1muF05Index_ = info.localIndex; + else if (muR == 1.0 && muF == 1.0) + centralIndex_ = info.localIndex; + else if (muR == 1.0 && muF == 2.0) + muR1muF2Index_ = info.localIndex; + else if (muR == 2.0 && muF == 0.5) + muR2muF05Index_ = info.localIndex; + else if (muR == 2.0 && muF == 1.0) + muR2muF1Index_ = info.localIndex; + else if (muR == 2.0 && muF == 2.0) + muR2muF2Index_ = info.localIndex; + else + isWellFormed_ = false; + } +} + + diff --git a/SimDataFormats/GeneratorProducts/src/UnknownWeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/UnknownWeightGroupInfo.cc new file mode 100644 index 0000000000000..2d56e8e177b31 --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/UnknownWeightGroupInfo.cc @@ -0,0 +1,8 @@ +#include "SimDataFormats/GeneratorProducts/interface/UnknownWeightGroupInfo.h" + +namespace gen { + UnknownWeightGroupInfo* UnknownWeightGroupInfo::clone() const { + return new UnknownWeightGroupInfo(*this); + } +} + diff --git a/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc b/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc new file mode 100644 index 0000000000000..2b8c5079a443f --- /dev/null +++ b/SimDataFormats/GeneratorProducts/src/WeightGroupInfo.cc @@ -0,0 +1,99 @@ +#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("LogicError", "WeightGroupInfo is abstract, so it's clone() method can't be implemented.\n"); + } + + WeightMetaInfo WeightGroupInfo::weightMetaInfo(int weightEntry) const { + return idsContained_.at(weightEntry); + } + + WeightMetaInfo WeightGroupInfo::weightMetaInfoByGlobalIndex(std::string wgtId, int weightEntry) const { + int entry = weightVectorEntry(wgtId, weightEntry); + if (entry < 0 || entry >= static_cast(idsContained_.size())) + throw std::range_error("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(const std::string& wgtId) const { + return weightVectorEntry(wgtId, 0); + } + + bool WeightGroupInfo::containsWeight(const std::string& wgtId, int weightEntry) const { + return weightVectorEntry(wgtId, weightEntry) != -1; + } + + int WeightGroupInfo::weightVectorEntry(const std::string& wgtId, int weightEntry) const { + // First try ordered search + if (indexInRange(weightEntry)) { + size_t orderedEntry = weightEntry - firstId_; + if (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 (firstId_ == -1 || weightEntry < firstId_) { + firstId_ = weightEntry; + // Reset to reflect that indices will be shifted + for (auto& entry : idsContained_) + entry.localIndex = entry.globalIndex - firstId_; + } + if (weightEntry > lastId_) + lastId_ = weightEntry; + + WeightMetaInfo info; + info.globalIndex = weightEntry; + info.localIndex = weightEntry - firstId_; + info.id = id; + info.label = label; + + if (idsContained_.size() < info.localIndex) { + idsContained_.resize(info.localIndex); + idsContained_.insert(idsContained_.begin()+info.localIndex, info); + } + else if (idsContained_.size() == info.localIndex) { + idsContained_.push_back(info); + } + else { + idsContained_.resize(info.localIndex+1); + idsContained_[info.localIndex] = info; + } + } + + std::vector WeightGroupInfo::containedIds() const { return idsContained_; } + + + bool WeightGroupInfo::indexInRange(int index) const { + return (index <= lastId_ && index >= firstId_); + } +} diff --git a/SimDataFormats/GeneratorProducts/src/classes.h b/SimDataFormats/GeneratorProducts/src/classes.h index 7039e45b3951c..c83887ba46e43 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 05bf110a49686..2ba9fe119446e 100644 --- a/SimDataFormats/GeneratorProducts/src/classes_def.xml +++ b/SimDataFormats/GeneratorProducts/src/classes_def.xml @@ -199,6 +199,8 @@ + + @@ -227,9 +229,33 @@ + + + + + + + + + + + + + + + + + + + + + + + +