Merge pull request #45241 from Dominic-Stafford/hw_lhe_match_fix_13_0

Herwig LHE matching fix and Matchbox relval for 13_0

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/DY_TuneCH3_13TeV_herwig_madgraph_matchbox_cff.py

@@ -0,0 +1,47 @@
+import FWCore.ParameterSet.Config as cms
+
+from Configuration.Generator.Herwig7Settings.Herwig7StableParticlesForDetector_cfi import *
+from Configuration.Generator.Herwig7Settings.Herwig7CH3TuneSettings_cfi import *
+
+generator = cms.EDFilter("Herwig7GeneratorFilter",
+    herwig7StableParticlesForDetectorBlock,
+    herwig7CH3SettingsBlock,
+    run = cms.string('InterfaceMatchboxTest'),
+    dumpConfig = cms.untracked.string('HerwigConfig.in'),
+    repository = cms.string('${HERWIGPATH}/HerwigDefaults.rpo'),
+    dataLocation = cms.string('${HERWIGPATH:-6}'),
+    generatorModule = cms.string('/Herwig/Generators/EventGenerator'),
+    eventHandlers = cms.string('/Herwig/EventHandlers'),
+    configFiles = cms.vstring(),
+    crossSection = cms.untracked.double(-1),
+    filterEfficiency = cms.untracked.double(1.0),
+    Matchbox = cms.vstring( 'read snippets/Matchbox.in',
+    'read snippets/PPCollider.in',
+    'cd /Herwig/EventHandlers',
+    'set EventHandler:LuminosityFunction:Energy 13000*GeV',
+    '## Model assumptions',
+    'read Matchbox/StandardModelLike.in',
+    'read Matchbox/DiagonalCKM.in',
+    '## Set the order of the couplings',
+    'cd /Herwig/MatrixElements/Matchbox',
+    'set Factory:OrderInAlphaS 0',
+    'set Factory:OrderInAlphaEW 2',
+    '## Select the process',
+    'do Factory:Process p p -> l+ l-',
+    '# read Matchbox/MadGraph-GoSam.in',
+    '# read Matchbox/MadGraph-MadGraph.in',
+    'read Matchbox/MadGraph-OpenLoops.in',
+    'set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 50*GeV',
+    'set /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 14000*GeV',
+    'cd /Herwig/MatrixElements/Matchbox',
+    'set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/LeptonPairMassScale',
+    'read Matchbox/MCatNLO-DefaultShower.in',
+    '# read Matchbox/NLO-NoShower.in',
+    '# read Matchbox/LO-NoShower.in',
+    'read Matchbox/FiveFlavourScheme.in',
+    'do /Herwig/MatrixElements/Matchbox/Factory:ProductionMode',
+    ),
+    parameterSets = cms.vstring('herwig7StableParticlesForDetector', 'Matchbox','herwig7CH3PDF', 'herwig7CH3AlphaS','herwig7CH3MPISettings'),
+)
+
+ProductionFilterSequence = cms.Sequence(generator)

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:

Configuration/Generator/python/concurrentLumisDisable.py

@@ -5,6 +5,7 @@
     "CosMuoGenProducer",
     "ExhumeGeneratorFilter",
     "Herwig7GeneratorFilter",
+    "Herwig7HadronizerFilter",
     "HydjetGeneratorFilter",
     "Hydjet2GeneratorFilter",
     "PyquenGeneratorFilter",

Configuration/PyReleaseValidation/python/relval_extendedgen.py

@@ -47,6 +47,8 @@
 workflows[535]=['',['TTbar_13TeV_Pow_herwig7','HARVESTGEN']]
 workflows[537]=['',['DYToLL012Jets_5FS_TuneCH3_13TeV_amcatnloFxFx_herwig7','HARVESTGEN']]
 workflows[538]=['',['DYToLL01234Jets_5FS_TuneCH3_13TeV_madgraphMLM_herwig7','HARVESTGEN']]
+workflows[539]=['',['DY_TuneCH3_13TeV_herwig_madgraph_matchbox','HARVESTGEN']]
+
 # External Decays
 
 workflows[521]=['WTolNuJets_LO_Mad_13TeV_py8_Ta',['WTolNu01234Jets_5f_LO_MLM_Madgraph_LHE_13TeV','Hadronizer_TuneCP5_13TeV_MLM_5f_max4j_LHE_pythia8_Tauola','HARVESTGEN2']]

Configuration/PyReleaseValidation/python/relval_steps.py

@@ -1716,6 +1716,7 @@ def gen2018RD(fragment,howMuch):
 steps['TTbar_13TeV_Pow_herwig7']=genvalid('Configuration/Generator/python/TT_13TeV_Pow_Herwig7_cff',step1LHEGenDQM)
 steps['DYToLL012Jets_5FS_TuneCH3_13TeV_amcatnloFxFx_herwig7']=genvalid('Configuration/Generator/python/DYToLL012Jets_5FS_TuneCH3_13TeV_amcatnloFxFx_herwig7_cff',merge([{'-n':'12'},step1LHEGenDQM]))
 steps['DYToLL01234Jets_5FS_TuneCH3_13TeV_madgraphMLM_herwig7']=genvalid('Configuration/Generator/python/DYToLL01234Jets_5FS_TuneCH3_13TeV_madgraphMLM_herwig7_cff',step1LHEGenDQM)
+steps['DY_TuneCH3_13TeV_herwig_madgraph_matchbox']=genvalid('Configuration/Generator/python/DY_TuneCH3_13TeV_herwig_madgraph_matchbox_cff',step1GenDefaults)
 
 
 # Heavy Ion

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