Merge pull request #8728 from Dominic-Stafford/hw_lhe_match_fix_10_6

[10.6.X] Add patches to propagate lhe numbering through herwig
cms-sw · Oct 10, 2023 · cc3ddc7 · cc3ddc7
2 parents 01d7f69 + 12d3b49
commit cc3ddc7
Show file tree

Hide file tree

Showing 5 changed files with 385 additions and 18 deletions.
diff --git a/LHEEventNum.patch b/LHEEventNum.patch
@@ -0,0 +1,199 @@
+diff -u ThePEG-2.2.2/LesHouches/LesHouchesEventHandler.cc leshoucheseventnumbers/ThePEG/LesHouchesEventHandler.cc
+--- ThePEG-2.2.2/LesHouches/LesHouchesEventHandler.cc	2021-01-20 21:18:33.000000001 +0100
++++ leshoucheseventnumbers/ThePEG/LesHouchesEventHandler.cc	2022-01-12 19:04:24.000000001 +0100
+@@ -244,9 +244,18 @@
+       const double fact = 
+       	theNormWeight ? 
+       	  double(selector().sum()/picobarn) : 1.;
+- 
+-      currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
++
++      //whether to use the LHE event number or not for the event identification
++      if(UseLHEEvent==0 || currentReader()->LHEEventNum() == -1) { 
++	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
+                                  generator()->currentEventNumber(), weight*fact )));
++      }
++      else if(UseLHEEvent==1 && currentReader()->LHEEventNum() != -1) {
++	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
++				  currentReader()->LHEEventNum(), weight*fact )));
++      }
++
++      
+       currentEvent()->optionalWeights() = currentReader()->optionalEventWeights();
+       // normalize the optional weights
+       for(map<string,double>::iterator it = currentEvent()->optionalWeights().begin();
+@@ -512,13 +521,13 @@
+ void LesHouchesEventHandler::persistentOutput(PersistentOStream & os) const {
+   os << stats << histStats << theReaders << theSelector
+      << oenum(theWeightOption) << theUnitTolerance << theCurrentReader << warnPNum
+-     << theNormWeight;
++     << theNormWeight << UseLHEEvent;
+ }
+
+ void LesHouchesEventHandler::persistentInput(PersistentIStream & is, int) {
+   is >> stats >> histStats >> theReaders >> theSelector
+      >> ienum(theWeightOption) >> theUnitTolerance >> theCurrentReader >> warnPNum
+-     >> theNormWeight;
++     >> theNormWeight >> UseLHEEvent;
+ }
+
+ ClassDescription<LesHouchesEventHandler>
+@@ -613,6 +622,22 @@
+      "Normalize the weights to the max cross section in pb",
+      1);
+
++  
++  static Switch<LesHouchesEventHandler,unsigned int> interfaceEventNumbering
++    ("EventNumbering",
++     "How to number the events",
++     &LesHouchesEventHandler::UseLHEEvent, 0, false, false);
++  static SwitchOption interfaceEventNumberingIncremental
++    (interfaceEventNumbering,
++     "Incremental",
++     "Standard incremental numbering (i.e. as they are generated)",
++     0);
++  static SwitchOption interfaceEventNumberingLHE
++    (interfaceEventNumbering,
++     "LHE",
++     "Corresponding to the LHE event number",
++     1);
++
+
+   interfaceLesHouchesReaders.rank(10);
+   interfaceWeightOption.rank(9);
+diff -u build/BUILD/slc7_amd64_gcc900/external/thepeg/2.2.2-cms2/ThePEG-2.2.2/LesHouches/LesHouchesEventHandler.h leshoucheseventnumbers/ThePEG/LesHouchesEventHandler.h
+--- ThePEG-2.2.2/LesHouches/LesHouchesEventHandler.h	2021-01-20 21:18:33.000000001 +0100
++++ leshoucheseventnumbers/ThePEG/LesHouchesEventHandler.h	2022-01-12 19:04:24.000000001 +0100
+@@ -71,7 +71,7 @@
+    */
+   LesHouchesEventHandler()
+     : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true),
+-      theNormWeight(0)
++      theNormWeight(0), UseLHEEvent(0)
+   {
+     selector().tolerance(unitTolerance());
+   }
+@@ -405,6 +405,11 @@
+    */
+   unsigned int theNormWeight;
+
++  /** 
++   * How to number the events
++   */
++  unsigned int UseLHEEvent;
++
+ public:
+
+   /** @cond EXCEPTIONCLASSES */
+diff -u build/BUILD/slc7_amd64_gcc900/external/thepeg/2.2.2-cms2/ThePEG-2.2.2/LesHouches/LesHouchesFileReader.cc leshoucheseventnumbers/ThePEG/LesHouchesFileReader.cc
+--- ThePEG-2.2.2/LesHouches/LesHouchesFileReader.cc	2021-07-21 00:36:31.000000001 +0200
++++ leshoucheseventnumbers/ThePEG/LesHouchesFileReader.cc	2022-01-12 19:04:24.000000001 +0100
+@@ -546,6 +546,14 @@
+        */
+       if(readingInitWeights_sc && !cfile.find("</weightgroup")) {
+ 	hs = cfile.getline();
++   //cout << "hs=" << hs << endl;
++   //cout << "weightinfo= " << weightinfo << endl;
++   //fix for potential new lines:
++   if(!cfile.find("<weight") and !cfile.find("</weightgroup")) {
++     weightinfo = weightinfo + hs;
++     //cout << "weightinfo fixed= " << weightinfo << endl;
++     continue;
++   }
+ 	istringstream isc(hs);
+ 	int ws = 0;
+ 	/* get the name that will be used to identify the scale 
+@@ -724,6 +732,8 @@
+     } 
+   }
+
++  LHEeventnum = -1; // set the LHEeventnum to -1, this will be the default if the tag <event_num> is not found
++  
+   // Now read any additional comments and named weights.
+   // read until the end of rwgt is found
+   bool readingWeights = false, readingaMCFast = false, readingMG5ClusInfo = false;
+@@ -797,6 +807,20 @@
+       optionalWeights[mg5scinfo.c_str()] = -333; //for the mg5 scale info weights we give them a weight -333 for future identification
+     }
+
++    //the event num tag
++    if(cfile.find("<event_num")) {
++      string hs = cfile.getline();
++      string stopDEL = "</event_num>"; //end delimiter
++      unsigned firstLim = hs.find(">") + 1;
++      string LHEeventnum_str = hs.substr(firstLim);
++      erase_substr(LHEeventnum_str,stopDEL);
++      LHEeventnum =  std::stol(LHEeventnum_str, nullptr, 10); 
++    }
++
++
++
++
++
+     //determine start of aMCFast weights
+     if(cfile.find("<applgrid")) { readingaMCFast = true;}
+     //determine start of optional weights
+diff -u build/BUILD/slc7_amd64_gcc900/external/thepeg/2.2.2-cms2/ThePEG-2.2.2/LesHouches/LesHouchesReader.cc leshoucheseventnumbers/ThePEG/LesHouchesReader.cc
+--- ThePEG-2.2.2/LesHouches/LesHouchesReader.cc	2021-01-20 21:18:33.000000001 +0100
++++ leshoucheseventnumbers/ThePEG/LesHouchesReader.cc	2022-01-12 19:04:24.000000001 +0100
+@@ -1179,6 +1179,7 @@
+   pos = mwrite(pos, hepeup.SPINUP[0], hepeup.NUP);
+   pos = mwrite(pos, lastweight);
+   pos = mwrite(pos, optionalWeights);
++  pos = mwrite(pos, LHEeventnum);
+   for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) {
+     pos = mwrite(pos, optionalWeightsNames[ff]);
+   }
+@@ -1225,6 +1226,7 @@
+   pos = mread(pos, optionalnpLO);
+   pos = mread(pos, optionalnpNLO);
+   pos = mread(pos, preweight);
++  pos = mread(pos, LHEeventnum);
+
+   // If we are skipping, we do not have to do anything else.
+   if ( skipping ) return true;
+@@ -1254,7 +1256,7 @@
+      << thePartonBinInstances
+      << theBeams << theIncoming << theOutgoing << theIntermediates
+      << reweights << preweights << preweight << reweightPDF << doInitPDFs
+-     << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO
++     << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO << LHEeventnum
+      << maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights
+      << theMomentumTreatment << useWeightWarnings << theReOpenAllowed
+      << theIncludeSpin;
+@@ -1275,7 +1277,7 @@
+      >> thePartonBinInstances
+      >> theBeams >> theIncoming >> theOutgoing >> theIntermediates
+      >> reweights >> preweights >> preweight >> reweightPDF >> doInitPDFs
+-     >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO
++     >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO >> LHEeventnum
+      >> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights
+      >> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed
+      >> theIncludeSpin;
+diff -u build/BUILD/slc7_amd64_gcc900/external/thepeg/2.2.2-cms2/ThePEG-2.2.2/LesHouches/LesHouchesReader.h leshoucheseventnumbers/ThePEG/LesHouchesReader.h
+--- ThePEG-2.2.2/LesHouches/LesHouchesReader.h	2021-01-20 21:18:33.000000001 +0100
++++ leshoucheseventnumbers/ThePEG/LesHouchesReader.h	2022-01-12 19:04:24.000000001 +0100
+@@ -286,6 +286,11 @@
+    */
+   const map<string,double>& optionalEventWeights() const { return optionalWeights; }
+
++  /** 
++   * Return the Les Houches event number associated with the current event
++   */
++  const long& LHEEventNum() const { return LHEeventnum; }
++
+   /**
+    * Return the optional npLO and npNLO
+    */
+@@ -831,6 +836,11 @@
+    */
+   map<string,double> optionalWeights;
+
++  /**
++   * The event number
++   */
++  long LHEeventnum;
++
+
+   /**
+    * If the maximum cross section of this reader has been increased
diff --git a/LHEEventNumFxFx.patch b/LHEEventNumFxFx.patch
@@ -0,0 +1,181 @@
+diff -u Herwig-7.2.2/MatrixElement/FxFx/FxFxEventHandler.cc ../../../../../../../leshoucheseventnumbers/FxFx/FxFxEventHandler.cc
+--- Herwig-7.2.2/MatrixElement/FxFx/FxFxEventHandler.cc	2021-01-20 20:46:46.000000001 +0100
++++ ../../../../../../../leshoucheseventnumbers/FxFx/FxFxEventHandler.cc	2022-01-12 19:04:24.000000001 +0100
+@@ -249,8 +249,15 @@
+
+       theLastXComb = currentReader()->getXComb();
+
+-      currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
+-				 generator()->currentEventNumber(), weight*fact)));
++      //whether to use the LHE event number or not for the event identification
++      if(UseLHEEvent==0 || currentReader()->LHEEventNum() == -1) { 
++	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
++                                 generator()->currentEventNumber(), weight*fact )));
++      }
++      else if(UseLHEEvent==1 && currentReader()->LHEEventNum() != -1) {
++	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
++				  currentReader()->LHEEventNum(), weight*fact )));
++      }
+       currentEvent()->optionalWeights() = currentReader()->optionalEventWeights();
+      // normalize the optional weights
+       for(map<string,double>::iterator it = currentEvent()->optionalWeights().begin();
+@@ -543,14 +550,14 @@
+ void FxFxEventHandler::persistentOutput(PersistentOStream & os) const {
+   os << stats << histStats << theReaders << theSelector
+      << oenum(theWeightOption) << theUnitTolerance << theCurrentReader << warnPNum
+-     << theNormWeight;
++     << theNormWeight << UseLHEEvent;
+
+ }
+
+ void FxFxEventHandler::persistentInput(PersistentIStream & is, int) {
+   is >> stats >> histStats >> theReaders >> theSelector
+      >> ienum(theWeightOption) >> theUnitTolerance >> theCurrentReader >> warnPNum
+-     >> theNormWeight;
++     >> theNormWeight >> UseLHEEvent;;
+ }
+
+ ClassDescription<FxFxEventHandler>
+@@ -644,6 +651,20 @@
+      "Normalize the weights to the max cross section in pb",
+      1);
+
++    static Switch<FxFxEventHandler,unsigned int> interfaceEventNumbering
++    ("EventNumbering",
++     "How to number the events",
++     &FxFxEventHandler::UseLHEEvent, 0, false, false);
++  static SwitchOption interfaceEventNumberingIncremental
++    (interfaceEventNumbering,
++     "Incremental",
++     "Standard incremental numbering (i.e. as they are generated)",
++     0);
++  static SwitchOption interfaceEventNumberingLHE
++    (interfaceEventNumbering,
++     "LHE",
++     "Corresponding to the LHE event number",
++     1);
+
+   interfaceFxFxReaders.rank(10);
+   interfaceWeightOption.rank(9);
+diff -u Herwig-7.2.2/MatrixElement/FxFx/FxFxEventHandler.h ../../../../../../../leshoucheseventnumbers/FxFx/FxFxEventHandler.h
+--- Herwig-7.2.2/MatrixElement/FxFx/FxFxEventHandler.h	2021-01-20 20:46:46.000000001 +0100
++++ ../../../../../../../leshoucheseventnumbers/FxFx/FxFxEventHandler.h	2022-01-12 19:04:24.000000001 +0100
+@@ -72,7 +72,7 @@
+    * The default constructor.
+    */
+   FxFxEventHandler()
+-    : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true), theNormWeight(0)
++    : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true), theNormWeight(0), UseLHEEvent(0)
+   {
+     selector().tolerance(unitTolerance());
+   }
+@@ -377,6 +377,11 @@
+    */
+   unsigned int theNormWeight;
+
++  /** 
++   * How to number the events
++   */
++  unsigned int UseLHEEvent;
++
+
+ public:
+
+diff -u Herwig-7.2.2/MatrixElement/FxFx/FxFxFileReader.cc ../../../../../../../leshoucheseventnumbers/FxFx/FxFxFileReader.cc
+--- Herwig-7.2.2/MatrixElement/FxFx/FxFxFileReader.cc	2021-01-20 20:46:46.000000001 +0100
++++ ../../../../../../../leshoucheseventnumbers/FxFx/FxFxFileReader.cc	2022-01-12 19:04:24.000000001 +0100
+@@ -707,6 +707,9 @@
+     } 
+   }
+
++  LHEeventnum = -1; // set the LHEeventnum to -1, this will be the default if the tag <event_num> is not found
++
++
+  // Now read any additional comments and named weights.
+   // read until the end of rwgt is found
+   bool readingWeights = false, readingaMCFast = false, readingMG5ClusInfo = false;
+@@ -767,6 +770,20 @@
+       erase_substr(mg5clusinfo,str_newline);
+       optionalWeights[mg5clusinfo.c_str()] = -222; //for the mg5 scale info weights we give them a weight -222 for future identification
+     }
++
++    //the event num tag
++    if(cfile.find("<event_num")) {
++      string hs = cfile.getline();
++      string stopDEL = "</event_num>"; //end delimiter
++      unsigned firstLim = hs.find(">") + 1;
++      string LHEeventnum_str = hs.substr(firstLim);
++      erase_substr(LHEeventnum_str,stopDEL);
++      LHEeventnum =  std::stol(LHEeventnum_str, nullptr, 10); 
++    }
++
++
++
++
+
+     //store MG5 clustering information 
+     if(cfile.find("<scales")) {
+diff -u Herwig-7.2.2/MatrixElement/FxFx/FxFxReader.cc ../../../../../../../leshoucheseventnumbers/FxFx/FxFxReader.cc
+--- Herwig-7.2.2/MatrixElement/FxFx/FxFxReader.cc	2021-01-20 20:46:46.000000001 +0100
++++ ../../../../../../../leshoucheseventnumbers/FxFx/FxFxReader.cc	2022-01-12 19:04:24.000000001 +0100
+@@ -1185,6 +1185,7 @@
+   pos = mwrite(pos, hepeup.SPINUP[0], hepeup.NUP);
+   pos = mwrite(pos, lastweight);
+   pos = mwrite(pos, optionalWeights);
++  pos = mwrite(pos, LHEeventnum);
+   for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) {
+     pos = mwrite(pos, optionalWeightsNames[ff]);
+   }
+@@ -1234,6 +1235,7 @@
+   pos = mread(pos, optionalnpLO);
+   pos = mread(pos, optionalnpNLO);
+   pos = mread(pos, preweight);
++  pos = mread(pos, LHEeventnum);
+
+   // If we are skipping, we do not have to do anything else.
+   if ( skipping ) return true;
+@@ -1263,7 +1265,7 @@
+      << thePartonBinInstances
+      << theBeams << theIncoming << theOutgoing << theIntermediates
+      << reweights << preweights << preweight << reweightPDF << doInitPDFs
+-     << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO
++     << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO << LHEeventnum
+      << maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights
+      << theMomentumTreatment << useWeightWarnings << theReOpenAllowed
+      << theIncludeSpin;
+@@ -1284,7 +1286,7 @@
+      >> thePartonBinInstances
+      >> theBeams >> theIncoming >> theOutgoing >> theIntermediates
+      >> reweights >> preweights >> preweight >> reweightPDF >> doInitPDFs
+-     >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO
++     >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO >> LHEeventnum
+      >> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights
+      >> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed
+      >> theIncludeSpin;
+diff -u Herwig-7.2.2/MatrixElement/FxFx/FxFxReader.h ../../../../../../../leshoucheseventnumbers/FxFx/FxFxReader.h
+--- Herwig-7.2.2/MatrixElement/FxFx/FxFxReader.h	2021-01-20 20:46:46.000000001 +0100
++++ ../../../../../../../leshoucheseventnumbers/FxFx/FxFxReader.h	2022-01-12 19:04:24.000000001 +0100
+@@ -290,6 +290,11 @@
+    */
+   const map<string,double>& optionalEventWeights() const { return optionalWeights; }
+
++  /** 
++   * Return the Les Houches event number associated with the current event
++   */
++  const long& LHEEventNum() const { return LHEeventnum; }
++
+   /**
+    * Return the optional npLO and npNLO
+    */
+@@ -834,6 +839,11 @@
+    * The optional weights associated to the last read events.
+    */
+   map<string,double> optionalWeights;
++  
++  /**
++   * The event number
++   */
++  long LHEeventnum;
+
+    /**
+    * If the maximum cross section of this reader has been increased
diff --git a/LesHouchesFileReader.patch b/LesHouchesFileReader.patch