LHE MadGraph reweighting implementation into DelphesPythia8_EDM4HEP.

converter/include/k4SimDelphes/DelphesEDM4HepConverter.h

@@ -106,7 +106,7 @@ class DelphesEDM4HepConverter {
   createExternalRecoMCLinks(const std::unordered_map<UInt_t, edm4hep::MCParticle>& mc_map);
 
 private:
-  void createEventHeader(const HepMCEvent* delphesEvent);
+  void createEventHeader(const HepMCEvent* delphesEvent, TClonesArray* lhefWeights);
 
   void processParticles(const TClonesArray* delphesCollection, std::string const& branch);
   void processTracks(const TClonesArray* delphesCollection, std::string const& branch);
@@ -191,4 +191,4 @@ CollectionT* DelphesEDM4HepConverter::createCollection(std::string const& name,
 
 } // namespace k4SimDelphes
 
-#endif
+#endif

converter/src/DelphesEDM4HepConverter.cc

@@ -132,13 +132,21 @@ void DelphesEDM4HepConverter::process(TTree* delphesTree) {
   // Make sure the shared collections are present
   registerGlobalCollections();
 
-  // filling the event header
+  // Retrieve event header from the "Event" branch
   auto* eventBranch = delphesTree->GetBranch("Event");
-
   if (eventBranch) {
     auto* delphesEvents = *(TClonesArray**)eventBranch->GetAddress();
     auto* delphesEvent = static_cast<HepMCEvent*>(delphesEvents->At(0));
-    createEventHeader(delphesEvent);
+
+    // Retrieve the weights array from the "WeightLHEF" branch
+    auto* weightBranch = delphesTree->GetBranch("WeightLHEF");
+    TClonesArray* lhefWeights = nullptr;
+    if (weightBranch) {
+      lhefWeights = *(TClonesArray**)weightBranch->GetAddress();
+    }
+
+    // call event header fn to store weights, doesn't store (in fn) if lhefweights is empty.
+    createEventHeader(delphesEvent, lhefWeights);
   }
 
   for (const auto& branch : m_branches) {
@@ -165,12 +173,23 @@ void DelphesEDM4HepConverter::process(TTree* delphesTree) {
 }
 
 // convert the eventHeader with metaData
-void DelphesEDM4HepConverter::createEventHeader(const HepMCEvent* delphesEvent) {
+void DelphesEDM4HepConverter::createEventHeader(const HepMCEvent* delphesEvent, TClonesArray* lhefWeights) {
   auto* collection = createCollection<edm4hep::EventHeaderCollection>(EVENTHEADER_NAME);
   auto cand = collection->create();
 
+  // Set basic event properties
   cand.setWeight(delphesEvent->Weight);
   cand.setEventNumber(delphesEvent->Number);
+
+  // if we have weights being read from LHE file, store them in _EventHeader_weights
+  if (lhefWeights) {
+    for (Int_t i = 0; i < lhefWeights->GetEntries(); ++i) {
+      // get entry in vector
+      LHEFWeight* weightEntry = static_cast<LHEFWeight*>(lhefWeights->At(i));
+      // append this weight to the event header
+      cand.addToWeights(weightEntry->Weight);
+    }
+  }
 }
 
 void DelphesEDM4HepConverter::processParticles(const TClonesArray* delphesCollection, std::string const& branch) {

standalone/src/DelphesInputReader.h

@@ -24,6 +24,10 @@ class DelphesInputReader {
                          TObjArray* stableParticleOutputArray, TObjArray* partonOutputArray) = 0;
 
   virtual TTree* converterTree() = 0;
+  virtual const std::vector<std::string>& getWeightNames() const {
+    static const std::vector<std::string> empty;
+    return empty;
+  }
 };
 
 #endif

standalone/src/DelphesMain.h

@@ -74,6 +74,16 @@ int doit(int argc, char* argv[], DelphesInputReader& inputReader) {
       for (auto& [name, coll] : edm4hepConverter.getCollections()) {
         frame.put(std::move(coll), name);
       }
+      // write out LHEF weights once (if available)
+      static bool weightNames_stored = false;
+      if (!weightNames_stored) {
+        const auto& weightNames = inputReader.getWeightNames();
+        if (!weightNames.empty()) {
+          std::cout << "Writing LHEF weight ID vector of strings with size " << weightNames.size() << std::endl;
+          frame.putParameter<std::vector<std::string>>(edm4hep::labels::EventWeightsNames, weightNames);
+          weightNames_stored = true;
+        }
+      }
       podioWriter.writeFrame(frame, "events");
 
       modularDelphes->Clear();

standalone/src/DelphesPythia8Common.h

@@ -190,4 +190,36 @@ void fillPartons(int id, double pMax, double etaMax, Pythia8::Event& event, Pyth
   }
 }
 
+std::vector<std::string> getWeightNames(const std::string& lheFilePath) {
+  std::ifstream lheFile(lheFilePath);
+  std::vector<std::string> weightIDs;
+
+  if (!lheFile.is_open()) {
+    std::cerr << "Failed to open LHE file: " << lheFilePath << std::endl;
+    return weightIDs;
+  }
+
+  std::string line;
+  while (std::getline(lheFile, line)) {
+    // Stop after the initrwgt block
+    if (line.find("</initrwgt>") != std::string::npos) {
+      break;
+    }
+
+    size_t idPos = line.find("weight id='");
+    if (idPos != std::string::npos) {
+      size_t start = line.find("'", idPos);
+      size_t end = line.find("'", start + 1);
+      if (start != std::string::npos && end != std::string::npos) {
+        std::string id = line.substr(start + 1, end - start - 1);
+        weightIDs.push_back(id);
+        std::cout << "Found weight ID: " << id << std::endl;
+      }
+    }
+  }
+
+  lheFile.close();
+  return weightIDs;
+}
+
 #endif

standalone/src/DelphesPythia8Reader.h

@@ -36,6 +36,9 @@ class DelphesPythia8Reader : public DelphesInputReader {
       PrintXS(m_pythia.get());
   }
 
+  const std::vector<std::string>& getWeightNames() const {
+    return m_weightNames;
+  }
   std::string init(Delphes* modularDelphes, int argc, char* argv[]) override {
     if (argc != 5) {
       return "";
@@ -105,6 +108,16 @@ class DelphesPythia8Reader : public DelphesInputReader {
         m_brancheEventLHEF = m_treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
         m_branchWeightLHEF = m_treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
 
+        if (!m_branchWeightLHEF) {
+          std::cerr << "Error: m_branchWeightLHEF failed to initialise." << std::endl;
+        } else {
+          m_weightNames = ::getWeightNames(m_pythia->word("Beams:LHEF"));
+          std::cout << "Found " << m_weightNames.size() << " weight names in initrwgt:\n";
+          for (const auto& id : m_weightNames) {
+            std::cout << "  " << id << "\n";
+          }
+        }
+
         m_allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
         m_stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
         m_partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
@@ -146,6 +159,16 @@ class DelphesPythia8Reader : public DelphesInputReader {
       }
     }
 
+    // reading weights
+    if (reader) {
+      if (m_branchWeightLHEF) {
+        m_branchWeightLHEF->Clear();
+        reader->AnalyzeWeight(m_branchWeightLHEF);
+      } else {
+        std::cerr << "Error: m_branchWeightLHEF is null." << std::endl;
+      }
+    }
+
     if (!m_pythia->next()) {
       // If failure because reached end of file then exit event loop
       if (m_pythia->info.atEndOfFile()) {
@@ -166,6 +189,9 @@ class DelphesPythia8Reader : public DelphesInputReader {
     m_procStopWatch.Start();
     ConvertInput(m_eventCounter, m_pythia.get(), m_branchEvent.get(), factory, allParticleOutputArray,
                  stableParticleOutputArray, partonOutputArray, &m_readStopWatch, &m_procStopWatch);
+
+    // fill branches not read in by pythia (LHEF reweighting branch)
+    m_treeWriter->Fill();
     ++m_eventCounter;
     return true;
   };
@@ -185,6 +211,8 @@ class DelphesPythia8Reader : public DelphesInputReader {
   std::unique_ptr<TTree> m_converterTree{nullptr};
 
   ExRootTreeBranch *m_brancheEventLHEF = 0, *m_branchWeightLHEF = 0;
+  // arrays to store weight names
+  std::vector<std::string> m_weightNames;
   TObjArray *m_stableParticleOutputArrayLHEF = 0, *m_allParticleOutputArrayLHEF = 0, *m_partonOutputArrayLHEF = 0;
   DelphesLHEFReader* reader = 0;
   Long64_t m_eventCounter{0}, m_errorCounter{0};