Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[13_0_X] Pythia8 photon flux #44258

Merged
merged 1 commit into from
Mar 7, 2024
Merged

[13_0_X] Pythia8 photon flux #44258

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 48 additions & 0 deletions GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,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 @@ -114,6 +151,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 @@ -184,6 +226,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 @@ -371,6 +414,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
Open in desktop
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
Open in desktop
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)