Skip to content

Commit

Permalink
Intoduce the photon flux into Pythia 8 for UPC studies.
Browse files Browse the repository at this point in the history
  • Loading branch information
Andre Govinda Stahl Leiton committed Feb 28, 2024
1 parent d862fb4 commit 4c0f63c
Show file tree
Hide file tree
Showing 3 changed files with 187 additions and 0 deletions.
48 changes: 48 additions & 0 deletions GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,43 @@ namespace CLHEP {

using namespace gen;

//Insert class for use w/ PDFPtr for proton-photon flux
//parameters hardcoded according to main70.cc in PYTHIA8 v3.10
class Nucleus2gamma2 : public Pythia8::PDF {
public:
// Constructor.
Nucleus2gamma2(int idBeamIn) : Pythia8::PDF(idBeamIn) {}

// Update the photon flux.
void xfUpdate(int, double x, double) override {
// lead
double radius = 0; // radius in [fm]
double z = 0;
if (idBeam == 1000822080) {
radius = 6.636;
z = 82;
}
// oxygen
else if (idBeam == 80160) {
radius = 3.02;
z = 8;
}

// Minimum impact parameter (~2*radius) [fm].
double bmin = 2 * radius;

// Per-nucleon mass for lead.
double m2 = pow2(0.9314);
double alphaEM = 0.007297353080;
double hbarc = 0.197;
double xi = x * sqrt(m2) * bmin / hbarc;
double bK0 = besselK0(xi);
double bK1 = besselK1(xi);
double intB = xi * bK1 * bK0 - 0.5 * pow2(xi) * (pow2(bK1) - pow2(bK0));
xgamma = 2. * alphaEM * pow2(z) / M_PI * intB;
}
};

class Pythia8Hadronizer : public Py8InterfaceBase {
public:
Pythia8Hadronizer(const edm::ParameterSet &params);
Expand Down Expand Up @@ -113,6 +150,11 @@ class Pythia8Hadronizer : public Py8InterfaceBase {
double fBeam1PZ;
double fBeam2PZ;

//PDFPtr for the photonFlux
//Following main70.cc example in PYTHIA8 v3.10
bool doProtonPhotonFlux = false;
Pythia8::PDFPtr photonFlux = nullptr;

//helper class to allow multiple user hooks simultaneously
std::shared_ptr<UserHooksVector> fUserHooksVector;
bool UserHooksSet;
Expand Down Expand Up @@ -182,6 +224,7 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params)
comEnergy(params.getParameter<double>("comEnergy")),
LHEInputFileName(params.getUntrackedParameter<std::string>("LHEInputFileName", "")),
fInitialState(PP),
doProtonPhotonFlux(params.getUntrackedParameter<bool>("doProtonPhotonFlux", false)),
UserHooksSet(false),
nME(-1),
nMEFiltered(-1),
Expand Down Expand Up @@ -369,6 +412,11 @@ bool Pythia8Hadronizer::initializeForInternalPartons() {
fMasterGen->settings.word("Beams:LHEF", lheFile_);
}

if (doProtonPhotonFlux) {
photonFlux = make_shared<Nucleus2gamma2>(1000822080);
fMasterGen->setPhotonFluxPtr(photonFlux, nullptr);
}

if (!fUserHooksVector.get())
fUserHooksVector.reset(new UserHooksVector);
(fUserHooksVector->hooks).clear();
Expand Down
88 changes: 88 additions & 0 deletions GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@
import FWCore.ParameterSet.Config as cms

process = cms.Process("PROD")

process.load("Configuration.StandardSequences.SimulationRandomNumberGeneratorSeeds_cff")

process.source = cms.Source("EmptySource")

process.generator = cms.EDFilter("Pythia8GeneratorFilter",
maxEventsToPrint = cms.untracked.int32(1),
pythiaPylistVerbosity = cms.untracked.int32(1),
filterEfficiency = cms.untracked.double(1.0),
pythiaHepMCVerbosity = cms.untracked.bool(False),
comEnergy = cms.double(5360.),
doProtonPhotonFlux = cms.untracked.bool(True),
#PPbarInitialState = cms.PSet(),
#SLHAFileForPythia8 = cms.string('Configuration/Generator/data/CSA07SUSYBSM_LM9p_sftsdkpyt_slha.out'),
#reweightGen = cms.PSet( # flat in pT
# pTRef = cms.double(15.0),
# power = cms.double(4.5)
#),
#reweightGenRap = cms.PSet( # flat in eta
# yLabSigmaFunc = cms.string("15.44/pow(x,0.0253)-12.56"),
# yLabPower = cms.double(2.),
# yCMSigmaFunc = cms.string("5.45/pow(x+64.84,0.34)"),
# yCMPower = cms.double(2.),
# pTHatMin = cms.double(15.),
# pTHatMax = cms.double(3000.)
#),
#reweightGenPtHatRap = cms.PSet( # flat in Pt and eta
# yLabSigmaFunc = cms.string("15.44/pow(x,0.0253)-12.56"),
# yLabPower = cms.double(2.),
# yCMSigmaFunc = cms.string("5.45/pow(x+64.84,0.34)"),
# yCMPower = cms.double(2.),
# pTHatMin = cms.double(15.),
# pTHatMax = cms.double(3000.)
#),
PythiaParameters = cms.PSet(
pythia8_example02 = cms.vstring('HardQCD:all = on',
'PhaseSpace:pTHatMin = 10.',#CM Edit 20->10
'PhotonParton:all = on',#Added from main70
'MultipartonInteractions:pT0Ref = 3.0',#Added from main70
'PDF:beamA2gamma = on',#Added from main70
#This option below crashes - debug
'PDF:proton2gammaSet = 0',#Added from main70
'PDF:useHardNPDFB = on',
'PDF:gammaFluxApprox2bMin = 13.272',
'PDF:beam2gammaApprox = 2',
'Photon:sampleQ2 = off'
),
parameterSets = cms.vstring('pythia8_example02')
)
)

# in order to use lhapdf PDF add a line like this to pythia8_example02:
# 'PDF:pSet = LHAPDF6:CT10'

process.load("FWCore.MessageLogger.MessageLogger_cfi")
process.MessageLogger = cms.Service("MessageLogger",
cerr = cms.untracked.PSet(
enable = cms.untracked.bool(False)
),
cout = cms.untracked.PSet(
default = cms.untracked.PSet(
limit = cms.untracked.int32(2)
),
enable = cms.untracked.bool(True)
)
)

process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
generator = cms.PSet(
initialSeed = cms.untracked.uint32(123456789),
)
)

process.maxEvents = cms.untracked.PSet(
input = cms.untracked.int32(10)
)

process.GEN = cms.OutputModule("PoolOutputModule",
fileName = cms.untracked.string('pythia8ex2.root')
)

process.p = cms.Path(process.generator)
process.outpath = cms.EndPath(process.GEN)

process.schedule = cms.Schedule(process.p, process.outpath)
51 changes: 51 additions & 0 deletions GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg_Fragment.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
import FWCore.ParameterSet.Config as cms

_generator = cms.EDFilter("Pythia8GeneratorFilter",
maxEventsToPrint = cms.untracked.int32(1),
pythiaPylistVerbosity = cms.untracked.int32(1),
filterEfficiency = cms.untracked.double(1.0),
pythiaHepMCVerbosity = cms.untracked.bool(False),
comEnergy = cms.double(5360.),
doProtonPhotonFlux = cms.untracked.bool(True),
#PPbarInitialState = cms.PSet(),
#SLHAFileForPythia8 = cms.string('Configuration/Generator/data/CSA07SUSYBSM_LM9p_sftsdkpyt_slha.out'),
#reweightGen = cms.PSet( # flat in pT
# pTRef = cms.double(15.0),
# power = cms.double(4.5)
#),
#reweightGenRap = cms.PSet( # flat in eta
# yLabSigmaFunc = cms.string("15.44/pow(x,0.0253)-12.56"),
# yLabPower = cms.double(2.),
# yCMSigmaFunc = cms.string("5.45/pow(x+64.84,0.34)"),
# yCMPower = cms.double(2.),
# pTHatMin = cms.double(15.),
# pTHatMax = cms.double(3000.)
#),
#reweightGenPtHatRap = cms.PSet( # flat in Pt and eta
# yLabSigmaFunc = cms.string("15.44/pow(x,0.0253)-12.56"),
# yLabPower = cms.double(2.),
# yCMSigmaFunc = cms.string("5.45/pow(x+64.84,0.34)"),
# yCMPower = cms.double(2.),
# pTHatMin = cms.double(15.),
# pTHatMax = cms.double(3000.)
#),
PythiaParameters = cms.PSet(
pythia8_example02 = cms.vstring('HardQCD:all = on',
'PhaseSpace:pTHatMin = 10.',#CM Edit 20->10
'PhotonParton:all = on',#Added from main70
'MultipartonInteractions:pT0Ref = 3.0',#Added from main70
'PDF:beamA2gamma = on',#Added from main70
#This option below crashes - debug
'PDF:proton2gammaSet = 0',#Added from main70
'PDF:useHardNPDFB = on',
'PDF:gammaFluxApprox2bMin = 13.272',
'PDF:beam2gammaApprox = 2',
'Photon:sampleQ2 = off'
),
parameterSets = cms.vstring('pythia8_example02')
)
)


from GeneratorInterface.Core.ExternalGeneratorFilter import ExternalGeneratorFilter
generator = ExternalGeneratorFilter(_generator)

0 comments on commit 4c0f63c

Please sign in to comment.