diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 943bbdd2b3cf67d2e26e0427f896ab3f58fcba2a..6529b5153565a74ddc39fc16b3e7c90f882f7aac 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -34,7 +34,6 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
- ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
@@ -45,6 +44,24 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example)
+add_executable (boundary_example boundary_example.cc)
+target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
+target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ CORSIKAcascade
+ ProcessTrackWriter
+ ProcessTrackingLine
+ CORSIKAprocesses
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+
+install (TARGETS boundary_example DESTINATION share/examples)
+CORSIKA_ADD_TEST (boundary_example)
+
add_executable (cascade_proton_example cascade_proton_example.cc)
target_compile_options(cascade_proton_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
@@ -63,6 +80,7 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
CORSIKAenvironment
CORSIKAprocesssequence
)
+
install (TARGETS cascade_proton_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_proton_example)
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
new file mode 100644
index 0000000000000000000000000000000000000000..98e01b89f0aaeaa556e8d0364b2be3fff4498ecc
--- /dev/null
+++ b/Documentation/Examples/boundary_example.cc
@@ -0,0 +1,343 @@
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/cascade/Cascade.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+#include <corsika/geometry/Sphere.h>
+
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/process/track_writer/TrackWriter.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::particles;
+using namespace corsika::random;
+using namespace corsika::setup;
+using namespace corsika::geometry;
+using namespace corsika::environment;
+
+using namespace std;
+using namespace corsika::units::si;
+
+class ProcessCut : public process::ContinuousProcess<ProcessCut> {
+
+ HEPEnergyType fECut;
+
+ HEPEnergyType fEnergy = 0_GeV;
+ HEPEnergyType fEmEnergy = 0_GeV;
+ int fEmCount = 0;
+ HEPEnergyType fInvEnergy = 0_GeV;
+ int fInvCount = 0;
+
+public:
+ ProcessCut(const HEPEnergyType cut)
+ : fECut(cut) {}
+
+ template <typename Particle>
+ bool isBelowEnergyCut(Particle& p) const {
+ // nuclei
+ if (p.GetPID() == particles::Code::Nucleus) {
+ auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
+ auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
+ if (ElabNuc < fECut || EcmNN < 10_GeV)
+ return true;
+ else
+ return false;
+ } else {
+ // TODO: center-of-mass energy hard coded
+ const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ if (p.GetEnergy() < fECut || Ecm < 10_GeV)
+ return true;
+ else
+ return false;
+ }
+ }
+
+ bool isEmParticle(Code pCode) const {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::Electron:
+ is_em = true;
+ break;
+ case Code::Positron:
+ is_em = true;
+ break;
+ case Code::Gamma:
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible(Code pCode) const {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::NuE:
+ is_inv = true;
+ break;
+ case Code::NuEBar:
+ is_inv = true;
+ break;
+ case Code::NuMu:
+ is_inv = true;
+ break;
+ case Code::NuMuBar:
+ is_inv = true;
+ break;
+ case Code::MuPlus:
+ is_inv = true;
+ break;
+ case Code::MuMinus:
+ is_inv = true;
+ break;
+
+ case Code::Neutron:
+ is_inv = true;
+ break;
+
+ case Code::AntiNeutron:
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+ template <typename Particle>
+ LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
+ cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
+ cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
+ const Code pid = p.GetPID();
+ if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0_m;
+ } else {
+ LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
+ const Code pid = p.GetPID();
+ HEPEnergyType energy = p.GetEnergy();
+ cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
+ << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
+ EProcessReturn ret = EProcessReturn::eOk;
+ if (isEmParticle(pid)) {
+ cout << "removing em. particle..." << endl;
+ fEmEnergy += energy;
+ fEmCount += 1;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ } else if (isInvisible(pid)) {
+ cout << "removing inv. particle..." << endl;
+ fInvEnergy += energy;
+ fInvCount += 1;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ } else if (isBelowEnergyCut(p)) {
+ cout << "removing low en. particle..." << endl;
+ fEnergy += energy;
+ // p.Delete();
+ ret = EProcessReturn::eParticleAbsorbed;
+ }
+ return ret;
+ }
+
+ void Init() {
+ fEmEnergy = 0. * 1_GeV;
+ fEmCount = 0;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
+ }
+
+ void ShowResults() {
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
+ << " ******************************" << endl;
+ }
+
+ HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+ HEPEnergyType GetCutEnergy() const { return fEnergy; }
+ HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+};
+
+template <bool deleteParticle>
+struct MyBoundaryCrossingProcess
+ : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
+
+ MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
+
+ template <typename Particle>
+ EProcessReturn DoBoundaryCrossing(Particle& p,
+ typename Particle::BaseNodeType const& from,
+ typename Particle::BaseNodeType const& to) {
+ std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
+
+ auto const& name = particles::GetName(p.GetPID());
+ auto const start = p.GetPosition().GetCoordinates();
+
+ fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
+ << start[2] / 1_m << '\n';
+
+ if constexpr (deleteParticle) { p.Delete(); }
+
+ return EProcessReturn::eOk;
+ }
+
+ void Init() {}
+
+private:
+ std::ofstream fFile;
+};
+
+//
+// The example main program for a particle cascade
+//
+int main() {
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+
+ // setup environment, geometry
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ auto& universe = *(env.GetUniverse());
+
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ auto outerMedium = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ // fraction of oxygen
+ auto const props =
+ outerMedium
+ ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Oxygen},
+ std::vector<float>{1.f}));
+
+ auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
+
+ innerMedium->SetModelProperties(props);
+
+ outerMedium->AddChild(std::move(innerMedium));
+
+ universe.AddChild(std::move(outerMedium));
+
+ // setup processes, decays and interactions
+ tracking_line::TrackingLine tracking;
+
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+ process::sibyll::Interaction sibyll;
+ process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
+ particles::Code::KPlus, particles::Code::KMinus,
+ particles::Code::K0Long, particles::Code::K0Short}};
+ ProcessCut cut(20_GeV);
+
+ process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+ MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
+
+ // assemble all processes into an ordered process list
+ auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
+
+ // setup particle stack, and add primary particles
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Proton;
+ const HEPMassType mass = particles::GetMass(Code::Proton);
+ const HEPEnergyType E0 = 50_TeV;
+
+ std::uniform_real_distribution distTheta(0., 180.);
+ std::uniform_real_distribution distPhi(0., 360.);
+ std::mt19937 rng;
+
+ for (int i = 0; i < 100; ++i) {
+ auto const theta = distTheta(rng);
+ auto const phi = distPhi(rng);
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+ Point pos(rootCS, 0_m, 0_m, 0_m);
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, pos, 0_ns});
+ }
+
+ // define air shower object, run simulation
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+ EAS.Init();
+ EAS.Run();
+
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
+ cut.ShowResults();
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+}
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c154db988f183799d3c18f3dc8033bd8f600c2d3..9855dfcb8bfa227c1b5da8487d73fe0d64b4941f 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -12,10 +12,10 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -39,9 +39,6 @@
#include <corsika/utl/CorsikaFenv.h>
-#include <boost/type_index.hpp>
-using boost::typeindex::type_id_with_cvr;
-
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -229,28 +226,36 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- environment::Environment env;
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ auto outerMedium = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Nitrogen, particles::Code::Oxygen},
- std::vector<float>{(float)1.-fox, fox}));
+ auto const props =
+ outerMedium
+ ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Nitrogen,
+ particles::Code::Oxygen},
+ std::vector<float>{1.f - fox, fox}));
- universe.AddChild(std::move(theMedium));
+ auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
- const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ innerMedium->SetModelProperties(props);
+
+ outerMedium->AddChild(std::move(innerMedium));
+
+ universe.AddChild(std::move(outerMedium));
// setup processes, decays and interactions
- tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
+ tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
@@ -259,29 +264,16 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- process::sibyll::Interaction sibyll(env);
- process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+ process::sibyll::Interaction sibyll;
+ process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
- // random::RNGManager::GetInstance().RegisterRandomStream("pythia");
- // process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
- // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
- // process::HadronicElasticModel::HadronicElasticInteraction
- // hadronicElastic(env);
-
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
- // auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut <<
- // trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss <<
- // cut << trackWriter;
- // auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
- auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
-
- // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
- // << "\n";
+ auto sequence = sibyll << sibyllNuc << decay << cut << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index b8851ffe622e1583a2ce55e60e7dd4a112065959..97e6c40ac5f4b80fd2a64430eec8b19fe81867ed 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -229,12 +229,13 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- environment::Environment env;
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ auto theMedium =
+ EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
@@ -248,18 +249,17 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
- tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
- stack_inspector::StackInspector<setup::Stack> p0(true);
+ tracking_line::TrackingLine tracking;
+ stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
- process::pythia::Interaction pythia(env);
+ process::pythia::Interaction pythia;
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay(trackedHadrons);
process::pythia::Decay decay(trackedHadrons);
@@ -274,8 +274,8 @@ int main() {
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
// auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
-
- auto sequence = p0 << pythia << decay << cut << trackWriter;
+
+ auto sequence = stackInspect << pythia << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
@@ -285,7 +285,7 @@ int main() {
stack.Clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::Proton::GetMass();
- const HEPEnergyType E0 = 100_GeV;
+ const HEPEnergyType E0 = 100_GeV;
double theta = 0.;
double phi = 0.;
@@ -305,10 +305,10 @@ int main() {
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
- stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
- beamCode, E0, plab, pos, 0_ns});
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, pos, 0_ns});
}
// define air shower object, run simulation
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 47b973cb138cbc7a52ba11e36b8347ec8497e883..0a94c410aeebbf353d225e9e7ea2822d814b1728 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -9,6 +9,7 @@ set (
LinearApproximationIntegrator.h
DensityFunction.h
Environment.h
+ NameModel.h
FlatExponential.h
)
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 97931aeffe2fe557b0328f97bbea8712a5dffc7e..d7ac8d14a98e19d852c7151c0069dc928d018e39 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -9,20 +9,17 @@
* the license.
*/
-#ifndef _include_Environment_h
-#define _include_Environment_h
+#ifndef _include_environment_Environment_h
+#define _include_environment_Environment_h
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
-#include <corsika/setup/SetupEnvironment.h>
#include <limits>
namespace corsika::environment {
- using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
-
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
@@ -34,16 +31,18 @@ namespace corsika::environment {
bool Contains(corsika::geometry::Point const&) const override { return true; }
};
- // template <typename IEnvironmentModel>
+ template <typename IEnvironmentModel>
class Environment {
public:
+ using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
+
Environment()
: fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem()}
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
- using IEnvironmentModel = corsika::setup::IEnvironmentModel;
+ // using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
@@ -63,9 +62,12 @@ namespace corsika::environment {
private:
corsika::geometry::CoordinateSystem const& fCoordinateSystem;
- BaseNodeType::VTNUPtr fUniverse;
+ typename BaseNodeType::VTNUPtr fUniverse;
};
+ // using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
+ // using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
+
} // namespace corsika::environment
#endif
diff --git a/Environment/NameModel.h b/Environment/NameModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..458c07819e75c6553c6fda4235bd445605218ec6
--- /dev/null
+++ b/Environment/NameModel.h
@@ -0,0 +1,27 @@
+/*
+ * (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_NameModel_h
+#define _include_NameModel_h
+
+#include <string>
+#include <utility>
+
+namespace corsika::environment {
+
+ template <typename T>
+ struct NameModel : public T {
+ virtual std::string const& GetName() const = 0;
+ virtual ~NameModel() = default;
+ };
+
+} // namespace corsika::environment
+
+#endif
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 9cd3be1a9e534c5f6b74bb4e584dac1c78e0d9bf..faadc94c6a070876d34653fcf5cdebfb9482b7af 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -12,6 +11,7 @@
#ifndef _include_VolumeTreeNode_H
#define _include_VolumeTreeNode_H
+#include <corsika/environment/IMediumModel.h>
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
@@ -20,9 +20,10 @@ namespace corsika::environment {
class Empty {}; //<! intended for usage as default template argument
- template <typename IModelProperties = Empty>
+ template <typename TModelProperties = Empty>
class VolumeTreeNode {
public:
+ using IModelProperties = TModelProperties;
using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
@@ -30,6 +31,7 @@ namespace corsika::environment {
VolumeTreeNode(VolUPtr pVolume = nullptr)
: fGeoVolume(std::move(pVolume)) {}
+ //! convenience function equivalent to Volume::Contains
bool Contains(corsika::geometry::Point const& p) const {
return fGeoVolume->Contains(p);
}
@@ -44,7 +46,7 @@ namespace corsika::environment {
}
/** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
- * \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
+ * \class Point \p p, or nullptr iff \p p is not contained in this volume.
*/
VolumeTreeNode<IModelProperties> const* GetContainingNode(
corsika::geometry::Point const& p) const {
@@ -106,14 +108,14 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
- IMPSharedPtr GetModelPropertiesPtr() const { return fModelProperties; }
+ auto HasModelProperties() const { return fModelProperties != nullptr; }
- template <typename TModelProperties, typename... Args>
+ template <typename ModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
- static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
+ static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
"unusable type provided");
- fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
+ fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
return fModelProperties;
}
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index 646692248b001e645af842ec70147a906e753571..3e294b214dd325ee9c082c32e2e0d46cdd70de3b 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -9,6 +9,7 @@ set (
set (
CORSIKAcascade_HEADERS
Cascade.h
+ testCascade.h
)
add_library (CORSIKAcascade INTERFACE)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 1a559a210eb16c60f3d8165b0b18de5589dbdbd3..c7a187447a5d3f6034600a05bffe09e3543cd7f0 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -20,9 +20,11 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
+#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
+#include <type_traits>
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
@@ -38,14 +40,12 @@ namespace corsika::cascade {
* plugged into the cascade simulation.
*
* <b>Tracking</b> must be a class according to the
- * TrackingInterface providing the functions: <code>void
- * Init();</code> and <code>auto GetTrack(Particle const& p)</auto>,
- * where the latter has a return type of <code>
- * geometry::Trajectory<corsika::geometry::Line or Helix> </code>
- *
- * <b>ProcessList</b> must be a ProcessSequence.
- * TimeOfIntersection(corsika::geometry::Line const& line,
+ * TrackingInterface providing the functions:
+ * <code>auto GetTrack(Particle const& p)</auto>,
+ * with the return type <code>geometry::Trajectory<corsika::geometry::Line>
+ * </code>
*
+ * <b>ProcessList</b> must be a ProcessSequence. *
* <b>Stack</b> is the storage object for particle data, i.e. with
* Particle class type <code>Stack::ParticleType</code>
*
@@ -56,6 +56,9 @@ namespace corsika::cascade {
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
using Particle = typename Stack::ParticleType;
+ using VolumeTreeNode =
+ std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
+ using MediumInterface = typename VolumeTreeNode::IModelProperties;
// we only want fully configured objects
Cascade() = delete;
@@ -65,8 +68,8 @@ namespace corsika::cascade {
* Cascade class cannot be default constructed, but needs a valid
* list of physics processes for configuration at construct time.
*/
- Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
- Stack& stack)
+ Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr,
+ ProcessList& pl, Stack& stack)
: fEnvironment(env)
, fTracking(tr)
, fProcessSequence(pl)
@@ -77,16 +80,29 @@ namespace corsika::cascade {
* All components of the Cascade simulation must be configured here.
*/
void Init() {
- fTracking.Init();
fProcessSequence.Init();
fStack.Init();
}
+ /**
+ * set the nodes for all particles on the stack according to their numerical
+ * position
+ */
+ void SetNodes() {
+ std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
+ auto const* numericalNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ p.SetNode(numericalNode);
+ });
+ }
+
/**
* The Run function is the main simulation loop, which processes
* particles from the Stack until the Stack is empty.
*/
void Run() {
+ SetNodes();
+
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
auto pNext = fStack.GetNextParticle();
@@ -113,7 +129,7 @@ namespace corsika::cascade {
using namespace corsika::units::si;
// determine geometric tracking
- corsika::setup::Trajectory step = fTracking.GetTrack(particle);
+ auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
@@ -126,19 +142,17 @@ namespace corsika::cascade {
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
- // convert next_step from grammage to length
- auto const* currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+ auto const* currentLogicalNode = particle.GetNode();
- if (currentNode == &*fEnvironment.GetUniverse()) {
- throw std::runtime_error("particle entered void universe");
- }
+ // assert that particle stays outside void Universe if it has no
+ // model properties set
+ assert(currentLogicalNode != &*fEnvironment.GetUniverse() ||
+ fEnvironment.GetUniverse()->HasModelProperties());
+ // convert next_step from grammage to length
LengthType const distance_interact =
- (next_interact == std::numeric_limits<double>::infinity() * 1_g / (1_m * 1_m))
- ? std::numeric_limits<double>::infinity() * 1_m
- : currentNode->GetModelProperties().ArclengthFromGrammage(step,
- next_interact);
+ currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
+ next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
@@ -161,7 +175,7 @@ namespace corsika::cascade {
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
- std::min({distance_interact, distance_decay, distance_max});
+ std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
std::cout << " move particle by : " << min_distance << std::endl;
@@ -172,13 +186,6 @@ namespace corsika::cascade {
step.LimitEndTo(min_distance);
- // particle.GetNode(); // previous VolumeNode
- particle.SetNode(
- currentNode); // NOTE @Max : here we need to distinguish: IF particle step is
- // limited by tracking (via fTracking.GetTrack()), THEN we need
- // to check/update VolumeNodes. In all other cases it is
- // guaranteed that we are still in the same volume
-
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcessSequence.DoContinuous(particle, step, fStack);
@@ -191,11 +198,9 @@ namespace corsika::cascade {
}
std::cout << "sth. happening before geometric limit ? "
- << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
-
- if (min_distance < distance_max) { // interaction to happen within geometric limit
- // check whether decay or interaction limits this step
+ << ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
+ if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
@@ -208,7 +213,7 @@ namespace corsika::cascade {
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
inv_lambda_count);
- } else {
+ } else if (min_distance == distance_decay) {
std::cout << "decay" << std::endl;
InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
@@ -218,18 +223,33 @@ namespace corsika::cascade {
const auto sample_process = uniDist(fRNG);
InverseTimeType inv_decay_count = 0 / second;
fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
+ } else { // step-length limitation within volume
+ std::cout << "step-length limitation" << std::endl;
}
+
+ auto const assertion = [&] {
+ auto const* numericalNodeAfterStep =
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+ return numericalNodeAfterStep == currentLogicalNode;
+ };
+
+ assert(assertion()); // numerical and logical nodes don't match
+ } else { // boundary crossing, step is limited by volume boundary
+ std::cout << "boundary crossing! next node = " << nextVol << std::endl;
+ particle.SetNode(nextVol);
+ // DoBoundary may delete the particle (or not)
+ fProcessSequence.DoBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
}
}
private:
- corsika::environment::Environment const& fEnvironment;
+ corsika::environment::Environment<MediumInterface> const& fEnvironment;
Tracking& fTracking;
ProcessList& fProcessSequence;
Stack& fStack;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
- };
+ }; // namespace corsika::cascade
} // namespace corsika::cascade
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 1e8a7a4deef47ccd8bd0f692b95951496afb743f..5ec522a670de7ea048d33397e868590ea9e6faed 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -11,7 +11,7 @@
#include <limits>
-#include <corsika/environment/Environment.h>
+#include <corsika/cascade/testCascade.h>
#include <corsika/cascade/Cascade.h>
@@ -25,7 +25,6 @@
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
-#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
@@ -46,11 +45,11 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
-environment::Environment MakeDummyEnv() {
- environment::Environment env; // dummy environment
+auto MakeDummyEnv() {
+ TestEnvironmentType env; // dummy environment
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<Sphere>(
+ auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
@@ -113,14 +112,14 @@ TEST_CASE("Cascade", "[Cascade]") {
rmng.RegisterRandomStream("cascade");
auto env = MakeDummyEnv();
- tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
+ tracking_line::TrackingLine tracking;
- stack_inspector::StackInspector<setup::Stack> p0(true);
+ stack_inspector::StackInspector<TestCascadeStack> p0(true);
const HEPEnergyType Ecrit = 85_MeV;
ProcessSplit p1(Ecrit);
auto sequence = p0 << p1;
- setup::Stack stack;
+ TestCascadeStack stack;
cascade::Cascade EAS(env, tracking, sequence, stack);
CoordinateSystem const& rootCS =
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
new file mode 100644
index 0000000000000000000000000000000000000000..c84d1e14bff303274e648a4d018e711c4ff9dfd5
--- /dev/null
+++ b/Framework/Cascade/testCascade.h
@@ -0,0 +1,33 @@
+
+/*
+ * (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 _Framework_Cascade_testCascade_h
+#define _Framework_Cascade_testCascade_h
+
+#include <corsika/environment/Environment.h>
+#include <corsika/setup/SetupStack.h>
+
+using TestEnvironmentType =
+ corsika::environment::Environment<corsika::environment::IMediumModel>;
+
+template <typename T>
+using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
+
+// combine particle data stack with geometry information for tracking
+template <typename StackIter>
+using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<
+ corsika::setup::detail::ParticleDataStack::PIType, SetupGeometryDataInterface,
+ StackIter>;
+using TestCascadeStack = corsika::stack::CombinedStack<
+ typename corsika::setup::detail::ParticleDataStack::StackImpl,
+ GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
+
+#endif
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index ab57a994db6cba3cce60ab3acdb028b69de4c3b1..2712c188d8daf8c65f354bf416711c7d4d71e32c 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -118,7 +118,7 @@ namespace corsika::geometry {
auto& operator-() const { return QuantityVector<dim>(-eVector); }
- auto normalized() const { return (*this) * (1 / norm()); }
+ auto normalized() const { return QuantityVector<dim>(eVector.normalized()); }
auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
};
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index e8a1a58bf0d7cdaa8a6fab5836a592c674ad3a7d..559dbd562694a6a38cc150fab71bee52059b9062 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -12,6 +11,7 @@
#ifndef _include_TRAJECTORY_H
#define _include_TRAJECTORY_H
+#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
diff --git a/Framework/ProcessSequence/BoundaryCrossingProcess.h b/Framework/ProcessSequence/BoundaryCrossingProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..7ca73cb3644fb70ecdb09bb356d3773feba14762
--- /dev/null
+++ b/Framework/ProcessSequence/BoundaryCrossingProcess.h
@@ -0,0 +1,35 @@
+/*
+ * (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_corsika_processes_BoundaryCrossingProcess_h_
+#define _include_corsika_processes_BoundaryCrossingProcess_h_
+
+#include <corsika/environment/Environment.h>
+#include <corsika/process/ProcessReturn.h>
+
+namespace corsika::process {
+ template <typename TDerived>
+ struct BoundaryCrossingProcess {
+ auto& GetRef() { return static_cast<TDerived&>(*this); }
+ auto const& GetRef() const { return static_cast<const TDerived&>(*this); }
+
+ /**
+ * This method is called when a particle crosses the boundary between the nodes
+ * \p from and \p to.
+ */
+ template <typename Particle, typename VTNType>
+ EProcessReturn DoBoundaryCrossing(Particle&, VTNType const& from, VTNType const& to);
+ };
+
+ template <class T>
+ std::true_type is_process_impl(BoundaryCrossingProcess<T> const* impl);
+} // namespace corsika::process
+
+#endif
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index d45d28d6c208ef80302130bbe3e2e1c837b1dfcf..9ec9072dabb1dab6ecb17213ddc84a45bb73cd5f 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -11,6 +11,7 @@ set (
set (
CORSIKAprocesssequence_HEADERS
BaseProcess.h
+ BoundaryCrossingProcess.h
ContinuousProcess.h
InteractionProcess.h
DecayProcess.h
@@ -33,6 +34,12 @@ install (
FILES ${CORSIKAprocesssequence_HEADERS}
DESTINATION include/${CORSIKAprocesssequence_NAMESPACE}
)
+
+target_link_libraries (
+ CORSIKAprocesssequence
+ INTERFACE
+ CORSIKAenvironment
+)
#-- -- -- -- -- -- -- --
#code unit testing
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
index 96472755376c4cf630f8ac92e92c5c08eb1c2f48..723cbf02fcc7ce443d968300293ccc44d1b40811 100644
--- a/Framework/ProcessSequence/InteractionProcess.h
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -36,7 +36,7 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoInteraction...
template <typename P, typename S>
- inline EProcessReturn DoInteraction(P&, S&);
+ EProcessReturn DoInteraction(P&, S&);
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t);
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 1056ccd9e4cdb4ec80fca29b7c08857de1e28fc9..c697d09f416d49025a0d4763f00851c047cedb62 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -13,6 +13,7 @@
#define _include_ProcessSequence_h_
#include <corsika/process/BaseProcess.h>
+#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
@@ -70,6 +71,24 @@ namespace corsika::process {
// example for a trait-based call:
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
+ template <typename Particle, typename VTNType>
+ EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from,
+ VTNType const& to) {
+ EProcessReturn ret = EProcessReturn::eOk;
+
+ if constexpr (std::is_base_of<BoundaryCrossingProcess<T1type>, T1type>::value ||
+ is_process_sequence<T1>::value) {
+ ret |= A.DoBoundaryCrossing(p, from, to);
+ }
+
+ if constexpr (std::is_base_of<BoundaryCrossingProcess<T2type>, T2type>::value ||
+ is_process_sequence<T2>::value) {
+ ret |= B.DoBoundaryCrossing(p, from, to);
+ }
+
+ return ret;
+ }
+
template <typename Particle, typename Track, typename Stack>
EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) {
EProcessReturn ret = EProcessReturn::eOk;
diff --git a/Framework/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt
index 84479e1230a92ef5a1773ff75ebdf72347c3f7db..bbbe8d32d484373f6ed62495893508c5996dd59f 100644
--- a/Framework/Random/CMakeLists.txt
+++ b/Framework/Random/CMakeLists.txt
@@ -19,6 +19,13 @@ set (
add_library (CORSIKArandom STATIC ${CORSIKArandom_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKArandom ${CORSIKArandom_NAMESPACE} ${CORSIKArandom_HEADERS})
+target_link_libraries (
+ CORSIKArandom
+ INTERFACE
+ CORSIKAutilities
+ CORSIKAunits
+ )
+
target_include_directories (
CORSIKArandom
PUBLIC
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index d145c1e88610579c8b8ef2ddc219a261b48b695d..b5704e762a667b8ff7d4c7442a068f561d2dee8a 100644
--- a/Framework/StackInterface/ParticleBase.h
+++ b/Framework/StackInterface/ParticleBase.h
@@ -90,7 +90,7 @@ namespace corsika::stack {
// protected:
/**
- @name Access to underlying stack dfata, these are service
+ @name Access to underlying stack fData, these are service
function for user classes. User code can only rely on GetIndex
and GetStackData to retrieve data
@{
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index 9d259e1c282b278ab9bcec4f37a337685e29d641..72f3a72fe2d9c5078f893b3f815ed880fd515118 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -33,8 +33,8 @@ namespace corsika::stack {
<b>Important:</b> ParticleInterface must inherit from ParticleBase !
*/
- template <typename> //, bool>
- class ParticleInterface; // forward decl
+ template <typename>
+ class ParticleInterface;
/**
The Stack class provides (and connects) the main particle data storage machinery.
@@ -98,8 +98,8 @@ namespace corsika::stack {
typedef StackDataType
StackImpl; ///< this is the type of the user-provided data structure
- template <typename SI> //, bool IsBase>
- using PIType = ParticleInterface<SI>; //, IsBase>;
+ template <typename SI>
+ using PIType = ParticleInterface<SI>;
/**
* Via the StackIteratorInterface and ConstStackIteratorInterface
@@ -209,7 +209,6 @@ namespace corsika::stack {
}
if (p.GetIndex() < GetSize() - 1) fData.Copy(GetSize() - 1, p.GetIndex());
DeleteLast();
- // p.SetInvalid();
}
/**
* delete this particle
diff --git a/Framework/Utilities/Singleton.h b/Framework/Utilities/Singleton.h
index 3780cc99a1398833a77729e26eb11f9b089ba753..80c272d58bff6cea52fbd2083fdfb7e4238b9b60 100644
--- a/Framework/Utilities/Singleton.h
+++ b/Framework/Utilities/Singleton.h
@@ -15,9 +15,6 @@
//#define OFFLINE_USE_GAMMA_SINGLETON
namespace corsika::utl {
-
-#ifndef OFFLINE_USE_GAMMA_SINGLETON
-
/**
* \class Singleton Singleton.h utl/Singleton.h
*
@@ -48,102 +45,18 @@ namespace corsika::utl {
template <typename T>
class Singleton {
public:
- static T& GetInstance()
-#ifdef __MAKECINT__
- ;
-#else
- {
+ static T& GetInstance() {
static T instance;
return instance;
}
-#endif
- protected:
- // derived class can call ctor and dtor
- Singleton() {}
- ~Singleton() {}
-
- private:
- // no one should do copies
- Singleton(const Singleton&);
- Singleton& operator=(const Singleton&);
- };
-
-#else
-
- /// classical Gamma singleton
- template <typename T>
- class Singleton {
- public:
- static T& GetInstance() {
- if (!fgInstance) fgInstance = new T;
- return *fgInstance;
- }
+ Singleton(const Singleton&) = delete;
+ Singleton& operator=(const Singleton&) = delete;
protected:
// derived class can call ctor and dtor
Singleton() {}
~Singleton() {}
-
- private:
- // no one should do copies
- Singleton(const Singleton&);
- Singleton& operator=(const Singleton&);
-
- static T* fgInstance = 0;
- };
-
-#endif
-
- /**
- * \class LeakingSingleton Singleton.h utl/Singleton.h
- *
- * \brief CRTP for leaking singleton
- *
- * This type of creation (Gamma singleton) leaks the object at
- * the end of the run, i.e. class T destructor does not get called
- * in at_exit().
- *
- * This singleton can be implemented as follows
- * \code
- * #include <utl/Singleton.h>
- *
- * class SomeClass : public utl::LeakingSingleton<SomeClass> {
- * ...
- * private:
- * // prevent creation, destruction
- * SomeClass() { }
- * ~SomeClass() { }
- *
- * friend class utl::LeakingSingleton<SomeClass>;
- * };
- * \endcode
- * LeakingSingleton automatically prevents copying of the derived
- * class.
- *
- * \author Darko Veberic
- * \date 9 Aug 2006
- * \version $Id: Singleton.h 25091 2014-01-30 09:49:57Z darko $
- * \ingroup stl
- */
-
- template <class T>
- class LeakingSingleton {
- public:
- static T& GetInstance() {
- static T* const instance = new T;
- return *instance;
- }
-
- protected:
- // derived class can call ctor and dtor
- LeakingSingleton() {}
- ~LeakingSingleton() {}
-
- private:
- // no one should do copies
- LeakingSingleton(const LeakingSingleton&);
- LeakingSingleton& operator=(const LeakingSingleton&);
};
} // namespace corsika::utl
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index a20158bb96623895c1af1454eae79d33740b6da5..960bbd3235e0de25b4898d5a68dd84813ca27ea4 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -125,8 +125,8 @@ namespace corsika::process::EnergyLoss {
// Barkas correction O(Z3) higher-order Born approximation
// see Appl. Phys. 85 (1999) 1249
- double A = 1;
- if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
+ // double A = 1;
+ // if (p.GetPID() == particles::Code::Nucleus) A = p.GetNuclearA();
// double const Erel = (p.GetEnergy()-p.GetMass()) / A / 1_keV;
// double const Llow = 0.01 * Erel;
// double const Lhigh = 1.5/pow(Erel, 0.4) + 45000./Zmat * pow(Erel, 1.6);
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index d04b8dcd60e7fa916aa092ce8d92253ac906d7a1..289968645c4c758d10f7230dbcf99babba40838c 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -24,27 +24,24 @@
#include <iostream>
using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
+using Particle = corsika::setup::Stack::ParticleType;
+using Track = corsika::setup::Trajectory;
namespace corsika::process::HadronicElasticModel {
void HadronicElasticInteraction::Init() {}
- HadronicElasticInteraction::HadronicElasticInteraction(
- environment::Environment const& env, units::si::CrossSectionType x,
- units::si::CrossSectionType y)
+ HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x,
+ units::si::CrossSectionType y)
: fX(x)
- , fY(y)
- , fEnvironment(env) {}
+ , fY(y) {}
template <>
units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
Particle const& p, Track&) {
using namespace units::si;
if (p.GetPID() == particles::Code::Proton) {
- auto const* currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ auto const* currentNode = p.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
@@ -65,8 +62,7 @@ namespace corsika::process::HadronicElasticModel {
avgCrossSection += CrossSection(s) * fractions[i];
}
- std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb"
- << std::endl;
+ std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl;
return avgCrossSection;
}();
@@ -89,8 +85,7 @@ namespace corsika::process::HadronicElasticModel {
using namespace units::si;
using namespace units::constants;
- const auto* currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto* currentNode = p.GetNode();
const auto& composition = currentNode->GetModelProperties().GetNuclearComposition();
const auto& components = composition.GetComponents();
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index 7f37905de1672f740f1ea7ec44051c5a730d1037..fae8c97242bac32426eca9089ac8a8ff611b631e 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -19,10 +19,6 @@
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
-namespace corsika::environment {
- class Environment;
-}
-
namespace corsika::process::HadronicElasticModel {
/**
* A simple model for elastic hadronic interactions based on the formulas
@@ -50,16 +46,14 @@ namespace corsika::process::HadronicElasticModel {
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream(
"HadronicElasticModel");
- corsika::environment::Environment const& fEnvironment;
inveV2 B(eV2 s) const;
corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const;
public:
- HadronicElasticInteraction(corsika::environment::Environment const&,
- // x & y values taken from DL for pp collisions
- units::si::CrossSectionType x = 0.0217 * units::si::barn,
- units::si::CrossSectionType y = 0.05608 * units::si::barn);
+ HadronicElasticInteraction( // x & y values taken from DL for pp collisions
+ units::si::CrossSectionType x = 0.0217 * units::si::barn,
+ units::si::CrossSectionType y = 0.05608 * units::si::barn);
void Init();
template <typename Particle, typename Track>
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index a57d967c04f065054e4f656652aca03feec15c60..092563e52bb66cf8130604100ed6a3f7e5b1c18d 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -32,13 +32,8 @@ using Track = Trajectory;
namespace corsika::process::pythia {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
-
- Interaction::Interaction(environment::Environment const& env)
- : fEnvironment(env) {}
- Interaction::~Interaction() {
- cout << "Pythia::Interaction n=" << fCount << endl;
- }
+ Interaction::~Interaction() { cout << "Pythia::Interaction n=" << fCount << endl; }
void Interaction::Init() {
@@ -49,58 +44,58 @@ namespace corsika::process::pythia {
fPythia.readString("Print:quiet = on");
// TODO: proper process initialization for MinBias needed
- fPythia.readString("HardQCD:all = on");
+ fPythia.readString("HardQCD:all = on");
fPythia.readString("ProcessLevel:resonanceDecays = off");
fPythia.init();
-
+
// any decays in pythia? if yes need to define which particles
- if(fInternalDecays){
- // define which particles are passed to corsika, i.e. which particles make it into history
- // even very shortlived particles like charm or pi0 are of interest here
- const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Pi0,
- particles::Code::KMinus, particles::Code::KPlus,
- particles::Code::K0Long, particles::Code::K0Short,
- particles::Code::SigmaPlus, particles::Code::SigmaMinus,
- particles::Code::Lambda0,
- particles::Code::Xi0, particles::Code::XiMinus,
- particles::Code::OmegaMinus,
- particles::Code::DPlus, particles::Code::DMinus, particles::Code::D0, particles::Code::D0Bar};
-
- Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
+ if (fInternalDecays) {
+ // define which particles are passed to corsika, i.e. which particles make it into
+ // history even very shortlived particles like charm or pi0 are of interest here
+ const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
+ particles::Code::PiPlus, particles::Code::PiMinus,
+ particles::Code::Pi0, particles::Code::KMinus,
+ particles::Code::KPlus, particles::Code::K0Long,
+ particles::Code::K0Short, particles::Code::SigmaPlus,
+ particles::Code::SigmaMinus, particles::Code::Lambda0,
+ particles::Code::Xi0, particles::Code::XiMinus,
+ particles::Code::OmegaMinus, particles::Code::DPlus,
+ particles::Code::DMinus, particles::Code::D0,
+ particles::Code::D0Bar};
+
+ Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
}
// basic initialization of cross section routines
- fSigma.init( &fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm );
+ fSigma.init(&fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm);
fInitialized = true;
}
}
- void Interaction::SetParticleListStable(const std::vector<particles::Code> particleList) {
- for (auto p : particleList)
- Interaction::SetStable( p );
+ void Interaction::SetParticleListStable(
+ std::vector<particles::Code> const& particleList) {
+ for (auto p : particleList) Interaction::SetStable(p);
}
void Interaction::SetUnstable(const particles::Code pCode) {
cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
- fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , true);
+ fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
}
void Interaction::SetStable(const particles::Code pCode) {
cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
- fPythia.particleData.mayDecay( static_cast<int>( particles::GetPDG(pCode) ) , false);
+ fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
}
-
- void Interaction::ConfigureLabFrameCollision( const particles::Code BeamId, const particles::Code TargetId,
- const units::si::HEPEnergyType BeamEnergy )
- {
+ void Interaction::ConfigureLabFrameCollision(
+ const particles::Code BeamId, const particles::Code TargetId,
+ const units::si::HEPEnergyType BeamEnergy) {
using namespace units::si;
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
-
+
// set beam
// beam id for pythia
auto const pdgBeam = static_cast<int>(particles::GetPDG(BeamId));
@@ -110,8 +105,8 @@ namespace corsika::process::pythia {
// set target
auto pdgTarget = static_cast<int>(particles::GetPDG(TargetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
- if(TargetId == particles::Code::Hydrogen)
- pdgTarget = static_cast<int>( particles::GetPDG(particles::Code::Proton));
+ if (TargetId == particles::Code::Hydrogen)
+ pdgTarget = static_cast<int>(particles::GetPDG(particles::Code::Proton));
std::stringstream stTarget;
stTarget << "Beams:idB = " << pdgTarget;
fPythia.readString(stTarget.str());
@@ -128,14 +123,15 @@ namespace corsika::process::pythia {
fPythia.init();
}
-
- bool Interaction::CanInteract(const corsika::particles::Code pCode)
- {
- return pCode == corsika::particles::Code::Proton || pCode == corsika::particles::Code::Neutron
- || pCode == corsika::particles::Code::AntiProton || pCode == corsika::particles::Code::AntiNeutron
- || pCode == corsika::particles::Code::PiMinus || pCode == corsika::particles::Code::PiPlus;
- }
-
+ bool Interaction::CanInteract(const corsika::particles::Code pCode) {
+ return pCode == corsika::particles::Code::Proton ||
+ pCode == corsika::particles::Code::Neutron ||
+ pCode == corsika::particles::Code::AntiProton ||
+ pCode == corsika::particles::Code::AntiNeutron ||
+ pCode == corsika::particles::Code::PiMinus ||
+ pCode == corsika::particles::Code::PiPlus;
+ }
+
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
@@ -143,27 +139,27 @@ namespace corsika::process::pythia {
using namespace units::si;
// interaction possible in pythia?
- if( TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen ){
- if( CanInteract(BeamId) && ValidCoMEnergy(CoMenergy) ){
- // input particle PDG
- auto const pdgCodeBeam = static_cast<int>( particles::GetPDG( BeamId ));
- auto const pdgCodeTarget = static_cast<int>( particles::GetPDG( TargetId ));
- const double ecm = CoMenergy / 1_GeV;
-
- // calculate cross section
- fSigma.calc( pdgCodeBeam, pdgCodeTarget, ecm);
- if(fSigma.hasSigmaTot()){
- const double sigEla = fSigma.sigmaEl();
- const double sigProd = fSigma.sigmaTot() - sigEla;
-
- return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
-
- } else
- throw std::runtime_error("pythia cross section init failed");
+ if (TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen) {
+ if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
+ // input particle PDG
+ auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(BeamId));
+ auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(TargetId));
+ const double ecm = CoMenergy / 1_GeV;
+
+ // calculate cross section
+ fSigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
+ if (fSigma.hasSigmaTot()) {
+ const double sigEla = fSigma.sigmaEl();
+ const double sigProd = fSigma.sigmaTot() - sigEla;
+
+ return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
+
+ } else
+ throw std::runtime_error("pythia cross section init failed");
} else {
- return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
- std::numeric_limits<double>::infinity() * 1_mb);
+ return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
+ std::numeric_limits<double>::infinity() * 1_mb);
}
} else {
throw std::runtime_error("invalid target for pythia");
@@ -206,7 +202,7 @@ namespace corsika::process::pythia {
<< " beam pid:" << p.GetPID() << endl;
// TODO: move limits into variables
- if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM) ) {
+ if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
// get target from environment
/*
@@ -214,8 +210,7 @@ namespace corsika::process::pythia {
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* currentNode = p.GetNode();
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
@@ -223,31 +218,29 @@ namespace corsika::process::pythia {
int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium
- const auto w = mediumComposition.GetFractions();
+ 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] =
+ [[maybe_unused]] auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
- [[maybe_unused]] auto elaCrossSectionCopy =
- elaCrossSection; // ONLY TO AVOID COMPILER WARNING
cout << "Interaction: IntLength: pythia return (mb): "
- << productionCrossSection / 1_mb << endl
- << "Interaction: IntLength: weight : " << w[i] << endl;
-
+ << productionCrossSection / 1_mb << endl
+ << "Interaction: IntLength: weight : " << w[i] << endl;
+
weightedProdCrossSection += w[i] * productionCrossSection;
}
cout << "Interaction: IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mb << endl
- << "Interaction: IntLength: average mass number: "
- << mediumComposition.GetAverageMassNumber() << endl;
+ << "Interaction: IntLength: average mass number: "
+ << mediumComposition.GetAverageMassNumber() << endl;
// calculate interaction length in medium
- GrammageType const int_length =
- mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection;
+ GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+ units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl;
@@ -257,7 +250,7 @@ namespace corsika::process::pythia {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
-
+
/**
In this function PYTHIA is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
@@ -271,7 +264,6 @@ namespace corsika::process::pythia {
using namespace units::si;
using namespace geometry;
-
const auto corsikaBeamId = p.GetPID();
cout << "Pythia::Interaction: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
@@ -283,7 +275,7 @@ namespace corsika::process::pythia {
}
if (process::pythia::Interaction::CanInteract(corsikaBeamId)) {
-
+
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -339,7 +331,7 @@ namespace corsika::process::pythia {
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ const auto* currentNode = p.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -353,66 +345,68 @@ namespace corsika::process::pythia {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, sigEla] =
+ [[maybe_unused]] const auto [sigProd, sigEla] =
GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
- [[maybe_unused]] auto sigElaCopy =
- sigEla; // to avoid not used warning in array binding
}
const auto corsikaTargetId =
mediumComposition.SampleTarget(cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << corsikaTargetId << endl;
-
- if (corsikaTargetId != particles::Code::Hydrogen && corsikaTargetId != particles::Code::Neutron
- && corsikaTargetId != particles::Code::Proton )
- throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
+
+ if (corsikaTargetId != particles::Code::Hydrogen &&
+ corsikaTargetId != particles::Code::Neutron &&
+ corsikaTargetId != particles::Code::Proton)
+ throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
-
+
if (eProjectileLab < 8.5_GeV || !ValidCoMEnergy(Ecm)) {
cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for PYTHIA");
-
- } else {
+
+ } else {
fCount++;
- ConfigureLabFrameCollision( corsikaBeamId, corsikaTargetId, eProjectileLab );
-
- // create event in pytia
- if(!fPythia.next())
- throw std::runtime_error("Pythia::DoInteraction: failed!");
+ ConfigureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
+
+ // create event in pytia
+ if (!fPythia.next()) throw std::runtime_error("Pythia::DoInteraction: failed!");
- // link to pythia stack
- Pythia8::Event& event = fPythia.event;
+ // link to pythia stack
+ Pythia8::Event& event = fPythia.event;
// print final state
- event.list();
+ event.list();
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
- for (int i = 0; i < event.size(); ++i){
- Pythia8::Particle &pp = event[i];
+ for (int i = 0; i < event.size(); ++i) {
+ Pythia8::Particle& pp = event[i];
// skip particles that have decayed in pythia
- if (!pp.isFinal()) continue;
+ if (!pp.isFinal()) continue;
- auto const pyId = particles::ConvertFromPDG(static_cast<particles::PDGCode>(pp.id()));
+ auto const pyId =
+ particles::ConvertFromPDG(static_cast<particles::PDGCode>(pp.id()));
- const MomentumVector pyPlab(rootCS, {pp.px()*1_GeV, pp.py()*1_GeV, pp.pz()*1_GeV});
+ const MomentumVector pyPlab(
+ rootCS, {pp.px() * 1_GeV, pp.py() * 1_GeV, pp.pz() * 1_GeV});
HEPEnergyType const pyEn = pp.e() * 1_GeV;
// add to corsika stack
auto pnew = p.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig, tOrig});
+ geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
+ tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
}
- cout << "conservation (all GeV): " << "Elab_final=" << Elab_final / 1_GeV
+ cout << "conservation (all GeV): "
+ << "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
// delete current particle
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
index f384da655322cac969588308c85094358c11cd80..3e1d8813b22a55ffef913ea8298414a8d441ec7d 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -19,10 +19,6 @@
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
-namespace corsika::environment {
- class Environment;
-}
-
namespace corsika::process::pythia {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
@@ -31,14 +27,14 @@ namespace corsika::process::pythia {
bool fInitialized = false;
public:
- Interaction(corsika::environment::Environment const& env);
+ Interaction() {}
~Interaction();
void Init();
- void SetParticleListStable(const std::vector<particles::Code>);
- void SetUnstable(const corsika::particles::Code );
- void SetStable(const corsika::particles::Code );
+ void SetParticleListStable(std::vector<particles::Code> const&);
+ void SetUnstable(const corsika::particles::Code);
+ void SetStable(const corsika::particles::Code);
bool WasInitialized() { return fInitialized; }
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
@@ -47,8 +43,9 @@ namespace corsika::process::pythia {
}
bool CanInteract(const corsika::particles::Code);
- void ConfigureLabFrameCollision(const corsika::particles::Code, const corsika::particles::Code,
- const corsika::units::si::HEPEnergyType);
+ void ConfigureLabFrameCollision(const corsika::particles::Code,
+ const corsika::particles::Code,
+ const corsika::units::si::HEPEnergyType);
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
@@ -66,9 +63,8 @@ namespace corsika::process::pythia {
corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
private:
- corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
+ corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
Pythia8::Pythia fPythia;
Pythia8::SigmaTotal fSigma;
const bool fInternalDecays = true;
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 418ff6973b6c51cc0c8a9cc6fa96c06e9c8988a3..1e6c57cca49e7fb53017776a8d4fb69790c9b278 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -9,10 +9,9 @@
* the license.
*/
-#include <corsika/process/pythia/Decay.h>
-#include <corsika/process/pythia/Interaction.h>
#include <Pythia8/Pythia.h>
#include <corsika/process/pythia/Decay.h>
+#include <corsika/process/pythia/Interaction.h>
#include <corsika/random/RNGManager.h>
@@ -96,26 +95,29 @@ TEST_CASE("Pythia", "[processes]") {
using namespace corsika;
using namespace corsika::units::si;
-TEST_CASE("pythia process"){
+TEST_CASE("pythia process") {
// setup environment, geometry
- environment::Environment env;
+ environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
+
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{cs, 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Hydrogen}, std::vector<float>{1.}));
+ std::vector<particles::Code>{particles::Code::Hydrogen},
+ std::vector<float>{1.}));
+ auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
universe.AddChild(std::move(theMedium));
- const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
-
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
1_m / second);
@@ -152,7 +154,7 @@ TEST_CASE("pythia process"){
}
SECTION("pythia interaction") {
-
+
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
@@ -163,15 +165,13 @@ TEST_CASE("pythia process"){
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::PiPlus, E0, plab, pos, 0_ns});
+ particle.SetNode(nodePtr);
+ process::pythia::Interaction model;
-
- process::pythia::Interaction model(env);
-
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoInteraction(particle,
- stack);
+ stack);
[[maybe_unused]] const GrammageType length =
- model.GetInteractionLength(particle, track);
+ model.GetInteractionLength(particle, track);
}
-
}
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index d54b974a5843857c9a3f580a9d40196a257c5c31..333b6c9894a9f1e7e84a77fecaa283afe9afe047 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -27,15 +27,12 @@ using std::cout;
using std::endl;
using std::tuple;
-using namespace corsika;
-using namespace corsika::setup;
-using Particle = Stack::StackIterator; // ParticleType;
-using Track = Trajectory;
+using SetupParticle = corsika::setup::Stack::StackIterator; // SetupParticleType;
+using Track = corsika::setup::Trajectory;
namespace corsika::process::sibyll {
- Interaction::Interaction(environment::Environment const& env)
- : fEnvironment(env) {}
+ Interaction::Interaction() {}
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
@@ -53,7 +50,7 @@ namespace corsika::process::sibyll {
if (fInternalDecays) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
- const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
+ const std::vector<particles::Code> hadronsWeWantTrackedByCorsika = {
particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Pi0, particles::Code::KMinus,
particles::Code::KPlus, particles::Code::K0Long,
@@ -64,7 +61,7 @@ namespace corsika::process::sibyll {
particles::Code::DMinus, particles::Code::D0,
particles::Code::D0Bar};
- Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
+ Interaction::SetParticleListStable(hadronsWeWantTrackedByCorsika);
}
fInitialized = true;
@@ -72,7 +69,7 @@ namespace corsika::process::sibyll {
}
void Interaction::SetParticleListStable(
- const std::vector<particles::Code> particleList) {
+ std::vector<particles::Code> const& particleList) {
for (auto p : particleList) Interaction::SetStable(p);
}
@@ -91,7 +88,7 @@ namespace corsika::process::sibyll {
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
- const units::si::HEPEnergyType CoMenergy) {
+ const units::si::HEPEnergyType CoMenergy) const {
using namespace units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
@@ -118,7 +115,8 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
+ units::si::GrammageType Interaction::GetInteractionLength(SetupParticle& p,
+ Track&) const {
using namespace units;
using namespace units::si;
@@ -162,8 +160,7 @@ namespace corsika::process::sibyll {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
- const auto currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ auto const* currentNode = p.GetNode();
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
@@ -209,7 +206,7 @@ namespace corsika::process::sibyll {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(Particle& p, Stack&) {
+ process::EProcessReturn Interaction::DoInteraction(SetupParticle& p, setup::Stack&) {
using namespace units;
using namespace utl;
@@ -283,8 +280,8 @@ namespace corsika::process::sibyll {
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
- const auto& mediumComposition =
+ auto const* currentNode = p.GetNode();
+ auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
/*
@@ -297,10 +294,9 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ [[maybe_unsused]] const auto [sigProd, sigEla] =
+ GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
- [[maybe_unused]] auto sigElaCopy =
- sigEla; // to avoid not used warning in array binding
}
const auto targetCode =
@@ -348,9 +344,6 @@ namespace corsika::process::sibyll {
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;
@@ -384,6 +377,8 @@ namespace corsika::process::sibyll {
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
}
+ // delete current particle
+ p.Delete();
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 21d4675b98f6d1733646aaf7184ff615243556ba..d2274d428e7819bb929b5a199e20ae2743da5a02 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -18,10 +18,6 @@
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
-namespace corsika::environment {
- class Environment;
-}
-
namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
@@ -31,23 +27,23 @@ namespace corsika::process::sibyll {
bool fInitialized = false;
public:
- Interaction(corsika::environment::Environment const& env);
+ Interaction();
~Interaction();
void Init();
- void SetParticleListStable(const std::vector<particles::Code>);
+ void SetParticleListStable(std::vector<particles::Code> const&);
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
bool WasInitialized() { return fInitialized; }
- bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+ bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
return (fMinEnergyCoM <= ecm) && (ecm <= fMaxEnergyCoM);
}
- int GetMaxTargetMassNumber() { return fMaxTargetMassNumber; }
- corsika::units::si::HEPEnergyType GetMinEnergyCoM() { return fMinEnergyCoM; }
- corsika::units::si::HEPEnergyType GetMaxEnergyCoM() { return fMaxEnergyCoM; }
- bool IsValidTarget(corsika::particles::Code TargetId) {
+ int GetMaxTargetMassNumber() const { return fMaxTargetMassNumber; }
+ corsika::units::si::HEPEnergyType GetMinEnergyCoM() const { return fMinEnergyCoM; }
+ corsika::units::si::HEPEnergyType GetMaxEnergyCoM() const { return fMaxEnergyCoM; }
+ bool IsValidTarget(corsika::particles::Code TargetId) const {
return (corsika::particles::GetNucleusA(TargetId) < fMaxTargetMassNumber) &&
corsika::particles::IsNucleus(TargetId);
}
@@ -55,10 +51,10 @@ namespace corsika::process::sibyll {
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
- const corsika::units::si::HEPEnergyType CoMenergy);
+ const corsika::units::si::HEPEnergyType CoMenergy) const;
template <typename Particle, typename Track>
- corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const;
/**
In this function SIBYLL is called to produce one event. The
@@ -69,7 +65,6 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
private:
- corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index d0db3df4a2e603fc403909e4f725d270df34f4e1..52c585916e5b274f5700d2e8cde726a4c664117f 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -29,46 +29,49 @@ using std::endl;
using std::tuple;
using std::vector;
-using namespace corsika;
-using namespace corsika::setup;
-using Particle = Stack::StackIterator; // ParticleType;
-using Track = Trajectory;
+using corsika::setup::SetupEnvironment;
+using SetupStack = corsika::setup::Stack;
+using Particle = SetupStack::StackIterator; // ParticleType;
+using Track = corsika::setup::Trajectory;
namespace corsika::process::sibyll {
- NuclearInteraction::NuclearInteraction(environment::Environment const& env,
- process::sibyll::Interaction& hadint)
+ template <>
+ NuclearInteraction<SetupEnvironment>::NuclearInteraction(
+ process::sibyll::Interaction& hadint, SetupEnvironment const& env)
: fEnvironment(env)
, fHadronicInteraction(hadint) {}
- NuclearInteraction::~NuclearInteraction() {
+ template <>
+ NuclearInteraction<SetupEnvironment>::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
- fTargetComponentsIndex.clear();
}
- void NuclearInteraction::Init() {
-
- using random::RNGManager;
-
- // initialize hadronic interaction module
- // TODO: safe to run multiple initializations?
- if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
-
- // check compatibility of energy ranges, someone could try to use low-energy model..
- if (!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
- !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
- throw std::runtime_error(
- "NuclearInteraction: hadronic interaction model incompatible!");
-
- // initialize nuclib
- // TODO: make sure this does not overlap with sibyll
- nuc_nuc_ini_();
+ template <>
+ void NuclearInteraction<SetupEnvironment>::PrintCrossSectionTable(
+ corsika::particles::Code pCode) {
+ using namespace corsika::particles;
+ const int k = fTargetComponentsIndex.at(pCode);
+ Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
+ Code::Neon, Code::Argon, Code::Iron};
+ cout << "en/A ";
+ for (auto& j : pNuclei) cout << std::setw(9) << j;
+ cout << endl;
- // initialize cross sections
- InitializeNuclearCrossSections();
+ // loop over energy bins
+ for (int i = 0; i < GetNEnergyBins(); ++i) {
+ cout << " " << i << " ";
+ for (auto& n : pNuclei) {
+ auto const j = GetNucleusA(n);
+ cout << " " << std::setprecision(5) << std::setw(8)
+ << cnucsignuc_.sigma[j - 1][k][i];
+ }
+ cout << endl;
+ }
}
- void NuclearInteraction::InitializeNuclearCrossSections() {
+ template <>
+ void NuclearInteraction<SetupEnvironment>::InitializeNuclearCrossSections() {
using namespace corsika::particles;
using namespace units::si;
@@ -77,12 +80,10 @@ namespace corsika::process::sibyll {
auto const allElementsInUniverse = std::invoke([&]() {
std::set<particles::Code> allElementsInUniverse;
auto collectElements = [&](auto& vtn) {
- if (auto const mp = vtn.GetModelPropertiesPtr();
- mp != nullptr) { // do not query Universe it self, it has no ModelProperties
- auto const& comp = mp->GetNuclearComposition().GetComponents();
+ if (vtn.HasModelProperties()) {
+ auto const& comp =
+ vtn.GetModelProperties().GetNuclearComposition().GetComponents();
for (auto const c : comp) allElementsInUniverse.insert(c);
- // std::for_each(comp.cbegin(), comp.cend(),
- // [&](particles::Code c) { allElementsInUniverse.insert(c); });
}
};
universe.walk(collectElements);
@@ -116,10 +117,10 @@ namespace corsika::process::sibyll {
const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
- for (int j = 1; j < fMaxNucleusAProjectile; ++j) {
+ for (int j = 1; j < gMaxNucleusAProjectile; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
- sigma_mc_(jj, ib, dsig, dsigela, fNSample, sig_out, dsig_out, sigqe_out,
+ sigma_mc_(jj, ib, dsig, dsigela, gNSample, sig_out, dsig_out, sigqe_out,
dsigqe_out);
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
@@ -135,28 +136,28 @@ namespace corsika::process::sibyll {
}
}
- void NuclearInteraction::PrintCrossSectionTable(corsika::particles::Code pCode) {
- using namespace corsika::particles;
- const int k = fTargetComponentsIndex.at(pCode);
- Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
- Code::Neon, Code::Argon, Code::Iron};
- cout << "en/A ";
- for (auto& j : pNuclei) cout << std::setw(9) << j;
- cout << endl;
+ template <>
+ void NuclearInteraction<SetupEnvironment>::Init() {
+ // initialize hadronic interaction module
+ // TODO: safe to run multiple initializations?
+ if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
- // loop over energy bins
- for (int i = 0; i < GetNEnergyBins(); ++i) {
- cout << " " << i << " ";
- for (auto& n : pNuclei) {
- auto const j = GetNucleusA(n);
- cout << " " << std::setprecision(5) << std::setw(8)
- << cnucsignuc_.sigma[j - 1][k][i];
- }
- cout << endl;
- }
+ // check compatibility of energy ranges, someone could try to use low-energy model..
+ if (!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
+ !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
+ throw std::runtime_error(
+ "NuclearInteraction: hadronic interaction model incompatible!");
+
+ // initialize nuclib
+ // TODO: make sure this does not overlap with sibyll
+ nuc_nuc_ini_();
+
+ // initialize cross sections
+ InitializeNuclearCrossSections();
}
- units::si::CrossSectionType NuclearInteraction::ReadCrossSectionTable(
+ template <>
+ units::si::CrossSectionType NuclearInteraction<SetupEnvironment>::ReadCrossSectionTable(
const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
@@ -174,8 +175,10 @@ namespace corsika::process::sibyll {
// TODO: remove elastic cross section?
template <>
+ template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) {
+ NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle& p,
+ const particles::Code TargetId) {
using namespace units::si;
if (p.GetPID() != particles::Code::Nucleus)
throw std::runtime_error(
@@ -211,7 +214,9 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType NuclearInteraction::GetInteractionLength(Particle& p, Track&) {
+ template <>
+ units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
+ Particle& p, Track&) {
using namespace units;
using namespace units::si;
@@ -261,8 +266,8 @@ namespace corsika::process::sibyll {
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
- if (ElabNuc >= 8.5_GeV && ECoMNN >= fMinEnergyPerNucleonCoM &&
- ECoMNN < fMaxEnergyPerNucleonCoM) {
+ if (ElabNuc >= 8.5_GeV && ECoMNN >= gMinEnergyPerNucleonCoM &&
+ ECoMNN < gMaxEnergyPerNucleonCoM) {
// get target from environment
/*
@@ -270,24 +275,22 @@ namespace corsika::process::sibyll {
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 =
+ auto const* const currentNode = p.GetNode();
+ auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
// weighted sum
int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium
- const auto w = mediumComposition.GetFractions();
+ 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] =
+ [[maybe_unused]] auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(p, targetId);
- [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): " << productionCrossSection / 1_mb
@@ -295,14 +298,14 @@ namespace corsika::process::sibyll {
weightedProdCrossSection += w[i] * productionCrossSection;
}
cout << "NuclearInteraction: "
- << "IntLength: weighted CrossSection (mb): "
- << weightedProdCrossSection / 1_mb << endl;
+ << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
+ << 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
+ << "interaction length (g/cm2): " << int_length * (1_cm * 1_cm / (0.001_kg))
<< endl;
return int_length;
@@ -312,7 +315,9 @@ namespace corsika::process::sibyll {
}
template <>
- process::EProcessReturn NuclearInteraction::DoInteraction(Particle& p, Stack& s) {
+ template <>
+ process::EProcessReturn NuclearInteraction<SetupEnvironment>::DoInteraction(
+ Particle& p, SetupStack& s) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
@@ -426,7 +431,7 @@ namespace corsika::process::sibyll {
//
// proton stand-in for nucleon
const auto beamId = particles::Proton::GetCode();
- const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
+ auto const* const currentNode = p.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
cout << "get nucleon-nucleus cross sections for target materials.." << endl;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index c62249bbefb0193d603714cc960743e921ee1a5d..b04bba538ec9cf2472f8d102fe4ec3bce070578d 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -16,10 +16,6 @@
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
-namespace corsika::environment {
- class Environment;
-}
-
namespace corsika::process::sibyll {
class Interaction; // fwd-decl
@@ -28,16 +24,15 @@ namespace corsika::process::sibyll {
*
*
**/
-
+ template <class TEnvironment>
class NuclearInteraction
- : public corsika::process::InteractionProcess<NuclearInteraction> {
+ : public corsika::process::InteractionProcess<NuclearInteraction<TEnvironment>> {
int fCount = 0;
int fNucCount = 0;
public:
- NuclearInteraction(corsika::environment::Environment const& env,
- corsika::process::sibyll::Interaction& hadint);
+ NuclearInteraction(corsika::process::sibyll::Interaction&, TEnvironment const&);
~NuclearInteraction();
void Init();
void InitializeNuclearCrossSections();
@@ -45,14 +40,15 @@ namespace corsika::process::sibyll {
corsika::units::si::CrossSectionType ReadCrossSectionTable(
const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() {
- return fMinEnergyPerNucleonCoM;
+ return gMinEnergyPerNucleonCoM;
}
corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() {
- return fMaxEnergyPerNucleonCoM;
+ return gMaxEnergyPerNucleonCoM;
}
- int GetMaxNucleusAProjectile() { return fMaxNucleusAProjectile; }
- int GetMaxNFragments() { return fMaxNFragments; }
- int GetNEnergyBins() { return fNEnBins; }
+ int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile; }
+ int constexpr GetMaxNFragments() { return gMaxNFragments; }
+ int constexpr GetNEnergyBins() { return gNEnBins; }
+
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
@@ -64,20 +60,21 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s);
private:
- corsika::environment::Environment const& fEnvironment;
+ TEnvironment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
std::map<corsika::particles::Code, int> fTargetComponentsIndex;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- const int fNSample = 500; // number of samples in MC estimation of cross section
- const int fMaxNucleusAProjectile = 56;
- const int fNEnBins = 6;
- const int fMaxNFragments = 60;
+ static constexpr int gNSample =
+ 500; // number of samples in MC estimation of cross section
+ static constexpr int gMaxNucleusAProjectile = 56;
+ static constexpr int gNEnBins = 6;
+ static constexpr int gMaxNFragments = 60;
// energy limits defined by table used for cross section in signuc.f
// 10**1 GeV to 10**6 GeV
- const corsika::units::si::HEPEnergyType fMinEnergyPerNucleonCoM =
+ static constexpr corsika::units::si::HEPEnergyType gMinEnergyPerNucleonCoM =
10. * 1e9 * corsika::units::si::electronvolt;
- const corsika::units::si::HEPEnergyType fMaxEnergyPerNucleonCoM =
+ static constexpr corsika::units::si::HEPEnergyType gMaxEnergyPerNucleonCoM =
1.e6 * 1e9 * corsika::units::si::electronvolt;
};
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index d5b71ae6db7562f6a69c32562dcef4197b4c3129..6acb2b4938282453a94997f69f9cdc01ad99e475 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -79,12 +79,13 @@ using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
- environment::Environment env;
+ environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
@@ -92,6 +93,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+ auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
@@ -116,8 +118,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
+ particle.SetNode(nodePtr);
- Interaction model(env);
+ Interaction model;
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
@@ -140,9 +143,10 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
+ particle.SetNode(nodePtr);
- Interaction hmodel(env);
- NuclearInteraction model(env, hmodel);
+ Interaction hmodel;
+ NuclearInteraction model(hmodel, env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
diff --git a/Processes/StackInspector/CMakeLists.txt b/Processes/StackInspector/CMakeLists.txt
index 96b502f60981b891a6c45c94b18013c1f7910477..f00f44c2ab6b5513630f5344545f988b8e78dc3f 100644
--- a/Processes/StackInspector/CMakeLists.txt
+++ b/Processes/StackInspector/CMakeLists.txt
@@ -28,6 +28,7 @@ set_target_properties (
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessStackInspector
+ CORSIKAcascade
CORSIKAunits
CORSIKAgeometry
CORSIKAsetup
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index fb749c0b8c7e6430c7283f079055979d099ee454..a22742165bec0a6f4a9bfe3b06e33e25ca2e9c3c 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -16,6 +16,7 @@
#include <corsika/logging/Logger.h>
+#include <corsika/cascade/testCascade.h>
#include <corsika/setup/SetupTrajectory.h>
#include <iostream>
@@ -51,7 +52,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
- << " pos=" << pos;
+ << " pos=" << pos << " node = " << iterP.GetNode();
// if (iterP.GetPID()==Code::Nucleus)
// cout << " nuc_ref=" << iterP.GetNucleusRef();
cout << endl;
@@ -76,3 +77,4 @@ void StackInspector<Stack>::Init() {
#include <corsika/setup/SetupStack.h>
template class process::stack_inspector::StackInspector<setup::Stack>;
+template class process::stack_inspector::StackInspector<TestCascadeStack>;
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 9021bf2796ce50768ca8f7c9c13a7d971f193c9c..319030f91d5347d47af58cddbdea607684a16677 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -9,13 +9,12 @@
* the license.
*/
-#include <corsika/process/tracking_line/TrackingLine.h>
-
#include <corsika/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
#include <algorithm>
#include <iostream>
@@ -26,120 +25,29 @@ using namespace corsika;
namespace corsika::process::tracking_line {
- template <class Stack, class Trajectory>
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TrackingLine<Stack, Trajectory>::TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere) {
- using namespace corsika::units::si;
- auto const& cs = fEnvironment.GetCoordinateSystem();
- geometry::Point const origin(cs, 0_m, 0_m, 0_m);
-
- auto const r0 = (line.GetR0() - origin);
- auto const v0 = line.GetV0();
- auto const c0 = (sphere.GetCenter() - origin);
+ TimeOfIntersection(corsika::geometry::Line const& line,
+ geometry::Sphere const& sphere) {
+ auto const delta = line.GetR0() - sphere.GetCenter();
+ auto const v = line.GetV0();
+ auto const vSqNorm = v.squaredNorm();
+ auto const R = sphere.GetRadius();
- auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
- auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
- sphere.GetRadius() * sphere.GetRadius();
-
- auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
+ auto const vDotDelta = v.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
//~ std::cout << "discriminant: " << discriminant << std::endl;
//~ std::cout << "alpha: " << alpha << std::endl;
//~ std::cout << "beta: " << beta << std::endl;
if (discriminant.magnitude() > 0) {
- (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
- return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
- (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+ return std::make_pair((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
} else {
return {};
}
}
-
- template <class Stack, class Trajectory>
- TrackingLine<Stack, Trajectory>::TrackingLine(
- corsika::environment::Environment const& pEnv)
- : fEnvironment(pEnv) {}
-
- template <class Stack, class Trajectory>
- void TrackingLine<Stack, Trajectory>::Init() {}
-
- template <class Stack, class Trajectory>
- Trajectory TrackingLine<Stack, Trajectory>::GetTrack(Particle const& p) {
- using std::cout;
- using std::endl;
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- geometry::Vector<SpeedType::dimension_type> const velocity =
- p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
-
- auto const currentPosition = p.GetPosition();
- std::cout << "TrackingLine pid: " << p.GetPID() << " , E = " << p.GetEnergy() / 1_GeV
- << " GeV" << std::endl;
- std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates() << std::endl;
- std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
- std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
-
- // to do: include effect of magnetic field
- geometry::Line line(currentPosition, velocity);
-
- auto const* currentVolumeNode =
- fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
- auto const& children = currentVolumeNode->GetChildNodes();
- auto const& excluded = currentVolumeNode->GetExcludedNodes();
-
- std::vector<TimeType> intersectionTimes;
-
- auto addIfIntersects = [&](auto& vtn) {
- auto const& volume = vtn.GetVolume();
- auto const& sphere = dynamic_cast<geometry::Sphere const&>(
- volume); // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
-
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
- auto const [t1, t2] = *opt;
- if (t1.magnitude() >= 0)
- intersectionTimes.push_back(t1);
- else if (t2.magnitude() >= 0)
- intersectionTimes.push_back(t2);
- }
- };
-
- for (auto const& child : children) { addIfIntersects(*child); }
-
- for (auto const* child : excluded) { addIfIntersects(*child); }
-
- addIfIntersects(*currentVolumeNode);
-
- auto const minIter =
- std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
-
- TimeType min;
-
- if (minIter == intersectionTimes.cend()) {
- min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
- // to handle the numerics properly
- //~ throw std::runtime_error("no intersection with anything!");
- } else {
- min = *minIter;
- }
-
- std::cout << " t-intersect: " << min << std::endl;
-
- return Trajectory(line, min);
- }
-
} // namespace corsika::process::tracking_line
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-using namespace corsika::setup;
-template class corsika::process::tracking_line::TrackingLine<setup::Stack,
- setup::Trajectory>;
-
-#include "testTrackingLineStack.h"
-template class corsika::process::tracking_line::TrackingLine<DummyStack,
- setup::Trajectory>;
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index cfc3568079bf7ad6433f2f28aee416230d9d5866..e847a3fd7ebef23583c632cc5b71af474c883428 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,39 +12,125 @@
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
-#include <map> // for std::pair
#include <optional>
+#include <type_traits>
+#include <utility>
namespace corsika::environment {
+ template <typename IEnvironmentModel>
class Environment;
-}
-namespace corsika::geometry {
- class Line;
- class Sphere;
-} // namespace corsika::geometry
+ template <typename IEnvironmentModel>
+ class VolumeTreeNode;
+} // namespace corsika::environment
namespace corsika::process {
namespace tracking_line {
- template <typename Stack, typename Trajectory>
- class TrackingLine { //
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+ TimeOfIntersection(corsika::geometry::Line const& line,
+ geometry::Sphere const& sphere);
- using Particle = typename Stack::StackIterator;
-
- corsika::environment::Environment const& fEnvironment;
+ class TrackingLine {
public:
- std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere);
-
- TrackingLine(corsika::environment::Environment const& pEnv);
-
- void Init();
-
- Trajectory GetTrack(Particle const& p);
+ TrackingLine(){};
+
+ template <typename Particle> // was Stack previously, and argument was
+ // Stack::StackIterator
+ auto GetTrack(Particle const& p) {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> const velocity =
+ p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
+
+ auto const currentPosition = p.GetPosition();
+ std::cout << "TrackingLine pid: " << p.GetPID()
+ << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine pos: "
+ << currentPosition.GetCoordinates()
+ // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
+ << std::endl;
+ std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
+ << " GeV " << std::endl;
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
+ // to do: include effect of magnetic field
+ geometry::Line line(currentPosition, velocity);
+
+ auto const* currentLogicalVolumeNode = p.GetNode();
+ //~ auto const* currentNumericalVolumeNode =
+ //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+ auto const numericallyInside =
+ currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
+
+ std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
+
+ auto const& children = currentLogicalVolumeNode->GetChildNodes();
+ auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
+
+ std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
+
+ // for entering from outside
+ auto addIfIntersects = [&](auto const& vtn) {
+ auto const& volume = vtn.GetVolume();
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ volume); // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+
+ if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+ auto const [t1, t2] = *opt;
+ std::cout << "intersection times: " << t1 / 1_s << "; "
+ << t2 / 1_s
+ // << " " << vtn.GetModelProperties().GetName()
+ << std::endl;
+ if (t1.magnitude() > 0)
+ intersections.emplace_back(t1, &vtn);
+ else if (t2.magnitude() > 0)
+ std::cout << "inside other volume" << std::endl;
+ }
+ };
+
+ for (auto const& child : children) { addIfIntersects(*child); }
+ for (auto const* ex : excluded) { addIfIntersects(*ex); }
+
+ {
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ currentLogicalVolumeNode->GetVolume());
+ // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+ auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+ intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
+ }
+
+ auto const minIter = std::min_element(
+ intersections.cbegin(), intersections.cend(),
+ [](auto const& a, auto const& b) { return a.first < b.first; });
+
+ TimeType min;
+
+ if (minIter == intersections.cend()) {
+ min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+ // to handle the numerics properly
+ throw std::runtime_error("no intersection with anything!");
+ } else {
+ min = minIter->first;
+ }
+
+ std::cout << " t-intersect: "
+ << min
+ // << " " << minIter->second->GetModelProperties().GetName()
+ << std::endl;
+
+ return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
+ velocity.norm() * min, minIter->second);
+ }
};
} // namespace tracking_line
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 00ad3e8e52e2d8c62f89fb3a22d3e4a31d956ab5..fb213f9d8e3729015a8ad2f19c092e3d476bf9dc 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -36,13 +35,13 @@ using namespace std;
using namespace corsika::units::si;
TEST_CASE("TrackingLine") {
- corsika::environment::Environment env; // dummy environment
+ environment::Environment<environment::Empty> env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
- tracking_line::TrackingLine<DummyStack, setup::Trajectory> tracking(env);
+ tracking_line::TrackingLine tracking;
SECTION("intersection with sphere") {
- Point const origin(cs, {0_m, 0_m, 0_m});
+ Point const origin(cs, {0_m, 0_m, -5_m});
Point const center(cs, {0_m, 0_m, 10_m});
Sphere const sphere(center, 1_m);
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
@@ -52,38 +51,46 @@ TEST_CASE("TrackingLine") {
setup::Trajectory traj(line, 12345_s);
auto const opt =
- tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+ tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
auto [t1, t2] = opt.value();
- REQUIRE(t1 / 9_s == Approx(1));
- REQUIRE(t2 / 11_s == Approx(1));
+ REQUIRE(t1 / 14_s == Approx(1));
+ REQUIRE(t2 / 16_s == Approx(1));
auto const optNoIntersection =
- tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+ tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
}
SECTION("maximally possible propagation") {
auto& universe = *(env.GetUniverse());
- //~ std::cout << env.GetUniverse().get() << std::endl;
-
- DummyParticle p(1_GeV, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV),
- Point(cs, 0_m, 0_m, 0_m));
-
auto const radius = 20_m;
- auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+ auto const* theMediumPtr = theMedium.get();
universe.AddChild(std::move(theMedium));
+ TestTrackingLineStack stack;
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::MuPlus,
+ 1_GeV,
+ {cs, {0_GeV, 0_GeV, 1_GeV}},
+ {cs, {0_m, 0_m, 0_km}},
+ 0_ns});
+ auto p = stack.GetNextParticle();
+ p.SetNode(theMediumPtr);
+
Point const origin(cs, {0_m, 0_m, 0_m});
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
- auto const traj = tracking.GetTrack(p);
+ auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
.GetComponents(cs)
diff --git a/Processes/TrackingLine/testTrackingLineStack.cc b/Processes/TrackingLine/testTrackingLineStack.cc
deleted file mode 100644
index 4c21923884ff69a76ec9acdc7428a313b1559734..0000000000000000000000000000000000000000
--- a/Processes/TrackingLine/testTrackingLineStack.cc
+++ /dev/null
@@ -1,37 +0,0 @@
-
-/*
- * (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_
-
-typedef corsika::units::si::hepmomentum_d MOMENTUM;
-
-struct DummyParticle {
- HEPEnergyType fEnergy;
- Vector<MOMENTUM> fMomentum;
- Point fPosition;
-
- DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
- : fEnergy(pEnergy)
- , fMomentum(pMomentum)
- , fPosition(pPosition) {}
-
- auto GetEnergy() const { return fEnergy; }
- auto GetMomentum() const { return fMomentum; }
- auto GetPosition() const { return fPosition; }
- auto GetPID() const { return corsika::particles::Code::Unknown; }
-};
-
-struct DummyStack {
- using ParticleType = DummyParticle;
-};
-
-#endif
diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h
index b825bc2f17b684e5fd753557296015d748181e57..8bdf766268ac90a41c52e99aadd30c1679f7d96f 100644
--- a/Processes/TrackingLine/testTrackingLineStack.h
+++ b/Processes/TrackingLine/testTrackingLineStack.h
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -12,34 +11,26 @@
#ifndef _include_process_trackinling_teststack_h_
#define _include_process_trackinling_teststack_h_
+#include <corsika/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
-typedef corsika::units::si::hepmomentum_d MOMENTUM;
-
-struct DummyParticle {
- corsika::units::si::HEPEnergyType fEnergy;
- corsika::geometry::Vector<MOMENTUM> fMomentum;
- corsika::geometry::Point fPosition;
-
- DummyParticle(corsika::units::si::HEPEnergyType pEnergy,
- corsika::geometry::Vector<MOMENTUM> pMomentum,
- corsika::geometry::Point pPosition)
- : fEnergy(pEnergy)
- , fMomentum(pMomentum)
- , fPosition(pPosition) {}
+using TestEnvironmentType =
+ corsika::environment::Environment<corsika::environment::Empty>;
- auto GetEnergy() const { return fEnergy; }
- auto GetMomentum() const { return fMomentum; }
- auto GetPosition() const { return fPosition; }
- auto GetPID() const { return corsika::particles::Code::Unknown; }
-};
+template <typename T>
+using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
-struct DummyStack {
- using ParticleType = DummyParticle;
- using StackIterator = DummyParticle;
-};
+// combine particle data stack with geometry information for tracking
+template <typename StackIter>
+using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<
+ corsika::setup::detail::ParticleDataStack::PIType, SetupGeometryDataInterface,
+ StackIter>;
+using TestTrackingLineStack = corsika::stack::CombinedStack<
+ typename corsika::setup::detail::ParticleDataStack::StackImpl,
+ GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
#endif
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index 08780b73b38311412a11500adb9d37a67a6d4228..83397da951f2d7594fea3e2c5a3dfa1d2b88da19 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -12,10 +12,13 @@
#ifndef _include_corsika_setup_environment_h_
#define _include_corsika_setup_environment_h_
+#include <corsika/environment/Environment.h>
#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NameModel.h>
namespace corsika::setup {
- using IEnvironmentModel = corsika::environment::IMediumModel;
-}
+ using IEnvironmentModel = environment::IMediumModel;
+ using SetupEnvironment = environment::Environment<IEnvironmentModel>;
+} // namespace corsika::setup
#endif
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index 9f42d76c92e960a0a71fcd01854aa2e37fb85de3..dfcfa757dc2d0ad183a1ae9919bd603fbe8bb250 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -20,15 +20,20 @@
// extension with geometry information for tracking
#include <corsika/environment/Environment.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/stack/CombinedStack.h>
#include <tuple>
+#include <utility>
#include <vector>
// definition of stack-data object to store geometry information
+template <typename TEnvType>
class GeometryData {
public:
+ using BaseNodeType = typename TEnvType::BaseNodeType;
+
// these functions are needed for the Stack interface
void Init() {}
void Clear() { fNode.clear(); }
@@ -38,10 +43,8 @@ public:
void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); }
// custom data access function
- void SetNode(const int i, const corsika::environment::BaseNodeType* v) { fNode[i] = v; }
- const corsika::environment::BaseNodeType* GetNode(const int i) const {
- return fNode[i];
- }
+ void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; }
+ auto const* GetNode(const int i) const { return fNode[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fNode.push_back(nullptr); }
@@ -51,37 +54,34 @@ public:
// custom private data section
private:
- std::vector<const corsika::environment::BaseNodeType*> fNode;
+ std::vector<const BaseNodeType*> fNode;
};
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
-template <typename T>
+template <typename T, typename TEnvType>
class GeometryDataInterface : public T {
public:
using T::GetIndex;
using T::GetStackData;
using T::SetParticleData;
+ using BaseNodeType = typename TEnvType::BaseNodeType;
// default version for particle-creation from input data
- void SetParticleData(const std::tuple<const corsika::environment::BaseNodeType*> v) {
+ void SetParticleData(const std::tuple<BaseNodeType const*> v) {
SetNode(std::get<0>(v));
}
- void SetParticleData(GeometryDataInterface<T>& parent,
- const std::tuple<const corsika::environment::BaseNodeType*>) {
+ void SetParticleData(GeometryDataInterface& parent,
+ const std::tuple<BaseNodeType const*>) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
void SetParticleData() { SetNode(nullptr); }
- void SetParticleData(GeometryDataInterface<T>& parent) {
+ void SetParticleData(GeometryDataInterface& parent) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
- void SetNode(const corsika::environment::BaseNodeType* v) {
- GetStackData().SetNode(GetIndex(), v);
- }
- const corsika::environment::BaseNodeType* GetNode() const {
- return GetStackData().GetNode(GetIndex());
- }
+ void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); }
+ auto const* GetNode() const { return GetStackData().GetNode(GetIndex()); }
};
namespace corsika::setup {
@@ -101,14 +101,18 @@ namespace corsika::setup {
using ParticleDataStack = corsika::stack::nuclear_extension::NuclearStackExtension<
corsika::stack::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>;
+ template <typename T>
+ using SetupGeometryDataInterface = GeometryDataInterface<T, setup::SetupEnvironment>;
+
// combine particle data stack with geometry information for tracking
template <typename StackIter>
using StackWithGeometryInterface =
corsika::stack::CombinedParticleInterface<ParticleDataStack::PIType,
- GeometryDataInterface, StackIter>;
+ SetupGeometryDataInterface, StackIter>;
using StackWithGeometry =
- corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl, GeometryData,
+ corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl,
+ GeometryData<setup::SetupEnvironment>,
StackWithGeometryInterface>;
} // namespace detail
diff --git a/ThirdParty/pythia8235.tgz b/ThirdParty/pythia8235.tgz
new file mode 100644
index 0000000000000000000000000000000000000000..552398c41da544421f187f0a5cf78f3739a1ee1d
Binary files /dev/null and b/ThirdParty/pythia8235.tgz differ
diff --git a/Tools/CMakeLists.txt b/Tools/CMakeLists.txt
index c0e4c12b9088e5384462b527158d1315fb47322d..c1347108d1d451651bb8a6a593c4acd7c8aab66d 100644
--- a/Tools/CMakeLists.txt
+++ b/Tools/CMakeLists.txt
@@ -1,4 +1,4 @@
-set (TOOLS_FILES plot_tracks.sh)
+set (TOOLS_FILES plot_tracks.sh plot_crossings.sh)
install (
FILES ${TOOLS_FILES}
diff --git a/Tools/plot_crossings.sh b/Tools/plot_crossings.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d4a318e20ab7fe013ed89afae57b5a3c9220ea98
--- /dev/null
+++ b/Tools/plot_crossings.sh
@@ -0,0 +1,39 @@
+#!/bin/sh
+
+# (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.
+
+# with this script you can plot an animation of output of TrackWriter
+
+crossings_dat=$1
+if [ -z "$crossings_dat" ]; then
+ echo "usage: $0 <crossings.dat> [output.gif]" >&2
+ exit 1
+fi
+
+output=$2
+if [ -z "$output" ]; then
+ output="$crossings_dat.gif"
+fi
+
+cat <<EOF | gnuplot
+set term gif animate size 600,600
+set output "$output"
+
+set xlabel "x / m"
+set ylabel "y / m"
+set zlabel "z / m"
+set title "CORSIKA 8 preliminary"
+
+do for [t=0:360:1] {
+ set view 60, t
+ splot "$crossings_dat" u 2:3:4 w points t ""
+}
+EOF
+
+exit $?