Store normalized weights from the start, cleanup

GeneratorInterface/Core/interface/WeightHelper.h

@@ -21,8 +21,8 @@ namespace gen {
         edm::OwnVector<gen::WeightGroupInfo> weightGroups() {
             return weightGroups_;
         }
-        std::unique_ptr<GenWeightProduct> weightProduct(std::vector<gen::WeightsInfo>);
-        std::unique_ptr<GenWeightProduct> weightProduct(std::vector<double>);
+        std::unique_ptr<GenWeightProduct> weightProduct(std::vector<gen::WeightsInfo>, float w0);
+        std::unique_ptr<GenWeightProduct> weightProduct(std::vector<double>, float w0);
         void setGroupInfo();
         bool currentGroupIsScale();
         bool currentGroupIsMEParam();

GeneratorInterface/Core/plugins/GenWeightProductProducer.cc

@@ -43,10 +43,8 @@ class GenWeightProductProducer : public edm::one::EDProducer<edm::BeginLuminosit
 // constructors and destructor
 //
 GenWeightProductProducer::GenWeightProductProducer(const edm::ParameterSet& iConfig) :
-    genLumiInfoToken_(consumes<GenLumiInfoHeader, edm::InLumi>(edm::InputTag("generator"))),
-        //iConfig.getUntrackedParameter<edm::InputTag>("lheSource", edm::InputTag("externalLHEProducer")))),
-    genEventToken_(consumes<GenEventInfoProduct>(edm::InputTag("generator")))
-        //iConfig.getUntrackedParameter<edm::InputTag>("lheSource", edm::InputTag("externalLHEProducer"))))
+    genLumiInfoToken_(consumes<GenLumiInfoHeader, edm::InLumi>(iConfig.getParameter<edm::InputTag>("genInfo"))),
+    genEventToken_(consumes<GenEventInfoProduct>(iConfig.getParameter<edm::InputTag>("genInfo")))
 {
   produces<GenWeightProduct>();
   produces<GenWeightInfoProduct, edm::Transition::BeginLuminosityBlock>();
@@ -63,14 +61,12 @@ void
 GenWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
   edm::Handle<GenEventInfoProduct> genEventInfo;
   iEvent.getByToken(genEventToken_, genEventInfo);
-  // Read weights from LHEEventProduct
-  auto weightProduct = weightHelper_.weightProduct(genEventInfo->weights());
+
+  float centralWeight = genEventInfo->weights().size() > 0 ? genEventInfo->weights().at(0) : 1.;
+  auto weightProduct = weightHelper_.weightProduct(genEventInfo->weights(), centralWeight);
   iEvent.put(std::move(weightProduct));
 }
 
-//void 
-//GenWeightProductProducer::endLuminosityBlockProduce(edm::LuminosityBlock& iLumi, edm::EventSetup const& iSetup) {}
-
 void
 GenWeightProductProducer::beginLuminosityBlockProduce(edm::LuminosityBlock& iLumi, edm::EventSetup const& iSetup) {
     if (weightNames_.size() == 0) {

GeneratorInterface/Core/plugins/LHEWeightProductProducer.cc

@@ -70,7 +70,7 @@ LHEWeightProductProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe
   edm::Handle<LHEEventProduct> lheEventInfo;
   iEvent.getByToken(lheEventToken_, lheEventInfo);
   // Read weights from LHEEventProduct
-  auto weightProduct = weightHelper_.weightProduct(lheEventInfo->weights());
+  auto weightProduct = weightHelper_.weightProduct(lheEventInfo->weights(), lheEventInfo->originalXWGTUP());
   iEvent.put(std::move(weightProduct));
 }
 

GeneratorInterface/Core/src/LHEWeightHelper.cc

@@ -136,6 +136,8 @@ namespace gen {
                 else {
                     weightGroups_.push_back(std::make_unique<gen::UnknownWeightGroupInfo>(name));
                 }
+                auto& group = weightGroups_.back();
+                group.setDescription("Test description");
             }
             else if (std::regex_match(headerLine, std::regex(".*<weight.*>.*\n*"))) {
                 std::string fullTag = headerLine;

GeneratorInterface/Core/src/WeightHelper.cc

@@ -85,8 +85,8 @@ namespace gen {
     }
 
     // TODO: Could probably recycle this code better
-    std::unique_ptr<GenWeightProduct> WeightHelper::weightProduct(std::vector<double> weights) {
-        auto weightProduct = std::make_unique<GenWeightProduct>();
+    std::unique_ptr<GenWeightProduct> WeightHelper::weightProduct(std::vector<double> weights, float w0) {
+        auto weightProduct = std::make_unique<GenWeightProduct>(w0);
         weightProduct->setNumWeightSets(weightGroups_.size());
         int weightGroupIndex = 0;
         for (unsigned int i = 0; i < weights.size(); i++) {
@@ -95,8 +95,8 @@ namespace gen {
         return std::move(weightProduct);
     }
 
-    std::unique_ptr<GenWeightProduct> WeightHelper::weightProduct(std::vector<gen::WeightsInfo> weights) {
-        auto weightProduct = std::make_unique<GenWeightProduct>();
+    std::unique_ptr<GenWeightProduct> WeightHelper::weightProduct(std::vector<gen::WeightsInfo> weights, float w0) {
+        auto weightProduct = std::make_unique<GenWeightProduct>(w0);
         weightProduct->setNumWeightSets(weightGroups_.size());
         int weightGroupIndex = 0;
         int i = 0;

GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc

@@ -197,7 +197,8 @@ void ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe
                             product.get(), _1));
 
   // Should also zero out the weights in the GenInfoProduct
-  auto weightProduct = weightHelper_.weightProduct(partonLevel->weights());
+  auto weightProduct = weightHelper_.weightProduct(partonLevel->weights(), 
+                                            partonLevel->originalXWGTUP());
   iEvent.put(std::move(weightProduct));
 
   product->setScales(partonLevel->scales());

PhysicsTools/NanoAOD/plugins/LHEWeightsTableProducer.cc

@@ -21,13 +21,13 @@
 #include <tinyxml2.h>
 
 namespace {
-  typedef std::pair<std::vector<gen::WeightGroupData>, std::vector<gen::WeightGroupData>> WeightGroupsToStore;
+  typedef std::vector<gen::WeightGroupData> WeightGroupDataContainer;
+  typedef std::array<std::vector<gen::WeightGroupData>, 2> WeightGroupsToStore;
 }  // namespace
 
 class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBlockCache<WeightGroupsToStore>> {
 public:
   LHEWeightsTableProducer(edm::ParameterSet const& params) :
-        lheToken_(consumes<LHEEventProduct>(params.getParameter<edm::InputTag>("lheInfo"))),
         lheWeightToken_(consumes<GenWeightProduct>(params.getParameter<edm::InputTag>("lheWeights"))),
         lheWeightInfoToken_(consumes<GenWeightInfoProduct, edm::InLumi>(params.getParameter<edm::InputTag>("lheWeights"))),
         genWeightToken_(consumes<GenWeightProduct>(params.getParameter<edm::InputTag>("genWeights"))),
@@ -44,8 +44,6 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
 
   void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
     // tables for LHE weights, may not be filled
-    auto const& lheInfo = iEvent.get(lheToken_);
-
     edm::Handle<GenWeightProduct> lheWeightHandle;
     iEvent.getByToken(lheWeightToken_, lheWeightHandle);
     const GenWeightProduct* lheWeightProduct = lheWeightHandle.product();
@@ -57,17 +55,16 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
     WeightsContainer genWeights = genWeightProduct->weights();
 
     auto lheWeightTables = std::make_unique<std::vector<nanoaod::FlatTable>>();
-    double w0 = lheInfo.originalXWGTUP();
     auto const& weightInfos = *luminosityBlockCache(iEvent.getLuminosityBlock().index());
 
-    addWeightGroupToTable(lheWeightTables, "LHE", weightInfos.first, lheWeights, w0);
-    addWeightGroupToTable(lheWeightTables, "GEN", weightInfos.second, genWeights, w0);
+    addWeightGroupToTable(lheWeightTables, "LHE", weightInfos.at(inLHE), lheWeights);
+    addWeightGroupToTable(lheWeightTables, "Gen", weightInfos.at(inGen), genWeights);
 
     iEvent.put(std::move(lheWeightTables));
   }
 
   void addWeightGroupToTable(std::unique_ptr<std::vector<nanoaod::FlatTable>>& lheWeightTables, const char* typeName,
-          const std::vector<gen::WeightGroupData>& weightInfos, WeightsContainer& lheWeights, double w0) const {
+          const WeightGroupDataContainer& weightInfos, WeightsContainer& allWeights) const {
     size_t typeCount = 0;
     gen::WeightType previousType = gen::WeightType::kUnknownWeights;
 
@@ -76,9 +73,23 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
         gen::WeightType weightType = groupInfo.group->weightType();
         if (previousType != weightType)
             typeCount = 0;
+
         std::string name = weightTypeNames_.at(weightType);
-        auto weights = weightType != gen::WeightType::kScaleWeights ? normalizedWeights(lheWeights, groupInfo, w0) :
-                            orderedScaleWeights(lheWeights, groupInfo, w0);
+        std::string label = groupInfo.group->name();
+
+        auto& weights = allWeights.at(groupInfo.index);
+        label.append("; ");
+        if (weightType == gen::WeightType::kScaleWeights && groupInfo.group->isWellFormed() 
+                && groupInfo.group->nIdsContained() < 10) {
+            weights = orderedScaleWeights(weights, 
+                            dynamic_cast<const gen::ScaleWeightGroupInfo*>(groupInfo.group));
+            label.append("[1] is mur=0.5 muf=1; [2] is mur=0.5 muf=2; [3] is mur=1 muf=0.5 ;"
+                        " [4] is mur=1 muf=1; [5] is mur=1 muf=2; [6] is mur=2 muf=0.5;"
+                        " [7] is mur=2 muf=1 ; [8] is mur=2 muf=2)");
+        }
+        else
+            label.append(groupInfo.group->description());
+
         entryName.append(name);
         entryName.append("Weight");
         if (typeCount > 0) {
@@ -88,7 +99,7 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
 
         lheWeightTables->emplace_back(weights.size(), entryName, false);
         lheWeightTables->back().addColumn<float>(
-            "", weights, groupInfo.group->name(), nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
+            "", weights, label, nanoaod::FlatTable::FloatColumn, lheWeightPrecision_);
 
         typeCount++;
         previousType = weightType;
@@ -112,33 +123,25 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
         storePerType[weightgroups_.at(i)] = maxGroupsPerType_.at(i);
 
     WeightGroupsToStore weightsToStore;
-    for (auto weightType : gen::allGenWeightTypes) {
+    for (auto weightType : gen::allWeightTypes) {
         auto lheWeights = weightDataPerType(lheWeightInfoHandle, weightType, storePerType[weightType]);
-        weightsToStore.first.insert(weightsToStore.first.end(), lheWeights.begin(), lheWeights.end());
+        weightsToStore.at(inLHE).insert(weightsToStore.at(inLHE).end(), lheWeights.begin(), lheWeights.end());
 
         auto genWeights = weightDataPerType(genWeightInfoHandle, weightType, storePerType[weightType]);
-        weightsToStore.second.insert(weightsToStore.second.end(), genWeights.begin(), genWeights.end());
+        weightsToStore.at(inGen).insert(weightsToStore.at(inGen).end(), genWeights.begin(), genWeights.end());
     }
 
     return std::make_shared<WeightGroupsToStore>(weightsToStore);
   }
 
-  std::vector<gen::WeightGroupData> weightDataPerType(edm::Handle<GenWeightInfoProduct>& weightsHandle, 
+  WeightGroupDataContainer weightDataPerType(edm::Handle<GenWeightInfoProduct>& weightsHandle, 
                                                 gen::WeightType weightType, int& maxStore) const {
-      std::vector<gen::WeightGroupData> group;
-      if (weightType == gen::WeightType::kPdfWeights) {
-        for (auto lhaid : pdfIds_) {
-            if (auto pdfGroup = weightsHandle->pdfGroupWithIndexByLHAID(lhaid)) {
-                group.push_back(pdfGroup.value());
-                maxStore--;
-                if (maxStore == 0)
-                    break;
-            }
-        }
-        return group;
+      WeightGroupDataContainer group;
+      if (weightType == gen::WeightType::kPdfWeights && pdfIds_.size() > 0) {
+          group = weightsHandle->pdfGroupsWithIndicesByLHAIDs(pdfIds_);
       }
-
-      group = weightsHandle->weightGroupsAndIndicesByType(weightType);
+      else
+          group = weightsHandle->weightGroupsAndIndicesByType(weightType);
 
       if (maxStore < 0 || static_cast<int>(group.size()) <= maxStore) {
           // Modify size in case one type of weight is present in multiple products
@@ -149,42 +152,27 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
   }
 
 
-  std::vector<float> normalizedWeights(WeightsContainer& lheWeights, const gen::WeightGroupData& meGroupInfo, double w0) const {
-    std::vector<float> normalizedWeights;
-    for (const auto& weight : lheWeights.at(meGroupInfo.index))
-        normalizedWeights.push_back(weight/w0);
-    return normalizedWeights;
-  }
+  std::vector<double> orderedScaleWeights(const std::vector<double>& scaleWeights, const gen::ScaleWeightGroupInfo* scaleGroup) const {
+
+    std::vector<double> weights;
+    weights.push_back(scaleWeights.at(scaleGroup->muR05muF05Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR05muF1Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR05muF2Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR1muF05Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->centralIndex()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR1muF2Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR2muF05Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR2muF1Index()));
+    weights.push_back(scaleWeights.at(scaleGroup->muR2muF2Index()));
 
-  std::vector<float> orderedScaleWeights(WeightsContainer& lheWeights, const gen::WeightGroupData& scaleGroupInfo, double w0) const {
-    const gen::ScaleWeightGroupInfo* scaleGroup = static_cast<const gen::ScaleWeightGroupInfo*>(scaleGroupInfo.group);
-    size_t scaleGroupIndex = scaleGroupInfo.index;
-
-    std::vector<float> normalizedWeights;
-    std::vector<double>& scaleWeights = lheWeights.at(scaleGroupIndex);
-
-    // nano ordering of mur=0.5 muf=0.5 ; [1] is mur=0.5 muf=1 ; [2] is mur=0.5 muf=2 ; [3] is mur=1 muf=0.5 ; 
-    // [4] is mur=1 muf=1 ; [5] is mur=1 muf=2 ; [6] is mur=2 muf=0.5 ; [7] is mur=2 muf=1 ; [8] is mur=2 muf=2 *
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR05muF05Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR05muF1Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR05muF2Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR1muF05Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->centralIndex())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR1muF2Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR2muF05Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR2muF1Index())/w0);
-    normalizedWeights.push_back(scaleWeights.at(scaleGroup->muR2muF2Index())/w0);
-
-    return normalizedWeights;
+    return weights;
   }
 
   // nothing to do here
   virtual void globalEndLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override {}
 
   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
     edm::ParameterSetDescription desc;
-    desc.add<edm::InputTag>("lheInfo", {"externalLHEProducer"})
-        ->setComment("tag(s) for the LHE information (LHEEventProduct and LHERunInfoProduct)");
     desc.add<edm::InputTag>("lheWeights");
     desc.add<edm::InputTag>("genWeights");
     desc.add<std::vector<std::string>>("weightgroups");
@@ -213,6 +201,7 @@ class LHEWeightsTableProducer : public edm::global::EDProducer<edm::LuminosityBl
 
   //std::unordered_map<std::string, int> weightGroupIndices_;
   int lheWeightPrecision_;
+  enum {inLHE, inGen};
 };
 
 #include "FWCore/Framework/interface/MakerMacros.h"

PhysicsTools/NanoAOD/python/nanogen_cff.py

@@ -6,18 +6,19 @@
 from PhysicsTools.NanoAOD.lheInfoTable_cfi import *
 from PhysicsTools.NanoAOD.genWeightsTable_cfi import *
 
-genWeights = cms.EDProducer("GenWeightProductProducer")
+genWeights = cms.EDProducer("GenWeightProductProducer",
+    genInfo = cms.InputTag("generator"))
 
 lheWeightsTable = cms.EDProducer(
     "LHEWeightsTableProducer",
-    lheInfo = cms.InputTag("externalLHEProducer"),
     lheWeights = cms.InputTag("externalLHEProducer"),
     genWeights = cms.InputTag("genWeights"),
     # Warning: you can use a full string, but only the first character is read.
     # 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']),
     maxGroupsPerType = cms.vint32([1, -1, 1, 2, 1]),
+    # If empty or not specified, all pdf weights are stored
     pdfIds = cms.vint32([91400, 306000, 260000]),
     lheWeightPrecision = cms.int32(14),
 )

SimDataFormats/GeneratorProducts/interface/GenWeightInfoProduct.h

@@ -36,7 +36,8 @@ class GenWeightInfoProduct {
         std::vector<gen::WeightGroupInfo*> weightGroupsByType(gen::WeightType type) const;
         std::vector<int> weightGroupIndicesByType(gen::WeightType type) const;
         std::vector<gen::WeightGroupData> weightGroupsAndIndicesByType(gen::WeightType type) const;
-        std::optional<gen::WeightGroupData> pdfGroupWithIndexByLHAID(size_t lhaid) const;
+        std::optional<gen::WeightGroupData> pdfGroupWithIndexByLHAID(int lhaid) const;
+        std::vector<gen::WeightGroupData> pdfGroupsWithIndicesByLHAIDs(const std::vector<int>& lhaids) const;
         void addWeightGroupInfo(gen::WeightGroupInfo* info);
         const int numberOfGroups() const { return weightGroupsInfo_.size(); }
 

SimDataFormats/GeneratorProducts/interface/GenWeightProduct.h

@@ -13,10 +13,12 @@ typedef std::vector<std::vector<double>> WeightsContainer;
 
 class GenWeightProduct {
     public:
-        GenWeightProduct() { weightsVector_ = {}; }
+        GenWeightProduct() { weightsVector_ = {}; centralWeight_ = 1.; }
+        GenWeightProduct(double w0) { weightsVector_ = {}; centralWeight_ = w0; }
         GenWeightProduct& operator=(GenWeightProduct&& other) {
-        weightsVector_ = std::move(other.weightsVector_);
-        return *this;
+            weightsVector_ = std::move(other.weightsVector_);
+            centralWeight_ = other.centralWeight_;
+            return *this;
         }
         ~GenWeightProduct() {}
 
@@ -31,12 +33,14 @@ class GenWeightProduct {
             if (static_cast<int>(weights.size()) <= weightNum) {
                 weights.resize(weightNum+1);
             }
-            weights[weightNum] = weight;
+            weights[weightNum] = weight/centralWeight_;
         }
         const WeightsContainer& weights() const { return weightsVector_; }
+        double centralWeight() const { return centralWeight_; }
 
     private:
         WeightsContainer weightsVector_;
+        double centralWeight_;
 };
 
 #endif // GeneratorEvent_LHEInterface_GenWeightProduct_h