diff --git a/Processes/Pythia/CMakeLists.txt b/Processes/Pythia/CMakeLists.txt
index 3af574c892dd423cdef3b7f018265c9f6f3912af..de1630fd928355c70bf23683408c83757891220f 100644
--- a/Processes/Pythia/CMakeLists.txt
+++ b/Processes/Pythia/CMakeLists.txt
@@ -5,12 +5,14 @@ set (
MODEL_SOURCES
Decay.cc
Random.cc
+ Interaction.cc
)
set (
MODEL_HEADERS
Decay.h
Random.h
+ Interaction.h
)
set (
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8c7d6d7527c6b0545dba8bb3872501ea3509ed71
--- /dev/null
+++ b/Processes/Pythia/Interaction.cc
@@ -0,0 +1,424 @@
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/pythia/Interaction.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/utl/COMBoost.h>
+
+#include <tuple>
+
+using std::cout;
+using std::endl;
+using std::tuple;
+
+using namespace corsika;
+using namespace corsika::setup;
+using Particle = Stack::StackIterator; // ParticleType;
+using Track = Trajectory;
+
+namespace corsika::process::pythia {
+
+ typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+
+ Interaction::Interaction(environment::Environment const& env)
+ : fEnvironment(env) {}
+
+ Interaction::~Interaction() {
+ cout << "Pythia::Interaction n=" << fCount << endl;
+ }
+
+ void Interaction::Init() {
+
+ using random::RNGManager;
+
+ // initialize Pythia
+ if (!fInitialized) {
+
+ fPythia.readString("Print:quiet = on");
+ // TODO: proper process initialization for MinBias needed
+ fPythia.readString("HardQCD:all = on");
+ fPythia.readString("ProcessLevel:resonanceDecays = off");
+
+ fPythia.init();
+
+ // any decays in pythia? if yes need to define which particles
+ if(fInternalDecays){
+ // 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
+ const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Pi0,
+ particles::Code::KMinus, particles::Code::KPlus,
+ particles::Code::K0Long, particles::Code::K0Short,
+ particles::Code::SigmaPlus, particles::Code::SigmaMinus,
+ particles::Code::Lambda0,
+ particles::Code::Xi0, particles::Code::XiMinus,
+ particles::Code::OmegaMinus,
+ particles::Code::DPlus, particles::Code::DMinus, particles::Code::D0, particles::Code::D0Bar};
+
+ Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
+ }
+
+ // basic initialization of cross section routines
+ fSigma.init( &fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm );
+
+ fInitialized = true;
+ }
+ }
+
+ void Interaction::SetParticleListStable(const std::vector<particles::Code> particleList) {
+ for (auto p : particleList)
+ Interaction::SetStable( p );
+ }
+
+ void Interaction::SetUnstable(const particles::Code pCode) {
+ cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
+ fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , true);
+ }
+
+ void Interaction::SetStable(const particles::Code pCode) {
+ cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
+ fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , false);
+ }
+
+
+ void Interaction::ConfigureLabFrameCollision( const particles::Code BeamId, const particles::Code TargetId,
+ const units::si::HEPEnergyType BeamEnergy )
+ {
+ using namespace units::si;
+ // 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>(particles::GetPDG(BeamId));
+ std::stringstream stBeam;
+ stBeam << "Beams:idA = " << pdgBeam;
+ fPythia.readString(stBeam.str());
+ // set target
+ auto pdgTarget = static_cast<int>(particles::GetPDG(TargetId));
+ // replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
+ if(TargetId == particles::Code::Hydrogen)
+ pdgTarget = static_cast<int>( particles::GetPDG(particles::Code::Proton));
+ std::stringstream stTarget;
+ stTarget << "Beams:idB = " << pdgTarget;
+ fPythia.readString(stTarget.str());
+ // set frame to lab. frame
+ fPythia.readString("Beams:frameType = 2");
+ // set beam energy
+ const double Elab = BeamEnergy / 1_GeV;
+ std::stringstream stEnergy;
+ stEnergy << "Beams:eA = " << Elab;
+ fPythia.readString(stEnergy.str());
+ // target at rest
+ fPythia.readString("Beams:eB = 0.");
+ // initialize this config
+ fPythia.init();
+ }
+
+
+ bool Interaction::CanInteract(const corsika::particles::Code pCode)
+ {
+ return pCode == corsika::particles::Code::Proton || pCode == corsika::particles::Code::Neutron
+ || pCode == corsika::particles::Code::AntiProton || pCode == corsika::particles::Code::AntiNeutron
+ || pCode == corsika::particles::Code::PiMinus || pCode == corsika::particles::Code::PiPlus;
+ }
+
+ tuple<units::si::CrossSectionType, units::si::CrossSectionType>
+ Interaction::GetCrossSection(const particles::Code BeamId,
+ const particles::Code TargetId,
+ const units::si::HEPEnergyType CoMenergy) {
+ using namespace units::si;
+
+ // interaction possible in pythia?
+ if( TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen ){
+ if( CanInteract(BeamId) && ValidCoMEnergy(CoMenergy) ){
+ // input particle PDG
+ auto const pdgCodeBeam = static_cast<int>( particles::GetPDG( BeamId ));
+ auto const pdgCodeTarget = static_cast<int>( particles::GetPDG( TargetId ));
+ const double ecm = CoMenergy / 1_GeV;
+
+ // calculate cross section
+ fSigma.calc( pdgCodeBeam, pdgCodeTarget, ecm);
+ if(fSigma.hasSigmaTot()){
+ const double sigEla = fSigma.sigmaEl();
+ const double sigProd = fSigma.sigmaTot() - sigEla;
+
+ return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn);
+
+ } else
+ throw std::runtime_error("pythia cross section init failed");
+
+ } else {
+ return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
+ std::numeric_limits<double>::infinity() * 1_mbarn);
+ }
+ } else {
+ throw std::runtime_error("invalid target for pythia");
+ }
+ }
+
+ template <>
+ units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
+
+ using namespace units;
+ using namespace units::si;
+ using namespace geometry;
+
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ const particles::Code corsikaBeamId = p.GetPID();
+
+ // beam particles for pythia : 1, 2, 3 for p, pi, k
+ // read from cross section code table
+ const bool kInteraction = CanInteract(corsikaBeamId);
+
+ // FOR NOW: assume target is at rest
+ process::pythia::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
+
+ // total momentum and energy
+ HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
+ process::pythia::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
+ pTotLab += p.GetMomentum();
+ pTotLab += pTarget;
+ auto const pTotLabNorm = pTotLab.norm();
+ // calculate cm. energy
+ const HEPEnergyType ECoM = sqrt(
+ (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
+
+ cout << "Interaction: LambdaInt: \n"
+ << " input energy: " << p.GetEnergy() / 1_GeV << endl
+ << " beam can interact:" << kInteraction << endl
+ << " beam pid:" << p.GetPID() << endl;
+
+ // TODO: move limits into variables
+ if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM) ) {
+
+ // get target from environment
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+ */
+ const auto currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // determine average interaction length
+ // weighted sum
+ int i = -1;
+ si::CrossSectionType weightedProdCrossSection = 0_mbarn;
+ // get weights of components from environment/medium
+ const auto w = mediumComposition.GetFractions();
+ // loop over components in medium
+ for (auto const targetId : mediumComposition.GetComponents()) {
+ i++;
+ cout << "Interaction: get interaction length for target: " << targetId << endl;
+
+ auto const [productionCrossSection, elaCrossSection] =
+ GetCrossSection(corsikaBeamId, targetId, ECoM);
+ [[maybe_unused]] auto elaCrossSectionCopy =
+ elaCrossSection; // ONLY TO AVOID COMPILER WARNING
+
+ cout << "Interaction: IntLength: pythia return (mb): "
+ << productionCrossSection / 1_mbarn << endl
+ << "Interaction: IntLength: weight : " << w[i] << endl;
+
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ }
+ cout << "Interaction: IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mbarn << endl
+ << "Interaction: IntLength: average mass number: "
+ << mediumComposition.GetAverageMassNumber() << endl;
+
+ // calculate interaction length in medium
+ GrammageType const int_length =
+ mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection;
+ cout << "Interaction: "
+ << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
+ << endl;
+
+ return int_length;
+ }
+
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+
+ /**
+ In this function PYTHIA is called to produce one event. The
+ event is copied (and boosted) into the shower lab frame.
+ */
+
+ template <>
+ process::EProcessReturn Interaction::DoInteraction(Particle& p, Stack&) {
+
+ using namespace units;
+ using namespace utl;
+ using namespace units::si;
+ using namespace geometry;
+
+
+ const auto corsikaBeamId = p.GetPID();
+ cout << "Pythia::Interaction: "
+ << "DoInteraction: " << corsikaBeamId << " interaction? "
+ << process::pythia::Interaction::CanInteract(corsikaBeamId) << endl;
+
+ if (particles::IsNucleus(corsikaBeamId)) {
+ // nuclei handled by different process, this should not happen
+ throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
+ }
+
+ if (process::pythia::Interaction::CanInteract(corsikaBeamId)) {
+
+ const CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ // position and time of interaction, not used in Sibyll
+ Point pOrig = p.GetPosition();
+ TimeType tOrig = p.GetTime();
+
+ // define target
+ // FOR NOW: target is always at rest
+ const auto eTargetLab = 0_GeV + constants::nucleonMass;
+ const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+ const FourVector PtargLab(eTargetLab, pTargetLab);
+
+ // define projectile
+ HEPEnergyType const eProjectileLab = p.GetEnergy();
+ auto const pProjectileLab = p.GetMomentum();
+
+ cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
+ << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
+ << endl;
+ cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
+ << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV << endl;
+
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Pythia projectile: +z
+ COMBoost const boost(PprojLab, constants::nucleonMass);
+
+ // just for show:
+ // boost projecticle
+ auto const PprojCoM = boost.toCoM(PprojLab);
+
+ // boost target
+ auto const PtargCoM = boost.toCoM(PtargLab);
+
+ cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: pbeam CoM: "
+ << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: ptarget CoM: "
+ << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+
+ cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
+ cout << "Interaction: time: " << tOrig << endl;
+
+ HEPEnergyType Etot = eProjectileLab + eTargetLab;
+ MomentumVector Ptot = p.GetMomentum();
+ // invariant mass, i.e. cm. energy
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+
+ // sample target mass number
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ const auto& mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // get cross sections for target materials
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from GetInteractionLength if possible
+ */
+ //#warning reading interaction cross section again, should not be necessary
+ auto const& compVec = mediumComposition.GetComponents();
+ std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
+
+ for (size_t i = 0; i < compVec.size(); ++i) {
+ auto const targetId = compVec[i];
+ const auto [sigProd, sigEla] =
+ GetCrossSection(corsikaBeamId, targetId, Ecm);
+ cross_section_of_components[i] = sigProd;
+ [[maybe_unused]] auto sigElaCopy =
+ sigEla; // to avoid not used warning in array binding
+ }
+
+ const auto corsikaTargetId =
+ mediumComposition.SampleTarget(cross_section_of_components, fRNG);
+ cout << "Interaction: target selected: " << corsikaTargetId << endl;
+
+ if (corsikaTargetId != particles::Code::Hydrogen && corsikaTargetId != particles::Code::Neutron
+ && corsikaTargetId != particles::Code::Proton )
+ throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
+
+ cout << "Interaction: "
+ << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
+ << " Ecm(GeV): " << Ecm / 1_GeV << endl;
+
+ if (eProjectileLab < 8.5_GeV || !ValidCoMEnergy(Ecm)) {
+ cout << "Interaction: "
+ << " DoInteraction: should have dropped particle.. "
+ << "THIS IS AN ERROR" << endl;
+ throw std::runtime_error("energy too low for PYTHIA");
+
+ } else {
+ fCount++;
+
+ ConfigureLabFrameCollision( corsikaBeamId, corsikaTargetId, eProjectileLab );
+
+ // create event in pytia
+ if(!fPythia.next())
+ throw std::runtime_error("Pythia::DoInteraction: failed!");
+
+ // link to pythia stack
+ Pythia8::Event& event = fPythia.event;
+ // print final state
+ event.list();
+
+ MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
+ for (int i = 0; i < event.size(); ++i){
+ Pythia8::Particle &pp = event[i];
+ // skip particles that have decayed in pythia
+ if (!pp.isFinal()) continue;
+
+ auto const pyId = particles::ConvertFromPDG(static_cast<particles::PDGCode>(pp.id()));
+
+ const MomentumVector pyPlab(rootCS, {pp.px()*1_GeV, pp.py()*1_GeV, pp.pz()*1_GeV});
+ HEPEnergyType const pyEn = pp.e() * 1_GeV;
+
+ // add to corsika stack
+ auto pnew = p.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig, tOrig});
+
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
+ }
+ cout << "conservation (all GeV): " << "Elab_final=" << Elab_final / 1_GeV
+ << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
+ }
+ // delete current particle
+ p.Delete();
+ }
+ return process::EProcessReturn::eOk;
+ }
+
+} // namespace corsika::process::pythia
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
new file mode 100644
index 0000000000000000000000000000000000000000..f384da655322cac969588308c85094358c11cd80
--- /dev/null
+++ b/Processes/Pythia/Interaction.h
@@ -0,0 +1,79 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _corsika_process_pythia_interaction_h_
+#define _corsika_process_pythia_interaction_h_
+
+#include <Pythia8/Pythia.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <tuple>
+
+namespace corsika::environment {
+ class Environment;
+}
+
+namespace corsika::process::pythia {
+
+ class Interaction : public corsika::process::InteractionProcess<Interaction> {
+
+ int fCount = 0;
+ bool fInitialized = false;
+
+ public:
+ Interaction(corsika::environment::Environment const& env);
+ ~Interaction();
+
+ void Init();
+
+ void SetParticleListStable(const std::vector<particles::Code>);
+ void SetUnstable(const corsika::particles::Code );
+ void SetStable(const corsika::particles::Code );
+
+ bool WasInitialized() { return fInitialized; }
+ bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+ using namespace corsika::units::si;
+ return (10_GeV < ecm) && (ecm < 1_PeV);
+ }
+
+ bool CanInteract(const corsika::particles::Code);
+ void ConfigureLabFrameCollision(const corsika::particles::Code, const corsika::particles::Code,
+ const corsika::units::si::HEPEnergyType);
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
+ GetCrossSection(const corsika::particles::Code BeamId,
+ const corsika::particles::Code TargetId,
+ const corsika::units::si::HEPEnergyType CoMenergy);
+
+ template <typename Particle, typename Track>
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
+
+ /**
+ In this function PYTHIA is called to produce one event. The
+ event is copied (and boosted) into the shower lab frame.
+ */
+
+ template <typename Particle, typename Stack>
+ corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
+
+ private:
+ corsika::environment::Environment const& fEnvironment;
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
+ Pythia8::Pythia fPythia;
+ Pythia8::SigmaTotal fSigma;
+ const bool fInternalDecays = true;
+ };
+
+} // namespace corsika::process::pythia
+
+#endif
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 05895edaca71bc5a54ad2dccebd19fc7cdd883dc..9fb2d757a33f1922a8567b1c4b8e8df4dbc2fb7d 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -10,6 +10,7 @@
*/
#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/pythia/Interaction.h>
#include <Pythia8/Pythia.h>
#include <corsika/random/RNGManager.h>
@@ -97,7 +98,7 @@ TEST_CASE("Pythia", "[processes]") {
using namespace corsika;
using namespace corsika::units::si;
-TEST_CASE("pytia decay"){
+TEST_CASE("pythia process"){
// setup environment, geometry
environment::Environment env;
@@ -111,7 +112,7 @@ TEST_CASE("pytia decay"){
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+ std::vector<particles::Code>{particles::Code::Hydrogen}, std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
@@ -154,4 +155,27 @@ TEST_CASE("pytia decay"){
}
+ SECTION("pythia interaction") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 100_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
+ auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ geometry::Point pos(cs, 0_m, 0_m, 0_m);
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::PiPlus, E0, plab, pos, 0_ns});
+
+
+ process::pythia::Interaction model(env);
+
+ model.Init();
+ /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoInteraction(particle,
+ stack);
+ [[maybe_unused]] const GrammageType length =
+ model.GetInteractionLength(particle, track);
+ }
+
}