Documentation/Examples/CMakeLists.txt

-
-CORSIKA_ADD_EXAMPLE (helix_example)
-target_link_libraries (helix_example CORSIKAgeometry CORSIKAunits)
-
-CORSIKA_ADD_EXAMPLE (particle_list_example)
-target_link_libraries (particle_list_example CORSIKAparticles CORSIKAunits CORSIKAprocesses ProcessSibyll ProcessQGSJetII)
-
-CORSIKA_ADD_EXAMPLE (geometry_example)
-target_link_libraries (geometry_example CORSIKAgeometry CORSIKAunits)
-
-CORSIKA_ADD_EXAMPLE (stack_example)
-target_link_libraries (stack_example SuperStupidStack CORSIKAunits)
-
-# address sanitizer is making this example too slow, so we only do "undefined"
-CORSIKA_ADD_EXAMPLE (boundary_example)
-target_link_libraries (boundary_example
-  SuperStupidStack
-  CORSIKAunits
-  CORSIKAlogging
-  CORSIKArandom
-  ProcessSibyll
-  CORSIKAcascade
-  ProcessTrackWriter
-  ProcessParticleCut
-  ProcessTrackingLine
-  ProcessPythia8
-  ProcessProposal
-  CORSIKAprocesses
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  )
-
-CORSIKA_ADD_EXAMPLE (cascade_example)
-target_link_libraries (cascade_example
-  SuperStupidStack
-  CORSIKAunits
-  CORSIKAlogging
-  CORSIKArandom
-  ProcessSibyll
-  ProcessProposal
-  CORSIKAcascade
-  ProcessEnergyLoss
-  ProcessTrackWriter
-  ProcessStackInspector
-  ProcessTrackingLine
-  ProcessParticleCut
-  ProcessHadronicElasticModel
-  ProcessStackInspector
-  CORSIKAprocesses
-  CORSIKAcascade
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  )
-
-if (Pythia8_FOUND)
-  CORSIKA_ADD_EXAMPLE (cascade_proton_example)
-  target_link_libraries (cascade_proton_example
-    SuperStupidStack
-    CORSIKAunits
-    CORSIKAlogging
-    CORSIKArandom
-    ProcessSibyll
-    ProcessPythia8
-    ProcessUrQMD
-    ProcessSwitch
-    CORSIKAcascade
-    ProcessCONEXSourceCut
-    ProcessEnergyLoss
-    ProcessTrackWriter
-    ProcessStackInspector
-    ProcessTrackingLine
-    ProcessParticleCut
-    ProcessOnShellCheck
-    ProcessHadronicElasticModel
-    ProcessStackInspector
-    CORSIKAprocesses
-    CORSIKAcascade
-    CORSIKAparticles
-    CORSIKAgeometry
-    CORSIKAenvironment
-    CORSIKAprocesssequence
-    )
-
-  CORSIKA_ADD_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
-  target_link_libraries (vertical_EAS
-    SuperStupidStack
-    CORSIKAunits
-    CORSIKAlogging
-    CORSIKArandom
-    CORSIKAcascade
-    ProcessProposal
-    ProcessPythia8
-    ProcessObservationPlane
-    ProcessInteractionCounter
-    ProcessTrackWriter
-    ProcessEnergyLoss
-    ProcessTrackingLine
-    ProcessParticleCut
-    ProcessOnShellCheck
-    ProcessStackInspector
-    ProcessLongitudinalProfile
-    ProcessCONEXSourceCut
-    CORSIKAprocesses
-    CORSIKAcascade
-    CORSIKAparticles
-    CORSIKAgeometry
-    CORSIKAenvironment
-    CORSIKAprocesssequence
-    )
-endif()
-
-CORSIKA_ADD_EXAMPLE (stopping_power stopping_power)
-target_link_libraries (stopping_power
-  SuperStupidStack
-  CORSIKAunits
-  ProcessEnergyLoss
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  )
-
-CORSIKA_ADD_EXAMPLE (staticsequence_example)
-target_link_libraries (staticsequence_example
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAgeometry
-  CORSIKAlogging)
-
-
-CORSIKA_ADD_EXAMPLE (em_shower RUN_OPTIONS 100.)
-target_link_libraries (em_shower
-  SuperStupidStack
-  CORSIKAunits
-  CORSIKAlogging
-  CORSIKArandom
-  ProcessSibyll
-  ProcessPythia8
-  ProcessUrQMD
-  ProcessSwitch
-  CORSIKAcascade
-  ProcessCONEXSourceCut
-  ProcessEnergyLoss
-  ProcessObservationPlane
-  ProcessInteractionCounter
-  ProcessTrackWriter
-  ProcessProposal
-  ProcessTrackingLine
-  ProcessParticleCut
-  ProcessOnShellCheck
-  ProcessStackInspector
-  ProcessLongitudinalProfile
-  CORSIKAprocesses
-  CORSIKAcascade
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  )

Documentation/Examples/boundary_example.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/process/particle_cut/ParticleCut.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;
-
-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;
-  }
-
-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<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());
-
-  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::Proton},
-                  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("sibyll");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::Decay decay;
-
-  process::particle_cut::ParticleCut cut(20_GeV, true, true);
-
-  process::track_writer::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.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;
-}

Documentation/Examples/cascade_example.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/energy_loss/EnergyLoss.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>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/ShowerAxis.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/particle_cut/ParticleCut.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>
-
-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;
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-
-  const LengthType height_atmosphere = 112.8_km;
-
-  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
-  const float fox = 0.20946;
-  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}));
-
-  auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
-
-  innerMedium->SetModelProperties(props);
-
-  outerMedium->AddChild(std::move(innerMedium));
-
-  universe.AddChild(std::move(outerMedium));
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  const Code beamCode = Code::Nucleus;
-  const int nuclA = 4;
-  const int nuclZ = int(nuclA / 2.15 + 0.7);
-  const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
-  const HEPEnergyType E0 = nuclA * 1_TeV;
-  double theta = 0.;
-  double phi = 0.;
-
-  Point const injectionPos(
-      rootCS, 0_m, 0_m,
-      height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
-
-  ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -5000_km}, env};
-
-  {
-    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;
-    stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                 corsika::stack::MomentumVector, geometry::Point,
-                                 units::si::TimeType, unsigned short, unsigned short>{
-        beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ});
-  }
-
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
-
-  random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
-  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay;
-  // cascade with only HE model ==> HE cut
-  process::particle_cut::ParticleCut cut(80_GeV, true, true);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-  process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
-
-  // assemble all processes into an ordered process list
-  auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
-                               << trackWriter;
-
-  // define air shower object, run simulation
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-
-  EAS.Run();
-
-  eLoss.PrintProfile(); // print longitudinal profile
-
-  cut.ShowResults();
-  const HEPEnergyType Efinal =
-      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
-  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-  cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
-       << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
-  cut.Reset();  
-}

Documentation/Examples/cascade_proton_example.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/cascade/Cascade.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#include <corsika/geometry/Sphere.h>
-
-#include <corsika/process/sibyll/Decay.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-
-#include <corsika/process/pythia/Decay.h>
-#include <corsika/process/pythia/Interaction.h>
-
-#include <corsika/process/track_writer/TrackWriter.h>
-
-#include <corsika/process/particle_cut/ParticleCut.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;
-
-//
-// 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<setup::IEnvironmentModel>;
-  EnvType env;
-  auto& universe = *(env.GetUniverse());
-
-  auto theMedium =
-      EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
-                                  1_km * std::numeric_limits<double>::infinity());
-
-  using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_kg / (1_m * 1_m * 1_m),
-      NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
-                         std::vector<float>{(float)1.}));
-
-  universe.AddChild(std::move(theMedium));
-
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  const Code beamCode = Code::Proton;
-  const HEPMassType mass = particles::Proton::GetMass();
-  const HEPEnergyType E0 = 100_GeV;
-  double theta = 0.;
-  double phi = 0.;
-
-  {
-    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-      return sqrt(Elab * Elab - m * m);
-    };
-    HEPMomentumType P0 = elab2plab(E0, mass);
-    auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
-      return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-                             -ptot * cos(theta));
-    };
-    auto const [px, py, pz] =
-        momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
-    auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
-    cout << "input particle: " << beamCode << endl;
-    cout << "input angles: theta=" << theta << " phi=" << phi << endl;
-    cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
-    Point pos(rootCS, 0_m, 0_m, 0_m);
-    stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            beamCode, E0, plab, pos, 0_ns});
-  }
-
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
-
-  random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
-  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  //  process::sibyll::Interaction sibyll(env);
-  process::pythia::Interaction pythia;
-  //  process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
-  //  process::sibyll::Decay decay;
-  process::pythia::Decay decay;
-  process::particle_cut::ParticleCut cut(20_GeV, true, true);
-
-  // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
-  // process::HadronicElasticModel::HadronicElasticInteraction
-  // hadronicElastic(env);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-
-  // assemble all processes into an ordered process list
-  // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
-  auto sequence = pythia << decay << cut << trackWriter << stackInspect;
-
-  // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
-  // << "\n";
-
-  // define air shower object, run simulation
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-
-  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;
-}

Documentation/Examples/em_shower.cc

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/environment/Environment.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/ShowerAxis.h>
-#include <corsika/geometry/Plane.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/StackProcess.h>
-#include <corsika/process/interaction_counter/InteractionCounter.h>
-#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
-#include <corsika/process/observation_plane/ObservationPlane.h>
-#include <corsika/process/particle_cut/ParticleCut.h>
-#include <corsika/process/proposal/ContinuousProcess.h>
-#include <corsika/process/proposal/Interaction.h>
-#include <corsika/process/track_writer/TrackWriter.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <iomanip>
-#include <iostream>
-#include <limits>
-#include <string>
-#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;
-
-void registerRandomStreams() {
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-  random::RNGManager::GetInstance().RegisterRandomStream("proposal");
-  random::RNGManager::GetInstance().SeedAll();
-}
-
-int main(int argc, char** argv) {
-  if (argc != 2) {
-    std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
-    return 1;
-  }
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  registerRandomStreams();
-
-  // setup environment, geometry
-  using EnvType = Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-  Point const center{rootCS, 0_m, 0_m, 0_m};
-  environment::LayeredSphericalAtmosphereBuilder builder{center};
-  builder.setNuclearComposition(
-      {{particles::Code::Nitrogen, particles::Code::Oxygen},
-       {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
-
-  builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
-  builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
-  builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
-  builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
-  builder.addLinearLayer(1e9_cm, 112.8_km);
-
-  builder.assemble(env);
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  const Code beamCode = Code::Electron;
-  auto const mass = particles::GetMass(beamCode);
-  const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
-  double theta = 0.;
-  auto const thetaRad = theta / 180. * M_PI;
-
-  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-    return sqrt((Elab - m) * (Elab + m));
-  };
-  HEPMomentumType P0 = elab2plab(E0, mass);
-  auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
-    return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
-  };
-
-  auto const [px, py, pz] = momentumComponents(thetaRad, P0);
-  auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
-  cout << "input particle: " << beamCode << endl;
-  cout << "input angles: theta=" << theta << endl;
-  cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
-       << endl;
-
-  auto const observationHeight = 1.4_km + builder.getEarthRadius();
-  auto const injectionHeight = 112.75_km + builder.getEarthRadius();
-  auto const t = -observationHeight * cos(thetaRad) +
-                 sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
-                      si::detail::static_pow<2>(injectionHeight));
-  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
-  Point const injectionPos =
-      showerCore +
-      Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
-
-  std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
-
-  stack.AddParticle(
-      std::tuple<particles::Code, units::si::HEPEnergyType,
-                 corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-          beamCode, E0, plab, injectionPos, 0_ns});
-
-  std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
-            << std::endl;
-
-  environment::ShowerAxis const showerAxis{injectionPos,
-                                           (showerCore - injectionPos) * 1.02, env};
-
-  // setup processes, decays and interactions
-
-  // PROPOSAL processs proposal{...};
-  process::particle_cut::ParticleCut cut(10_GeV, false, true);
-  process::proposal::Interaction proposal(env, cut.GetECut());
-  process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
-  process::interaction_counter::InteractionCounter proposalCounted(proposal);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-
-  // long. profile; columns for gamma, e+, e- still need to be added
-  process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
-
-  Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-  process::observation_plane::ObservationPlane observationLevel(obsPlane,
-                                                                "particles.dat");
-
-  auto sequence = proposalCounted << em_continuous << longprof << cut << observationLevel
-                                  << trackWriter;
-  // define air shower object, run simulation
-  tracking_line::TrackingLine tracking;
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-
-  // to fix the point of first interaction, uncomment the following two lines:
-  //  EAS.SetNodes();
-  //  EAS.forceInteraction();
-
-  EAS.Run();
-
-  cut.ShowResults();
-  em_continuous.ShowResults();
-  observationLevel.ShowResults();
-  cout << "Cascade energy cut: " << EAS.GetEnergyCut()/1_GeV << " GeV" << endl;
-  const HEPEnergyType Efinal =
-    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround() + EAS.GetEnergyCut();
-  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-  observationLevel.Reset();
-  cut.Reset();
-  em_continuous.Reset();
-  
-  auto const hists = proposalCounted.GetHistogram();
-  hists.saveLab("inthist_lab.txt");
-  hists.saveCMS("inthist_cms.txt");
-
-  longprof.save("longprof.txt");
-
-  std::ofstream finish("finished");
-  finish << "run completed without error" << std::endl;
-}

Documentation/Examples/particle_list_example.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/particles/ParticleProperties.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/qgsjetII/ParticleConversion.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <iomanip>
-#include <iostream>
-#include <string>
-using namespace corsika::units;
-using namespace corsika::units::si;
-using namespace corsika::particles;
-using namespace std;
-
-//
-// The example main program for a particle list
-//
-int main() {
-  cout << "------------------------------------------"
-       << "particles in CORSIKA"
-       << "------------------------------------------" << endl;
-  cout << std::setw(20) << "Name"
-       << " | " << std::setw(10) << "PDG-id"
-       << " | " << std::setw(10) << "SIBYLL-id"
-       << " | " << std::setw(10) << "QGSJETII-id"
-       << " | " << std::setw(18) << "PDG-mass (GeV)"
-       << " | " << std::setw(18) << "SIBYLL-mass (GeV)"
-       << " | " << endl;
-  cout << std::setw(104) << std::setfill('-') << "-" << endl;
-  for (auto p : getAllParticles()) {
-    if (!IsNucleus(p)) {
-      corsika::process::sibyll::SibyllCode sib_id =
-          corsika::process::sibyll::ConvertToSibyll(p);
-      auto const sib_mass =
-          (sib_id != corsika::process::sibyll::SibyllCode::Unknown
-               ? to_string(corsika::process::sibyll::GetSibyllMass(p) / 1_GeV)
-               : "--");
-      auto const qgs_id = corsika::process::qgsjetII::ConvertToQgsjetII(p);
-      cout << std::setw(20) << std::setfill(' ') << p << " | " << std::setw(10)
-           << static_cast<int>(GetPDG(p)) << " | " << std::setw(10)
-           << (sib_id != corsika::process::sibyll::SibyllCode::Unknown
-                   ? to_string(static_cast<int>(sib_id))
-                   : "--")
-           << " | " << std::setw(10)
-           << (qgs_id != corsika::process::qgsjetII::QgsjetIICode::Unknown
-                   ? to_string(static_cast<int>(qgs_id))
-                   : "--")
-           << " | " << std::setw(18) << std::setprecision(5) << GetMass(p) / 1_GeV
-           << " | " << std::setw(18) << std::setprecision(5) << sib_mass << " | " << endl;
-    }
-  }
-  cout << std::setw(104) << std::setfill('-') << "-" << endl;
-}

Documentation/Examples/stack_example.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/particles/ParticleProperties.h>
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-
-#include <cassert>
-#include <iomanip>
-#include <iostream>
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::stack;
-using namespace corsika::geometry;
-using namespace std;
-
-void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
-  const geometry::CoordinateSystem& rootCS =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  for (int i = 0; i < 11; ++i) {
-    s.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::Electron, 1.5_GeV * i,
-            corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
-            geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-  }
-}
-
-void read(corsika::stack::super_stupid::SuperStupidStack& s) {
-  assert(s.GetSize() == 11); // stack has 11 particles
-
-  HEPEnergyType total_energy;
-  int i = 0;
-  for (auto& p : s) {
-    total_energy += p.GetEnergy();
-    // particles are electrons with 1.5 GeV energy times i
-    assert(p.GetPID() == particles::Code::Electron);
-    assert(p.GetEnergy() == 1.5_GeV * (i++));
-  }
-}
-
-int main() {
-  corsika::stack::super_stupid::SuperStupidStack s;
-  fill(s);
-  read(s);
-  return 0;
-}

Documentation/Examples/staticsequence_example.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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 <array>
-#include <iomanip>
-#include <iostream>
-
-#include <corsika/process/ProcessSequence.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::process;
-using namespace std;
-
-const int nData = 10;
-
-class Process1 : public BaseProcess<Process1> {
-public:
-  Process1() {}
-  template <typename D, typename T, typename S>
-  EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] += 1;
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process2 : public BaseProcess<Process2> {
-public:
-  Process2() {}
-
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process3 : public BaseProcess<Process3> {
-public:
-  Process3() {}
-
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D&, T&, S&) const {
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process4 : public BaseProcess<Process4> {
-public:
-  Process4(const double v)
-      : fV(v) {}
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] *= fV;
-    return EProcessReturn::eOk;
-  }
-
-private:
-  double fV;
-};
-
-struct DummyData {
-  double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-};
-struct DummyStack {
-  void clear() {}
-};
-struct DummyTrajectory {};
-
-void modular() {
-
-  Process1 m1;
-  Process2 m2;
-  Process3 m3;
-  Process4 m4(0.9);
-
-  auto sequence = m1 << m2 << m3 << m4;
-
-  DummyData particle;
-  DummyTrajectory track;
-
-  const int n = 1000;
-  for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
-
-  for (int i = 0; i < nData; ++i) {
-    // cout << p.p[i] << endl;
-    // assert(p.p[i] == n-i*100);
-  }
-
-  cout << " done (nothing...) " << endl;
-}
-
-int main() {
-  modular();
-  return 0;
-}

Documentation/Examples/stopping_power.cc

-/*
- * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/ShowerAxis.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <fstream>
-#include <iostream>
-#include <limits>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::particles;
-using namespace corsika::geometry;
-using namespace corsika::environment;
-
-using namespace std;
-using namespace corsika::units::si;
-
-//
-// This example demonstrates the energy loss of muons as function of beta*gamma (=p/m)
-//
-int main() {
-  feenableexcept(FE_INVALID);
-
-  // setup environment, geometry
-  using EnvType = Environment<IMediumModel>;
-  EnvType env;
-  env.GetUniverse()->SetModelProperties<HomogeneousMedium<IMediumModel>>(
-      1_g / cube(1_cm), NuclearComposition{{particles::Code::Unknown}, {1.f}});
-
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  Point const injectionPos(
-      rootCS, 0_m, 0_m,
-      112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
-
-  environment::ShowerAxis showerAxis{injectionPos,
-                                     Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
-  process::energy_loss::EnergyLoss eLoss{showerAxis, 300_MeV};
-
-  setup::Stack stack;
-
-  std::ofstream file("dEdX.dat");
-  file << "# beta*gamma, dE/dX / eV/(g/cm²)" << std::endl;
-
-  for (HEPEnergyType E0 = 300_MeV; E0 < 1_PeV; E0 *= 1.05) {
-    stack.Clear();
-    const Code beamCode = Code::MuPlus;
-    const HEPMassType mass = GetMass(beamCode);
-    double theta = 0.;
-    double phi = 0.;
-
-    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;
-
-    stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            beamCode, E0, plab, injectionPos, 0_ns});
-
-    auto const p = stack.GetNextParticle();
-    HEPEnergyType dE = eLoss.TotalEnergyLoss(p, 1_g / square(1_cm));
-    file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
-  }
-}

Documentation/Examples/vertical_EAS.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/environment/Environment.h>
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/ShowerAxis.h>
-#include <corsika/geometry/Plane.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/StackProcess.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/process/interaction_counter/InteractionCounter.h>
-#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
-#include <corsika/process/observation_plane/ObservationPlane.h>
-#include <corsika/process/on_shell_check/OnShellCheck.h>
-#include <corsika/process/particle_cut/ParticleCut.h>
-#include <corsika/process/pythia/Decay.h>
-#include <corsika/process/sibyll/Decay.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/proposal/ContinuousProcess.h>
-#include <corsika/process/proposal/Interaction.h>
-#include <corsika/process/switch_process/SwitchProcess.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-#include <corsika/process/urqmd/UrQMD.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <iomanip>
-#include <iostream>
-#include <limits>
-#include <string>
-#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;
-
-void registerRandomStreams() {
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-  random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
-  random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
-  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  random::RNGManager::GetInstance().RegisterRandomStream("urqmd");
-  random::RNGManager::GetInstance().RegisterRandomStream("proposal");
-
-  random::RNGManager::GetInstance().SeedAll();
-}
-
-int main(int argc, char** argv) {
-  if (argc != 4) {
-    std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV>" << std::endl;
-    return 1;
-  }
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  registerRandomStreams();
-
-  // setup environment, geometry
-  using EnvType = Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-  Point const center{rootCS, 0_m, 0_m, 0_m};
-  environment::LayeredSphericalAtmosphereBuilder builder{center};
-  builder.setNuclearComposition(
-      {{particles::Code::Nitrogen, particles::Code::Oxygen},
-       {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
-
-  builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
-  builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
-  builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
-  builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
-  builder.addLinearLayer(1e9_cm, 112.8_km);
-
-  builder.assemble(env);
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  const Code beamCode = Code::Nucleus;
-  unsigned short const A = std::stoi(std::string(argv[1]));
-  unsigned short Z = std::stoi(std::string(argv[2]));
-  auto const mass = particles::GetNucleusMass(A, Z);
-  const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
-  double theta = 0.;
-  auto const thetaRad = theta / 180. * M_PI;
-
-  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-    return sqrt((Elab - m) * (Elab + m));
-  };
-  HEPMomentumType P0 = elab2plab(E0, mass);
-  auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
-    return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
-  };
-
-  auto const [px, py, pz] = momentumComponents(thetaRad, P0);
-  auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
-  cout << "input particle: " << beamCode << endl;
-  cout << "input angles: theta=" << theta << endl;
-  cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
-       << endl;
-
-  auto const observationHeight = 0_km + builder.getEarthRadius();
-  auto const injectionHeight = 112.75_km + builder.getEarthRadius();
-  auto const t = -observationHeight * cos(thetaRad) +
-                 sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
-                      si::detail::static_pow<2>(injectionHeight));
-  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
-  Point const injectionPos =
-      showerCore +
-      Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
-
-  std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
-
-  if (A != 1) {
-    stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                 corsika::stack::MomentumVector, geometry::Point,
-                                 units::si::TimeType, unsigned short, unsigned short>{
-        beamCode, E0, plab, injectionPos, 0_ns, A, Z});
-
-  } else {
-    stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::Proton, E0, plab, injectionPos, 0_ns});
-  }
-
-  std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
-            << std::endl;
-
-  environment::ShowerAxis const showerAxis{injectionPos,
-                                           (showerCore - injectionPos) * 1.02, env};
-
-  // setup processes, decays and interactions
-
-  process::sibyll::Interaction sibyll;
-  process::interaction_counter::InteractionCounter sibyllCounted(sibyll);
-
-  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::interaction_counter::InteractionCounter sibyllNucCounted(sibyllNuc);
-
-  process::pythia::Decay decayPythia;
-
-  // use sibyll decay routine for decays of particles unknown to pythia
-  process::sibyll::Decay decaySibyll{{
-      Code::N1440Plus,
-      Code::N1440MinusBar,
-      Code::N1440_0,
-      Code::N1440_0Bar,
-      Code::N1710Plus,
-      Code::N1710MinusBar,
-      Code::N1710_0,
-      Code::N1710_0Bar,
-
-      Code::Pi1300Plus,
-      Code::Pi1300Minus,
-      Code::Pi1300_0,
-
-      Code::KStar0_1430_0,
-      Code::KStar0_1430_0Bar,
-      Code::KStar0_1430_Plus,
-      Code::KStar0_1430_MinusBar,
-  }};
-
-  decaySibyll.PrintDecayConfig();
-
-  // PROPOSAL processs proposal{...};
-  process::particle_cut::ParticleCut cut{60_GeV, false, true};
-  process::proposal::Interaction proposal(env, cut.GetECut());
-  process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
-  process::interaction_counter::InteractionCounter proposalCounted(proposal);
-
-  process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
-
-  process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
-
-  Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-  process::observation_plane::ObservationPlane observationLevel(obsPlane,
-                                                                "particles.dat");
-
-  process::UrQMD::UrQMD urqmd;
-  process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
-
-  // assemble all processes into an ordered process list
-
-  auto sibyllSequence = sibyllNucCounted << sibyllCounted;
-  process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
-                                                       55_GeV);
-  auto decaySequence = decayPythia << decaySibyll;
-
-  auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << em_continuous
-				<< cut << longprof << observationLevel;
-
-  // define air shower object, run simulation
-  tracking_line::TrackingLine tracking;
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-
-  // to fix the point of first interaction, uncomment the following two lines:
-  //  EAS.SetNodes();
-  //  EAS.forceInteraction();
-
-  EAS.Run();
-
-
-  cut.ShowResults();
-  em_continuous.ShowResults();
-  observationLevel.ShowResults();
-  cout << "Cascade energy cut: " << EAS.GetEnergyCut()/1_GeV << " GeV" << endl;
-  const HEPEnergyType Efinal =
-    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround() + EAS.GetEnergyCut();
-  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-  observationLevel.Reset();
-  cut.Reset();
-  em_continuous.Reset();
-
-  auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
-    urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
-
-  hists.saveLab("inthist_lab.txt");
-  hists.saveCMS("inthist_cms.txt");
-
-  hists.saveLab("inthist_lab.txt");
-  hists.saveCMS("inthist_cms.txt");
-
-  longprof.save("longprof.txt");
-
-  std::ofstream finish("finished");
-  finish << "run completed without error" << std::endl;
-}

Environment/CMakeLists.txt

-set (
-  ENVIRONMENT_SOURCES
-  LayeredSphericalAtmosphereBuilder.cc
-  ShowerAxis.cc
-)
-
-set (
-  ENVIRONMENT_HEADERS
-  VolumeTreeNode.h
-  IMediumModel.h
-  NuclearComposition.h
-  HomogeneousMedium.h
-  InhomogeneousMedium.h
-  HomogeneousMedium.h
-  LinearApproximationIntegrator.h
-  DensityFunction.h
-  Environment.h
-  NameModel.h
-  BaseExponential.h
-  FlatExponential.h
-  SlidingPlanarExponential.h
-  LayeredSphericalAtmosphereBuilder.h
-  ShowerAxis.h
-  IMagneticFieldModel.h
-  UniformMagneticField.h
-  IRefractiveIndexModel.h
-  UniformRefractiveIndex.h
-  )
-
-set (
-  ENVIRONMENT_NAMESPACE
-  corsika/environment
-  )
-
-add_library (CORSIKAenvironment STATIC ${ENVIRONMENT_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
-
-set_target_properties (
-  CORSIKAenvironment
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${ENVIRONMENT_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAenvironment
-  CORSIKAgeometry
-  CORSIKAparticles
-  CORSIKAunits
-  CORSIKArandom
-  )
-
-target_include_directories (
-  CORSIKAenvironment
-  PUBLIC
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  TARGETS CORSIKAenvironment
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testEnvironment)
-target_link_libraries (
-  testEnvironment
-  CORSIKAsetup
-  CORSIKAenvironment
-  CORSIKAtesting
-  )
-
-CORSIKA_ADD_TEST(testShowerAxis)
-target_link_libraries (
-  testShowerAxis
-  CORSIKAsetup
-  CORSIKAenvironment
-  CORSIKAtesting
-  )

Environment/Environment.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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 <limits>
-
-namespace corsika::environment {
-  struct Universe : public corsika::geometry::Sphere {
-    Universe(corsika::geometry::CoordinateSystem const& pCS)
-        : corsika::geometry::Sphere(
-              corsika::geometry::Point{pCS, 0 * corsika::units::si::meter,
-                                       0 * corsika::units::si::meter,
-                                       0 * corsika::units::si::meter},
-              corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
-
-    bool Contains(corsika::geometry::Point const&) const override { return true; }
-  };
-
-  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;
-
-    auto& GetUniverse() { return fUniverse; }
-    auto const& GetUniverse() const { return fUniverse; }
-
-    auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
-
-    // factory method for creation of VolumeTreeNodes
-    template <class TVolumeType, typename... TVolumeArgs>
-    static auto CreateNode(TVolumeArgs&&... args) {
-      static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
-                    "unusable type provided, needs to be derived from "
-                    "\"corsika::geometry::Volume\"");
-
-      return std::make_unique<BaseNodeType>(
-          std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
-    }
-
-  private:
-    corsika::geometry::CoordinateSystem const& fCoordinateSystem;
-    typename BaseNodeType::VTNUPtr fUniverse;
-  };
-
-  // using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
-  // using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
-
-} // namespace corsika::environment

Environment/FlatExponential.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/environment/BaseExponential.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  //clang-format off
-  /**
-   * flat exponential density distribution with
-   * \f[
-   *  \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
-   *    \vec{a} \right).
-   * \f]
-   * \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
-   * with the scale parameter \f$ \lambda \f$.
-   */
-  //clang-format on
-  template <class T>
-  class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
-    geometry::Vector<units::si::dimensionless_d> const fAxis;
-    NuclearComposition const fNuclComp;
-
-    using Base = BaseExponential<FlatExponential<T>>;
-
-  public:
-    FlatExponential(geometry::Point const& vP0,
-                    geometry::Vector<units::si::dimensionless_d> const& vAxis,
-                    units::si::MassDensityType vRho, units::si::LengthType vLambda,
-                    NuclearComposition vNuclComp)
-        : Base(vP0, vRho, vLambda)
-        , fAxis(vAxis)
-        , fNuclComp(vNuclComp) {}
-
-    units::si::MassDensityType GetMassDensity(geometry::Point const& vP) const override {
-      return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
-    }
-
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::LengthType vTo) const override {
-      return Base::IntegratedGrammage(vLine, vTo, fAxis);
-    }
-
-    units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::GrammageType vGrammage) const override {
-      return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
-    }
-  };
-} // namespace corsika::environment

Environment/HomogeneousMedium.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cassert>
-
-/**
- * a homogeneous medium
- */
-
-namespace corsika::environment {
-
-  template <class T>
-  class HomogeneousMedium : public T {
-    corsika::units::si::MassDensityType const fDensity;
-    NuclearComposition const fNuclComp;
-
-  public:
-    HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
-                      NuclearComposition pNuclComp)
-        : fDensity(pDensity)
-        , fNuclComp(pNuclComp) {}
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const override {
-      return fDensity;
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::LengthType pTo) const override {
-      using namespace corsika::units::si;
-      return pTo * fDensity;
-    }
-
-    corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return pGrammage / fDensity;
-    }
-  };
-
-} // namespace corsika::environment

Environment/IMediumModel.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  class IMediumModel {
-  public:
-    virtual ~IMediumModel() = default; // LCOV_EXCL_LINE
-
-    virtual corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const = 0;
-
-    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
-    // approach; for now, only lines are supported
-    virtual corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::LengthType) const = 0;
-
-    virtual corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType) const = 0;
-
-    virtual NuclearComposition const& GetNuclearComposition() const = 0;
-  };
-
-} // namespace corsika::environment

Environment/InhomogeneousMedium.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-
-/**
- * A general inhomogeneous medium. The mass density distribution TDensityFunction must be
- * a \f$C^2\f$-function.
- */
-
-namespace corsika::environment {
-
-  template <class T, class TDensityFunction>
-  class InhomogeneousMedium : public T {
-    NuclearComposition const fNuclComp;
-    TDensityFunction const fDensityFunction;
-
-  public:
-    template <typename... Args>
-    InhomogeneousMedium(NuclearComposition pNuclComp, Args&&... rhoArgs)
-        : fNuclComp(pNuclComp)
-        , fDensityFunction(rhoArgs...) {}
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const& p) const override {
-      return fDensityFunction.EvaluateAt(p);
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
-        corsika::units::si::LengthType pTo) const override {
-      return fDensityFunction.IntegrateGrammage(pLine, pTo);
-    }
-
-    corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
-    }
-  };
-
-} // namespace corsika::environment

Environment/LayeredSphericalAtmosphereBuilder.cc

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/environment/FlatExponential.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
-#include <corsika/environment/SlidingPlanarExponential.h>
-
-using namespace corsika::environment;
-
-void LayeredSphericalAtmosphereBuilder::checkRadius(units::si::LengthType r) const {
-  if (r <= previousRadius_) {
-    throw std::runtime_error("radius must be greater than previous");
-  }
-}
-
-void LayeredSphericalAtmosphereBuilder::setNuclearComposition(
-    NuclearComposition composition) {
-  composition_ = std::make_unique<NuclearComposition>(composition);
-}
-
-void LayeredSphericalAtmosphereBuilder::addExponentialLayer(
-    units::si::GrammageType b, units::si::LengthType c,
-    units::si::LengthType upperBoundary) {
-  auto const radius = earthRadius_ + upperBoundary;
-  checkRadius(radius);
-  previousRadius_ = radius;
-
-  auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
-      std::make_unique<geometry::Sphere>(center_, radius));
-
-  auto const rho0 = b / c;
-  std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
-
-  node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(
-      center_, rho0, -c, *composition_, earthRadius_);
-
-  layers_.push(std::move(node));
-}
-
-void LayeredSphericalAtmosphereBuilder::addLinearLayer(
-    units::si::LengthType c, units::si::LengthType upperBoundary) {
-  using namespace units::si;
-
-  auto const radius = earthRadius_ + upperBoundary;
-  checkRadius(radius);
-  previousRadius_ = radius;
-
-  units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
-  auto const rho0 = b / c;
-
-  std::cout << "rho0 = " << rho0;
-
-  auto node = std::make_unique<VolumeTreeNode<environment::IMediumModel>>(
-      std::make_unique<geometry::Sphere>(center_, radius));
-  node->SetModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
-
-  layers_.push(std::move(node));
-}
-
-Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
-  Environment<IMediumModel> env;
-  assemble(env);
-  return env;
-}
-
-void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
-  auto& universe = env.GetUniverse();
-  auto* outmost = universe.get();
-
-  while (!layers_.empty()) {
-    auto l = std::move(layers_.top());
-    auto* tmp = l.get();
-    outmost->AddChild(std::move(l));
-    layers_.pop();
-    outmost = tmp;
-  }
-}

Environment/LayeredSphericalAtmosphereBuilder.h

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/environment/Environment.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/units/PhysicalConstants.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <memory>
-#include <stack>
-
-namespace corsika::environment {
-
-  class LayeredSphericalAtmosphereBuilder {
-    std::unique_ptr<NuclearComposition> composition_;
-    geometry::Point center_;
-    units::si::LengthType previousRadius_{units::si::LengthType::zero()};
-    units::si::LengthType earthRadius_;
-
-    std::stack<VolumeTreeNode<environment::IMediumModel>::VTNUPtr>
-        layers_; // innermost layer first
-
-    void checkRadius(units::si::LengthType) const;
-
-  public:
-    LayeredSphericalAtmosphereBuilder(
-        corsika::geometry::Point center,
-        units::si::LengthType earthRadius = units::constants::EarthRadius::Mean)
-        : center_(center)
-        , earthRadius_(earthRadius) {}
-
-    void setNuclearComposition(NuclearComposition);
-
-    void addExponentialLayer(units::si::GrammageType, units::si::LengthType,
-                             units::si::LengthType);
-
-    auto size() const { return layers_.size(); }
-
-    void addLinearLayer(units::si::LengthType, units::si::LengthType);
-
-    void assemble(Environment<IMediumModel>&);
-    Environment<IMediumModel> assemble();
-
-    /**
-     * Get the current Earth radius.
-     */
-    units::si::LengthType getEarthRadius() const { return earthRadius_; };
-  };
-
-} // namespace corsika::environment

Environment/LinearApproximationIntegrator.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <limits>
-
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-
-namespace corsika::environment {
-  template <class TDerived>
-  class LinearApproximationIntegrator {
-    auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
-
-  public:
-    auto IntegrateGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-        corsika::units::si::LengthType length) const {
-      auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-      return (c0 + 0.5 * c1 * length) * length;
-    }
-
-    auto ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-        corsika::units::si::GrammageType grammage) const {
-      auto const c0 = GetImplementation().fRho(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
-    }
-
-    auto MaximumLength(corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-                       [[maybe_unused]] double relError) const {
-      using namespace corsika::units::si;
-      [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      // todo: provide a real, working implementation
-      return 1_m * std::numeric_limits<double>::infinity();
-    }
-  };
-} // namespace corsika::environment

Environment/NuclearComposition.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cassert>
-#include <functional>
-#include <numeric>
-#include <random>
-#include <stdexcept>
-#include <vector>
-
-
-    
-namespace corsika::environment {
-  class NuclearComposition {
-    std::vector<float> const fNumberFractions; //!< relative fractions of number density
-    std::vector<corsika::particles::Code> const
-        fComponents; //!< particle codes of consitutents
-
-    double const fAvgMassNumber;
-
-    std::size_t hash_;
-    
-    template <class AConstIterator, class BConstIterator>
-    class WeightProviderIterator {
-      AConstIterator fAIter;
-      BConstIterator fBIter;
-
-    public:
-      using value_type = double;
-      using iterator_category = std::input_iterator_tag;
-      using pointer = value_type*;
-      using reference = value_type&;
-      using difference_type = ptrdiff_t;
-
-      WeightProviderIterator(AConstIterator a, BConstIterator b)
-          : fAIter(a)
-          , fBIter(b) {}
-
-      value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
-
-      WeightProviderIterator& operator++() { // prefix ++
-        ++fAIter;
-        ++fBIter;
-        return *this;
-      }
-
-      auto operator==(WeightProviderIterator other) { return fAIter == other.fAIter; }
-
-      auto operator!=(WeightProviderIterator other) { return !(*this == other); }
-    };
-
-  public:
-    NuclearComposition(std::vector<corsika::particles::Code> pComponents,
-                       std::vector<float> pFractions)
-        : fNumberFractions(pFractions)
-        , fComponents(pComponents)
-        , fAvgMassNumber(std::inner_product(
-              pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
-              std::plus<double>(), [](auto const compID, auto const fraction) -> double {
-                if (particles::IsNucleus(compID)) {
-                  return particles::GetNucleusA(compID) * fraction;
-                } else {
-                  return particles::GetMass(compID) /
-                         units::si::ConvertSIToHEP(units::constants::u) * fraction;
-                }
-              })) {
-      assert(pComponents.size() == pFractions.size());
-      auto const sumFractions =
-          std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-      if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
-        throw std::runtime_error("element fractions do not add up to 1");
-      }
-      updateHash();
-    }
-
-    template <typename TFunction>
-    auto WeightedSum(TFunction func) const {
-      using ResultQuantity = decltype(func(*fComponents.cbegin()));
-
-      auto const prod = [&](auto const compID, auto const fraction) {
-        return func(compID) * fraction;
-      };
-
-      if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
-        return std::inner_product(
-            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
-            ResultQuantity::zero(), // .zero() is defined for quantity types only
-            std::plus<ResultQuantity>(), prod);
-      } else {
-        return std::inner_product(
-            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
-            ResultQuantity(0), // in other cases we have to use a bare 0
-            std::plus<ResultQuantity>(), prod);
-      }
-    }
-
-    auto size() const { return fNumberFractions.size(); }
-
-    auto const& GetFractions() const { return fNumberFractions; }
-    auto const& GetComponents() const { return fComponents; }
-    auto const GetAverageMassNumber() const { return fAvgMassNumber; }
-
-    template <class TRNG>
-    corsika::particles::Code SampleTarget(
-        std::vector<corsika::units::si::CrossSectionType> const& sigma,
-        TRNG& randomStream) const {
-      using namespace corsika::units::si;
-
-      assert(sigma.size() == fNumberFractions.size());
-
-      std::discrete_distribution channelDist(
-          WeightProviderIterator<decltype(fNumberFractions.begin()),
-                                 decltype(sigma.begin())>(fNumberFractions.begin(),
-                                                          sigma.begin()),
-          WeightProviderIterator<decltype(fNumberFractions.begin()),
-                                 decltype(sigma.end())>(fNumberFractions.end(),
-                                                        sigma.end()));
-
-      auto const iChannel = channelDist(randomStream);
-      return fComponents[iChannel];
-    }
-
-    // Note: when this class ever modifies its internal data, the hash
-    // must be updated, too!
-    size_t hash() const { return hash_; }
-    
-  private:
-    void updateHash() const {
-      std::vector<std::size_t> hashes;
-      for (float ifrac : GetFractions())
-	hashes.push_back(std::hash<float>{}(ifrac));
-      for (corsika::particles::Code icode : GetComponents())
-	hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
-      std::size_t h = std::hash<double>{}(GetAverageMassNumber());
-      for (std::size_t ih : hashes)
-	h = h ^ (ih<<1);
-      hash_ = h;
-    }
-
-  };
-
-} // namespace corsika::environment
-