Herwig lhe matching fix 10_6

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

@@ -65,6 +65,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_;
@@ -110,8 +111,20 @@ 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)
@@ -155,9 +168,38 @@ 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.";
-	return false;
+    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;
+    } else 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;
 }
 
 std::unique_ptr<GenLumiInfoHeader> Herwig7Hadronizer::getGenLumiInfoHeader() const {

GeneratorInterface/LHEInterface/interface/LHEEvent.h

@@ -55,9 +55,11 @@ class LHEEvent {
 
 	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); }
 
@@ -95,6 +97,7 @@ class LHEEvent {
         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

@@ -188,7 +188,7 @@ ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
   if (!partonLevel_) {
     throw edm::Exception(edm::errors::EventGenerationFailure) << "No lhe event found in ExternalLHEProducer::produce().  "
     << "The likely cause is that the lhe file contains fewer events than were requested, which is possible "
-    << "in case of phase space integration or uneweighting efficiency problems.";
+    << "in case of phase space integration or unweighting efficiency problems.";
   }
 
   std::unique_ptr<LHEEventProduct> product(
@@ -203,6 +203,7 @@ ExternalLHEProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
                 boost::bind(&LHEEventProduct::addWeight,
                             product.get(), _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 python
+
+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,10 @@ 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 +417,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

@@ -36,7 +36,7 @@ namespace lhef {
 LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
                    std::istream &in) :
   runInfo(runInfo), weights_(0), counted(false), 
-  readAttemptCounter(0), npLO_(-99), npNLO_(-99)
+  readAttemptCounter(0), npLO_(-99), npNLO_(-99), evtnum_(-1)
 
 {
 	hepeup.NUP = 0;
@@ -116,7 +116,7 @@ LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
 LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
                    const HEPEUP &hepeup) :
 	runInfo(runInfo), hepeup(hepeup), counted(false), readAttemptCounter(0),
-        npLO_(-99), npNLO_(-99)
+        npLO_(-99), npNLO_(-99), evtnum_(-1)
 {
 }
 
@@ -126,7 +126,7 @@ LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
                    const std::vector<std::string> &comments) :
 	runInfo(runInfo), hepeup(hepeup), pdf(pdf ? new PDF(*pdf) : nullptr),
 	comments(comments), counted(false), readAttemptCounter(0),
-	npLO_(-99), npNLO_(-99)
+	npLO_(-99), npNLO_(-99), evtnum_(-1)
 {
 }
 
@@ -138,7 +138,7 @@ LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
 	comments(product.comments_begin(), product.comments_end()),
 	counted(false), readAttemptCounter(0),
 	originalXWGTUP_(product.originalXWGTUP()), scales_(product.scales()),
-	npLO_(product.npLO()), npNLO_(product.npNLO())
+	npLO_(product.npLO()), npNLO_(product.npNLO()), evtnum_(product.evtnum())
 {
 }
 

GeneratorInterface/LHEInterface/src/LHEReader.cc

@@ -153,6 +153,7 @@ class LHEReader::XMLHandler : public XMLDocument::Handler {
         int                             npLO;
         int                             npNLO;
         std::vector<float>              scales;
+        int                             evtnum=-1;
 };
 
 static void attributesToDom(DOMElement *dom, const Attributes &attributes)
@@ -234,7 +235,11 @@ void LHEReader::XMLHandler::startElement(const XMLCh *const uri,
 
         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;
   } else if (mode == kInit) {
@@ -576,12 +581,14 @@ LHEReader::~LHEReader()
 	  lheevent->addWeight(gen::WeightsInfo(info[i].first,num));
 	}
 	lheevent->setNpLO(handler->npLO);
-        lheevent->setNpNLO(handler->npNLO);
-        //fill scales
-        if (!handler->scales.empty()) {
-          lheevent->setScales(handler->scales);
-        }
-        return lheevent;
+    lheevent->setNpNLO(handler->npNLO);
+    lheevent->setEvtNum(handler->evtnum);
+    handler->evtnum = -1;
+    //fill scales
+    if (!handler->scales.empty()) {
+      lheevent->setScales(handler->scales);
+    }
+    return lheevent;
 	}
       }
     }

SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h

@@ -33,6 +33,7 @@ class LHEEventProduct {
 	  scales_ = std::move(other.scales_);
 	  npLO_ = std::move(other.npLO_);
 	  npNLO_ = std::move(other.npNLO_);
+	  evtnum_ = std::move(other.evtnum_);
 	}
 	LHEEventProduct& operator=(LHEEventProduct&& other) {
 	  hepeup_ = std::move(other.hepeup_);
@@ -43,6 +44,7 @@ class LHEEventProduct {
 	  scales_ = std::move(other.scales_);
 	  npLO_ = std::move(other.npLO_);
 	  npNLO_ = std::move(other.npNLO_);
+	  evtnum_ = std::move(other.evtnum_);
 	  return *this;
 	}
 	~LHEEventProduct() {}
@@ -61,9 +63,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 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(); }
@@ -125,6 +129,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