Herwig lhe matching fix (including multithreading fixes)

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,

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')
 )
 
 

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'),
 )

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',

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',

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,

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:

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;
 }
 

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

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());

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)

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()

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() {}
 

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);

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

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();]]>