Top recoil hook cleaned, backport to CMSSW_13_2_X

GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc

@@ -19,6 +19,7 @@ using namespace Pythia8;
 
 #include "GeneratorInterface/Pythia8Interface/plugins/ReweightUserHooks.h"
 #include "GeneratorInterface/Pythia8Interface/interface/CustomHook.h"
+#include "GeneratorInterface/Pythia8Interface/plugins/TopRecoilHook.h"
 
 // PS matchning prototype
 //
@@ -152,6 +153,9 @@ class Pythia8Hadronizer : public Py8InterfaceBase {
   //Generic customized hooks vector
   std::shared_ptr<UserHooksVector> fCustomHooksVector;
 
+  //RecoilToTop userhook
+  std::shared_ptr<TopRecoilHook> fTopRecoilHook;
+
   int EV1_nFinal;
   bool EV1_vetoOn;
   int EV1_maxVetoCount;
@@ -425,6 +429,14 @@ bool Pythia8Hadronizer::initializeForInternalPartons() {
     (fUserHooksVector->hooks).push_back(fPowhegHooksBB4L);
   }
 
+  bool TopRecoilHook1 = fMasterGen->settings.flag("TopRecoilHook:doTopRecoilIn");
+  if (TopRecoilHook1) {
+    edm::LogInfo("Pythia8Interface") << "Turning on RecoilToTop hook from Pythia8Interface";
+    if (!fTopRecoilHook.get())
+      fTopRecoilHook.reset(new TopRecoilHook());
+    (fUserHooksVector->hooks).push_back(fTopRecoilHook);
+  }
+
   //adapted from main89.cc in pythia8 examples
   bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
   bool internalMerging = !(fMasterGen->settings.word("Merging:Process") == "void");
@@ -594,6 +606,13 @@ bool Pythia8Hadronizer::initializeForExternalPartons() {
     (fUserHooksVector->hooks).push_back(fPowhegHooksBB4L);
   }
 
+  bool TopRecoilHook1 = fMasterGen->settings.flag("TopRecoilHook:doTopRecoilIn");
+  if (TopRecoilHook1) {
+    edm::LogInfo("Pythia8Interface") << "Turning on RecoilToTop hook from Pythia8Interface";
+    if (!fTopRecoilHook.get())
+      fTopRecoilHook.reset(new TopRecoilHook());
+    (fUserHooksVector->hooks).push_back(fTopRecoilHook);
+  }
   //adapted from main89.cc in pythia8 examples
   bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
   bool internalMerging = !(fMasterGen->settings.word("Merging:Process") == "void");

GeneratorInterface/Pythia8Interface/plugins/TopRecoilHook.h

@@ -0,0 +1,131 @@
+// TopRecoilHook.h is a part of the PYTHIA event generator.
+// Copyright (C) 2020 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// Includes a user hook that corrects emission in top decay for dipole
+// from gluon to W, to instead be from gluon to top.
+
+// Important: the top mass shift analysis encoded here is very primitive,
+// does not perform well at all, and should not be taken seriously.
+// The important part is that you see how the different scenarios
+// should be set up to operate as intended.
+
+#include "Pythia8/Pythia.h"
+namespace Pythia8 {
+
+  //==========================================================================
+
+  // Write own derived UserHooks class for modified emission in top decay.
+
+  class TopRecoilHook : public UserHooks {
+  public:
+    // Constructor.
+    //  doTopRecoil : eikonal correction in GW dipole on/off when no MEC applied.
+    //  useOldDipole  : in GW dipole, use partons before or after branching.
+    //  doList        : diagnostic output; set false for production runs.
+    TopRecoilHook(bool doTopRecoilIn = true, bool useOldDipoleIn = false, bool doListIn = false) {
+      doTopRecoil = doTopRecoilIn;
+      useOldDipole = useOldDipoleIn;
+      doList = doListIn;
+      // Constructor also creates some histograms for analysis inside User Hook.
+      wtCorr = new Hist("corrective weight", 100, 0., 2.);
+    }
+
+    // Destructor prints histogram.
+    ~TopRecoilHook() override { delete wtCorr; }
+
+    // Initialise. Only use hook for simple showers with recoilToColoured = off.
+    bool initAfterBeams() override {
+      // Switch off if recoilToColoured = on.
+      bool recoilToColoured = settingsPtr->flag("TimeShower:recoilToColoured");
+      if (recoilToColoured)
+        doTopRecoil = false;
+      // Flag if W mass term is already accounted for (true) or not (false).
+      recoilDeadCone = settingsPtr->flag("TimeShower:recoilDeadCone");
+      // All ok.
+      return true;
+    }
+
+    // Allow a veto after an FSR emission
+    bool canVetoFSREmission() override { return doTopRecoil; }
+
+    // Access the event after an FSR emission, specifically inside top decay.
+    bool doVetoFSREmission(int sizeOld, const Event& event, int iSys, bool inResonance) override {
+      // Check that we are inside a resonance decay.
+      if (!inResonance)
+        return false;
+
+      // Check that it is a top decay.
+      int iTop = partonSystemsPtr->getInRes(iSys);
+      if (iTop == 0 || event[iTop].idAbs() != 6)
+        return false;
+
+      // Skip first emission, where ME corrections are already made.
+      int sizeOut = partonSystemsPtr->sizeOut(iSys);
+      if (sizeOut == 2)
+        return false;
+
+      // Location of trial new particles: radiator, emitted, recoiler.
+      int iRad = sizeOld;
+      int iEmt = sizeOld + 1;
+      int iRec = sizeOld + 2;
+
+      // The above partons are after emission;
+      // alternatively use the ones before.
+      if (useOldDipole) {
+        iRad = event[iRad].mother1();
+        iRec = event[iRec].mother1();
+      }
+
+      // Check if newly emitted gluon matches (anti)top colour line.
+      if (event[iEmt].id() != 21)
+        return false;
+      if (event[iTop].id() == 6) {
+        if (event[iEmt].col() != event[iTop].col())
+          return false;
+      } else {
+        if (event[iEmt].acol() != event[iTop].acol())
+          return false;
+      }
+
+      // Recoiler should now be a W, else something is wrong.
+      if (event[iRec].idAbs() != 24) {
+        return false;
+      }
+
+      // Denominator: eikonal weight with W as recoiler.
+      double pRadRec = event[iRad].p() * event[iRec].p();
+      double pRadEmt = event[iRad].p() * event[iEmt].p();
+      double pRecEmt = event[iRec].p() * event[iEmt].p();
+      double wtW = 2. * pRadRec / (pRadEmt * pRecEmt) - pow2(event[iRad].m() / pRadEmt);
+      // If recoilDeadCone = on, include W mass term in denominator.
+      if (recoilDeadCone)
+        wtW -= pow2(event[iRec].m() / pRecEmt);
+
+      // Numerator: eikonal weight with top as recoiler.
+      double pRadTop = event[iRad].p() * event[iTop].p();
+      double pTopEmt = event[iTop].p() * event[iEmt].p();
+      double wtT =
+          2. * pRadTop / (pRadEmt * pTopEmt) - pow2(event[iRad].m() / pRadEmt) - pow2(event[iTop].m() / pTopEmt);
+
+      // Histogram weight ratio.
+      wtCorr->fill(wtT / wtW);
+
+      // List relevant properties.
+      if (doList) {
+        partonSystemsPtr->list();
+        event.list();
+      }
+
+      // Accept/reject emission. Smooth suppression or step function.
+      return (wtT < wtW * rndmPtr->flat());
+    }
+
+  private:
+    // Options and Histograms.
+    bool doTopRecoil, useOldDipole, doList, recoilDeadCone;
+    Hist* wtCorr;
+  };
+
+}  // namespace Pythia8

GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc

@@ -104,6 +104,11 @@ namespace gen {
     fMasterGen->settings.addParm("PTFilter:quarkRapidity", 10.0, true, true, 0.0, 10.);
     fMasterGen->settings.addParm("PTFilter:quarkPt", -.1, true, true, -.1, 100.);
 
+    //add settings for RecoilToTop tool
+    fMasterGen->settings.addFlag("TopRecoilHook:doTopRecoilIn", false);
+    fMasterGen->settings.addFlag("TopRecoilHook:useOldDipoleIn", false);
+    fMasterGen->settings.addFlag("TopRecoilHook:doListIn", false);
+
     //add settings for powheg resonance scale calculation
     fMasterGen->settings.addFlag("POWHEGres:calcScales", false);
     fMasterGen->settings.addFlag("POWHEG:bb4l", false);

GeneratorInterface/Pythia8Interface/test/TopRecoilHook_test_cfg.py

@@ -0,0 +1,92 @@
+import sys
+import FWCore.ParameterSet.Config as cms
+import FWCore.ParameterSet.VarParsing as VarParsing
+#import pythia8
+
+options = VarParsing.VarParsing ('standard')
+options.register('runOnly', '', VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "Run only specified analysis")
+options.register('yodafile', 'test_toprecoil.yoda', VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "Name of yoda output file")
+options.register('recoiltoTop','on',VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "TopRecoilHook setting")
+options.register('recoiltoBottom','off',VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "RecoilToColor setting")
+options.setDefault('maxEvents', 1000)
+options.register('alphaSvalue', 0.1365,VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "AlphaS value setting")
+options.register('topmass', 172.5,VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "Top mass setting")
+if(hasattr(sys, "argv")):
+    options.parseArguments()
+print(options)
+
+process = cms.Process("runRivetAnalysis")
+
+from Configuration.Generator.Pythia8CommonSettings_cfi import *
+from Configuration.Generator.Pythia8CUEP8M1Settings_cfi import *
+from Configuration.Generator.Pythia8aMCatNLOSettings_cfi import *
+from IOMC.RandomEngine.RandomServiceHelper import RandomNumberServiceHelper
+#randSvc = RandomNumberServiceHelper(process.RandomNumberGeneratorService)
+#randSvc.populate()
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+
+randSvc = RandomNumberServiceHelper(process.RandomNumberGeneratorService)
+randSvc.populate()
+
+process.load("FWCore.MessageLogger.MessageLogger_cfi")
+process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(1000)
+process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(options.maxEvents))
+
+process.source = cms.Source("EmptySource")
+
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.generator = cms.EDFilter("Pythia8GeneratorFilter",
+     comEnergy = cms.double(13000.0),
+     crossSection = cms.untracked.double(421.1),
+     filterEfficiency = cms.untracked.double(1),
+     maxEventsToPrint = cms.untracked.int32(0),
+     pythiaHepMCVerbosity = cms.untracked.bool(False),
+     pythiaPylistVerbosity = cms.untracked.int32(1),
+     PythiaParameters = cms.PSet(
+       pythia8CommonSettingsBlock,
+       pythia8CUEP8M1SettingsBlock,
+       pythia8aMCatNLOSettingsBlock,
+       processParameters = cms.vstring(
+             'Main:timesAllowErrors = 10000',
+             'ParticleDecays:limitTau0 = on',
+             'ParticleDecays:tauMax = 10',
+             'Tune:ee 7',
+             'Tune:pp 14',      # Monash tune
+             'Top:gg2ttbar    = on',
+             'Top:qqbar2ttbar = on',
+             '6:m0 =  %s'%options.topmass,    # top mass'
+             'TopRecoilHook:doTopRecoilIn = %s'%options.recoiltoTop,
+             'TimeShower:recoilToColoured = %s'%options.recoiltoBottom,
+             'TimeShower:alphaSvalue = %s' %options.alphaSvalue
+       ),
+         parameterSets = cms.vstring('processParameters',)
+     )
+)
+
+process.load("GeneratorInterface.RivetInterface.rivetAnalyzer_cfi")
+
+#pythia.Settings.listAll()
+
+if options.runOnly:
+    process.rivetAnalyzer.AnalysisNames = cms.vstring(options.runOnly)
+else:
+    process.rivetAnalyzer.AnalysisNames = cms.vstring(
+        'CMS_2018_I1663958', # diff xs lepton+jets
+        'CMS_2018_I1662081', # event variables lepton+jets
+        'MC_TOPMASS_LJETS', # MC analysis for lepton+jets top mass
+        'CMS_2018_I1690148',  # jet substructure
+        'MC_BFRAG_LJETS',  # b fragmentation
+        'CMS_2018_I1703993', # diff xs dilepton
+        'CMS_2019_I1764472', #boosted top mass
+        'MC_TOPMASS_LSMT', #dilepton invariant mass, soft muon +primary
+        'ATLAS_2019_I1724098', #missing momentum analysis
+    )
+process.rivetAnalyzer.OutputFile      = options.yodafile
+process.rivetAnalyzer.HepMCCollection = cms.InputTag("generator:unsmeared")
+process.rivetAnalyzer.CrossSection    = 831.76 # NNLO (arXiv:1303.6254)
+
+process.p = cms.Path(process.generator*process.rivetAnalyzer)
+
+