IAP GITLAB

Skip to content
Snippets Groups Projects
Commit d586c7e1 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

preliminary, compiling version of Pythia 8.307 with simplified nuclear 

target
parent fad486a2
No related branches found
No related tags found
1 merge request!427Resolve "upgrade pythia to version 8.3xx"
Hide whitespace changes
Inline Side-by-side
corsika/detail/modules/pythia8/Interaction.inl
+ 261
− 125
View file @ d586c7e1
......@@ -10,6 +10,7 @@
#include <tuple>
#include <boost/filesystem/path.hpp>
#include <fmt/core.h>
#include <corsika/modules/pythia8/Interaction.hpp>
......@@ -25,100 +26,109 @@ namespace corsika::pythia8 {
CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_);
}
inline Interaction::Interaction(bool const print_listing)
: Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
, print_listing_(print_listing) {
CORSIKA_LOG_INFO("Configuring Pythia8 from: {}", CORSIKA_Pythia8_XML_DIR);
// initialize Pythia
// reduce output from pythia if set to "on"
Pythia8::Pythia::readString("Print:quiet = on");
// check if data in particle data file is minimally consistent. Very verbose! set to
// "off"! we do not change the basic file provided by pythia.
Pythia8::Pythia::readString("Check:particleData = off");
Pythia8::Pythia::readString("Check:event = on"); // default: on
Pythia8::Pythia::readString("Check:levelParticleData = 12"); // 1 is default
readString("SoftQCD:all = on");
readString("LowEnergyQCD:all = on");
readString("MultipartonInteraction:reuseInit = 3");
readString(fmt::format("MultipartonInteraction:initFile = {}", corsika_data("pythia8/MPI_init_file")));
readString("Beam:allowVariableEnergy = on");
readString("Beam:allowIDASwitch = on");
readString("Beam:frameType = 3"); // no special frame, need to specify full 3-mom
readString("Check:epTolErr = 0.1"); // relax energy-monetum conservation, copied from Pythia 8 main183.cc
// we can't test this block, LCOV_EXCL_START
if (!Pythia8::Pythia::init()) {
throw std::runtime_error("Pythia::Interaction: Initialization failed!");
}
// LCOV_EXCL_STOP
// any decays in pythia? if yes need to define which particles
if (internalDecays_) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
std::vector<Code> const HadronsWeWantTrackedByCorsika = {
Code::PiPlus, Code::PiMinus, Code::Pi0, Code::KMinus, Code::KPlus,
Code::K0Long, Code::K0Short, Code::SigmaPlus, Code::SigmaMinus, Code::Lambda0,
Code::Xi0, Code::XiMinus, Code::OmegaMinus, Code::DPlus, Code::DMinus,
Code::D0, Code::D0Bar};
Interaction::setStable(HadronsWeWantTrackedByCorsika);
}
inline Interaction::Interaction(boost::filesystem::path const& mpiInitFile, bool const print_listing)
: print_listing_(print_listing) {
// Main Pythia object for managing the cascade evolution.
// Can also do decays, but no hard processes.
pythiaMain_.readString("ProcessLevel:all = off");
pythiaMain_.readString("211:mayDecay = on"); // TODO: probably not necessary, but ask Torjbörn maybe
pythiaMain_.readString("13:mayDecay = on");
pythiaMain_.readString("321:mayDecay = on");
pythiaMain_.readString("130:mayDecay = on");
// Redure statistics printout to relevant ones.
pythiaMain_.readString("Stat:showProcessLevel = off");
pythiaMain_.readString("Stat:showPartonLevel = off");
// Add Argon, since not in Particle data. id:all = name antiName
// spinType chargeType colType m0 mWidth mMin mMax tau0.
pythiaMain_.readString("1000180400:all = 40Ar 40Arbar 1 54 0 "
"37.22474 0. 0. 0. 0.");
if (!pythiaMain_.init())
throw std::runtime_error("Pythia::Interaction: Initialization failed!");
// TODO: do we need this?
//~ CORSIKA_LOG_INFO("Configuring Pythia8 from: {}", CORSIKA_Pythia8_XML_DIR);
// Secondary Pythia object for performing individual collisions.
// Variable incoming beam type and energy.
pythiaColl_.readString("Beams:allowVariableEnergy = on");
pythiaColl_.readString("Beams:allowIDAswitch = on");
// Set up for fixed-target collisions.
pythiaColl_.readString("Beams:frameType = 3"); // arbitrary frame, need to define full 4-momenta
pythiaColl_.settings.parm("Beams:pzA", eMaxLab_ / 1_GeV);
pythiaColl_.readString("Beams:pzB = 0.");
// Must use the soft and low-energy QCD processes.
pythiaColl_.readString("SoftQCD:all = on");
pythiaColl_.readString("LowEnergyQCD:all = on");
// Decays to be done by pythiaMain_.
pythiaColl_.readString("HadronLevel:Decay = off");
// Reduce printout and relax energy-momentum conservation.
pythiaColl_.readString("Print:quiet = on");
pythiaColl_.readString("Check:epTolErr = 0.1");
// Redure statistics printout to relevant ones.
pythiaColl_.readString("Stat:showProcessLevel = off");
pythiaColl_.readString("Stat:showPartonLevel = off");
bool const reuse = true; // TODO: make this more flexible
// Reuse MPI initialization file if it exists; else create a new one.
if (reuse) pythiaColl_.readString("MultipartonInteractions:reuseInit = 3");
else pythiaColl_.readString("MultipartonInteractions:reuseInit = 1");
pythiaColl_.settings.word("MultipartonInteractions:initFile", mpiInitFile.native());
// initialize
// we can't test this block, LCOV_EXCL_START
if (!pythiaColl_.init())
throw std::runtime_error("Pythia::Interaction: Initialization failed!");
// LCOV_EXCL_STOP
//~ // any decays in pythia? if yes need to define which particles
//~ if (internalDecays_) {
//~ // define which particles are passed to corsika, i.e. which particles make it into
//~ // history even very shortlived particles like charm or pi0 are of interest here
//~ std::vector<Code> const HadronsWeWantTrackedByCorsika = {
//~ Code::PiPlus, Code::PiMinus, Code::Pi0, Code::KMinus, Code::KPlus,
//~ Code::K0Long, Code::K0Short, Code::SigmaPlus, Code::SigmaMinus, Code::Lambda0,
//~ Code::Xi0, Code::XiMinus, Code::OmegaMinus, Code::DPlus, Code::DMinus,
//~ Code::D0, Code::D0Bar};
//~ Interaction::setStable(HadronsWeWantTrackedByCorsika);
//~ }
}
inline void Interaction::setStable(std::vector<Code> const& particleList) {
for (auto p : particleList) Interaction::setStable(p);
}
//~ inline void Interaction::setStable(std::vector<Code> const& particleList) {
//~ for (auto p : particleList) Interaction::setStable(p);
//~ }
inline void Interaction::setUnstable(Code const pCode) {
CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} unstable..", pCode);
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
}
//~ inline void Interaction::setUnstable(Code const pCode) {
//~ CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} unstable..", pCode);
//~ Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
//~ }
inline void Interaction::setStable(Code const pCode) {
CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} stable..", pCode);
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
//~ inline void Interaction::setStable(Code const pCode) {
//~ CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} stable..", pCode);
//~ Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
//~ }
inline bool Interaction::isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtS) const {
// TODO: rethink this function. Pythia can do much more now.
return true;
}
inline void Interaction::configureLabFrameCollision(Code const projectileId,
Code const targetId,
MomentumVector const& projectileMom,
MomentumVector const& targetMom) {
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
// set beam
// beam id for pythia
auto const pdgBeam = static_cast<int>(get_PDG(projectileId));
// set target
auto pdgTarget = static_cast<int>(get_PDG(targetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
if (targetId == Code::Hydrogen)
pdgTarget = static_cast<int>(get_PDG(Code::Proton));
if (is_nucleus(projectileId)) // not yet possible with Pythia
return false;
// TODO: check sqrtS for validity
setBeamIDs(pdgBeam, pdgTarget);
auto const projMomGeV = projectileMom.getComponents() * (1/1_GeV);
auto const tarMomGeV = targetMom.getComponents() * (1/1_GeV);
setKinematics(projMomGeV[0], projMomGeV[1], projMomGeV[2],
tarMomGeV[0], tarMomGeV[1], tarMomGeV[2]);
// initialize this config
// we can't test this block, LCOV_EXCL_START
if (!Pythia8::Pythia::init())
throw std::runtime_error("Pythia::Interaction: Initialization failed!");
// LCOV_EXCL_STOP
return std::find(validTargets_.begin(), validTargets_.end(), targetId) != validTargets_.end();
}
inline bool Interaction::canInteract(Code const pCode) const {
......@@ -128,12 +138,14 @@ namespace corsika::pythia8 {
CrossSectionType
Interaction::getCrossSection(Code const projectileId, Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
FourMomentum const& targetP4) const
{
HEPEnergyType const CoMenergy = (projectileP4 + targetP4).getNorm();
HEPEnergyType const CoMenergy = is_nucleus(targetId) ?
(projectileP4 + targetP4 / get_nucleus_A(targetId)).getNorm() : (projectileP4 + targetP4).getNorm();
if (!isValid(projectileId, targetId, CoMenergy)) {
return {CrossSectionType::zero(), CrossSectionType::zero()};
return CrossSectionType::zero();
}
// input particle PDG
......@@ -141,21 +153,36 @@ namespace corsika::pythia8 {
auto const pdgCodeTarget = static_cast<int>(get_PDG(targetId));
double const ecm_GeV = CoMenergy * (1 / 1_GeV);
auto const sigTot = pythia.getSigmaTotal(pdgCodeBeam, 2212, ecm_GeV) * 1_mb;
auto const sigEla = pythia.getSigmaPartial(pdgCodeBeam, 2212, evm_GeV, 2) * 1_mb;
auto const sigProd = sigTot - sigEla;
if (is_nucleus(targetId)) {
// avg. no. of sub-collisions, from 2108.03481, for Nitrogen target...
// TODO: generalize
double const nCollAvg = (sigTot < 31._mb) ? 1. + (0.017 / 1_mb) * sigTot
: 1.2 + (0.0105 / 1_mb) * sigTot;
// TODO: is this valid also for production XS? Paper says total XS...
return get_nucleus_A(targetId) * sigProd / nCollAvg;
auto const sigTot = pythiaColl_.getSigmaTotal(pdgCodeBeam, 2212, ecm_GeV) * 1_mb;
auto const sigEla = pythiaColl_.getSigmaPartial(pdgCodeBeam, 2212, ecm_GeV, 2) * 1_mb;
//~ auto const sigProd = sigTot - sigEla;
// TODO: handle hydrogen as proton
if (!is_nucleus(targetId)) {
return sigTot;
} else {
return sigProd;
return sigTot * get_nucleus_A(targetId) / getAverageSubcollisions(targetId, sigTot);
}
}
double Interaction::getAverageSubcollisions(Code targetId, CrossSectionType sigTot) const {
auto const Z = get_nucleus_Z(targetId);
auto const A = get_nucleus_A(targetId);
double nCollAvg;
if (Z == 7 && A == 14) nCollAvg = (sigTot < 31._mb) // this is from the paper
? 1. + 0.017/1_mb * sigTot : 1.2 + 0.0105/1_mb * sigTot;
else if (Z == 8 && A == 16) nCollAvg = (sigTot < 16._mb) // provided by Torbjörn Sjöstrand
? 1. + 0.0245/1_mb * sigTot : 1.2 + 0.012/1_mb * sigTot;
else if (Z == 18 && A == 40) nCollAvg = (sigTot < 10._mb)
? 1. + 0.050/1_mb * sigTot : 1.28 + 0.022/1_mb * sigTot;
else {
CORSIKA_LOG_ERROR(fmt::format("Pythia8 cross-sections not defined for ({},{}) nucleus", A, Z));
nCollAvg = 0;
}
return nCollAvg;
}
template <class TView>
......@@ -187,7 +214,6 @@ namespace corsika::pythia8 {
TimeType const tOrig = projectile.getTime();
CORSIKA_LOG_DEBUG("Interaction: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
CORSIKA_LOG_DEBUG("Interaction: position of interaction: ", pOrig.getCoordinates());
CORSIKA_LOG_DEBUG("Interaction: time: {}", tOrig);
......@@ -198,36 +224,144 @@ namespace corsika::pythia8 {
eProjectileLab / 1_GeV, sqrtS / 1_GeV);
count_++;
configureLabFrameCollision(projectileId, targetId, projectileP4.getSpacelikeComponents(), targetP4.getSpacelikeComponents());
CrossSectionType const sigmaHp = pythiaColl.getSigmaTotal(idNow, 2212, eCMNow);
// create event in pytia. LCOV_EXCL_START: we don't validate pythia8 internals
if (!Pythia8::Pythia::next())
throw std::runtime_error("Pythia::DoInteraction: failed!");
// LCOV_EXCL_STOP
// LCOV_EXCL_START, we don't validate pythia8 internals
if (print_listing_) {
// print final state
event.list();
}
// LCOV_EXCL_STOP
MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
for (int i = 0; i < event.size(); ++i) {
Pythia8::Particle const& p8p = event[i];
int const idNow = static_cast<int>(get_PDG(projectileId));
// References to the two event records. Clear main event record.
Pythia8::Event& eventMain = pythiaMain_.event;
Pythia8::Event& eventColl = pythiaColl_.event;
eventMain.clear();
COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
auto const proj4MomLab = labFrameBoost.toCoM(projectileP4);
auto const& rotCS = labFrameBoost.getRotatedCS();
auto const pProjLab = proj4MomLab.getSpaceLikeComponents().getComponents(rotCS);
auto constexpr invGeV = 1 / 1_GeV;
Pythia8::Vec4 const pNow{pProjLab.getX() * invGeV, pProjLab.getY() * invGeV, pProjLab.getZ() * invGeV, proj4MomLab.getTimeLikeComponent() * invGeV};
// Insert incoming particle in cleared main event record.
eventMain.append(90, -11, 0, 0, 1, 1, 0, 0, pNow, pNow.mCalc());
int const iHad = eventMain.append(idNow, 12, 0, 0, 0, 0, 0, 0, pNow, projectile.getMass() * (1/1_GeV));
Pythia8::Vec4 const vNow{}; // production vertex; useless but necessary (?) TODO: ask TS
eventMain[iHad].vProd(vNow);
auto const [Anow, Znow] = std::invoke([targetId]() {
if (targetId == Code::Proton) {
return std::make_pair(1, 1);
} else if (targetId == Code::Neutron) {
return std::make_pair(1, 0);
} else if (is_nucleus(targetId)) {
return std::make_pair<int, int>(get_nucleus_A(targetId), get_nucleus_Z(targetId));
} else {
return std::make_pair(0, 0); //TODO: what to do in this case?
}
});
// Set up for collisions on a nucleus.
int np = Znow;
int nn = Anow - Znow;
int sizeOld = 0;
int sizeNew = 0;
Pythia8::Vec4 const dirNow = pNow / pNow.pAbs();
Pythia8::Rndm& rndm = pythiaMain_.rndm;
double constexpr mp = get_mass(Code::Proton) / 1_GeV;
double const sqrtSNN_GeV = (pNow + Pythia8::Vec4{0, 0, 0, mp}).mCalc();
auto const nCollAvg = getAverageSubcollisions(targetId, pythiaColl_.getSigmaTotal(idNow, 2212, sqrtSNN_GeV) * 1_mb);
double const probMore = 1. - 1. / nCollAvg;
// Loop over varying number of hit nucleons in target nucleus.
for (int iColl = 1; iColl <= Anow; ++iColl) {
if (iColl > 1 && rndm.flat() > probMore) break;
// Pick incoming projectile: trivial for first subcollision, else ...
int iProj = iHad;
int procType = 0;
// ... find highest-pLongitudinal particle from latest subcollision.
if (iColl > 1) {
iProj = 0;
double pMax = 0.;
for (int i = sizeOld; i < sizeNew; ++i)
if ( eventMain[i].isFinal() && eventMain[i].isHadron()) {
double const pp = Pythia8::dot3(dirNow, eventMain[i].p());
if (pp > pMax) {
iProj = i;
pMax = pp;
}
}
// No further subcollision if no particle with enough energy.
if ( iProj == 0 || eventMain[iProj].e() - eventMain[iProj].m()
< eKinMinLab_/1_GeV) break;
// Choose process; only SD or ND at perturbative energies.
double const eCMSub = (eventMain[iProj].p() + Pythia8::Vec4{0, 0, 0, mp}).mCalc();
if (eCMSub > 10.) procType = (rndm.flat() < probSD_) ? 4 : 1;
}
// Pick one p or n from target.
int const idProj = eventMain[iProj].id();
bool const doProton = rndm.flat() < (np / double(np + nn));
if (doProton) np--;
else nn--;
int const idNuc = doProton ? 2212 : 2112;
// Perform the projectile-nucleon subcollision.
pythiaColl_.setBeamIDs(idProj, idNuc);
pythiaColl_.setKinematics(eventMain[iProj].p(), Pythia8::Vec4());
pythiaColl_.next(procType);
// Insert target nucleon. Mothers are (0,iProj) to mark who it
// interacted with. Always use proton mass for simplicity.
int const statusNuc = (iColl == 1) ? -181 : -182;
int const iNuc = eventMain.append( idNuc, statusNuc, 0, iProj, 0, 0, 0, 0,
0., 0., 0., mp, mp);
eventMain[iNuc].vProdAdd(vNow);
// Insert secondary produced particles (but skip intermediate partons)
// into main event record and shift to correct production vertex.
sizeOld = eventMain.size();
for (int iSub = 3; iSub < eventColl.size(); ++iSub) {
if (!eventColl[iSub].isFinal()) continue;
int iNew = eventMain.append(eventColl[iSub]);
eventMain[iNew].mothers(iNuc, iProj);
eventMain[iNew].vProdAdd(vNow);
}
sizeNew = eventMain.size();
// Update daughters of colliding hadrons and other history.
eventMain[iProj].daughters(sizeOld, sizeNew - 1);
eventMain[iNuc].daughters(sizeOld, sizeNew - 1);
eventMain[iProj].statusNeg();
double dTau = (iColl == 1) ? (vNow.e() - eventMain[iHad].tProd())
* eventMain[iHad].m() / eventMain[iHad].e() : 0.;
eventMain[iProj].tau(dTau);
// End of loop over interactions in a nucleus.
}
MomentumVector Plab_final{labFrameBoost.getOriginalCS()};
auto Elab_final = HEPEnergyType::zero();
for (Pythia8::Particle const& p8p: eventColl) {
// skip particles that have decayed in pythia
if (!p8p.isFinal()) continue;
auto const pyId = convert_from_PDG(static_cast<PDGCode>(p8p.id()));
MomentumVector const pyPlab(labCS,
MomentumVector const pyPlab(rotCS,
{p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
auto const pyP = labFrameBoost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPlab});
HEPEnergyType const mass = get_mass(pyId);
HEPEnergyType const Ekin = sqrt(pyPlab.getSquaredNorm() + mass * mass) - mass;
HEPEnergyType const Ekin = sqrt(pyP.getSpaceLikeComponents().getSquaredNorm() + mass * mass) - mass;
// add to corsika stack
auto pnew =
......@@ -236,6 +370,8 @@ namespace corsika::pythia8 {
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
eventMain.clear();
CORSIKA_LOG_DEBUG(
"conservation (all GeV): "
......
corsika/modules/pythia8/Interaction.hpp
+ 51
− 30
View file @ d586c7e1
......@@ -8,54 +8,57 @@
#pragma once
#include <tuple>
#include <boost/filesystem/path.hpp>
#include <corsika/framework/utility/CorsikaData.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/modules/pythia8/Pythia8.hpp>
#include <tuple>
namespace corsika::pythia8 {
class Interaction : public InteractionProcess<Interaction>, public Pythia8::Pythia {
class Interaction : public InteractionProcess<Interaction> {
public:
Interaction(bool const print_listing = false);
Interaction(boost::filesystem::path const& mpiInitFile = corsika_data("Pythia/main184.mpi"),
bool const print_listing = false);
~Interaction();
void setStable(std::vector<Code> const&);
void setUnstable(Code const);
void setStable(Code const);
//~ void setStable(std::vector<Code> const&);
//~ void setUnstable(Code const);
//~ void setStable(Code const);
bool isValidCoMEnergy(HEPEnergyType const ecm) const {
return (10_GeV < ecm) && (ecm < 1_PeV);
}
bool canInteract(Code const) const;
void configureLabFrameCollision(Code const, Code const, HEPEnergyType const);
//~ void configureLabFrameCollision(Code const, Code const, HEPEnergyType const);
bool isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtS) const;
/**
* Returns inelastic AND elastic cross sections.
*
* These cross sections must correspond to the process described in doInteraction
* AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
* nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
* Sibyll must be useful inside the NuclearInteraction model, which requires that.
*
* @param projectile is the Code of the projectile
* @param target is the Code of the target
* @param sqrtSnn is the center-of-mass energy (per nucleon pair)
* @param Aprojectil is the mass number of the projectils, if it is a nucleus
* @param Atarget is the mass number of the target, if it is a nucleus
*
* @return a tuple of: inelastic cross section, elastic cross section
*/
std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
Code const projectile, Code const target, FourMomentum const& projectileP4,
FourMomentum const& targetP4) const;
//~ /**
//~ * Returns inelastic AND elastic cross sections.
//~ *
//~ * These cross sections must correspond to the process described in doInteraction
//~ * AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
//~ * nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
//~ * Sibyll must be useful inside the NuclearInteraction model, which requires that.
//~ *
//~ * @param projectile is the Code of the projectile
//~ * @param target is the Code of the target
//~ * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
//~ * @param Aprojectil is the mass number of the projectils, if it is a nucleus
//~ * @param Atarget is the mass number of the target, if it is a nucleus
//~ *
//~ * @return a tuple of: inelastic cross section, elastic cross section
//~ */
//~ std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
//~ Code const projectile, Code const target, FourMomentum const& projectileP4,
//~ FourMomentum const& targetP4) const;
/**
* Returns inelastic (production) cross section.
......@@ -74,7 +77,7 @@ namespace corsika::pythia8 {
*/
CrossSectionType getCrossSection(Code const projectile, Code const target,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const;
FourMomentum const& targetP4) const; // non-const for now
/**
* In this function PYTHIA is called to produce one event. The
......@@ -83,12 +86,30 @@ namespace corsika::pythia8 {
template <typename TView>
void doInteraction(TView& output, Code const projectileId, Code const targetId,
FourMomentum const& projectileP4, FourMomentum const& targetP4);
/**
* return average number of sub-collisions in a nucleus, using the parameterizations of
* Sjöstrand and Utheim, EPJ C 82 (2022) 1, 21, arXiv:2108.03481 [hep-ph]
*
* @param targetId target (nucleus
* @param sigTot projectile-nucleon (=proton) cross-secion
*/
double getAverageSubcollisions(Code targetId, CrossSectionType sigTot) const;
static std::array constexpr validTargets_{Code::Oxygen, Code::Nitrogen, Code::Argon, Code::Hydrogen, Code::Proton, Code::Neutron};
private:
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
bool const internalDecays_ = true;
//~ bool const internalDecays_ = true;
int count_ = 0;
bool print_listing_ = false;
bool const print_listing_ = false;
Pythia8::Pythia pythiaMain_;
Pythia8::Pythia mutable pythiaColl_; // Pythia's getSigma...() are not marked const...
double const probSD_ = 0.3; // Fraction of single diffractive events beyond first collision in nucleus.
HEPEnergyType const eMaxLab_ = 1e18_eV;
HEPEnergyType const eKinMinLab_ = 0.2_GeV;
};
} // namespace corsika::pythia8
......
modules/pythia8/CMakeLists.txt
+ 4
− 4
View file @ d586c7e1
......@@ -24,7 +24,7 @@ file (MAKE_DIRECTORY ${CORSIKA_Pythia8_MODULE_DIR})
add_library (C8::ext::pythia8 STATIC IMPORTED GLOBAL)
if ("x_${USE_Pythia8_C8}" STREQUAL "x_SYSTEM")
find_package (Pythia8 8245 EXACT REQUIRED)
find_package (Pythia8 8307 EXACT REQUIRED)
message (STATUS "Using system-level Pythia8 version ${Pythia8_VERSION} at ${Pythia8_PREFIX}")
set (Pythia8_INCLUDE_DIRS ${Pythia8_INCLUDE_DIR})
set_target_properties (
......@@ -50,14 +50,14 @@ if ("x_${USE_Pythia8_C8}" STREQUAL "x_SYSTEM")
else ()
set (_C8_Pythia8_VERSION "8245")
set (_C8_Pythia8_VERSION "8307")
message (STATUS "Building modules/pythia8 using pythia${_C8_Pythia8_VERSION}-stripped.tar.bz2")
message (STATUS "This will take a bit.....")
# this is not a full Pythia8 install, it is a bit simplified, e.g. no pythia8-config, no examples
# and also, it is fitted into the normal CORSIKA8 install "include/corsika_modules/Pythia8, lib/corsika"
ExternalProject_Add (pythia8
URL ${CMAKE_CURRENT_SOURCE_DIR}/pythia${_C8_Pythia8_VERSION}-stripped.tar.bz2
URL_MD5 d3e951a2f101e8cfec26405cb61db83b
URL ${CMAKE_CURRENT_SOURCE_DIR}/pythia${_C8_Pythia8_VERSION}.tgz
URL_MD5 4ee7aea85cc66db13e92722f150c51c5
SOURCE_DIR "${CMAKE_CURRENT_BINARY_DIR}/pythia8/source"
INSTALL_DIR "${CMAKE_CURRENT_BINARY_DIR}/pythia8/install"
CONFIGURE_COMMAND ./configure --prefix=${CMAKE_CURRENT_BINARY_DIR}/pythia8/install
......
modules/pythia8/pythia8245-stripped.tar.bz2 deleted 100644 → 0
+ 0
− 0
View file @ fad486a2
File deleted
modules/pythia8/pythia8307.tgz 0 → 120000
+ 1
− 0
View file @ d586c7e1
/cr/users/reininghaus/Downloads/pythia8307.tgz
\ No newline at end of file
tests/modules/testPythia8.cpp
+ 14
− 14
View file @ d586c7e1
......@@ -89,7 +89,7 @@ using namespace corsika;
template <typename TStackView>
auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
MomentumVector sum{vCS, 0_eV, 0_eV, 0_eV};
MomentumVector sum{vCS};
for (auto const& p : view) { sum += p.getMomentum(); }
return sum;
}
......@@ -182,21 +182,21 @@ TEST_CASE("Pythia8Interface", "modules") {
CHECK(collision.canInteract(Code::PiPlus));
CHECK_FALSE(collision.canInteract(Code::Electron));
// nuclei not supported
std::tuple<CrossSectionType, CrossSectionType> xs_test =
collision.getCrossSectionInelEla(
Code::Proton, Code::Hydrogen,
{sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
{rootCS, {0_eV, 0_eV, 100_GeV}}},
{Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
CHECK(std::get<0>(xs_test) / 1_mb == Approx(314).margin(2));
CHECK(std::get<1>(xs_test) / 1_mb == Approx(69).margin(2));
collision.doInteraction(view, Code::Proton, Code::Hydrogen,
//~ // nuclei not supported
//~ std::tuple<CrossSectionType, CrossSectionType> xs_test =
//~ collision.getCrossSectionInelEla(
//~ Code::Proton, Code::Hydrogen,
//~ {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
//~ {rootCS, {0_eV, 0_eV, 100_GeV}}},
//~ {Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
//~ CHECK(std::get<0>(xs_test) / 1_mb == Approx(314).margin(2));
//~ CHECK(std::get<1>(xs_test) / 1_mb == Approx(69).margin(2));
collision.doInteraction(view, Code::Proton, Code::Nitrogen,
{sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
{rootCS, {0_eV, 0_eV, 100_GeV}}},
{Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
CHECK(view.getSize() == 12);
{Nitrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
CHECK(view.getSize() > 0);
}
SECTION("pythia too low energy") {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum