diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 03204dde72a73c0fed4733e968b89266b6f76a6f..c55ac0b43b86e1956750618ec65e868a0d1334e0 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -201,7 +201,7 @@ namespace corsika::cascade {
actual_inv_length);
const auto sample_process = uniDist(fRNG);
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
- fProcessSequence.SelectInteraction(particle, fStack, sample_process,
+ fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
inv_lambda_count);
} else {
std::cout << "decay" << std::endl;
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 4d4d3fbfdaff3d075f6cc613da78ece584ed7142..1056ccd9e4cdb4ec80fca29b7c08857de1e28fc9 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -132,24 +132,24 @@ namespace corsika::process {
return tot;
}
- template <typename Particle, typename Stack>
+ template <typename Particle, typename Track, typename Stack>
EProcessReturn SelectInteraction(
- Particle& p, Stack& s,
+ Particle& vP, Track& vT, Stack& vS,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (is_process_sequence<T1type>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- A.SelectInteraction(p, s, lambda_select, lambda_inv_count);
+ A.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_count += A.GetInverseInteractionLength(p, s);
+ lambda_inv_count += A.GetInverseInteractionLength(vP, vT);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
- A.DoInteraction(p, s);
+ A.DoInteraction(vP, vS);
return EProcessReturn::eInteracted;
}
} // end branch A
@@ -157,15 +157,15 @@ namespace corsika::process {
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- B.SelectInteraction(p, s, lambda_select, lambda_inv_count);
+ B.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_count += B.GetInverseInteractionLength(p, s);
+ lambda_inv_count += B.GetInverseInteractionLength(vP, vT);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
- B.DoInteraction(p, s);
+ B.DoInteraction(vP, vS);
return EProcessReturn::eInteracted;
}
} // end branch A
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index e2f277e379c7ac0b393fb898ec7e6fd896dfaaf0..ef85de7af7f61c829e66f04058b89d48706ddc03 100644
--- a/Framework/StackInterface/ParticleBase.h
+++ b/Framework/StackInterface/ParticleBase.h
@@ -48,6 +48,7 @@ namespace corsika::stack {
class ParticleBase {
public:
+ using StackIteratorType = StackIterator;
ParticleBase() = default;
private:
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index bfb5acfbd10f5910ede79649df123bdc8aa296c5..e917f0a721913666da4eea3df37e5c3162c9ec7d 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -19,22 +19,6 @@
#include <corsika/stack/SecondaryView.h>
-// SFINAE test
-template <typename T>
-class HasGetIndexFromIterator {
-private:
- typedef char YesType[1];
- typedef char NoType[2];
-
- template <typename C>
- static YesType& test(decltype(&C::GetIndexFromIterator));
- template <typename C>
- static NoType& test(...);
-
-public:
- enum { value = sizeof(test<T>(0)) == sizeof(YesType) };
-};
-
/**
All classes around management of particles on a stack.
*/
@@ -123,9 +107,14 @@ namespace corsika::stack {
ParticleInterface, StackType>;
/**
- * this is the full type of the declared ParticleInterface: typedef typename
+ * this is the full type of the user-declared ParticleInterface
+ */
+ using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
+ /**
+ * In all programming context, the object to access, copy, and
+ * transport particle data is via the StackIterator
*/
- typedef typename StackIterator::ParticleInterfaceType ParticleType;
+ using ParticleType = StackIterator;
// friends are needed since they need access to protected members
friend class StackIteratorInterface<
@@ -216,7 +205,7 @@ namespace corsika::stack {
/**
* delete this particle
*/
- void Delete(ParticleType p) { Delete(p.GetIterator()); }
+ void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
/**
* delete last particle on stack by decrementing stack size
diff --git a/Framework/StackInterface/testCombinedStack.cc b/Framework/StackInterface/testCombinedStack.cc
index 63e984b1b95a9583156df1ce0bad72e5f85573fc..11a3314a2e60a517b9c1b197f913ffec1e7cc421 100644
--- a/Framework/StackInterface/testCombinedStack.cc
+++ b/Framework/StackInterface/testCombinedStack.cc
@@ -89,7 +89,6 @@ using CombinedTestInterfaceType =
TestParticleInterface2, StackIter>;
using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInterfaceType>;
-typedef StackTest::ParticleType Particle;
TEST_CASE("Combined Stack", "[stack]") {
@@ -161,10 +160,9 @@ TEST_CASE("Combined Stack", "[stack]") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto iter = s.AddParticle(std::tuple{9.9});
- Particle& p = *iter; // also this is valid to access particle data
- p.SetData2(2);
+ iter.SetData2(2);
REQUIRE(s.GetSize() == 1);
- p.AddSecondary(std::tuple{4.4});
+ iter.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 2);
// p.AddSecondary(3.3, 2.2, 1.);
// REQUIRE(s.GetSize() == 3);
@@ -259,7 +257,6 @@ using CombinedTestInterfaceType2 =
using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
CombinedTestInterfaceType2>;
-typedef StackTest2::ParticleType Particle2;
TEST_CASE("Combined Stack - multi", "[stack]") {
diff --git a/Framework/StackInterface/testSecondaryView.cc b/Framework/StackInterface/testSecondaryView.cc
index df8521386b3b4f49214ad92f2a69906b1c804b55..08d27c169bfa6144abd203603eae74dc1e134c98 100644
--- a/Framework/StackInterface/testSecondaryView.cc
+++ b/Framework/StackInterface/testSecondaryView.cc
@@ -32,7 +32,6 @@ using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
-typedef StackTest::ParticleType Particle;
TEST_CASE("SecondaryStack", "[stack]") {
diff --git a/Framework/StackInterface/testStackInterface.cc b/Framework/StackInterface/testStackInterface.cc
index 448aab92731f41c7722aaff369b1517ab76c7c25..5d4de9dae7b1cb49e6d4915eb565e2a9bf1c8416 100644
--- a/Framework/StackInterface/testStackInterface.cc
+++ b/Framework/StackInterface/testStackInterface.cc
@@ -33,7 +33,7 @@ using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
-typedef StackTest::ParticleType Particle;
+typedef StackTest::ParticleInterfaceType Particle;
TEST_CASE("Stack", "[Stack]") {
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index 880496a219b8674a2cca984c6ab3e8615ecd1470..357964b39c77d5d1ddcdaa9d3ae290debd61d911 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -18,6 +18,9 @@ add_custom_command (
set (
MODEL_SOURCES
ParticleConversion.cc
+ Interaction.cc
+ Decay.cc
+ NuclearInteraction.cc
sibyll2.3c.f
nuclib.f
signuc.f
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1e0d1ca2b64be8db9444c94894fc6983a01e62ba
--- /dev/null
+++ b/Processes/Sibyll/Decay.cc
@@ -0,0 +1,189 @@
+
+/*
+ * (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/sibyll/Decay.h>
+
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/SibStack.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using std::cout;
+using std::endl;
+using std::tuple;
+using std::vector;
+
+using namespace corsika;
+using namespace corsika::setup;
+using Particle = Stack::StackIterator; // ParticleType;
+using Track = Trajectory;
+
+namespace corsika::process::sibyll {
+
+ Decay::Decay() {}
+ Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
+ void Decay::Init() {
+ setHadronsUnstable();
+ setTrackedParticlesStable();
+ }
+
+ void Decay::setTrackedParticlesStable() {
+ /*
+ Sibyll is hadronic generator
+ only hadrons decay
+ */
+ // set particles unstable
+ setHadronsUnstable();
+ // make tracked particles stable
+ cout << "Interaction: setting tracked hadrons stable.." << endl;
+ const vector<particles::Code> particleList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
+ for (auto p : particleList) {
+ // set particle stable by setting table value negative
+ const int sibid = process::sibyll::ConvertToSibyllRaw(p);
+ s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
+ }
+ }
+
+ void Decay::setUnstable(const particles::Code pCode) {
+ int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
+ }
+
+ void Decay::setStable(const particles::Code pCode) {
+ int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ }
+
+ void Decay::setAllStable() {
+ // name? also makes EM particles stable
+
+ cout << "Decay: setting all particles stable.." << endl;
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. particles::ListOfParticles()
+ for (auto& p : corsika2sibyll) {
+ // cout << (int)p << endl;
+ const int sibCode = static_cast<int>(p);
+ // skip unknown and antiparticles
+ if (sibCode < 1) continue;
+ s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
+ }
+ }
+
+ void Decay::setHadronsUnstable() {
+
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. particles::ListOfParticles()
+ cout << "Sibyll: setting hadrons unstable.." << endl;
+ // make ALL particles unstable, then set EM stable
+ for (int sibCode : corsika2sibyll) {
+ if (sibCode < 1) continue;
+ s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
+ }
+ // set Leptons and Proton and Neutron stable
+ // use stack to loop over particles
+ constexpr particles::Code particleList[] = {
+ particles::Code::Proton, particles::Code::Neutron, particles::Code::Electron,
+ particles::Code::Positron, particles::Code::NuE, particles::Code::NuEBar,
+ particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::NuMu,
+ particles::Code::NuMuBar};
+
+ for (auto p : particleList) {
+ const int sibid = process::sibyll::ConvertToSibyllRaw(p);
+ s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
+ }
+ }
+
+ template <>
+ units::si::TimeType Decay::GetLifetime(Particle const& p) {
+ using namespace units::si;
+
+ HEPEnergyType E = p.GetEnergy();
+ HEPMassType m = particles::GetMass(p.GetPID());
+
+ const double gamma = E / m;
+
+ const TimeType t0 = particles::GetLifetime(p.GetPID());
+ auto const lifetime = gamma * t0;
+
+ const auto mkin =
+ (E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
+ cout << "Decay: code: " << p.GetPID() << endl;
+ cout << "Decay: MinStep: t0: " << t0 << endl;
+ cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
+ cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
+ << endl;
+ cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV << " "
+ << m / 1_GeV * m / 1_GeV << endl;
+ auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
+ cout << "Decay: MinStep: gamma: " << gamma << endl;
+ cout << "Decay: MinStep: tau: " << lifetime << endl;
+
+ return lifetime;
+ }
+
+ template <>
+ void Decay::DoDecay(Particle& p, Stack&) {
+ using geometry::Point;
+ using namespace units::si;
+
+ fCount++;
+ SibStack ss;
+ ss.Clear();
+ const particles::Code pCode = p.GetPID();
+ // copy particle to sibyll stack
+ ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
+ p.GetMomentum(),
+ // setting particle mass with Corsika values, may be inconsistent
+ // with sibyll internal values
+ particles::GetMass(pCode));
+ // remember position
+ Point const decayPoint = p.GetPosition();
+ TimeType const t0 = p.GetTime();
+ // remove original particle from corsika stack
+ p.Delete();
+ // set all particles/hadrons unstable
+ // setHadronsUnstable();
+ setUnstable(pCode);
+ // call sibyll decay
+ cout << "Decay: calling Sibyll decay routine.." << endl;
+ decsib_();
+ // reset to stable
+ setStable(pCode);
+ // print output
+ int print_unit = 6;
+ sib_list_(print_unit);
+
+ // copy particles from sibyll stack to corsika
+ for (auto& psib : ss) {
+ // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
+ if (psib.HasDecayed()) continue;
+ // add to corsika stack
+ p.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(),
+ psib.GetMomentum(), decayPoint, t0});
+ }
+ // empty sibyll stack
+ ss.Clear();
+ }
+
+} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 24ccffcb8766824fd7f5da6050096d6da897ead6..c078af7b8763c4d6df9015177ed764f6a663b460 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -12,14 +12,8 @@
#ifndef _include_corsika_process_sibyll_decay_h_
#define _include_corsika_process_sibyll_decay_h_
-#include <corsika/process/DecayProcess.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/process/sibyll/SibStack.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/DecayProcess.h>
namespace corsika::process {
@@ -29,191 +23,21 @@ namespace corsika::process {
int fCount = 0;
public:
- Decay() {}
- ~Decay() { std::cout << "Sibyll::Decay n=" << fCount << std::endl; }
- void Init() {
- setHadronsUnstable();
- setTrackedParticlesStable();
- }
-
- void setTrackedParticlesStable() {
- /*
- Sibyll is hadronic generator
- only hadrons decay
- */
- // set particles unstable
- setHadronsUnstable();
- // make tracked particles stable
- std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
- const std::vector<corsika::particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
- for (auto p : particleList) {
- // set particle stable by setting table value negative
- const int sibid = process::sibyll::ConvertToSibyllRaw(p);
- s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
- }
- }
-
- void setUnstable(const corsika::particles::Code pCode) {
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
- s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
- }
-
- void setStable(const corsika::particles::Code pCode) {
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- }
-
- void setAllStable() {
- // name? also makes EM particles stable
-
- using std::cout;
- using std::endl;
-
- std::cout << "Decay: setting all particles stable.." << std::endl;
-
- // loop over all particles in sibyll
- // should be changed to loop over human readable list
- // i.e. corsika::particles::ListOfParticles()
- for (auto& p : corsika2sibyll) {
- // std::cout << (int)p << std::endl;
- const int sibCode = static_cast<int>(p);
- // skip unknown and antiparticles
- if (sibCode < 1) continue;
- s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
- }
- }
+ Decay();
+ ~Decay();
+ void Init();
- void setHadronsUnstable() {
-
- using std::cout;
- using std::endl;
- using namespace corsika::units::si;
-
- // name? also makes EM particles stable
-
- // loop over all particles in sibyll
- // should be changed to loop over human readable list
- // i.e. corsika::particles::ListOfParticles()
- std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
- // make ALL particles unstable, then set EM stable
- for (int sibCode : corsika2sibyll) {
- // std::cout << (int)p << std::endl;
- // skip unknown and antiparticles
- if (sibCode < 1) continue;
- // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
- // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
- s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
- // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
- // std::endl;
- }
- // set Leptons and Proton and Neutron stable
- // use stack to loop over particles
- constexpr corsika::particles::Code particleList[] = {
- corsika::particles::Code::Proton, corsika::particles::Code::Neutron,
- corsika::particles::Code::Electron, corsika::particles::Code::Positron,
- corsika::particles::Code::NuE, corsika::particles::Code::NuEBar,
- corsika::particles::Code::MuMinus, corsika::particles::Code::MuPlus,
- corsika::particles::Code::NuMu, corsika::particles::Code::NuMuBar};
-
- for (auto p : particleList) {
- // set particle stable by setting table value negative
- // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
- // << " stable in Sibyll .." << endl;
- const int sibid = process::sibyll::ConvertToSibyllRaw(p);
- s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
- }
- }
+ void setTrackedParticlesStable();
+ void setUnstable(const corsika::particles::Code pCode);
+ void setStable(const corsika::particles::Code pCode);
+ void setAllStable();
+ void setHadronsUnstable();
template <typename Particle>
- corsika::units::si::TimeType GetLifetime(Particle const& p) {
- using std::cout;
- using std::endl;
- using namespace corsika::units::si;
-
- corsika::units::si::HEPEnergyType E = p.GetEnergy();
- corsika::units::si::HEPMassType m = corsika::particles::GetMass(p.GetPID());
-
- const double gamma = E / m;
-
- const corsika::units::si::TimeType t0 =
- corsika::particles::GetLifetime(p.GetPID());
- auto const lifetime = gamma * t0;
-
- const auto mkin =
- (E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
- cout << "Decay: code: " << p.GetPID() << endl;
- cout << "Decay: MinStep: t0: " << t0 << endl;
- cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
- cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
- << endl;
- cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV
- << " " << m / 1_GeV * m / 1_GeV << endl;
- // cout << "Decay: sib mass: " << s_mass1_.am2[
- // process::sibyll::ConvertToSibyllRaw(p.GetPID()) ] << endl;
- auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
- cout << "Decay: MinStep: gamma: " << gamma << endl;
- cout << "Decay: MinStep: tau: " << lifetime << endl;
-
- return lifetime;
- }
+ corsika::units::si::TimeType GetLifetime(Particle const& p);
template <typename Particle, typename Stack>
- void DoDecay(Particle& p, Stack&) {
- using corsika::geometry::Point;
- using namespace corsika::units::si;
-
- fCount++;
- SibStack ss;
- ss.Clear();
- const corsika::particles::Code pCode = p.GetPID();
- // copy particle to sibyll stack
- ss.AddParticle(
- // std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- // corsika::stack::MomentumVector, corsika::geometry::Point,
- // corsika::units::si::TimeType>{
- corsika::process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
- p.GetMomentum(),
- // setting particle mass with Corsika values, may be inconsistent
- // with sibyll internal values
- // TODO: #warning setting particle mass with Corsika values, may be
- // inconsistent with sibyll internal values
- corsika::particles::GetMass(pCode));
- // remember position
- Point const decayPoint = p.GetPosition();
- TimeType const t0 = p.GetTime();
- // remove original particle from corsika stack
- p.Delete();
- // set all particles/hadrons unstable
- // setHadronsUnstable();
- setUnstable(pCode);
- // call sibyll decay
- std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
- decsib_();
- // reset to stable
- setStable(pCode);
- // print output
- int print_unit = 6;
- sib_list_(print_unit);
-
- // copy particles from sibyll stack to corsika
- for (auto& psib : ss) {
- // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
- if (psib.HasDecayed()) continue;
- // add to corsika stack
- p.AddSecondary(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType>{
- corsika::process::sibyll::ConvertFromSibyll(psib.GetPID()),
- psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0});
- }
- // empty sibyll stack
- ss.Clear();
- }
+ void DoDecay(Particle& p, Stack&);
};
} // namespace sibyll
} // namespace corsika::process
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..893de4f51263b21b374d4292c9c6c18f58cbb879
--- /dev/null
+++ b/Processes/Sibyll/Interaction.cc
@@ -0,0 +1,362 @@
+
+/*
+ * (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/sibyll/Interaction.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/SibStack.h>
+#include <corsika/process/sibyll/sibyll2.3c.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::sibyll {
+
+ Interaction::Interaction(environment::Environment const& env)
+ : fEnvironment(env) {}
+
+ Interaction::~Interaction() {
+ cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
+ }
+
+ void Interaction::Init() {
+
+ using random::RNGManager;
+
+ // initialize Sibyll
+ if (!fInitialized) {
+ sibyll_ini_();
+
+ fInitialized = true;
+ }
+ }
+
+ tuple<units::si::CrossSectionType, units::si::CrossSectionType, int>
+ Interaction::GetCrossSection(const particles::Code BeamId,
+ const particles::Code TargetId,
+ const units::si::HEPEnergyType CoMenergy) {
+ using namespace units::si;
+ double sigProd, sigEla, dummy, dum1, dum3, dum4;
+ double dumdif[3];
+ const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
+ const double dEcm = CoMenergy / 1_GeV;
+ int iTarget = -1;
+ if (particles::IsNucleus(TargetId)) {
+ iTarget = particles::GetNucleusA(TargetId);
+ if (iTarget > 18 || iTarget == 0)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 are allowed.");
+ sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+ } else if (TargetId == particles::Proton::GetCode()) {
+ sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+ iTarget = 1;
+ } else {
+ // no interaction in sibyll possible, return infinite cross section? or throw?
+ sigProd = std::numeric_limits<double>::infinity();
+ sigEla = std::numeric_limits<double>::infinity();
+ }
+ return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
+ }
+
+ 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 sibyll : 1, 2, 3 for p, pi, k
+ // read from cross section code table
+ const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
+
+ const HEPMassType nucleon_mass =
+ 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
+
+ // FOR NOW: assume target is at rest
+ MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
+
+ // total momentum and energy
+ HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ 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 && ECoM >= 10_GeV) {
+
+ // 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;
+ double avgTargetMassNumber = 0.;
+ 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, numberOfNucleons] =
+ GetCrossSection(corsikaBeamId, targetId, ECoM);
+ [[maybe_unused]] auto elaCrossSectionCopy =
+ elaCrossSection; // ONLY TO AVOID COMPILER WARNING
+
+ cout << "Interaction: "
+ << " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn
+ << endl;
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ avgTargetMassNumber += w[i] * numberOfNucleons;
+ }
+ cout << "Interaction: "
+ << "IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mbarn << endl;
+
+ // calculate interaction length in medium
+ //#warning check interaction length units
+ GrammageType const int_length =
+ avgTargetMassNumber * 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 SIBYLL 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 << "ProcessSibyll: "
+ << "DoInteraction: " << corsikaBeamId << " interaction? "
+ << process::sibyll::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 SIBYLL!");
+ }
+
+ if (process::sibyll::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 Sibyll is always a single nucleon
+ auto constexpr nucleon_mass =
+ 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
+ // FOR NOW: target is always at rest
+ const auto eTargetLab = 0_GeV + nucleon_mass;
+ 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 Sibyll projectile: +z
+ COMBoost const boost(PprojLab, nucleon_mass);
+
+ // 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, nNuc] =
+ GetCrossSection(corsikaBeamId, targetId, Ecm);
+ cross_section_of_components[i] = sigProd;
+ [[maybe_unused]] int ideleteme =
+ nNuc; // to avoid not used warning in array binding
+ [[maybe_unused]] auto sigElaCopy =
+ sigEla; // to avoid not used warning in array binding
+ }
+
+ const auto targetCode = currentNode->GetModelProperties().SampleTarget(
+ cross_section_of_components, fRNG);
+ cout << "Interaction: target selected: " << targetCode << endl;
+ /*
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be
+ allowed air in atmosphere also contains some Argon.
+ */
+ int targetSibCode = -1;
+ if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
+ if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
+ cout << "Interaction: sibyll code: " << targetSibCode << endl;
+ if (targetSibCode > 18 || targetSibCode < 1)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 or protons are "
+ "allowed.");
+
+ // beam id for sibyll
+ const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
+
+ cout << "Interaction: "
+ << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
+ << " Ecm(GeV): " << Ecm / 1_GeV << endl;
+ if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
+ cout << "Interaction: "
+ << " DoInteraction: should have dropped particle.. "
+ << "THIS IS AN ERROR" << endl;
+ throw std::runtime_error("energy too low for SIBYLL");
+ // p.Delete(); delete later... different process
+ } else {
+ fCount++;
+ // Sibyll does not know about units..
+ const double sqs = Ecm / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_(kBeam, targetSibCode, sqs);
+ // running decays
+ // setTrackedParticlesStable();
+ decsib_();
+ // print final state
+ int print_unit = 6;
+ sib_list_(print_unit);
+ fNucCount += get_nwounded() - 1;
+
+ // delete current particle
+ p.Delete();
+
+ // add particles from sibyll to stack
+ // link to sibyll stack
+ SibStack ss;
+
+ MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
+ for (auto& psib : ss) {
+
+ // skip particles that have decayed in Sibyll
+ if (psib.HasDecayed()) continue;
+
+ // transform energy to lab. frame
+ auto const pCoM = psib.GetMomentum();
+ HEPEnergyType const eCoM = psib.GetEnergy();
+ auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
+
+ // add to corsika stack
+ auto pnew = p.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ process::sibyll::ConvertFromSibyll(psib.GetPID()),
+ Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+ tOrig});
+
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
+ Ecm_final += psib.GetEnergy();
+ }
+ cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << endl
+ << "Elab_final=" << Elab_final / 1_GeV
+ << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
+ }
+ }
+ return process::EProcessReturn::eOk;
+ }
+
+} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 6a85592fe5af8b62d7e75f34db80abe60e247c16..40463f6657c2b6e877be5a018b7cfcf67ae6f884 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -12,18 +12,15 @@
#ifndef _corsika_process_sibyll_interaction_h_
#define _corsika_process_sibyll_interaction_h_
-#include <corsika/process/InteractionProcess.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/FourVector.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/process/sibyll/SibStack.h>
-#include <corsika/process/sibyll/sibyll2.3c.h>
+#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/COMBoost.h>
+#include <tuple>
+
+namespace corsika::environment {
+ class Environment;
+}
namespace corsika::process::sibyll {
@@ -34,28 +31,12 @@ namespace corsika::process::sibyll {
bool fInitialized = false;
public:
- Interaction(corsika::environment::Environment const& env)
- : fEnvironment(env) {}
- ~Interaction() {
- std::cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount
- << std::endl;
- }
-
- void Init() {
+ Interaction(corsika::environment::Environment const& env);
+ ~Interaction();
- using corsika::random::RNGManager;
- using std::cout;
- using std::endl;
-
- // initialize Sibyll
- if (!fInitialized) {
- sibyll_ini_();
- fInitialized = true;
- }
- }
+ void Init();
bool WasInitialized() { return fInitialized; }
-
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
using namespace corsika::units::si;
return (10_GeV < ecm) && (ecm < 1_PeV);
@@ -65,123 +46,10 @@ namespace corsika::process::sibyll {
int>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
- const corsika::units::si::HEPEnergyType CoMenergy) {
- using namespace corsika::units::si;
- double sigProd, sigEla, dummy, dum1, dum3, dum4;
- double dumdif[3];
- const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
- const double dEcm = CoMenergy / 1_GeV;
- int iTarget = -1;
- if (corsika::particles::IsNucleus(TargetId)) {
- iTarget = corsika::particles::GetNucleusA(TargetId);
- if (iTarget > 18 || iTarget == 0)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 are allowed.");
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
- } else if (TargetId == corsika::particles::Proton::GetCode()) {
- sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
- iTarget = 1;
- } else {
- // no interaction in sibyll possible, return infinite cross section? or throw?
- sigProd = std::numeric_limits<double>::infinity();
- sigEla = std::numeric_limits<double>::infinity();
- }
- return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
- }
+ const corsika::units::si::HEPEnergyType CoMenergy);
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
-
- using namespace corsika::units;
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- using std::cout;
- using std::endl;
-
- // coordinate system, get global frame of reference
- CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
- const particles::Code corsikaBeamId = p.GetPID();
-
- // beam particles for sibyll : 1, 2, 3 for p, pi, k
- // read from cross section code table
- const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
- const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
-
- // FOR NOW: assume target is at rest
- MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
-
- // total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
- 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
-
- std::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 && ECoM >= 10_GeV) {
-
- // 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;
- double avgTargetMassNumber = 0.;
- 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, numberOfNucleons] =
- GetCrossSection(corsikaBeamId, targetId, ECoM);
- [[maybe_unused]] auto elaCrossSectionCopy =
- elaCrossSection; // ONLY TO AVOID COMPILER WARNING
-
- std::cout << "Interaction: "
- << " IntLength: sibyll return (mb): "
- << productionCrossSection / 1_mbarn << std::endl;
- weightedProdCrossSection += w[i] * productionCrossSection;
- avgTargetMassNumber += w[i] * numberOfNucleons;
- }
- cout << "Interaction: "
- << "IntLength: weighted CrossSection (mb): "
- << weightedProdCrossSection / 1_mbarn << endl;
-
- // calculate interaction length in medium
- //#warning check interaction length units
- GrammageType const int_length =
- avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
- std::cout << "Interaction: "
- << "interaction length (g/cm2): "
- << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
-
- return int_length;
- }
-
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
/**
In this function SIBYLL is called to produce one event. The
@@ -189,193 +57,7 @@ namespace corsika::process::sibyll {
*/
template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
-
- using namespace corsika::units;
- using namespace corsika::utl;
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- using std::cout;
- using std::endl;
-
- const auto corsikaBeamId = p.GetPID();
- cout << "ProcessSibyll: "
- << "DoInteraction: " << corsikaBeamId << " interaction? "
- << process::sibyll::CanInteract(corsikaBeamId) << endl;
-
- if (corsika::particles::IsNucleus(corsikaBeamId)) {
- // nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
- }
-
- if (process::sibyll::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 Sibyll is always a single nucleon
- auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
- // FOR NOW: target is always at rest
- const auto eTargetLab = 0_GeV + nucleon_mass;
- 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 Sibyll projectile: +z
- COMBoost const boost(PprojLab, nucleon_mass);
-
- // 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, nNuc] =
- GetCrossSection(corsikaBeamId, targetId, Ecm);
- cross_section_of_components[i] = sigProd;
- [[maybe_unused]] int ideleteme =
- nNuc; // to avoid not used warning in array binding
- [[maybe_unused]] auto sigElaCopy =
- sigEla; // to avoid not used warning in array binding
- }
-
- const auto targetCode = currentNode->GetModelProperties().SampleTarget(
- cross_section_of_components, fRNG);
- cout << "Interaction: target selected: " << targetCode << endl;
- /*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be
- allowed air in atmosphere also contains some Argon.
- */
- int targetSibCode = -1;
- if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
- if (targetCode == corsika::particles::Proton::GetCode()) targetSibCode = 1;
- cout << "Interaction: sibyll code: " << targetSibCode << endl;
- if (targetSibCode > 18 || targetSibCode < 1)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 or protons are "
- "allowed.");
-
- // beam id for sibyll
- const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
-
- std::cout << "Interaction: "
- << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
- << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
- if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
- std::cout << "Interaction: "
- << " DoInteraction: should have dropped particle.. "
- << "THIS IS AN ERROR" << std::endl;
- throw std::runtime_error("energy too low for SIBYLL");
- // p.Delete(); delete later... different process
- } else {
- fCount++;
- // Sibyll does not know about units..
- const double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, targetSibCode, sqs);
- // running decays
- // setTrackedParticlesStable();
- decsib_();
- // print final state
- int print_unit = 6;
- sib_list_(print_unit);
- fNucCount += get_nwounded() - 1;
-
- // delete current particle
- p.Delete();
-
- // add particles from sibyll to stack
- // link to sibyll stack
- SibStack ss;
-
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
- for (auto& psib : ss) {
-
- // skip particles that have decayed in Sibyll
- if (psib.HasDecayed()) continue;
-
- // transform energy to lab. frame
- auto const pCoM = psib.GetMomentum();
- HEPEnergyType const eCoM = psib.GetEnergy();
- auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
-
- // add to corsika stack
- auto pnew =
- p.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
- process::sibyll::ConvertFromSibyll(psib.GetPID()),
- Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
- tOrig});
-
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
- Ecm_final += psib.GetEnergy();
- }
- std::cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
- << std::endl
- << "Elab_final=" << Elab_final / 1_GeV
- << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
- << std::endl;
- }
- }
- return process::EProcessReturn::eOk;
- }
+ corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
private:
corsika::environment::Environment const& fEnvironment;
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b3602cb57198aeb57619c37476ef3050594431da
--- /dev/null
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -0,0 +1,523 @@
+
+/*
+ * (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/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/process/sibyll/nuclib.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using std::cout;
+using std::endl;
+using std::tuple;
+using std::vector;
+
+using namespace corsika;
+using namespace corsika::setup;
+using Particle = Stack::StackIterator; // ParticleType;
+using Track = Trajectory;
+
+namespace corsika::process::sibyll {
+
+ NuclearInteraction::NuclearInteraction(environment::Environment const& env,
+ process::sibyll::Interaction& hadint)
+ : fEnvironment(env)
+ , fHadronicInteraction(hadint) {}
+
+ NuclearInteraction::~NuclearInteraction() {
+ cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
+ }
+
+ void NuclearInteraction::Init() {
+
+ using random::RNGManager;
+
+ // initialize hadronic interaction module
+ // TODO: safe to run multiple initializations?
+ if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
+
+ // initialize nuclib
+ // TODO: make sure this does not overlap with sibyll
+ nuc_nuc_ini_();
+ }
+
+ // TODO: remove number of nucleons, avg target mass is available in environment
+ template <>
+ tuple<units::si::CrossSectionType, units::si::CrossSectionType>
+ NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) {
+ using namespace units::si;
+ double sigProd;
+ auto const pCode = p.GetPID();
+ if (pCode != particles::Code::Nucleus)
+ throw std::runtime_error(
+ "NuclearInteraction: GetCrossSection: particle not a nucleus!");
+
+ auto const iBeam = p.GetNuclearA();
+ HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam;
+ cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeam
+ << " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV
+ << endl;
+
+ // use nuclib to calc. nuclear cross sections
+ // TODO: for now assumes air with hard coded composition
+ // extend to arbitrary mixtures, requires smarter initialization
+ // get nuclib projectile code: nucleon number
+ if (iBeam > 56 || iBeam < 2) {
+ cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" << endl
+ << "A=" << iBeam << endl;
+ throw std::runtime_error(
+ "NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
+ "NUCLIB!");
+ }
+
+ const double dElabNuc = LabEnergyPerNuc / 1_GeV;
+ // TODO: these limitations are still sibyll specific.
+ // available target nuclei depends on the hadronic interaction model and the
+ // initialization
+ if (dElabNuc < 10.)
+ throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
+
+ // TODO: these limitations are still sibyll specific.
+ // available target nuclei depends on the hadronic interaction model and the
+ // initialization
+ if (particles::IsNucleus(TargetId)) {
+ const int iTarget = particles::GetNucleusA(TargetId);
+ if (iTarget > 18 || iTarget == 0)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 are allowed.");
+ cout << "NuclearInteraction: calling signuc.." << endl;
+ cout << "WARNING: using hard coded cross section for Nucleus-Air with "
+ "SIBYLL! (fix me!)"
+ << endl;
+ // TODO: target id is not used because cross section is still hard coded and fixed
+ // to air.
+ signuc_(iBeam, dElabNuc, sigProd);
+ cout << "cross section: " << sigProd << endl;
+ return std::make_tuple(sigProd * 1_mbarn, 0_mbarn);
+ }
+ return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
+ std::numeric_limits<double>::infinity() * 1_mbarn);
+ }
+
+ template <>
+ units::si::GrammageType NuclearInteraction::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();
+ if (!particles::IsNucleus(corsikaBeamId)) {
+ // no nuclear interaction
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+ // check if target-style nucleus (enum)
+ if (corsikaBeamId != particles::Code::Nucleus)
+ throw std::runtime_error(
+ "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
+ "projectiles should use NuclearStackExtension!");
+
+ // read from cross section code table
+ const HEPMassType nucleon_mass =
+ 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
+
+ // FOR NOW: assume target is at rest
+ corsika::stack::MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+
+ // total momentum and energy
+ HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ int const nuclA = p.GetNuclearA();
+ auto const ElabNuc = p.GetEnergy() / nuclA;
+
+ corsika::stack::MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.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
+ auto const ECoMNN = sqrt(2. * ElabNuc * nucleon_mass);
+ cout << "NuclearInteraction: LambdaInt: \n"
+ << " input energy: " << Elab / 1_GeV << endl
+ << " input energy CoM: " << ECoM / 1_GeV << endl
+ << " beam pid:" << corsikaBeamId << endl
+ << " beam A: " << nuclA << endl
+ << " input energy per nucleon: " << ElabNuc / 1_GeV << endl
+ << " input energy CoM per nucleon: " << ECoMNN / 1_GeV << endl;
+ // throw std::runtime_error("stop here");
+
+ // energy limits
+ // TODO: values depend on hadronic interaction model !! this is sibyll specific
+ if (ElabNuc >= 8.5_GeV && ECoMNN >= 10_GeV) {
+
+ // 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 << "NuclearInteraction: get interaction length for target: " << targetId
+ << endl;
+ auto const [productionCrossSection, elaCrossSection] =
+ GetCrossSection(p, targetId);
+ [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
+
+ cout << "NuclearInteraction: "
+ << "IntLength: nuclib return (mb): " << productionCrossSection / 1_mbarn
+ << endl;
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ }
+ cout << "NuclearInteraction: "
+ << "IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mbarn << endl;
+
+ // calculate interaction length in medium
+ GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+ units::constants::u / weightedProdCrossSection;
+ cout << "NuclearInteraction: "
+ << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
+ << endl;
+
+ return int_length;
+ } else {
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+ }
+
+ template <>
+ process::EProcessReturn NuclearInteraction::DoInteraction(Particle& p, Stack& s) {
+
+ // this routine superimposes different nucleon-nucleon interactions
+ // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
+
+ using namespace units;
+ using namespace utl;
+ using namespace units::si;
+ using namespace geometry;
+
+ const auto ProjId = p.GetPID();
+ // TODO: calculate projectile mass in nuclearStackExtension
+ // const auto ProjMass = p.GetMass();
+ cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
+
+ if (!IsNucleus(ProjId)) {
+ cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
+ // this should not happen
+ // throw std::runtime_error("Non nuclear projectile in NUCLIB!");
+ return process::EProcessReturn::eOk;
+ }
+
+ // check if target-style nucleus (enum)
+ if (ProjId != particles::Code::Nucleus)
+ throw std::runtime_error(
+ "NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
+ "should use NuclearStackExtension!");
+
+ auto const ProjMass =
+ p.GetNuclearZ() * particles::Proton::GetMass() +
+ (p.GetNuclearA() - p.GetNuclearZ()) * particles::Neutron::GetMass();
+ cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
+
+ fCount++;
+
+ const CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ // position and time of interaction, not used in NUCLIB
+ Point pOrig = p.GetPosition();
+ TimeType tOrig = p.GetTime();
+
+ cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
+ cout << "Interaction: time: " << tOrig << endl;
+
+ // projectile nucleon number
+ const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
+ if (kAProj > 56) throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
+
+ // kinematics
+ // define projectile nucleus
+ HEPEnergyType const eProjectileLab = p.GetEnergy();
+ auto const pProjectileLab = p.GetMomentum();
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
+ << "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV
+ << endl;
+
+ // define projectile nucleon
+ HEPEnergyType const eProjectileNucLab = p.GetEnergy() / kAProj;
+ auto const pProjectileNucLab = p.GetMomentum() / kAProj;
+ const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
+
+ cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV << endl
+ << "NuclearInteraction: pProjNucleon lab: "
+ << pProjectileNucLab.GetComponents() / 1_GeV << endl;
+
+ // define target
+ // always a nucleon
+ auto constexpr nucleon_mass =
+ 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
+ // target is always at rest
+ const auto eTargetNucLab = 0_GeV + nucleon_mass;
+ const auto pTargetNucLab =
+ corsika::stack::MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+ const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
+
+ cout << "NuclearInteraction: etarget lab: " << eTargetNucLab / 1_GeV << endl
+ << "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV
+ << endl;
+
+ // center-of-mass energy in nucleon-nucleon frame
+ auto const PtotNN4 = PtargNucLab + PprojNucLab;
+ HEPEnergyType EcmNN = PtotNN4.GetNorm();
+ cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
+
+ if (!fHadronicInteraction.ValidCoMEnergy(EcmNN)) {
+ cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
+ "interaction model!"
+ << endl;
+ throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
+ }
+
+ // define boost to NUCLEON-NUCLEON frame
+ COMBoost const boost(PprojNucLab, nucleon_mass);
+ // boost projecticle
+ auto const PprojNucCoM = boost.toCoM(PprojNucLab);
+
+ // boost target
+ auto const PtargNucCoM = boost.toCoM(PtargNucLab);
+
+ cout << "Interaction: ebeam CoM: " << PprojNucCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: pbeam CoM: "
+ << PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ cout << "Interaction: etarget CoM: " << PtargNucCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: ptarget CoM: "
+ << PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+
+ // sample target nucleon number
+ //
+ // proton stand-in for nucleon
+ const auto beamId = particles::Proton::GetCode();
+ const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ const auto& mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ cout << "get nucleon-nucleus cross sections for target materials.." << endl;
+ // get cross sections for target materials
+ // using nucleon-target-nucleus cross section!!!
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from GetInteractionLength if possible
+ */
+ auto const& compVec = mediumComposition.GetComponents();
+ vector<si::CrossSectionType> cross_section_of_components(compVec.size());
+
+ for (size_t i = 0; i < compVec.size(); ++i) {
+ auto const targetId = compVec[i];
+ cout << "target component: " << targetId << endl;
+ cout << "beam id: " << beamId << endl;
+ const auto [sigProd, sigEla, nNuc] =
+ fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
+ cross_section_of_components[i] = sigProd;
+ [[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
+ [[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
+ }
+
+ const auto targetCode =
+ currentNode->GetModelProperties().SampleTarget(cross_section_of_components, fRNG);
+ cout << "Interaction: target selected: " << targetCode << endl;
+ /*
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be
+ allowed air in atmosphere also contains some Argon.
+ */
+ int kATarget = -1;
+ if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
+ if (targetCode == particles::Proton::GetCode()) kATarget = 1;
+ cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
+ if (kATarget > 18 || kATarget < 1)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 or protons are "
+ "allowed.");
+ // end of target sampling
+
+ // superposition
+ cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. " << endl;
+ // get nucleon-nucleon cross section
+ // (needed to determine number of nucleon-nucleon scatterings)
+ const auto protonId = particles::Proton::GetCode();
+ const auto [prodCrossSection, elaCrossSection, dum] =
+ fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
+ [[maybe_unused]] auto dumCopy = dum; // ONLY TO AVOID COMPILER WARNING
+ const double sigProd = prodCrossSection / 1_mbarn;
+ const double sigEla = elaCrossSection / 1_mbarn;
+ // sample number of interactions (only input variables, output in common cnucms)
+ // nuclear multiple scattering according to glauber (r.i.p.)
+ int_nuc_(kATarget, kAProj, sigProd, sigEla);
+
+ cout << "number of nucleons in target : " << kATarget << endl
+ << "number of wounded nucleons in target : " << cnucms_.na << endl
+ << "number of nucleons in projectile : " << kAProj << endl
+ << "number of wounded nucleons in project. : " << cnucms_.nb << endl
+ << "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl
+ << "number of elastic nucleons in target : " << cnucms_.nael << endl
+ << "number of elastic nucleons in project. : " << cnucms_.nbel << endl
+ << "impact parameter: " << cnucms_.b << endl;
+
+ // calculate fragmentation
+ cout << "calculating nuclear fragments.." << endl;
+ // number of interactions
+ // include elastic
+ const int nElasticNucleons = cnucms_.nbel;
+ const int nInelNucleons = cnucms_.nb;
+ const int nIntProj = nInelNucleons + nElasticNucleons;
+ const double impactPar = cnucms_.b; // only needed to avoid passing common var.
+ int nFragments;
+ // number of fragments is limited to 60
+ int AFragments[60];
+ // call fragmentation routine
+ // input: target A, projectile A, number of int. nucleons in projectile, impact
+ // parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
+ // in pf in common fragments, neglected
+ fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
+
+ // this should not occur but well :)
+ if (nFragments > 60)
+ throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
+
+ cout << "number of fragments: " << nFragments << endl;
+ for (int j = 0; j < nFragments; ++j)
+ cout << "fragment: " << j << " A=" << AFragments[j]
+ << " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1]
+ << " pz=" << fragments_.ppp[j][2] << endl;
+
+ cout << "adding nuclear fragments to particle stack.." << endl;
+ // put nuclear fragments on corsika stack
+ for (int j = 0; j < nFragments; ++j) {
+ particles::Code specCode;
+ const auto nuclA = AFragments[j];
+ // get Z from stability line
+ const auto nuclZ = int(nuclA / 2.15 + 0.7);
+
+ // TODO: do we need to catch single nucleons??
+ if (nuclA == 1)
+ // TODO: sample neutron or proton
+ specCode = particles::Code::Proton;
+ else
+ specCode = particles::Code::Nucleus;
+
+ // TODO: mass of nuclei?
+ const HEPMassType mass =
+ particles::Proton::GetMass() * nuclZ +
+ (nuclA - nuclZ) * particles::Neutron::GetMass(); // this neglects binding energy
+
+ cout << "NuclearInteraction: adding fragment: " << specCode << endl;
+ cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << endl;
+ cout << "NuclearInteraction: mass: " << mass / 1_GeV << endl;
+
+ // CORSIKA 7 way
+ // spectators inherit momentum from original projectile
+ const double mass_ratio = mass / ProjMass;
+
+ cout << "NuclearInteraction: mass ratio " << mass_ratio << endl;
+
+ auto const Plab = PprojLab * mass_ratio;
+
+ cout << "NuclearInteraction: fragment momentum: "
+ << Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+
+ if (nuclA == 1)
+ // add nucleon
+ p.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
+ stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+ tOrig});
+ else
+ // add nucleus
+ p.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+ tOrig, nuclA, nuclZ});
+ }
+
+ // add elastic nucleons to corsika stack
+ // TODO: the elastic interaction could be external like the inelastic interaction,
+ // e.g. use existing ElasticModel
+ cout << "adding elastically scattered nucleons to particle stack.." << endl;
+ for (int j = 0; j < nElasticNucleons; ++j) {
+ // TODO: sample proton or neutron
+ auto const elaNucCode = particles::Code::Proton;
+
+ // CORSIKA 7 way
+ // elastic nucleons inherit momentum from original projectile
+ // neglecting momentum transfer in interaction
+ const double mass_ratio = particles::GetMass(elaNucCode) / ProjMass;
+ auto const Plab = PprojLab * mass_ratio;
+
+ p.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ elaNucCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
+ pOrig, tOrig});
+ }
+
+ // add inelastic interactions
+ cout << "calculate inelastic nucleon-nucleon interactions.." << endl;
+ for (int j = 0; j < nInelNucleons; ++j) {
+ // TODO: sample neutron or proton
+ auto pCode = particles::Proton::GetCode();
+ // temporarily add to stack, will be removed after interaction in DoInteraction
+ cout << "inelastic interaction no. " << j << endl;
+ auto inelasticNucleon = p.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
+ geometry::Point, units::si::TimeType>{
+ pCode, PprojNucLab.GetTimeLikeComponent(),
+ PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
+ // create inelastic interaction
+ cout << "calling HadronicInteraction..." << endl;
+ fHadronicInteraction.DoInteraction(inelasticNucleon, s);
+ }
+
+ // delete parent particle
+ p.Delete();
+
+ cout << "NuclearInteraction: DoInteraction: done" << endl;
+
+ return process::EProcessReturn::eOk;
+ }
+
+} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index cd2735029bd739405cce4b7ad52680e0f2016141..cf1f6593330fddc72e867dd7d40387747d6ee9be 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -12,18 +12,23 @@
#ifndef _corsika_process_sibyll_nuclearinteraction_h_
#define _corsika_process_sibyll_nuclearinteraction_h_
-#include <corsika/process/InteractionProcess.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/sibyll/nuclib.h>
+#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/COMBoost.h>
+
+namespace corsika::environment {
+ class Environment;
+}
namespace corsika::process::sibyll {
+ class Interaction; // fwd-decl
+
+ /**
+ *
+ *
+ **/
+
class NuclearInteraction
: public corsika::process::InteractionProcess<NuclearInteraction> {
@@ -32,506 +37,19 @@ namespace corsika::process::sibyll {
public:
NuclearInteraction(corsika::environment::Environment const& env,
- corsika::process::sibyll::Interaction& hadint)
- : fEnvironment(env)
- , fHadronicInteraction(hadint) {}
- ~NuclearInteraction() {
- std::cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
- << std::endl;
- }
-
- void Init() {
-
- using corsika::random::RNGManager;
- using std::cout;
- using std::endl;
-
- // initialize hadronic interaction module
- // TODO: safe to run multiple initializations?
- if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
+ corsika::process::sibyll::Interaction& hadint);
+ ~NuclearInteraction();
+ void Init();
- // initialize nuclib
- // TODO: make sure this does not overlap with sibyll
- nuc_nuc_ini_();
- }
-
- // TODO: remove number of nucleons, avg target mass is available in environment
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
- GetCrossSection(Particle const& p, const corsika::particles::Code TargetId) {
- using namespace corsika::units::si;
- double sigProd;
- auto const pCode = p.GetPID();
- if (pCode != corsika::particles::Code::Nucleus)
- throw std::runtime_error(
- "NuclearInteraction: GetCrossSection: particle not a nucleus!");
-
- auto const iBeam = p.GetNuclearA();
- HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam;
- std::cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= "
- << iBeam << " TargetId= " << TargetId
- << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV << std::endl;
-
- // use nuclib to calc. nuclear cross sections
- // TODO: for now assumes air with hard coded composition
- // extend to arbitrary mixtures, requires smarter initialization
- // get nuclib projectile code: nucleon number
- if (iBeam > 56 || iBeam < 2) {
- std::cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
- << std::endl
- << "A=" << iBeam << std::endl;
- throw std::runtime_error(
- "NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
- "NUCLIB!");
- }
-
- const double dElabNuc = LabEnergyPerNuc / 1_GeV;
- // TODO: these limitations are still sibyll specific.
- // available target nuclei depends on the hadronic interaction model and the
- // initialization
- if (dElabNuc < 10.)
- throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
-
- // TODO: these limitations are still sibyll specific.
- // available target nuclei depends on the hadronic interaction model and the
- // initialization
- if (corsika::particles::IsNucleus(TargetId)) {
- const int iTarget = corsika::particles::GetNucleusA(TargetId);
- if (iTarget > 18 || iTarget == 0)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 are allowed.");
- std::cout << "NuclearInteraction: calling signuc.." << std::endl;
- std::cout << "WARNING: using hard coded cross section for Nucleus-Air with "
- "SIBYLL! (fix me!)"
- << std::endl;
- // TODO: target id is not used because cross section is still hard coded and fixed
- // to air.
- signuc_(iBeam, dElabNuc, sigProd);
- std::cout << "cross section: " << sigProd << std::endl;
- return std::make_tuple(sigProd * 1_mbarn, 0_mbarn);
- }
- return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
- std::numeric_limits<double>::infinity() * 1_mbarn);
- }
+ GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
-
- using namespace corsika::units;
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- using std::cout;
- using std::endl;
-
- // coordinate system, get global frame of reference
- CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
- const particles::Code corsikaBeamId = p.GetPID();
- if (!corsika::particles::IsNucleus(corsikaBeamId)) {
- // no nuclear interaction
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
- // check if target-style nucleus (enum)
- if (corsikaBeamId != corsika::particles::Code::Nucleus)
- throw std::runtime_error(
- "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
- "projectiles should use NuclearStackExtension!");
-
- // read from cross section code table
- const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
-
- // FOR NOW: assume target is at rest
- MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-
- // total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
- int const nuclA = p.GetNuclearA();
- auto const ElabNuc = p.GetEnergy() / nuclA;
-
- MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.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
- auto const ECoMNN = sqrt(2. * ElabNuc * nucleon_mass);
- std::cout << "NuclearInteraction: LambdaInt: \n"
- << " input energy: " << Elab / 1_GeV << endl
- << " input energy CoM: " << ECoM / 1_GeV << endl
- << " beam pid:" << corsikaBeamId << endl
- << " beam A: " << nuclA << endl
- << " input energy per nucleon: " << ElabNuc / 1_GeV << endl
- << " input energy CoM per nucleon: " << ECoMNN / 1_GeV << endl;
- // throw std::runtime_error("stop here");
-
- // energy limits
- // TODO: values depend on hadronic interaction model !! this is sibyll specific
- if (ElabNuc >= 8.5_GeV && ECoMNN >= 10_GeV) {
-
- // 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 << "NuclearInteraction: get interaction length for target: " << targetId
- << endl;
- auto const [productionCrossSection, elaCrossSection] =
- GetCrossSection(p, targetId);
- [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
-
- std::cout << "NuclearInteraction: "
- << "IntLength: nuclib return (mb): "
- << productionCrossSection / 1_mbarn << std::endl;
- weightedProdCrossSection += w[i] * productionCrossSection;
- }
- cout << "NuclearInteraction: "
- << "IntLength: weighted CrossSection (mb): "
- << weightedProdCrossSection / 1_mbarn << endl;
-
- // calculate interaction length in medium
- GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
- corsika::units::constants::u /
- weightedProdCrossSection;
- std::cout << "NuclearInteraction: "
- << "interaction length (g/cm2): "
- << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
-
- return int_length;
- } else {
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
- }
+ corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&);
template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
-
- // this routine superimposes different nucleon-nucleon interactions
- // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
-
- using namespace corsika::units;
- using namespace corsika::utl;
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- using std::cout;
- using std::endl;
-
- const auto ProjId = p.GetPID();
- // TODO: calculate projectile mass in nuclearStackExtension
- // const auto ProjMass = p.GetMass();
- cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
-
- if (!IsNucleus(ProjId)) {
- cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
- // this should not happen
- // throw std::runtime_error("Non nuclear projectile in NUCLIB!");
- return process::EProcessReturn::eOk;
- }
-
- // check if target-style nucleus (enum)
- if (ProjId != corsika::particles::Code::Nucleus)
- throw std::runtime_error(
- "NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
- "should use NuclearStackExtension!");
-
- auto const ProjMass =
- p.GetNuclearZ() * corsika::particles::Proton::GetMass() +
- (p.GetNuclearA() - p.GetNuclearZ()) * corsika::particles::Neutron::GetMass();
- cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
-
- fCount++;
-
- const CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
- // position and time of interaction, not used in NUCLIB
- Point pOrig = p.GetPosition();
- TimeType tOrig = p.GetTime();
-
- cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
- cout << "Interaction: time: " << tOrig << endl;
-
- // projectile nucleon number
- const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
- if (kAProj > 56)
- throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
-
- // kinematics
- // define projectile nucleus
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- auto const pProjectileLab = p.GetMomentum();
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
-
- cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
- << "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV
- << endl;
-
- // define projectile nucleon
- HEPEnergyType const eProjectileNucLab = p.GetEnergy() / kAProj;
- auto const pProjectileNucLab = p.GetMomentum() / kAProj;
- const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
-
- cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV
- << endl
- << "NuclearInteraction: pProjNucleon lab: "
- << pProjectileNucLab.GetComponents() / 1_GeV << endl;
-
- // define target
- // always a nucleon
- auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
- // target is always at rest
- const auto eTargetNucLab = 0_GeV + nucleon_mass;
- const auto pTargetNucLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
-
- cout << "NuclearInteraction: etarget lab: " << eTargetNucLab / 1_GeV << endl
- << "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV
- << endl;
-
- // center-of-mass energy in nucleon-nucleon frame
- auto const PtotNN4 = PtargNucLab + PprojNucLab;
- HEPEnergyType EcmNN = PtotNN4.GetNorm();
- cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
-
- if (!fHadronicInteraction.ValidCoMEnergy(EcmNN)) {
- cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
- "interaction model!"
- << endl;
- throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
- }
-
- // define boost to NUCLEON-NUCLEON frame
- COMBoost const boost(PprojNucLab, nucleon_mass);
- // boost projecticle
- auto const PprojNucCoM = boost.toCoM(PprojNucLab);
-
- // boost target
- auto const PtargNucCoM = boost.toCoM(PtargNucLab);
-
- cout << "Interaction: ebeam CoM: " << PprojNucCoM.GetTimeLikeComponent() / 1_GeV
- << endl
- << "Interaction: pbeam CoM: "
- << PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
- cout << "Interaction: etarget CoM: " << PtargNucCoM.GetTimeLikeComponent() / 1_GeV
- << endl
- << "Interaction: ptarget CoM: "
- << PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
-
- // sample target nucleon number
- //
- // proton stand-in for nucleon
- const auto beamId = corsika::particles::Proton::GetCode();
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
- const auto& mediumComposition =
- currentNode->GetModelProperties().GetNuclearComposition();
- cout << "get nucleon-nucleus cross sections for target materials.." << endl;
- // get cross sections for target materials
- // using nucleon-target-nucleus cross section!!!
- /*
- Here we read the cross section from the interaction model again,
- should be passed from GetInteractionLength if possible
- */
- 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];
- cout << "target component: " << targetId << endl;
- cout << "beam id: " << beamId << endl;
- const auto [sigProd, sigEla, nNuc] =
- fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
- cross_section_of_components[i] = sigProd;
- [[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
- [[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
- }
-
- const auto targetCode = currentNode->GetModelProperties().SampleTarget(
- cross_section_of_components, fRNG);
- cout << "Interaction: target selected: " << targetCode << endl;
- /*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be
- allowed air in atmosphere also contains some Argon.
- */
- int kATarget = -1;
- if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
- if (targetCode == corsika::particles::Proton::GetCode()) kATarget = 1;
- cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
- if (kATarget > 18 || kATarget < 1)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 or protons are "
- "allowed.");
- // end of target sampling
-
- // superposition
- cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "
- << endl;
- // get nucleon-nucleon cross section
- // (needed to determine number of nucleon-nucleon scatterings)
- const auto protonId = corsika::particles::Proton::GetCode();
- const auto [prodCrossSection, elaCrossSection, dum] =
- fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
- [[maybe_unused]] auto dumCopy = dum; // ONLY TO AVOID COMPILER WARNING
- const double sigProd = prodCrossSection / 1_mbarn;
- const double sigEla = elaCrossSection / 1_mbarn;
- // sample number of interactions (only input variables, output in common cnucms)
- // nuclear multiple scattering according to glauber (r.i.p.)
- int_nuc_(kATarget, kAProj, sigProd, sigEla);
-
- cout << "number of nucleons in target : " << kATarget << endl
- << "number of wounded nucleons in target : " << cnucms_.na << endl
- << "number of nucleons in projectile : " << kAProj << endl
- << "number of wounded nucleons in project. : " << cnucms_.nb << endl
- << "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl
- << "number of elastic nucleons in target : " << cnucms_.nael << endl
- << "number of elastic nucleons in project. : " << cnucms_.nbel << endl
- << "impact parameter: " << cnucms_.b << endl;
-
- // calculate fragmentation
- cout << "calculating nuclear fragments.." << endl;
- // number of interactions
- // include elastic
- const int nElasticNucleons = cnucms_.nbel;
- const int nInelNucleons = cnucms_.nb;
- const int nIntProj = nInelNucleons + nElasticNucleons;
- const double impactPar = cnucms_.b; // only needed to avoid passing common var.
- int nFragments;
- // number of fragments is limited to 60
- int AFragments[60];
- // call fragmentation routine
- // input: target A, projectile A, number of int. nucleons in projectile, impact
- // parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
- // in pf in common fragments, neglected
- fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
-
- // this should not occur but well :)
- if (nFragments > 60)
- throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
-
- cout << "number of fragments: " << nFragments << endl;
- for (int j = 0; j < nFragments; ++j)
- cout << "fragment: " << j << " A=" << AFragments[j]
- << " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1]
- << " pz=" << fragments_.ppp[j][2] << endl;
-
- cout << "adding nuclear fragments to particle stack.." << endl;
- // put nuclear fragments on corsika stack
- for (int j = 0; j < nFragments; ++j) {
- corsika::particles::Code specCode;
- const auto nuclA = AFragments[j];
- // get Z from stability line
- const auto nuclZ = int(nuclA / 2.15 + 0.7);
-
- // TODO: do we need to catch single nucleons??
- if (nuclA == 1)
- // TODO: sample neutron or proton
- specCode = corsika::particles::Code::Proton;
- else
- specCode = corsika::particles::Code::Nucleus;
-
- // TODO: mass of nuclei?
- const HEPMassType mass =
- corsika::particles::Proton::GetMass() * nuclZ +
- (nuclA - nuclZ) *
- corsika::particles::Neutron::GetMass(); // this neglects binding energy
-
- cout << "NuclearInteraction: adding fragment: " << specCode << endl;
- cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << endl;
- cout << "NuclearInteraction: mass: " << mass / 1_GeV << endl;
-
- // CORSIKA 7 way
- // spectators inherit momentum from original projectile
- const double mass_ratio = mass / ProjMass;
-
- cout << "NuclearInteraction: mass ratio " << mass_ratio << endl;
-
- auto const Plab = PprojLab * mass_ratio;
-
- cout << "NuclearInteraction: fragment momentum: "
- << Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
-
- if (nuclA == 1)
- // add nucleon
- p.AddSecondary(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType>{
- specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
- pOrig, tOrig});
- else
- // add nucleus
- p.AddSecondary(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType, unsigned short, unsigned short>{
- specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
- pOrig, tOrig, nuclA, nuclZ});
- }
-
- // add elastic nucleons to corsika stack
- // TODO: the elastic interaction could be external like the inelastic interaction,
- // e.g. use existing ElasticModel
- cout << "adding elastically scattered nucleons to particle stack.." << endl;
- for (int j = 0; j < nElasticNucleons; ++j) {
- // TODO: sample proton or neutron
- auto const elaNucCode = corsika::particles::Code::Proton;
-
- // CORSIKA 7 way
- // elastic nucleons inherit momentum from original projectile
- // neglecting momentum transfer in interaction
- const double mass_ratio = corsika::particles::GetMass(elaNucCode) / ProjMass;
- auto const Plab = PprojLab * mass_ratio;
-
- p.AddSecondary(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType>{
- elaNucCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
- pOrig, tOrig});
- }
-
- // add inelastic interactions
- cout << "calculate inelastic nucleon-nucleon interactions.." << endl;
- for (int j = 0; j < nInelNucleons; ++j) {
- // TODO: sample neutron or proton
- auto pCode = corsika::particles::Proton::GetCode();
- // temporarily add to stack, will be removed after interaction in DoInteraction
- cout << "inelastic interaction no. " << j << endl;
- auto inelasticNucleon = p.AddSecondary(
- std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
- corsika::stack::MomentumVector, corsika::geometry::Point,
- corsika::units::si::TimeType>{
- pCode, PprojNucLab.GetTimeLikeComponent(),
- PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
- // create inelastic interaction
- cout << "calling HadronicInteraction..." << endl;
- fHadronicInteraction.DoInteraction(inelasticNucleon, s);
- }
-
- // delete parent particle
- p.Delete();
-
- cout << "NuclearInteraction: DoInteraction: done" << endl;
-
- return process::EProcessReturn::eOk;
- }
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s);
private:
corsika::environment::Environment const& fEnvironment;
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 4af036ccd7d48f84da06c731937b2275c25f6ef7..9021bf2796ce50768ca8f7c9c13a7d971f193c9c 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -1,3 +1,14 @@
+
+/*
+ * (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/tracking_line/TrackingLine.h>
#include <corsika/environment/Environment.h>
@@ -126,8 +137,6 @@ namespace corsika::process::tracking_line {
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-// using Track = Trajectory;
template class corsika::process::tracking_line::TrackingLine<setup::Stack,
setup::Trajectory>;
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index b5e715cd21dba7814ca4143e3e59b6f199d03d04..cfc3568079bf7ad6433f2f28aee416230d9d5866 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -31,7 +31,7 @@ namespace corsika::process {
template <typename Stack, typename Trajectory>
class TrackingLine { //
- using Particle = typename Stack::ParticleType;
+ using Particle = typename Stack::StackIterator;
corsika::environment::Environment const& fEnvironment;
diff --git a/Processes/TrackingLine/testTrackingLineStack.cc b/Processes/TrackingLine/testTrackingLineStack.cc
index 0a718634a7f15cb9ca5a53c6675ecda9fb94c4e3..4c21923884ff69a76ec9acdc7428a313b1559734 100644
--- a/Processes/TrackingLine/testTrackingLineStack.cc
+++ b/Processes/TrackingLine/testTrackingLineStack.cc
@@ -1,3 +1,14 @@
+
+/*
+ * (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 _include_process_trackinling_teststack_h_
#define _include_process_trackinling_teststack_h_
diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h
index 6d8a46c1a4ca6ac5bcec0bca04a5a767b7f7fbc5..b825bc2f17b684e5fd753557296015d748181e57 100644
--- a/Processes/TrackingLine/testTrackingLineStack.h
+++ b/Processes/TrackingLine/testTrackingLineStack.h
@@ -1,3 +1,14 @@
+
+/*
+ * (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 _include_process_trackinling_teststack_h_
#define _include_process_trackinling_teststack_h_
@@ -28,6 +39,7 @@ struct DummyParticle {
struct DummyStack {
using ParticleType = DummyParticle;
+ using StackIterator = DummyParticle;
};
#endif