diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 06b4467e4d89fc2576b5f6b031738e441d890a60..943bbdd2b3cf67d2e26e0427f896ab3f58fcba2a 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -45,6 +45,27 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example)
+add_executable (cascade_proton_example cascade_proton_example.cc)
+target_compile_options(cascade_proton_example PRIVATE -g) # do not skip asserts
+target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ ProcessPythia
+ CORSIKAcascade
+ ProcessStackInspector
+ ProcessTrackWriter
+ ProcessTrackingLine
+ ProcessHadronicElasticModel
+ CORSIKAprocesses
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+install (TARGETS cascade_proton_example DESTINATION share/examples)
+CORSIKA_ADD_TEST (cascade_proton_example)
+
add_executable (staticsequence_example staticsequence_example.cc)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
target_link_libraries (staticsequence_example
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 764eeeba0cd920bf2fa137d6b2f0afb99e7072b8..c154db988f183799d3c18f3dc8033bd8f600c2d3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -242,9 +242,8 @@ int main() {
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Nitrogen,
- particles::Code::Oxygen},
- std::vector<float>{(float)1. - fox, fox}));
+ std::vector<particles::Code>{particles::Code::Nitrogen, particles::Code::Oxygen},
+ std::vector<float>{(float)1.-fox, fox}));
universe.AddChild(std::move(theMedium));
@@ -259,6 +258,7 @@ int main() {
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
process::sibyll::Interaction sibyll(env);
process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
process::sibyll::Decay decay(trackedHadrons);
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b8851ffe622e1583a2ce55e60e7dd4a112065959
--- /dev/null
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -0,0 +1,325 @@
+
+/*
+ * (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/cascade/Cascade.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
+#include <corsika/process/stack_inspector/StackInspector.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+#include <corsika/geometry/Sphere.h>
+
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/pythia/Interaction.h>
+
+#include <corsika/process/track_writer/TrackWriter.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::setup;
+using namespace corsika::geometry;
+using namespace corsika::environment;
+
+using namespace std;
+using namespace corsika::units::si;
+
+class ProcessCut : public process::ContinuousProcess<ProcessCut> {
+
+ HEPEnergyType fECut;
+
+ HEPEnergyType fEnergy = 0_GeV;
+ HEPEnergyType fEmEnergy = 0_GeV;
+ int fEmCount = 0;
+ HEPEnergyType fInvEnergy = 0_GeV;
+ int fInvCount = 0;
+
+public:
+ ProcessCut(const HEPEnergyType cut)
+ : fECut(cut) {}
+
+ template <typename Particle>
+ bool isBelowEnergyCut(Particle& p) const {
+ // nuclei
+ if (p.GetPID() == particles::Code::Nucleus) {
+ auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
+ auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
+ if (ElabNuc < fECut || EcmNN < 10_GeV)
+ return true;
+ else
+ return false;
+ } else {
+ // TODO: center-of-mass energy hard coded
+ const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ if (p.GetEnergy() < fECut || Ecm < 10_GeV)
+ return true;
+ else
+ return false;
+ }
+ }
+
+ bool isEmParticle(Code pCode) const {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::Electron:
+ is_em = true;
+ break;
+ case Code::Positron:
+ is_em = true;
+ break;
+ case Code::Gamma:
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible(Code pCode) const {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::NuE:
+ is_inv = true;
+ break;
+ case Code::NuEBar:
+ is_inv = true;
+ break;
+ case Code::NuMu:
+ is_inv = true;
+ break;
+ case Code::NuMuBar:
+ is_inv = true;
+ break;
+ case Code::MuPlus:
+ is_inv = true;
+ break;
+ case Code::MuMinus:
+ is_inv = true;
+ break;
+
+ case Code::Neutron:
+ is_inv = true;
+ break;
+
+ case Code::AntiNeutron:
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+ template <typename Particle>
+ LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
+ cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
+ cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
+ const Code pid = p.GetPID();
+ if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0_m;
+ } else {
+ LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
+ const Code pid = p.GetPID();
+ HEPEnergyType energy = p.GetEnergy();
+ cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
+ EProcessReturn ret = EProcessReturn::eOk;
+ if (isEmParticle(pid)) {
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += energy;
+ fEmCount += 1;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ } else if (isInvisible(pid)) {
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += energy;
+ fInvCount += 1;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ } else if (isBelowEnergyCut(p)) {
+ cout << "removing low en. particle..." << endl;
+ fEnergy += energy;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ }
+ return ret;
+ }
+
+ void Init() {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
+ }
+
+ void ShowResults() {
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
+ << " ******************************" << endl;
+ }
+
+ HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+ HEPEnergyType GetCutEnergy() const { return fEnergy; }
+ HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+};
+
+//
+// The example main program for a particle cascade
+//
+int main() {
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+
+ // setup environment, geometry
+ environment::Environment env;
+ auto& universe = *(env.GetUniverse());
+
+ auto theMedium = environment::Environment::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Hydrogen},
+ std::vector<float>{(float)1.}));
+
+ universe.AddChild(std::move(theMedium));
+
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ // setup processes, decays and interactions
+ tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
+ stack_inspector::StackInspector<setup::Stack> p0(true);
+
+ const std::vector<particles::Code> trackedHadrons = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
+
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+ // process::sibyll::Interaction sibyll(env);
+ process::pythia::Interaction pythia(env);
+ // process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+ // process::sibyll::Decay decay(trackedHadrons);
+ process::pythia::Decay decay(trackedHadrons);
+ ProcessCut cut(20_GeV);
+
+ // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
+ // process::HadronicElasticModel::HadronicElasticInteraction
+ // hadronicElastic(env);
+
+ process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+
+ // assemble all processes into an ordered process list
+ // auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
+ // auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
+
+ auto sequence = p0 << pythia << decay << cut << trackWriter;
+
+ // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
+ // << "\n";
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Proton;
+ const HEPMassType mass = particles::Proton::GetMass();
+ const HEPEnergyType E0 = 100_GeV;
+ double theta = 0.;
+ double phi = 0.;
+
+ {
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt(Elab * Elab - m * m);
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+ Point pos(rootCS, 0_m, 0_m, 0_m);
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{
+ beamCode, E0, plab, pos, 0_ns});
+ }
+
+ // define air shower object, run simulation
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+ EAS.Init();
+ EAS.Run();
+
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
+ cut.ShowResults();
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+}
diff --git a/Processes/Pythia/CMakeLists.txt b/Processes/Pythia/CMakeLists.txt
index b23c6d005f272854ac976104de24262996e98f47..adc953090c31505acf4fd4c6c85dd85435ca136c 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 9e2d5dd39cd7a02297c22f298802cc79d5f040a9..418ff6973b6c51cc0c8a9cc6fa96c06e9c8988a3 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -9,6 +9,8 @@
* the license.
*/
+#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/pythia/Interaction.h>
#include <Pythia8/Pythia.h>
#include <corsika/process/pythia/Decay.h>
@@ -94,7 +96,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;
@@ -108,7 +110,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));
@@ -148,4 +150,28 @@ TEST_CASE("pytia decay") {
stack);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
+
+ 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);
+ }
+
}