Skip to content

Commit

Permalink
Merge pull request #42337 from vslokenb/OneMoreTopRecoil_backport
Browse files Browse the repository at this point in the history
Top recoil hook cleaned, backport to CMSSW_13_1_X
  • Loading branch information
cmsbuild authored Aug 23, 2023
2 parents a8c35e4 + b327196 commit 5153b29
Show file tree
Hide file tree
Showing 4 changed files with 248 additions and 0 deletions.
19 changes: 19 additions & 0 deletions GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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
//
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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");
Expand Down
132 changes: 132 additions & 0 deletions GeneratorInterface/Pythia8Interface/plugins/TopRecoilHook.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
// 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();
}
//std::cout << "BEGIN TOP TEST: " << (wtT / wtW) << " | END TOP TEST" << endl;

// 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
5 changes: 5 additions & 0 deletions GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
92 changes: 92 additions & 0 deletions GeneratorInterface/Pythia8Interface/test/TopRecoilHook_test_cfg.py
Original file line number Diff line number Diff line change
@@ -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)


0 comments on commit 5153b29

Please sign in to comment.