From 4c0f63ce89ba232f9ac891bf2c772d2eac7380f5 Mon Sep 17 00:00:00 2001 From: Andre Govinda Stahl Leiton Date: Thu, 29 Feb 2024 00:52:15 +0100 Subject: [PATCH] Intoduce the photon flux into Pythia 8 for UPC studies. --- .../plugins/Pythia8Hadronizer.cc | 48 ++++++++++ .../test/pythia8ex2_ConvertToMain70_cfg.py | 88 +++++++++++++++++++ ...pythia8ex2_ConvertToMain70_cfg_Fragment.py | 51 +++++++++++ 3 files changed, 187 insertions(+) create mode 100644 GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg.py create mode 100644 GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg_Fragment.py diff --git a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc index b274ecaeec290..9cc53a808b953 100644 --- a/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc +++ b/GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc @@ -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 ¶ms); @@ -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 fUserHooksVector; bool UserHooksSet; @@ -182,6 +224,7 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet ¶ms) comEnergy(params.getParameter("comEnergy")), LHEInputFileName(params.getUntrackedParameter("LHEInputFileName", "")), fInitialState(PP), + doProtonPhotonFlux(params.getUntrackedParameter("doProtonPhotonFlux", false)), UserHooksSet(false), nME(-1), nMEFiltered(-1), @@ -369,6 +412,11 @@ bool Pythia8Hadronizer::initializeForInternalPartons() { fMasterGen->settings.word("Beams:LHEF", lheFile_); } + if (doProtonPhotonFlux) { + photonFlux = make_shared(1000822080); + fMasterGen->setPhotonFluxPtr(photonFlux, nullptr); + } + if (!fUserHooksVector.get()) fUserHooksVector.reset(new UserHooksVector); (fUserHooksVector->hooks).clear(); diff --git a/GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg.py b/GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg.py new file mode 100644 index 0000000000000..579672d942cc2 --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg.py @@ -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) diff --git a/GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg_Fragment.py b/GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg_Fragment.py new file mode 100644 index 0000000000000..98d9c6105bdb7 --- /dev/null +++ b/GeneratorInterface/Pythia8Interface/test/pythia8ex2_ConvertToMain70_cfg_Fragment.py @@ -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)