Skip to content

Commit 2bf15f2

Browse files
cmsbuildcmsbuild
authored
Merge pull request #47428 from tomasra/pythia8hadronizer_hepmc2to3
Port recent changes from Pythia8Hadronizer to Pythia8HepMC3Hadronizer
2 parents 8984e39 + cd3cd05 commit 2bf15f2

File tree

2 files changed

+100
-8
lines changed

2 files changed

+100
-8
lines changed

‎GeneratorInterface/Pythia8Interface/plugins/Pythia8HepMC3Hadronizer.cc‎

Lines changed: 86 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,9 @@ using namespace Pythia8;
1717

1818
#include "GeneratorInterface/Pythia8Interface/interface/Py8HMC3InterfaceBase.h"
1919

20-
#include "GeneratorInterface/Pythia8Interface/plugins/ReweightUserHooks.h"
20+
#include "ReweightUserHooks.h"
2121
#include "GeneratorInterface/Pythia8Interface/interface/CustomHook.h"
22+
#include "TopRecoilHook.h"
2223

2324
// PS matchning prototype
2425
//
@@ -75,6 +76,49 @@ namespace CLHEP {
7576

7677
using namespace gen;
7778

79+
//Insert class for use w/ PDFPtr for proton-photon flux
80+
//parameters hardcoded according to main70.cc in PYTHIA8 v3.10
81+
class Nucleus2gamma2 : public Pythia8::PDF {
82+
private:
83+
double radius;
84+
int z;
85+
86+
public:
87+
// Constructor.
88+
Nucleus2gamma2(int idBeamIn, double R = -1.0, int Z = -1) : Pythia8::PDF(idBeamIn), radius(R), z(Z) {}
89+
90+
void xfUpdate(int, double x, double) override {
91+
if (z == -1) {
92+
// lead
93+
if (idBeam == 1000822080)
94+
z = 82;
95+
}
96+
if (radius == -1) {
97+
// lead
98+
if (idBeam == 1000822080)
99+
radius = 6.636;
100+
}
101+
102+
if (z < 0 || radius < 0)
103+
throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
104+
<< " Invalid photon flux input parameters: beam ID= " << idBeam << " , radius= " << radius << " , z= " << z
105+
<< "\n";
106+
107+
// Minimum impact parameter (~2*radius) [fm].
108+
double bmin = 2 * radius;
109+
110+
// Per-nucleon mass for lead.
111+
double m2 = pow2(0.9314);
112+
double alphaEM = 0.007297353080;
113+
double hbarc = 0.197;
114+
double xi = x * sqrt(m2) * bmin / hbarc;
115+
double bK0 = besselK0(xi);
116+
double bK1 = besselK1(xi);
117+
double intB = xi * bK1 * bK0 - 0.5 * pow2(xi) * (pow2(bK1) - pow2(bK0));
118+
xgamma = 2. * alphaEM * pow2(z) / M_PI * intB;
119+
}
120+
};
121+
78122
class Pythia8HepMC3Hadronizer : public Py8HMC3InterfaceBase {
79123
public:
80124
Pythia8HepMC3Hadronizer(const edm::ParameterSet &params);
@@ -112,6 +156,10 @@ class Pythia8HepMC3Hadronizer : public Py8HMC3InterfaceBase {
112156
double fBeam1PZ;
113157
double fBeam2PZ;
114158

159+
//PDFPtr for the photonFlux
160+
//Following main70.cc example in PYTHIA8 v3.10
161+
edm::ParameterSet photonFluxParams;
162+
115163
//helper class to allow multiple user hooks simultaneously
116164
std::shared_ptr<UserHooksVector> fUserHooksVector;
117165
bool UserHooksSet;
@@ -150,6 +198,9 @@ class Pythia8HepMC3Hadronizer : public Py8HMC3InterfaceBase {
150198
//Generic customized hooks vector
151199
std::shared_ptr<UserHooksVector> fCustomHooksVector;
152200

201+
//RecoilToTop userhook
202+
std::shared_ptr<TopRecoilHook> fTopRecoilHook;
203+
153204
int EV1_nFinal;
154205
bool EV1_vetoOn;
155206
int EV1_maxVetoCount;
@@ -213,6 +264,10 @@ Pythia8HepMC3Hadronizer::Pythia8HepMC3Hadronizer(const edm::ParameterSet &params
213264
// avoid filling weights twice (from v8.30x)
214265
toHepMC.set_store_weights(false);
215266

267+
if (params.exists("PhotonFlux")) {
268+
photonFluxParams = params.getParameter<edm::ParameterSet>("PhotonFlux");
269+
}
270+
216271
// Reweight user hook
217272
//
218273
if (params.exists("reweightGen")) {
@@ -364,6 +419,20 @@ bool Pythia8HepMC3Hadronizer::initializeForInternalPartons() {
364419
fMasterGen->settings.word("Beams:LHEF", lheFile_);
365420
}
366421

422+
if (!photonFluxParams.empty()) {
423+
const auto &beamTypeA = photonFluxParams.getParameter<int>("beamTypeA");
424+
const auto &beamTypeB = photonFluxParams.getParameter<int>("beamTypeB");
425+
const auto &radiusA = photonFluxParams.getUntrackedParameter<double>("radiusA", -1.0);
426+
const auto &radiusB = photonFluxParams.getUntrackedParameter<double>("radiusB", -1.0);
427+
const auto &zA = photonFluxParams.getUntrackedParameter<int>("zA", -1);
428+
const auto &zB = photonFluxParams.getUntrackedParameter<int>("zB", -1);
429+
Pythia8::PDFPtr photonFluxA =
430+
fMasterGen->settings.flag("PDF:beamA2gamma") ? make_shared<Nucleus2gamma2>(beamTypeA, radiusA, zA) : nullptr;
431+
Pythia8::PDFPtr photonFluxB =
432+
fMasterGen->settings.flag("PDF:beamB2gamma") ? make_shared<Nucleus2gamma2>(beamTypeB, radiusB, zB) : nullptr;
433+
fMasterGen->setPhotonFluxPtr(photonFluxA, photonFluxB);
434+
}
435+
367436
if (!fUserHooksVector.get())
368437
fUserHooksVector = std::make_shared<UserHooksVector>();
369438
(fUserHooksVector->hooks).clear();
@@ -412,6 +481,14 @@ bool Pythia8HepMC3Hadronizer::initializeForInternalPartons() {
412481
(fUserHooksVector->hooks).push_back(fPowhegHooksBB4L);
413482
}
414483

484+
bool TopRecoilHook1 = fMasterGen->settings.flag("TopRecoilHook:doTopRecoilIn");
485+
if (TopRecoilHook1) {
486+
edm::LogInfo("Pythia8Interface") << "Turning on RecoilToTop hook from Pythia8Interface";
487+
if (!fTopRecoilHook.get())
488+
fTopRecoilHook = std::make_shared<TopRecoilHook>();
489+
(fUserHooksVector->hooks).push_back(fTopRecoilHook);
490+
}
491+
415492
//adapted from main89.cc in pythia8 examples
416493
bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
417494
bool internalMerging = !(fMasterGen->settings.word("Merging:Process") == "void");
@@ -571,6 +648,14 @@ bool Pythia8HepMC3Hadronizer::initializeForExternalPartons() {
571648
(fUserHooksVector->hooks).push_back(fPowhegHooksBB4L);
572649
}
573650

651+
bool TopRecoilHook1 = fMasterGen->settings.flag("TopRecoilHook:doTopRecoilIn");
652+
if (TopRecoilHook1) {
653+
edm::LogInfo("Pythia8Interface") << "Turning on RecoilToTop hook from Pythia8Interface";
654+
if (!fTopRecoilHook.get())
655+
fTopRecoilHook = std::make_shared<TopRecoilHook>();
656+
(fUserHooksVector->hooks).push_back(fTopRecoilHook);
657+
}
658+
574659
//adapted from main89.cc in pythia8 examples
575660
bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
576661
bool internalMerging = !(fMasterGen->settings.word("Merging:Process") == "void");

‎GeneratorInterface/Pythia8Interface/src/Py8HMC3InterfaceBase.cc‎

Lines changed: 14 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -70,9 +70,15 @@ namespace gen {
7070
}
7171

7272
bool Py8HMC3InterfaceBase::readSettings(int) {
73+
//Pythia 8's default value for first argument to constructor
74+
const string xmlDir = "../share/Pythia8/xmldoc";
75+
bool printBanner = true;
76+
if (fParameters.exists("printBanner")) {
77+
printBanner = fParameters.getUntrackedParameter<bool>("printBanner");
78+
}
7379
if (!fMasterGen.get())
74-
fMasterGen = std::make_unique<Pythia>();
75-
fDecayer = std::make_unique<Pythia>();
80+
fMasterGen = std::make_unique<Pythia>(xmlDir, printBanner);
81+
fDecayer = std::make_unique<Pythia>(xmlDir, printBanner);
7682

7783
//add settings for resonance decay filter
7884
fMasterGen->settings.addFlag("BiasedTauDecayer:filter", false);
@@ -98,21 +104,22 @@ namespace gen {
98104
fMasterGen->settings.addParm("PTFilter:quarkRapidity", 10.0, true, true, 0.0, 10.);
99105
fMasterGen->settings.addParm("PTFilter:quarkPt", -.1, true, true, -.1, 100.);
100106

107+
//add settings for RecoilToTop tool
108+
fMasterGen->settings.addFlag("TopRecoilHook:doTopRecoilIn", false);
109+
fMasterGen->settings.addFlag("TopRecoilHook:useOldDipoleIn", false);
110+
fMasterGen->settings.addFlag("TopRecoilHook:doListIn", false);
111+
101112
//add settings for powheg resonance scale calculation
102113
fMasterGen->settings.addFlag("POWHEGres:calcScales", false);
103114
fMasterGen->settings.addFlag("POWHEG:bb4l", false);
104-
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:onlyDistance1", false);
105115
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:veto", false);
106-
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:dryRun", false);
107-
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoAtPL", false);
108116
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoQED", false);
109-
fMasterGen->settings.addFlag("POWHEG:bb4l:PartonLevel:veto", false);
110-
fMasterGen->settings.addFlag("POWHEG:bb4l:PartonLevel:excludeFSRConflicting", false);
111117
fMasterGen->settings.addFlag("POWHEG:bb4l:DEBUG", false);
112118
fMasterGen->settings.addFlag("POWHEG:bb4l:ScaleResonance:veto", false);
113119
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoDipoleFrame", false);
114120
fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:pTpythiaVeto", false);
115121
fMasterGen->settings.addParm("POWHEG:bb4l:pTminVeto", 10.0, true, true, 0.0, 10.);
122+
fMasterGen->settings.addFlag("POWHEG:bb4l:vetoAllRadtypes", false);
116123

117124
fMasterGen->setRndmEnginePtr(p8RndmEngine_);
118125
fDecayer->setRndmEnginePtr(p8RndmEngine_);

0 commit comments

Comments
 (0)