diff --git a/converter/include/k4SimDelphes/DelphesEDM4HepConverter.h b/converter/include/k4SimDelphes/DelphesEDM4HepConverter.h index 69cdc6cd..4eb5066b 100644 --- a/converter/include/k4SimDelphes/DelphesEDM4HepConverter.h +++ b/converter/include/k4SimDelphes/DelphesEDM4HepConverter.h @@ -106,7 +106,7 @@ class DelphesEDM4HepConverter { createExternalRecoMCLinks(const std::unordered_map& 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 \ No newline at end of file diff --git a/converter/src/DelphesEDM4HepConverter.cc b/converter/src/DelphesEDM4HepConverter.cc index a53bf1e2..9aabdfaf 100644 --- a/converter/src/DelphesEDM4HepConverter.cc +++ b/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(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(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(lhefWeights->At(i)); + // append this weight to the event header + cand.addToWeights(weightEntry->Weight); + } + } } void DelphesEDM4HepConverter::processParticles(const TClonesArray* delphesCollection, std::string const& branch) { diff --git a/standalone/src/DelphesInputReader.h b/standalone/src/DelphesInputReader.h index 2a94ecd0..686c169a 100644 --- a/standalone/src/DelphesInputReader.h +++ b/standalone/src/DelphesInputReader.h @@ -24,6 +24,10 @@ class DelphesInputReader { TObjArray* stableParticleOutputArray, TObjArray* partonOutputArray) = 0; virtual TTree* converterTree() = 0; + virtual const std::vector& getWeightNames() const { + static const std::vector empty; + return empty; + } }; #endif diff --git a/standalone/src/DelphesMain.h b/standalone/src/DelphesMain.h index fde2c890..44f8e4b1 100644 --- a/standalone/src/DelphesMain.h +++ b/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>(edm4hep::labels::EventWeightsNames, weightNames); + weightNames_stored = true; + } + } podioWriter.writeFrame(frame, "events"); modularDelphes->Clear(); diff --git a/standalone/src/DelphesPythia8Common.h b/standalone/src/DelphesPythia8Common.h index b7912b7d..aec0f2e2 100644 --- a/standalone/src/DelphesPythia8Common.h +++ b/standalone/src/DelphesPythia8Common.h @@ -190,4 +190,36 @@ void fillPartons(int id, double pMax, double etaMax, Pythia8::Event& event, Pyth } } +std::vector getWeightNames(const std::string& lheFilePath) { + std::ifstream lheFile(lheFilePath); + std::vector 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("") != 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 diff --git a/standalone/src/DelphesPythia8Reader.h b/standalone/src/DelphesPythia8Reader.h index 02b35fbb..1ae3208b 100644 --- a/standalone/src/DelphesPythia8Reader.h +++ b/standalone/src/DelphesPythia8Reader.h @@ -36,6 +36,9 @@ class DelphesPythia8Reader : public DelphesInputReader { PrintXS(m_pythia.get()); } + const std::vector& 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 m_converterTree{nullptr}; ExRootTreeBranch *m_brancheEventLHEF = 0, *m_branchWeightLHEF = 0; + // arrays to store weight names + std::vector m_weightNames; TObjArray *m_stableParticleOutputArrayLHEF = 0, *m_allParticleOutputArrayLHEF = 0, *m_partonOutputArrayLHEF = 0; DelphesLHEFReader* reader = 0; Long64_t m_eventCounter{0}, m_errorCounter{0};