diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 38ca6094df214d15d5c1da0853358314cdd10e5a..65f909b01fe06208a1344510426c5b018f0544f6 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -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): "
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index 6385430f9a0046dab9e4284b41ec7a3688f47545..c9ebe4965fd1a0cd170020e8aca4d0227f3a9233 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -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
diff --git a/modules/pythia8/CMakeLists.txt b/modules/pythia8/CMakeLists.txt
index 06edb532098c0d54e8c48e48793b63e4b29ab781..1f6449a81f29e45b55ff0d1d5e001d9e25a33481 100644
--- a/modules/pythia8/CMakeLists.txt
+++ b/modules/pythia8/CMakeLists.txt
@@ -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
diff --git a/modules/pythia8/pythia8245-stripped.tar.bz2 b/modules/pythia8/pythia8245-stripped.tar.bz2
deleted file mode 100644
index 0ad96435f488cac464393fee13b4d5f4a38794fd..0000000000000000000000000000000000000000
Binary files a/modules/pythia8/pythia8245-stripped.tar.bz2 and /dev/null differ
diff --git a/modules/pythia8/pythia8307.tgz b/modules/pythia8/pythia8307.tgz
new file mode 120000
index 0000000000000000000000000000000000000000..c0a689fb4c46b5ddc653e4d08bc6ac3245ff8c06
--- /dev/null
+++ b/modules/pythia8/pythia8307.tgz
@@ -0,0 +1 @@
+/cr/users/reininghaus/Downloads/pythia8307.tgz
\ No newline at end of file
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index f1b9892a51412ca963d0636726c97a449a1e00c7..6f43f003cc79ee342553162158d03092f9db098f 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -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") {