Add lhe numbering to allow get correct matching between lhe and herwig events

Dominic-Stafford · Dominic-Stafford · commit c7d31e677890 · 2023-08-28T18:48:26.000+02:00
Read lhe event numbers in, if present

Add HadroniserFilter to Herwig which correctly matches LHE numbers between CMSSW and herwig

Add option to Herwig input fragements to use LHE numbering

Change test examples to use new HadroniserFilter+ LHE numbering

code style

code style

Call addLHEnumbers from mergeLHE.py if asked to number events and not using the DefaultLHEMerger

Number events before merge (so numbers are also available to LHEReader)

Include evtnum in LHEEventProduct (necesary for multithreading)
diff --git a/Configuration/Generator/python/DYToLL01234Jets_5FS_TuneCH3_13TeV_madgraphMLM_herwig7_cff.py b/Configuration/Generator/python/DYToLL01234Jets_5FS_TuneCH3_13TeV_madgraphMLM_herwig7_cff.py
@@ -6,7 +6,7 @@
 from Configuration.Generator.Herwig7Settings.Herwig7MGMergingSettings_cfi import *
 
 
-generator = cms.EDFilter("Herwig7GeneratorFilter",
+generator = cms.EDFilter("Herwig7HadronizerFilter",
     herwig7CH3SettingsBlock,
     herwig7StableParticlesForDetectorBlock,
     herwig7MGMergingSettingsBlock,
diff --git a/Configuration/Generator/python/DYToLL012Jets_5FS_TuneCH3_13TeV_amcatnloFxFx_herwig7_cff.py b/Configuration/Generator/python/DYToLL012Jets_5FS_TuneCH3_13TeV_amcatnloFxFx_herwig7_cff.py
@@ -5,7 +5,7 @@
 from Configuration.Generator.Herwig7Settings.Herwig7MGMergingSettings_cfi import *
 
 
-generator = cms.EDFilter("Herwig7GeneratorFilter",
+generator = cms.EDFilter("Herwig7HadronizerFilter",
     herwig7CH3SettingsBlock,
     herwig7StableParticlesForDetectorBlock,
     herwig7MGMergingSettingsBlock,
@@ -40,7 +40,7 @@
     outputFile = cms.string('cmsgrid_final.lhe'),
     scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh'),
     generateConcurrently = cms.untracked.bool(True),
-    postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe')
+    postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-n', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe')
 )
 
 
diff --git a/Configuration/Generator/python/DYToll01234Jets_5f_LO_MLM_Madgraph_LHE_13TeV_cff.py b/Configuration/Generator/python/DYToll01234Jets_5f_LO_MLM_Madgraph_LHE_13TeV_cff.py
@@ -7,5 +7,5 @@
   args = cms.vstring('/cvmfs/cms.cern.ch/phys_generator/gridpacks/UL/13TeV/madgraph/V5_2.6.5/dyellell01234j_5f_LO_MLM_v2/DYJets_HT-incl_slc6_amd64_gcc630_CMSSW_9_3_16_tarball.tar.xz','false','slc6_amd64_gcc630','CMSSW_9_3_16'),
   nEvents = cms.untracked.uint32(10),
   generateConcurrently = cms.untracked.bool(True),
-  postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),
+  postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-n', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),
 )
diff --git a/Configuration/Generator/python/Herwig7Settings/Herwig7LHECommonSettings_cfi.py b/Configuration/Generator/python/Herwig7Settings/Herwig7LHECommonSettings_cfi.py
@@ -14,6 +14,7 @@
 
         # set the weight option (e.g. for MC@NLO)
         'set LesHouchesHandler:WeightOption VarNegWeight',
+        'set LesHouchesHandler:EventNumbering LHE',
 
         'set /Herwig/Generators/EventGenerator:EventHandler /Herwig/EventHandlers/LesHouchesHandler',
         'create ThePEG::Cuts /Herwig/Cuts/NoCuts',
diff --git a/Configuration/Generator/python/Herwig7Settings/Herwig7MGMergingSettings_cfi.py b/Configuration/Generator/python/Herwig7Settings/Herwig7MGMergingSettings_cfi.py
@@ -14,6 +14,7 @@
         'set LesHouchesHandler:HadronizationHandler /Herwig/Hadronization/ClusterHadHandler',
         'set LesHouchesHandler:DecayHandler /Herwig/Decays/DecayHandler',
         'set LesHouchesHandler:WeightOption VarNegWeight',
+        'set LesHouchesHandler:EventNumbering LHE',
         'set /Herwig/Generators/EventGenerator:EventHandler  /Herwig/EventHandlers/LesHouchesHandler',
         'create ThePEG::Cuts /Herwig/Cuts/NoCuts',
         'cd /Herwig/EventHandlers',
diff --git a/Configuration/Generator/python/TT_13TeV_Pow_Herwig7_cff.py b/Configuration/Generator/python/TT_13TeV_Pow_Herwig7_cff.py
@@ -7,7 +7,7 @@
 from Configuration.Generator.Herwig7Settings.Herwig7LHEPowhegSettings_cfi import *
 
 
-generator = cms.EDFilter("Herwig7GeneratorFilter",
+generator = cms.EDFilter("Herwig7HadronizerFilter",
     herwig7CH3SettingsBlock,
     herwig7StableParticlesForDetectorBlock,
     herwig7LHECommonSettingsBlock,
diff --git a/Configuration/Generator/python/TTbar_Pow_LHE_13TeV_cff.py b/Configuration/Generator/python/TTbar_Pow_LHE_13TeV_cff.py
@@ -7,7 +7,7 @@
     outputFile = cms.string('cmsgrid_final.lhe'),
     scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh'),
     generateConcurrently = cms.untracked.bool(True),
-    postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),
+    postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-n', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),
 )
 
 #Link to datacards:
diff --git a/GeneratorInterface/Herwig7Interface/plugins/Herwig7Hadronizer.cc b/GeneratorInterface/Herwig7Interface/plugins/Herwig7Hadronizer.cc
@@ -63,6 +63,7 @@ class Herwig7Hadronizer : public Herwig7Interface, public gen::BaseHadronizer {
   unsigned int eventsToPrint;
 
   ThePEG::EventPtr thepegEvent;
+  bool haveEvt = false;
 
   std::shared_ptr<lhef::LHEProxy> proxy_;
   const std::string handlerDirectory_;
@@ -101,9 +102,18 @@ bool Herwig7Hadronizer::initializeForInternalPartons() {
 }
 
 bool Herwig7Hadronizer::initializeForExternalPartons() {
-  edm::LogError("Herwig7 interface")
-      << "Read in of LHE files is not supported in this way. You can read them manually if necessary.";
-  return false;
+  if (currentLumiBlock == firstLumiBlock) {
+    std::ifstream runFile(runFileName + ".run");
+    if (runFile.fail())  //required for showering of LHE files
+    {
+      initRepository(paramSettings);
+    }
+    if (!initGenerator()) {
+      edm::LogInfo("Generator|Herwig7Hadronizer") << "No run step for Herwig chosen. Program will be aborted.";
+      exit(0);
+    }
+  }
+  return true;
 }
 
 bool Herwig7Hadronizer::declareStableParticles(const std::vector<int>& pdgIds) { return false; }
@@ -141,8 +151,41 @@ bool Herwig7Hadronizer::generatePartonsAndHadronize() {
 }
 
 bool Herwig7Hadronizer::hadronize() {
-  edm::LogError("Herwig7 interface")
-      << "Read in of LHE files is not supported in this way. You can read them manually if necessary.";
+  if (!haveEvt) {
+    try {
+      thepegEvent = eg_->shoot();
+      haveEvt = true;
+    } catch (std::exception& exc) {
+      edm::LogWarning("Generator|Herwig7Hadronizer")
+          << "EGPtr::shoot() thrown an exception, event skipped: " << exc.what();
+      return false;
+    }
+  }
+  int evtnum = lheEvent()->evtnum();
+  if (evtnum == -1) {
+    edm::LogError("Generator|Herwig7Hadronizer")
+        << "Event number not set in lhe file, needed for correctly aligning Herwig and LHE events!";
+    return false;
+  }
+  if (thepegEvent->number() < evtnum) {
+    edm::LogError("Herwig7 interface") << "Herwig does not seem to be generating events in order, did you set "
+                                          "/Herwig/EventHandlers/FxFxLHReader:AllowedToReOpen Yes?";
+    return false;
+  } else if (thepegEvent->number() == evtnum) {
+    haveEvt = false;
+    if (!thepegEvent) {
+      edm::LogWarning("Generator|Herwig7Hadronizer") << "thepegEvent not initialized";
+      return false;
+    }
+
+    event() = convert(thepegEvent);
+    if (!event().get()) {
+      edm::LogWarning("Generator|Herwig7Hadronizer") << "genEvent not initialized";
+      return false;
+    }
+    return true;
+  }
+  edm::LogWarning("Generator|Herwig7Hadronizer") << "Event " << evtnum << " not generated (likely skipped in merging)";
   return false;
 }
 
diff --git a/GeneratorInterface/LHEInterface/interface/LHEEvent.h b/GeneratorInterface/LHEInterface/interface/LHEEvent.h
@@ -52,9 +52,11 @@ namespace lhef {
 
     int npLO() const { return npLO_; }
     int npNLO() const { return npNLO_; }
+    int evtnum() const { return evtnum_; }
 
     void setNpLO(int n) { npLO_ = n; }
     void setNpNLO(int n) { npNLO_ = n; }
+    void setEvtNum(int n) { evtnum_ = n; }
 
     void addComment(const std::string &line) { comments.push_back(line); }
 
@@ -93,6 +95,7 @@ namespace lhef {
     std::vector<float> scales_;  //scale value used to exclude EWK-produced partons from matching
     int npLO_;                   //number of partons for LO process (used to steer matching/merging)
     int npNLO_;                  //number of partons for NLO process (used to steer matching/merging)
+    int evtnum_;  //The number of the event (needed to ensure the correct LHE events are saved for MG +Herwig)
   };
 
 }  // namespace lhef
diff --git a/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc b/GeneratorInterface/LHEInterface/plugins/ExternalLHEProducer.cc
@@ -192,6 +192,7 @@ void ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe
                 partonLevel_->weights().end(),
                 std::bind(&LHEEventProduct::addWeight, product.get(), std::placeholders::_1));
   product->setScales(partonLevel_->scales());
+  product->setEvtNum(partonLevel_->evtnum());
   if (nPartonMapping_.empty()) {
     product->setNpLO(partonLevel_->npLO());
     product->setNpNLO(partonLevel_->npNLO());
diff --git a/GeneratorInterface/LHEInterface/scripts/addLHEnumbers.py b/GeneratorInterface/LHEInterface/scripts/addLHEnumbers.py
@@ -0,0 +1,47 @@
+#!/usr/bin/env python3
+
+from __future__ import print_function
+import logging
+import argparse
+import sys
+import os
+import re
+
+
+def number_events(input_file, output_file=None, offset=0):
+    if output_file is None:
+        output_file = input_file
+    if not os.path.exists(os.path.dirname(os.path.realpath(output_file))):
+        os.makedirs(os.path.dirname(os.path.realpath(output_file)))
+
+    nevent = offset
+    with open('tmp.txt', 'w') as fw:
+        with open(input_file, 'r') as ftmp:
+            for line in ftmp:
+                if re.search('\s*</event>', line):
+                    nevent += 1
+                    fw.write('<event_num num="' + str(nevent) +  '"> ' + str(nevent) + '</event_num>\n')
+                fw.write(line)
+    if output_file is not None:
+        os.rename("tmp.txt", output_file)
+    else:
+        os.rename("tmp.txt", input_file)
+    return nevent
+
+
+if __name__=="__main__":
+
+    parser = argparse.ArgumentParser(
+        description="Add numbers to lhe")
+    parser.add_argument("input_file", type=str,
+                        help="Input LHE file path.")
+    parser.add_argument("-o", "--output-file", default=None, type=str,
+                        help="Output LHE file path. If not specified, output to input file")
+    args = parser.parse_args()
+
+    logging.info('>>> launch addLHEnumbers.py in %s' % os.path.abspath(os.getcwd()))
+
+    logging.info('>>> Input file: [%s]' % args.input_file)
+    logging.info('>>> Write to output: %s ' % args.output_file)
+
+    number_events(args.input_file, args.output_file)
diff --git a/GeneratorInterface/LHEInterface/scripts/mergeLHE.py b/GeneratorInterface/LHEInterface/scripts/mergeLHE.py
@@ -9,6 +9,9 @@
 import os
 import re
 
+import addLHEnumbers
+
+
 class BaseLHEMerger(object):
     """Base class of the LHE merge schemes"""
 
@@ -62,7 +65,7 @@ def check_header_compatibility(self):
             % ', '.join([str(len(lines)) for lines in self._header_lines]))
         assert all([
             len(self._header_lines[0]) == len(lines) for lines in self._header_lines]
-            ), inconsistent_error_info % "line number not matches"
+            ), inconsistent_error_info % "line number does not match"
         inconsistent_lines_set = [set() for _ in self._header_lines]
         for line_zip in zip(*self._header_lines):
             if any([k in line_zip[0] for k in allow_diff_keys]):
@@ -253,7 +256,7 @@ def merge(self):
                     sign = lambda x: -1 if x < 0 else 1
                     for line in ftmp:
                         event_line += 1
-                        if re.search('\s*<event.*>', line):
+                        if re.search('\s*<event.*>', line) and not re.search('\s*<event_num.*>', line):
                             event_line = 0
                         if event_line == 1:
                             # modify the XWGTUP appeared in the first line of the
@@ -366,6 +369,8 @@ def main(argv = None):
                         help=("Bypass the compatibility check for the headers. If true, the header and init block "
                              "will be just a duplicate from the first input file, and events are concatenated without "
                              "modification."))
+    parser.add_argument("-n", "--number-events", action='store_true',
+                        help=("Add a tag to number each lhe event. Needed for Herwig to find correct lhe events"))
     parser.add_argument("--debug", action='store_true',
                         help="Use the debug mode.")
     args = parser.parse_args(argv)
@@ -392,6 +397,11 @@ def main(argv = None):
     if not os.path.exists(os.path.dirname(os.path.realpath(args.output_file))):
         os.makedirs(os.path.dirname(os.path.realpath(args.output_file)))
 
+    if args.number_events:
+        offset = 0
+        for input_file in input_files:
+            offset += addLHEnumbers.number_events(input_file, offset=offset)
+
     # Check arguments
     assert len(input_files) > 0, 'Input LHE files should be more than 0.'
     if len(input_files) == 1:
@@ -408,7 +418,8 @@ def main(argv = None):
     elif args.force_cpp_merger:
         lhe_merger = ExternalCppLHEMerger(input_files, args.output_file)
     else:
-        lhe_merger = DefaultLHEMerger(input_files, args.output_file, bypass_check=args.bypass_check)
+        lhe_merger = DefaultLHEMerger(
+            input_files, args.output_file, bypass_check=args.bypass_check)
 
     # Do merging
     lhe_merger.merge()
diff --git a/GeneratorInterface/LHEInterface/src/LHEEvent.cc b/GeneratorInterface/LHEInterface/src/LHEEvent.cc
@@ -40,7 +40,8 @@ namespace lhef {
         counted(false),
         readAttemptCounter(0),
         npLO_(-99),
-        npNLO_(-99)
+        npNLO_(-99),
+        evtnum_(-1)
 
   {
     hepeup.NUP = 0;
@@ -106,7 +107,7 @@ namespace lhef {
   }
 
   LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo, const HEPEUP &hepeup)
-      : runInfo(runInfo), hepeup(hepeup), counted(false), readAttemptCounter(0), npLO_(-99), npNLO_(-99) {}
+      : runInfo(runInfo), hepeup(hepeup), counted(false), readAttemptCounter(0), npLO_(-99), npNLO_(-99), evtnum_(-1) {}
 
   LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
                      const HEPEUP &hepeup,
@@ -119,7 +120,8 @@ namespace lhef {
         counted(false),
         readAttemptCounter(0),
         npLO_(-99),
-        npNLO_(-99) {}
+        npNLO_(-99),
+        evtnum_(-1) {}
 
   LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo, const LHEEventProduct &product)
       : runInfo(runInfo),
@@ -132,7 +134,8 @@ namespace lhef {
         originalXWGTUP_(product.originalXWGTUP()),
         scales_(product.scales()),
         npLO_(product.npLO()),
-        npNLO_(product.npNLO()) {}
+        npNLO_(product.npNLO()),
+        evtnum_(product.evtnum()) {}
 
   LHEEvent::~LHEEvent() {}
 
diff --git a/GeneratorInterface/LHEInterface/src/LHEReader.cc b/GeneratorInterface/LHEInterface/src/LHEReader.cc
@@ -141,6 +141,7 @@ namespace lhef {
     int npLO;
     int npNLO;
     std::vector<float> scales;
+    int evtnum = -1;
   };
 
   static void attributesToDom(DOMElement *dom, const Attributes &attributes) {
@@ -215,6 +216,9 @@ namespace lhef {
 
           scales.push_back(scaleval);
         }
+      } else if (name == "event_num") {
+        const char *evtnumstr = XMLSimpleStr(attributes.getValue(XMLString::transcode("num")));
+        sscanf(evtnumstr, "%d", &evtnum);
       }
       xmlEventNodes.push_back(elem);
       return;
@@ -526,6 +530,8 @@ namespace lhef {
           }
           lheevent->setNpLO(handler->npLO);
           lheevent->setNpNLO(handler->npNLO);
+          lheevent->setEvtNum(handler->evtnum);
+          handler->evtnum = -1;
           //fill scales
           if (!handler->scales.empty()) {
             lheevent->setScales(handler->scales);
diff --git a/SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h b/SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h
@@ -39,9 +39,11 @@ class LHEEventProduct {
 
   int npLO() const { return npLO_; }
   int npNLO() const { return npNLO_; }
+  int evtnum() const { return evtnum_; }
 
   void setNpLO(int n) { npLO_ = n; }
   void setNpNLO(int n) { npNLO_ = n; }
+  void setEvtNum(int n) { evtnum_ = n; }
 
   const lhef::HEPEUP &hepeup() const { return hepeup_; }
   const PDF *pdf() const { return pdf_.get(); }
@@ -108,6 +110,7 @@ class LHEEventProduct {
   std::vector<float> scales_;  //scale value used to exclude EWK-produced partons from matching
   int npLO_;                   //number of partons for LO process (used to steer matching/merging)
   int npNLO_;                  //number of partons for NLO process (used to steer matching/merging)
+  int evtnum_;  //The number of the event (needed to ensure the correct LHE events are saved for MG +Herwig)
 };
 
 #endif  // GeneratorEvent_LHEInterface_LHEEventProduct_h
diff --git a/SimDataFormats/GeneratorProducts/src/classes_def.xml b/SimDataFormats/GeneratorProducts/src/classes_def.xml
@@ -209,12 +209,13 @@
 	<class name="LHERunInfoProduct" ClassVersion="10">
   <version ClassVersion="10" checksum="2098560362"/>
  </class>
-	<class name="LHEEventProduct" ClassVersion="14">
+	<class name="LHEEventProduct" ClassVersion="15">
   <version ClassVersion="10" checksum="1026978494"/>
   <version ClassVersion="11" checksum="619154414"/>
   <version ClassVersion="12" checksum="259352862"/>
   <version ClassVersion="13" checksum="3927731647"/>
   <version ClassVersion="14" checksum="905719452"/>
+  <version ClassVersion="15" checksum="898288635"/>
  </class>
  <ioread sourceClass = "LHEEventProduct" version="[13]" targetClass="LHEEventProduct" source = "std::auto_ptr<gen::PdfInfo> pdf_;" target="pdf_">
    <![CDATA[pdf_.reset(onfile.pdf_.get()); onfile.pdf_.release();]]>

Original file line number	Diff line number	Diff line change
`@@ -7,5 +7,5 @@`
`7`	`7`	`args = cms.vstring('/cvmfs/cms.cern.ch/phys_generator/gridpacks/UL/13TeV/madgraph/V5_2.6.5/dyellell01234j_5f_LO_MLM_v2/DYJets_HT-incl_slc6_amd64_gcc630_CMSSW_9_3_16_tarball.tar.xz','false','slc6_amd64_gcc630','CMSSW_9_3_16'),`
`8`	`8`	`nEvents = cms.untracked.uint32(10),`
`9`	`9`	`generateConcurrently = cms.untracked.bool(True),`
`10`		`- postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),`
	`10`	`+ postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-n', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),`
`11`	`11`	`)`
Original file line number	Diff line number	Diff line change
`@@ -7,7 +7,7 @@`
`7`	`7`	`outputFile = cms.string('cmsgrid_final.lhe'),`
`8`	`8`	`scriptName = cms.FileInPath('GeneratorInterface/LHEInterface/data/run_generic_tarball_cvmfs.sh'),`
`9`	`9`	`generateConcurrently = cms.untracked.bool(True),`
`10`		`- postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),`
	`10`	`+ postGenerationCommand = cms.untracked.vstring('mergeLHE.py', '-n', '-i', 'thread*/cmsgrid_final.lhe', '-o', 'cmsgrid_final.lhe'),`
`11`	`11`	`)`
`12`	`12`
`13`	`13`	`#Link to datacards:`