Processes/CONEXSourceCut/testCONEXSourceCut.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/Environment.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::environment;
-using namespace corsika::geometry;
-using namespace corsika::units::si;
-
-TEST_CASE("CONEXSourceCut") {
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-  random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
-
-  feenableexcept(FE_INVALID);
-
-  // 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, conex::earthRadius};
-  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);
-
-  const HEPEnergyType E0 = 1_PeV;
-  double thetaDeg = 60.;
-  auto const thetaRad = thetaDeg / 180. * M_PI;
-
-  auto const observationHeight = 1.4_km + conex::earthRadius;
-  auto const injectionHeight = 112.75_km + conex::earthRadius;
-  auto const t =
-      -observationHeight * cos(thetaRad) +
-      sqrt(-units::si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
-           units::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;
-
-  environment::ShowerAxis const showerAxis{injectionPos,
-                                           (showerCore - injectionPos) * 1.02, env};
-
-  // need to initialize Sibyll, done in constructor:
-  process::sibyll::Interaction sibyll;
-  [[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-
-  corsika::process::conex_source_cut::CONEXSourceCut conex(
-      center, showerAxis, t, injectionHeight, E0,
-      particles::GetPDG(particles::Code::Proton));
-
-  HEPEnergyType const Eem{1_PeV};
-  auto const momentum = showerAxis.GetDirection() * Eem;
-
-  auto const emPosition = showerCore + showerAxis.GetDirection() * (-20_km);
-
-  std::cout << "position injection: "
-            << injectionPos.GetCoordinates(conex.GetObserverCS()) << " "
-            << injectionPos.GetCoordinates(rootCS) << std::endl;
-  std::cout << "position core: " << showerCore.GetCoordinates(conex.GetObserverCS())
-            << " " << showerCore.GetCoordinates(rootCS) << std::endl;
-  std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " "
-            << emPosition.GetCoordinates(rootCS) << std::endl;
-
-  conex.addParticle(0, Eem, 0_eV, emPosition, momentum.normalized(), 0_s);
-
-  conex.SolveCE();
-}

Processes/EnergyLoss/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  EnergyLoss.cc
-  )
-
-set (
-  MODEL_HEADERS
-  EnergyLoss.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/energy_loss
-  )
-
-add_library (ProcessEnergyLoss STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessEnergyLoss ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessEnergyLoss
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessEnergyLoss
-  CORSIKAunits
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAsetup
-  )
-
-target_include_directories (
-  ProcessEnergyLoss 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessEnergyLoss
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-#CORSIKA_ADD_TEST (testNullModel testNullModel.cc)
-#target_link_libraries (
-#  testNullModel  ProcessNullModel
-#  CORSIKAsetup
-#  CORSIKAgeometry
-#  CORSIKAunits
-#  CORSIKAthirdparty # for catch2
-#  )
-

Processes/EnergyLoss/EnergyLoss.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/process/energy_loss/EnergyLoss.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/geometry/Line.h>
-
-#include <cmath>
-#include <fstream>
-#include <iostream>
-#include <limits>
-#include <numeric>
-
-using namespace std;
-
-using namespace corsika;
-using namespace corsika::units::si;
-using SetupParticle = corsika::setup::Stack::ParticleType;
-using SetupTrack = corsika::setup::Trajectory;
-
-using namespace corsika::process::energy_loss;
-
-EnergyLoss::EnergyLoss(environment::ShowerAxis const& shower_axis,
-                       corsika::units::si::HEPEnergyType emCut)
-    : shower_axis_(shower_axis)
-    , emCut_(emCut)
-    , profile_(int(shower_axis.maximumX() / dX_) + 1) {}
-
-auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-  return sqrt((Elab - m) * (Elab + m));
-};
-
-/**
- *   PDG2018, passage of particles through matter
- *
- * Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I =
- * \sum_i a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
- * (\f$\sim Z_i\f$) of the \f$i\f$-th element. This can also be used for shell
- * corrections or density effects.
- *
- * The \f$I_{\mathrm{eff}}\f$ of compounds is not better than a few percent, if not
- * measured explicitly.
- *
- * For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
- *
- */
-HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const dX) {
-
-  // all these are material constants and have to come through Environment
-  // right now: values for nitrogen_D
-  // 7 nitrogen_gas 82.0 0.49976 D E 0.0011653 0.0 1.7378 4.1323 0.15349 3.2125 10.54
-  auto Ieff = 82.0_eV;
-  [[maybe_unused]] auto Zmat = 7;
-  auto ZoverA = 0.49976_mol / 1_g;
-  const double x0 = 1.7378;
-  const double x1 = 4.1323;
-  const double Cbar = 10.54;
-  const double delta0 = 0.0;
-  const double aa = 0.15349;
-  const double sk = 3.2125;
-  // end of material constants
-
-  // this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
-  auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
-  HEPEnergyType const E = p.GetEnergy();
-  HEPMassType const m = p.GetMass();
-  double const gamma = E / m;
-  int const Z = p.GetChargeNumber();
-  int const Z2 = Z * Z;
-  HEPMassType constexpr me = particles::Electron::GetMass();
-  auto const m2 = m * m;
-  auto constexpr me2 = me * me;
-  double const gamma2 = gamma * gamma;
-
-  double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2    (1-1/gamma)*(1+1/gamma);
-                                              // (gamma_2-1)/gamma_2 = (1-1/gamma2);
-  double constexpr c2 = 1;                    // HEP convention here c=c2=1
-  cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
-  [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
-  HEPMassType const Wmax =
-      2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
-  // approx, but <<1%    HEPMassType const Wmax = 2*me*c2*beta2*gamma2;      for HEAVY
-  // PARTICLES Wmax ~ 2me v2 for non-relativistic particles
-  cout << "BetheBloch Wmax=" << Wmax << endl;
-
-  // Sternheimer parameterization, density corrections towards high energies
-  // NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
-  // MISSING
-  cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm() / m
-       << endl;
-  double const x = log10(p.GetMomentum().GetNorm() / m);
-  double delta = 0;
-  if (x >= x1) {
-    delta = 2 * (log(10)) * x - Cbar;
-  } else if (x < x1 && x >= x0) {
-    delta = 2 * (log(10)) * x - Cbar + aa * pow((x1 - x), sk);
-  } else if (x < x0) { // and IF conductor (otherwise, this is 0)
-    delta = delta0 * pow(100, 2 * (x - x0));
-  }
-  cout << "BetheBloch delta=" << delta << endl;
-
-  // with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
-
-  // shell correction, <~100MeV
-  // need more clarity about formulas and units
-  const double Cadj = 0;
-  /*
-  // https://www.nap.edu/read/20066/chapter/8#104
-  HEPEnergyType Iadj = 12_eV * Z + 7_eV;  // Iadj<163eV
-  if (Iadj>=163_eV)
-    Iadj = 9.76_eV * Z + 58.8_eV * pow(Z, -0.19);  // Iadj>=163eV
-  double const Cadj = (0.422377/eta2 + 0.0304043/(eta2*eta2) -
-  0.00038106/(eta2*eta2*eta2)) * 1e-6 * Iadj*Iadj + (3.858019/eta2 -
-  0.1667989/(eta2*eta2) + 0.00157955/(eta2*eta2*eta2)) * 1e-9 * Iadj*Iadj*Iadj;
-  */
-
-  // 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 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);
-  // double const barkas = Z * Llow*Lhigh/(Llow+Lhigh); // RU, I think the Z was
-  // missing...
-  double const barkas = 1; // does not work yet
-
-  // Bloch correction for O(Z4) higher-order Born approximation
-  // see Appl. Phys. 85 (1999) 1249
-  const double alpha = 1. / 137.035999173;
-  double const y2 = Z * Z * alpha * alpha / beta2;
-  double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
-
-  // cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
-  // endl;
-
-  double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
-  return -K * Z2 * ZoverA / beta2 *
-         (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
-}
-
-// radiation losses according to PDG 2018, ch. 33 ref. [5]
-HEPEnergyType EnergyLoss::RadiationLosses(SetupParticle const& vP,
-                                          GrammageType const vDX) {
-  // simple-minded hard-coded value for b(E) inspired by data from
-  // http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O.
-  auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g;
-  return -vP.GetEnergy() * b * vDX;
-}
-
-HEPEnergyType EnergyLoss::TotalEnergyLoss(SetupParticle const& vP,
-                                          GrammageType const vDX) {
-  return BetheBloch(vP, vDX) + RadiationLosses(vP, vDX);
-}
-
-process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack const& t) {
-  if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
-
-  GrammageType const dX =
-      p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
-  cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
-       << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << endl;
-  HEPEnergyType dE = TotalEnergyLoss(p, dX);
-  auto E = p.GetEnergy();
-  const auto Ekin = E - p.GetMass();
-  auto Enew = E + dE;
-  cout << "EnergyLoss  dE=" << dE / 1_MeV << "MeV, "
-       << " E=" << E / 1_GeV << "GeV,  Ekin=" << Ekin / 1_GeV << ", Enew=" << Enew / 1_GeV
-       << "GeV" << endl;
-  auto status = process::EProcessReturn::eOk;
-  if (E < emCut_) {
-    Enew = emCut_;
-    status = process::EProcessReturn::eParticleAbsorbed;
-  }
-  p.SetEnergy(Enew);
-  MomentumUpdate(p, Enew);
-  FillProfile(t, dE);
-  return status;
-}
-
-LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
-                                     SetupTrack const& vTrack) const {
-  if (vParticle.GetChargeNumber() == 0) {
-    return units::si::meter * std::numeric_limits<double>::infinity();
-  }
-
-  auto constexpr dX = 1_g / square(1_cm);
-  auto const dEdX = -TotalEnergyLoss(vParticle, dX) / dX; // dE > 0
-  //~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
-
-  // in any case: never go below 0.99*emCut_ This needs to be
-  // slightly smaller than emCut_ since, either this Step is limited
-  // by energy_lim, then the particle is stopped in a very short
-  // range (before doing anythin else) and is then removed
-  // instantly. The exact position where it reaches emCut is not
-  // important, the important fact is that its E_kin is zero
-  // afterwards.
-  //
-  const auto energy = vParticle.GetEnergy();
-  auto energy_lim = std::max(0.9 * energy, 0.99 * emCut_);
-
-  auto const maxGrammage = (energy - energy_lim) / dEdX;
-
-  return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack,
-                                                                         maxGrammage);
-}
-
-void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
-                                corsika::units::si::HEPEnergyType Enew) {
-  HEPMomentumType Pnew = elab2plab(Enew, vP.GetMass());
-  auto pnew = vP.GetMomentum();
-  vP.SetMomentum(pnew * Pnew / pnew.GetNorm());
-}
-
-void EnergyLoss::FillProfile(SetupTrack const& vTrack, const HEPEnergyType dE) {
-
-  GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0));
-  GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
-  const auto deltaX = grammageEnd - grammageStart;
-
-  const int binStart = grammageStart / dX_;
-  const int binEnd = grammageEnd / dX_;
-
-  std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
-            << grammageEnd << std::endl;
-
-  auto energyCount = HEPEnergyType::zero();
-
-  auto fill = [&](int bin, GrammageType weight) {
-    if (deltaX > dX_threshold_) {
-      auto const increment = -dE * weight / deltaX;
-      profile_[bin] += increment;
-      energyCount += increment;
-
-      std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
-                << std::endl;
-    }
-  };
-
-  // fill longitudinal profile
-  fill(binStart, (1 + binStart) * dX_ - grammageStart);
-  fill(binEnd, grammageEnd - binEnd * dX_);
-
-  if (binStart == binEnd) { fill(binStart, -dX_); }
-
-  for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
-
-  std::cout << "total energy added to histogram: " << energyCount << std::endl;
-}
-
-void EnergyLoss::PrintProfile() const {
-  std::ofstream file("EnergyLossProfile.dat");
-  file << "# EnergyLoss profile" << std::endl
-       << "# lower X bin edge [g/cm2]  dE/dX [GeV/g/cm2]" << endl;
-  double const deltaX = dX_ / 1_g * square(1_cm);
-  for (size_t i = 0; i < profile_.size(); ++i) {
-    file << std::scientific << std::setw(15) << i * deltaX << std::setw(15)
-         << profile_.at(i) / (deltaX * 1_GeV) << endl;
-  }
-}
-
-HEPEnergyType EnergyLoss::GetTotal() const {
-  return std::accumulate(profile_.cbegin(), profile_.cend(), HEPEnergyType::zero());
-}

Processes/EnergyLoss/EnergyLoss.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/ShowerAxis.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/process/ContinuousProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <map>
-
-namespace corsika::process::energy_loss {
-
-  class EnergyLoss : public corsika::process::ContinuousProcess<EnergyLoss> {
-
-    using MeVgcm2 = decltype(1e6 * units::si::electronvolt / units::si::gram *
-                             units::si::square(1e-2 * units::si::meter));
-
-    void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
-
-  public:
-    EnergyLoss(environment::ShowerAxis const& showerAxis,
-               corsika::units::si::HEPEnergyType emCut);
-
-    process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
-                                         setup::Trajectory const&);
-    units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
-                                        setup::Trajectory const&) const;
-    units::si::HEPEnergyType GetTotal() const;
-    void PrintProfile() const;
-    static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&,
-                                               const units::si::GrammageType);
-    static units::si::HEPEnergyType RadiationLosses(setup::Stack::ParticleType const&,
-                                                    const units::si::GrammageType);
-    static units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
-                                                    const units::si::GrammageType);
-
-  private:
-    void FillProfile(setup::Trajectory const&, units::si::HEPEnergyType);
-
-    units::si::GrammageType const dX_ = std::invoke([]() {
-      using namespace units::si;
-      return 10_g / square(1_cm);
-    }); // profile binning
-
-  private:
-    environment::ShowerAxis const& shower_axis_;
-    corsika::units::si::HEPEnergyType emCut_;
-    std::vector<units::si::HEPEnergyType> profile_; // longitudinal profile
-  };
-
-  units::si::GrammageType const dX_threshold_ = std::invoke([]() {
-    using namespace units::si;
-    return 0.0001_g / square(1_cm);
-  });
-} // namespace corsika::process::energy_loss

Processes/HadronicElasticModel/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  HadronicElasticModel.cc
-  )
-
-set (
-  MODEL_HEADERS
-  HadronicElasticModel.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/hadronic_elastic_model
-  )
-
-add_library (ProcessHadronicElasticModel STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessHadronicElasticModel ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessHadronicElasticModel
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessHadronicElasticModel
-  CORSIKAunits
-  CORSIKAgeometry
-  CORSIKAsetup
-  CORSIKAutilities
-  )
-
-target_include_directories (
-  ProcessHadronicElasticModel 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessHadronicElasticModel
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-# CORSIKA_ADD_TEST (testProcessHadronicElasticModel)
-# target_link_libraries (
-  # testProcessHadronicElasticModel
-  # ProcessHadronicElasticModel
-  # CORSIKAsetup
-  # CORSIKAgeometry
-  # CORSIKAunits
-  # CORSIKAthirdparty # for catch2
-  # )

Processes/HadronicElasticModel/HadronicElasticModel.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/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/FourVector.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/random/ExponentialDistribution.h>
-#include <corsika/utl/COMBoost.h>
-
-#include <corsika/setup/SetupStack.h>
-
-#include <iomanip>
-#include <iostream>
-
-using SetupParticle = corsika::setup::Stack::ParticleType;
-using SetupView = corsika::setup::StackView;
-
-using namespace corsika::process::HadronicElasticModel;
-using namespace corsika;
-
-HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x,
-                                                       units::si::CrossSectionType y)
-    : fX(x)
-    , fY(y) {}
-
-template <>
-units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
-    SetupParticle const& p) {
-  using namespace units::si;
-  if (p.GetPID() == particles::Code::Proton) {
-    auto const* currentNode = p.GetNode();
-    auto const& mediumComposition =
-        currentNode->GetModelProperties().GetNuclearComposition();
-
-    auto const& components = mediumComposition.GetComponents();
-    auto const& fractions = mediumComposition.GetFractions();
-
-    auto const projectileMomentum = p.GetMomentum();
-    auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
-    auto const projectileEnergy = p.GetEnergy();
-
-    auto const avgCrossSection = [&]() {
-      CrossSectionType avgCrossSection = 0_b;
-
-      for (size_t i = 0; i < fractions.size(); ++i) {
-        auto const targetMass = particles::GetMass(components[i]);
-        auto const s = detail::static_pow<2>(projectileEnergy + targetMass) -
-                       projectileMomentumSquaredNorm;
-        avgCrossSection += CrossSection(s) * fractions[i];
-      }
-
-      std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl;
-
-      return avgCrossSection;
-    }();
-
-    auto const avgTargetMassNumber = mediumComposition.GetAverageMassNumber();
-
-    GrammageType const interactionLength =
-        avgTargetMassNumber * units::constants::u / avgCrossSection;
-
-    return interactionLength;
-  } else {
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-}
-
-template <>
-process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupView& view) {
-  using namespace units::si;
-  using namespace units::constants;
-
-  auto p = view.GetProjectile();
-  if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
-
-  const auto* currentNode = p.GetNode();
-  const auto& composition = currentNode->GetModelProperties().GetNuclearComposition();
-  const auto& components = composition.GetComponents();
-
-  std::vector<units::si::CrossSectionType> cross_section_of_components(
-      composition.GetComponents().size());
-
-  auto const projectileMomentum = p.GetMomentum();
-  auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
-  auto const projectileEnergy = p.GetEnergy();
-
-  for (size_t i = 0; i < components.size(); ++i) {
-    auto const targetMass = particles::GetMass(components[i]);
-    auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) -
-                   projectileMomentumSquaredNorm;
-    cross_section_of_components[i] = CrossSection(s);
-  }
-
-  const auto targetCode = composition.SampleTarget(cross_section_of_components, fRNG);
-
-  auto const targetMass = particles::GetMass(targetCode);
-
-  std::uniform_real_distribution phiDist(0., 2 * M_PI);
-
-  geometry::FourVector const projectileLab(projectileEnergy, projectileMomentum);
-  geometry::FourVector const targetLab(
-      targetMass, geometry::Vector<units::si::hepmomentum_d>(
-                      projectileMomentum.GetCoordinateSystem(), {0_eV, 0_eV, 0_eV}));
-  utl::COMBoost const boost(projectileLab, targetMass);
-
-  auto const projectileCoM = boost.toCoM(projectileLab);
-  auto const targetCoM = boost.toCoM(targetLab);
-
-  auto const pProjectileCoMSqNorm = projectileCoM.GetSpaceLikeComponents().squaredNorm();
-  auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm);
-
-  auto const eProjectileCoM = projectileCoM.GetTimeLikeComponent();
-  auto const eTargetCoM = targetCoM.GetTimeLikeComponent();
-
-  auto const sqrtS = eProjectileCoM + eTargetCoM;
-  auto const s = units::si::detail::static_pow<2>(sqrtS);
-
-  auto const B = this->B(s);
-  std::cout << B << std::endl;
-
-  random::ExponentialDistribution tDist(1 / B);
-  auto const absT = [&]() {
-    decltype(tDist(fRNG)) absT;
-    auto const maxT = 4 * pProjectileCoMSqNorm;
-
-    do {
-      // |t| cannot become arbitrarily large, max. given by GER eq. (4.16), so we just
-      // throw again until we have an acceptable value. Note that the formula holds in
-      // any frame despite of what is stated in the book.
-      absT = tDist(fRNG);
-    } while (absT >= maxT);
-
-    return absT;
-  }();
-
-  std::cout << "HadronicElasticInteraction: s = " << s * invGeVsq
-            << " GeV²; absT = " << absT * invGeVsq
-            << " GeV² (max./GeV² = " << 4 * invGeVsq * projectileMomentumSquaredNorm
-            << ')' << std::endl;
-
-  auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm)));
-  auto const phi = phiDist(fRNG);
-
-  auto const projectileScatteredLab =
-      boost.fromCoM(geometry::FourVector<HEPEnergyType, geometry::Vector<hepmomentum_d>>(
-          eProjectileCoM,
-          geometry::Vector<hepmomentum_d>(projectileMomentum.GetCoordinateSystem(),
-                                          {pProjectileCoMNorm * sin(theta) * cos(phi),
-                                           pProjectileCoMNorm * sin(theta) * sin(phi),
-                                           pProjectileCoMNorm * cos(theta)})));
-
-  p.SetMomentum(projectileScatteredLab.GetSpaceLikeComponents());
-  p.SetEnergy(
-      sqrt(projectileScatteredLab.GetSpaceLikeComponents().squaredNorm() +
-           units::si::detail::static_pow<2>(particles::GetMass(
-               p.GetPID())))); // Don't use energy from boost. It can be smaller than
-                               // the momentum due to limited numerical accuracy.
-
-  return process::EProcessReturn::eOk;
-}
-
-HadronicElasticInteraction::inveV2 HadronicElasticInteraction::B(eV2 s) const {
-  using namespace units::constants;
-  auto constexpr b_p = 2.3;
-  auto const result =
-      (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq;
-  std::cout << "B(" << s << ") = " << result / invGeVsq << " GeV¯²" << std::endl;
-  return result;
-}
-
-units::si::CrossSectionType HadronicElasticInteraction::CrossSection(
-    SquaredHEPEnergyType s) const {
-  using namespace units::si;
-  using namespace units::constants;
-  // assuming every target behaves like a proton, fX and fY are universal
-  CrossSectionType const sigmaTotal =
-      fX * pow(s * invGeVsq, gfEpsilon) + fY * pow(s * invGeVsq, -gfEta);
-
-  // according to Schuler & Sjöstrand, PRD 49, 2257 (1994)
-  // (we ignore rho because rho^2 is just ~2 %)
-  auto const sigmaElastic =
-      units::si::detail::static_pow<2>(sigmaTotal) /
-      (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s)));
-
-  std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl;
-  std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mb << " mb" << std::endl;
-  return sigmaElastic;
-}

Processes/InteractionCounter/CMakeLists.txt

-set (
-  MODEL_HEADERS
-  InteractionCounter.h
-  InteractionHistogram.h
-  )
-
-set (
-  MODEL_SOURCES
-  InteractionHistogram.cc
-)
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/interaction_counter
-  )
-
-add_library (ProcessInteractionCounter STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessInteractionCounter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessInteractionCounter
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessInteractionCounter
-  CORSIKAunits
-  CORSIKAprocesssequence
-  CORSIKAthirdparty
-  C8::ext::boost
-  )
-
-target_include_directories (
-  ProcessInteractionCounter 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testInteractionCounter SOURCES
-  testInteractionCounter.cc
-  ${MODEL_HEADERS}
-)
-
-target_link_libraries (
-  testInteractionCounter
-  ProcessInteractionCounter
-  CORSIKAsetup
-  CORSIKAunits
-  CORSIKAenvironment
-  CORSIKAtesting
-  )

Processes/InteractionCounter/InteractionCounter.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.
- */
-
-#pragma once
-
-#include <corsika/process/InteractionProcess.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/interaction_counter/InteractionHistogram.h>
-#include <corsika/setup/SetupStack.h>
-
-namespace corsika::process::interaction_counter {
-  /*!
-   * Wrapper around an InteractionProcess that fills histograms of the number
-   * of calls to DoInteraction() binned in projectile energy (both in
-   * lab and center-of-mass frame) and species
-   */
-  template <class TCountedProcess>
-  class InteractionCounter
-      : public InteractionProcess<InteractionCounter<TCountedProcess>> {
-
-    TCountedProcess& process_;
-    InteractionHistogram histogram_;
-
-  public:
-    InteractionCounter(TCountedProcess& process)
-        : process_(process) {}
-
-    template <typename TSecondaryView>
-    auto DoInteraction(TSecondaryView& view) {
-      auto const projectile = view.GetProjectile();
-      auto const massNumber = projectile.GetNode()
-                                  ->GetModelProperties()
-                                  .GetNuclearComposition()
-                                  .GetAverageMassNumber();
-      auto const massTarget = massNumber * units::constants::nucleonMass;
-
-      if (auto const projectile_id = projectile.GetPID();
-          projectile_id == particles::Code::Nucleus) {
-        auto const A = projectile.GetNuclearA();
-        auto const Z = projectile.GetNuclearZ();
-        histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget, A, Z);
-      } else {
-        histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget);
-      }
-      return process_.DoInteraction(view);
-    }
-
-    template <typename TParticle>
-    auto GetInteractionLength(TParticle const& particle) const {
-      return process_.GetInteractionLength(particle);
-    }
-
-    auto const& GetHistogram() const { return histogram_; }
-  };
-
-} // namespace corsika::process::interaction_counter

Processes/InteractionCounter/InteractionHistogram.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/process/interaction_counter/InteractionHistogram.h>
-
-#include <fstream>
-#include <string>
-
-using namespace corsika::process::interaction_counter;
-
-void InteractionHistogram::fill(particles::Code projectile_id,
-                                units::si::HEPEnergyType lab_energy,
-                                units::si::HEPEnergyType mass_target, int A, int Z) {
-  using namespace units::si;
-  if (projectile_id == particles::Code::Nucleus) {
-    auto const sqrtS =
-        sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) +
-             mass_target * mass_target + 2 * lab_energy * mass_target);
-
-    int32_t const pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
-
-    if (nuclear_inthist_cms_.count(pdg) == 0) {
-      nuclear_inthist_lab_.emplace(
-          pdg,
-          boost::histogram::make_histogram(
-              boost::histogram::axis::regular<double,
-                                              boost::histogram::axis::transform::log>(
-                  num_bins_lab, lower_edge_lab, upper_edge_lab)));
-
-      nuclear_inthist_cms_.emplace(
-          pdg,
-          boost::histogram::make_histogram(
-              boost::histogram::axis::regular<double,
-                                              boost::histogram::axis::transform::log>(
-                  num_bins_cms, lower_edge_cms, upper_edge_cms)));
-    }
-
-    nuclear_inthist_lab_[pdg](lab_energy / 1_GeV);
-    nuclear_inthist_cms_[pdg](sqrtS / 1_GeV);
-  } else {
-    auto const projectile_mass = particles::GetMass(projectile_id);
-    auto const sqrtS = sqrt(projectile_mass * projectile_mass +
-                            mass_target * mass_target + 2 * lab_energy * mass_target);
-
-    interaction_histogram_cms_(
-        static_cast<corsika::particles::CodeIntType>(projectile_id), sqrtS / 1_GeV);
-    interaction_histogram_lab_(
-        static_cast<corsika::particles::CodeIntType>(projectile_id), lab_energy / 1_GeV);
-  }
-}
-
-void InteractionHistogram::saveLab(std::string const& filename) const {
-  save(interaction_histogram_lab_, nuclear_inthist_lab_, filename, "lab system");
-}
-
-void InteractionHistogram::saveCMS(std::string const& filename) const {
-  save(interaction_histogram_lab_, nuclear_inthist_cms_, filename,
-       "center-of-mass system");
-}
-
-InteractionHistogram& InteractionHistogram::operator+=(
-    InteractionHistogram const& other) {
-  interaction_histogram_lab_ += other.interaction_histogram_lab_;
-  interaction_histogram_cms_ += other.interaction_histogram_cms_;
-
-  nuclear_inthist_lab_ += other.nuclear_inthist_lab_;
-  nuclear_inthist_cms_ += other.nuclear_inthist_cms_;
-
-  return *this;
-}
-
-InteractionHistogram InteractionHistogram::operator+(InteractionHistogram other) const {
-  other.nuclear_inthist_lab_ += nuclear_inthist_lab_;
-  other.nuclear_inthist_cms_ += nuclear_inthist_cms_;
-
-  other.interaction_histogram_lab_ += interaction_histogram_lab_;
-  other.interaction_histogram_cms_ += interaction_histogram_cms_;
-
-  return other;
-}
-
-void InteractionHistogram::saveHist(hist_type const& hist, std::string const& filename,
-                                    std::string const& comment) {
-  std::ofstream file;
-  file.open(filename);
-  saveHist(hist, file, comment);
-}
-
-void InteractionHistogram::saveHist(hist_type const& hist, std::ofstream& file,
-                                    std::string const& comment) {
-  auto const& energy_axis = hist.axis(1);
-
-  file << "# interaction count histogram (" << comment << ")" << std::endl
-       << "# " << energy_axis.size() << " bins between " << energy_axis.bin(0).lower()
-       << " and " << energy_axis.bin(energy_axis.size() - 1).upper() << " GeV"
-       << std::endl;
-
-  for (particles::CodeIntType p = 0;
-       p < static_cast<particles::CodeIntType>(particles::Code::LastParticle); ++p) {
-
-    if (auto pdg = static_cast<particles::PDGCodeType>(
-            particles::GetPDG(static_cast<particles::Code>(p)));
-        pdg < 1'000'000'000l) {
-      file << "# " << static_cast<particles::Code>(p) << std::endl;
-      file << pdg;
-      for (int i = 0; i < energy_axis.size(); ++i) { file << ' ' << hist.at(p, i); }
-      file << std::endl;
-    }
-  }
-}
-
-void InteractionHistogram::saveHistMap(nuclear_hist_type const& histMap,
-                                       std::ofstream& file) {
-  file << "# nuclei" << std::endl;
-  for (auto const& [pdg, hist] : histMap) {
-    auto const num_ebins_nucl = hist.axis(0).size();
-
-    file << pdg << ' ';
-    for (int i = 0; i < num_ebins_nucl; ++i) { file << ' ' << hist.at(i); }
-    file << std::endl;
-  }
-}
-
-void InteractionHistogram::save(hist_type const& hist, nuclear_hist_type const& histMap,
-                                std::string const& filename, std::string const& comment) {
-  std::ofstream file;
-  file.open(filename);
-
-  saveHist(hist, file, comment);
-  saveHistMap(histMap, file);
-}
-
-InteractionHistogram::InteractionHistogram()
-    : interaction_histogram_cms_{boost::histogram::make_histogram(
-          boost::histogram::axis::integer<short>(0, numIds),
-          boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
-              num_bins_cms, lower_edge_cms, upper_edge_cms))}
-    , interaction_histogram_lab_{boost::histogram::make_histogram(
-          boost::histogram::axis::integer<short>(0, numIds),
-          boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
-              num_bins_lab, lower_edge_lab, upper_edge_lab))} {}

Processes/InteractionCounter/InteractionHistogram.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/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <boost/histogram.hpp>
-#include <boost/histogram/ostream.hpp>
-#include <fstream>
-#include <functional>
-#include <map>
-#include <utility>
-
-namespace corsika::process::interaction_counter {
-  using hist_type = decltype(boost::histogram::make_histogram(
-      std::declval<boost::histogram::axis::integer<particles::CodeIntType>>(),
-      std::declval<boost::histogram::axis::regular<
-          double, boost::histogram::axis::transform::log>>()));
-
-  using nucl_hist_type = decltype(boost::histogram::make_histogram(
-      std::declval<boost::histogram::axis::regular<
-          double, boost::histogram::axis::transform::log>>()));
-
-  using nuclear_hist_type = std::map<int32_t, nucl_hist_type>;
-
-  template <typename TKey, typename TVal>
-  auto operator+=(std::map<TKey, TVal>& a, std::map<TKey, TVal> b) {
-    a.merge(b);
-    for (auto const& [id, hist] : b) { a[id] += hist; }
-    return a;
-  }
-
-  class InteractionHistogram {
-    static auto constexpr lower_edge_cms = 1., upper_edge_cms = 1e8;  // GeV sqrt s
-    static auto constexpr lower_edge_lab = 1., upper_edge_lab = 1e12; // GeV lab
-    static auto constexpr num_bins_lab = 120, num_bins_cms = 80;
-    static auto constexpr numIds =
-        static_cast<particles::CodeIntType>(particles::Code::LastParticle);
-
-    hist_type interaction_histogram_cms_;
-    hist_type interaction_histogram_lab_;
-
-    /*!
-     * These maps map PDG nuclei codes to their corresponding interaction histograms
-     */
-    nuclear_hist_type nuclear_inthist_lab_, nuclear_inthist_cms_;
-
-  public:
-    InteractionHistogram();
-
-    //! fill both CMS and lab histograms at the same time
-    void fill(particles::Code projectile_id, units::si::HEPEnergyType lab_energy,
-              units::si::HEPEnergyType mass_target, int A = 0, int Z = 0);
-
-    auto CMSHists() const {
-      return std::pair<decltype(interaction_histogram_cms_) const&,
-                       decltype(nuclear_inthist_cms_) const&>{interaction_histogram_cms_,
-                                                              nuclear_inthist_cms_};
-    }
-
-    auto labHists() const {
-      return std::pair<decltype(interaction_histogram_lab_) const&,
-                       decltype(nuclear_inthist_lab_) const&>(interaction_histogram_lab_,
-                                                              nuclear_inthist_lab_);
-    }
-
-    void saveLab(std::string const& filename) const;
-
-    void saveCMS(std::string const& filename) const;
-
-    InteractionHistogram& operator+=(InteractionHistogram const& other);
-
-    InteractionHistogram operator+(InteractionHistogram other) const;
-
-  private:
-    /*!
-     * Save a histogram into a text file. The \arg comment string is written
-     * into the header of the text file.
-     * This method is static so that you can sum the histograms of multiple
-     * InteractionProcesses before saving them into the same file.
-     */
-
-    static void saveHist(hist_type const& hist, std::string const& filename,
-                         std::string const& comment = "");
-
-    static void saveHist(hist_type const& hist, std::ofstream& file,
-                         std::string const& comment = "");
-
-    static void saveHistMap(nuclear_hist_type const& histMap, std::ofstream& file);
-
-    /*!
-     * save both the "normal" particle histograms as well as the "nuclear" histograms
-     * into the same file
-     */
-
-    static void save(hist_type const& hist, nuclear_hist_type const& histMap,
-                     std::string const& filename, std::string const& comment = "");
-  };
-
-} // namespace corsika::process::interaction_counter

Processes/InteractionCounter/testInteractionCounter.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/process/interaction_counter/InteractionCounter.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/setup/SetupStack.h>
-
-#include <catch2/catch.hpp>
-
-#include <numeric>
-
-using namespace corsika;
-using namespace corsika::process::interaction_counter;
-using namespace corsika::units;
-using namespace corsika::units::si;
-
-auto setupEnvironment(particles::Code target_code) {
-  // setup environment, geometry
-  auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
-  auto& universe = *(env->GetUniverse());
-  const geometry::CoordinateSystem& 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>{target_code},
-                                      std::vector<float>{1.}));
-
-  auto const* nodePtr = theMedium.get();
-  universe.AddChild(std::move(theMedium));
-
-  return std::make_tuple(std::move(env), &cs, nodePtr);
-}
-
-template <typename TNodeType>
-auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
-                geometry::CoordinateSystem const& cs) {
-  auto stack = std::make_unique<setup::Stack>();
-  auto constexpr mN = corsika::units::constants::nucleonMass;
-
-  geometry::Point const origin(cs, {0_m, 0_m, 0_m});
-  corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
-
-  HEPEnergyType const E0 =
-      sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
-  auto particle =
-      stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                    corsika::stack::MomentumVector, geometry::Point,
-                                    units::si::TimeType, unsigned short, unsigned short>{
-          particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ});
-
-  particle.SetNode(vNodePtr);
-  return std::make_tuple(
-      std::move(stack),
-      std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
-}
-
-template <typename TNodeType>
-auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
-                TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
-  auto stack = std::make_unique<setup::Stack>();
-
-  geometry::Point const origin(cs, {0_m, 0_m, 0_m});
-  corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
-
-  HEPEnergyType const E0 =
-      sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
-           pLab.squaredNorm());
-  auto particle = stack->AddParticle(
-      std::tuple<particles::Code, units::si::HEPEnergyType,
-                 corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-          vProjectileType, E0, pLab, origin, 0_ns});
-
-  particle.SetNode(vNodePtr);
-  return std::make_tuple(
-      std::move(stack),
-      std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
-}
-
-struct DummyProcess {
-  template <typename TParticle>
-  GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
-    return 100_g / 1_cm / 1_cm;
-  }
-
-  template <typename TParticle>
-  auto DoInteraction([[maybe_unused]] TParticle& projectile) {
-    return nullptr;
-  }
-};
-
-TEST_CASE("InteractionCounter") {
-  DummyProcess d;
-  InteractionCounter countedProcess(d);
-
-  SECTION("GetInteractionLength") {
-    REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
-  }
-
-  auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
-  [[maybe_unused]] auto& env_dummy = env;
-
-  SECTION("DoInteraction nucleus") {
-    unsigned short constexpr A = 14, Z = 7;
-    auto [stackPtr, secViewPtr] = setupStack(A, Z, 105_TeV, nodePtr, *csPtr);
-    REQUIRE(stackPtr->getEntries() == 1);
-    REQUIRE(secViewPtr->getEntries() == 0);
-
-    auto const ret = countedProcess.DoInteraction(*secViewPtr);
-    REQUIRE(ret == nullptr);
-
-    auto const& h = countedProcess.GetHistogram().labHists().second.at(1'000'070'140);
-    REQUIRE(h.at(50) == 1);
-    REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
-    auto const& h2 = countedProcess.GetHistogram().CMSHists().second.at(1'000'070'140);
-    REQUIRE(h2.at(32) == 1); // bin 1.584 .. 1.995 TeV √s
-    REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-  }
-
-  SECTION("DoInteraction Lambda") {
-    auto constexpr code = particles::Code::Lambda0;
-    auto constexpr codeInt = static_cast<particles::CodeIntType>(code);
-    auto [stackPtr, secViewPtr] = setupStack(code, 105_TeV, nodePtr, *csPtr);
-    REQUIRE(stackPtr->getEntries() == 1);
-    REQUIRE(secViewPtr->getEntries() == 0);
-
-    auto const ret = countedProcess.DoInteraction(*secViewPtr);
-    REQUIRE(ret == nullptr);
-
-    auto const& h = countedProcess.GetHistogram().labHists().first;
-    REQUIRE(h.at(codeInt, 50) == 1);
-    REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
-    auto const& h2 = countedProcess.GetHistogram().CMSHists().first;
-    REQUIRE(h2.at(codeInt, 32) == 1); // bin 1.584 .. 1.995 TeV √s
-    REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-  }
-}

Processes/LongitudinalProfile/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  LongitudinalProfile.cc
-  )
-
-set (
-  MODEL_HEADERS
-  LongitudinalProfile.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/longitudinal_profile
-  )
-
-add_library (ProcessLongitudinalProfile STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessLongitudinalProfile ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessLongitudinalProfile
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessLongitudinalProfile
-  CORSIKAenvironment
-  CORSIKAunits
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAsetup
-  )
-
-target_include_directories (
-  ProcessLongitudinalProfile 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessLongitudinalProfile
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-# CORSIKA_ADD_TEST(testNullModel)
-#target_link_libraries (
-#  testNullModel  ProcessNullModel
-#  CORSIKAsetup
-#  CORSIKAgeometry
-#  CORSIKAunits
-#  CORSIKAthirdparty # for catch2
-#  )
-

Processes/LongitudinalProfile/LongitudinalProfile.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/process/longitudinal_profile/LongitudinalProfile.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <cmath>
-#include <iomanip>
-#include <limits>
-
-using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
-
-using namespace corsika::process::longitudinal_profile;
-using namespace corsika::units::si;
-
-LongitudinalProfile::LongitudinalProfile(environment::ShowerAxis const& shower_axis)
-    : shower_axis_{shower_axis}
-    , profiles_{static_cast<unsigned int>(shower_axis.maximumX() / dX_) + 1} {}
-
-template <>
-corsika::process::EProcessReturn LongitudinalProfile::DoContinuous(Particle const& vP,
-                                                                   Track const& vTrack) {
-  auto const pid = vP.GetPID();
-
-  GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0));
-  GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
-
-  const int binStart = std::ceil(grammageStart / dX_);
-  const int binEnd = std::floor(grammageEnd / dX_);
-
-  for (int b = binStart; b <= binEnd; ++b) {
-    if (pid == particles::Code::Gamma) {
-      profiles_.at(b)[ProfileIndex::Gamma]++;
-    } else if (pid == particles::Code::Positron) {
-      profiles_.at(b)[ProfileIndex::Positron]++;
-    } else if (pid == particles::Code::Electron) {
-      profiles_.at(b)[ProfileIndex::Electron]++;
-    } else if (pid == particles::Code::MuPlus) {
-      profiles_.at(b)[ProfileIndex::MuPlus]++;
-    } else if (pid == particles::Code::MuMinus) {
-      profiles_.at(b)[ProfileIndex::MuMinus]++;
-    } else if (particles::IsHadron(pid)) {
-      profiles_.at(b)[ProfileIndex::Hadron]++;
-    }
-  }
-
-  return corsika::process::EProcessReturn::eOk;
-}
-
-void LongitudinalProfile::save(std::string const& filename) {
-  std::ofstream f{filename};
-  f << "# X / g·cm¯², gamma, e+, e-, mu+, mu-, all hadrons" << std::endl;
-  for (size_t b = 0; b < profiles_.size(); ++b) {
-    f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm));
-    for (auto const& N : profiles_.at(b)) {
-      f << std::setw(width_) << std::setprecision(precision_) << std::scientific << N;
-    }
-    f << std::endl;
-  }
-}

Processes/LongitudinalProfile/LongitudinalProfile.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.
- */
-
-#pragma once
-
-#include <corsika/environment/ShowerAxis.h>
-#include <corsika/process/ContinuousProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <array>
-#include <fstream>
-#include <limits>
-#include <string>
-
-namespace corsika::process::longitudinal_profile {
-
-  class LongitudinalProfile
-      : public corsika::process::ContinuousProcess<LongitudinalProfile> {
-
-  public:
-    LongitudinalProfile(environment::ShowerAxis const&);
-
-    template <typename Particle, typename Track>
-    corsika::process::EProcessReturn DoContinuous(Particle const&, Track const&);
-
-    template <typename Particle, typename Track>
-    corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&) {
-      return units::si::meter * std::numeric_limits<double>::infinity();
-    }
-
-    void save(std::string const&);
-
-  private:
-    units::si::GrammageType const dX_ = std::invoke([]() {
-      using namespace units::si;
-      return 10_g / square(1_cm);
-    }); // profile binning
-
-    environment::ShowerAxis const& shower_axis_;
-    using ProfileEntry = std::array<uint32_t, 6>;
-    enum ProfileIndex {
-      Gamma = 0,
-      Positron = 1,
-      Electron = 2,
-      MuPlus = 3,
-      MuMinus = 4,
-      Hadron = 5
-    };
-    std::vector<ProfileEntry> profiles_; // longitudinal profile
-
-    static int const width_ = 14;
-    static int const precision_ = 6;
-  };
-
-} // namespace corsika::process::longitudinal_profile

Processes/NullModel/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  NullModel.cc
-  )
-
-set (
-  MODEL_HEADERS
-  NullModel.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/null_model
-  )
-
-add_library (ProcessNullModel STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessNullModel ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessNullModel
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessNullModel
-  CORSIKAunits
-  CORSIKAgeometry
-  CORSIKAsetup
-  )
-
-target_include_directories (
-  ProcessNullModel 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessNullModel
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST (testNullModel)
-target_link_libraries (
-  testNullModel
-  ProcessNullModel
-  CORSIKAsetup
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAtesting
-  )

Processes/NullModel/NullModel.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/process/null_model/NullModel.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-using namespace corsika;
-namespace corsika::process::null_model {
-
-  NullModel::NullModel(units::si::LengthType maxStepLength)
-      : fMaxStepLength(maxStepLength) {}
-
-  template <>
-  process::EProcessReturn NullModel::DoContinuous(setup::Stack::ParticleType&,
-                                                  setup::Trajectory&) const {
-    return process::EProcessReturn::eOk;
-  }
-
-  template <>
-  units::si::LengthType NullModel::MaxStepLength(setup::Stack::ParticleType&,
-                                                 setup::Trajectory&) const {
-    return fMaxStepLength;
-  }
-
-} // namespace corsika::process::null_model

Processes/NullModel/NullModel.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/process/BaseProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process::null_model {
-
-  class NullModel : public corsika::process::BaseProcess<NullModel> {
-    corsika::units::si::LengthType const fMaxStepLength;
-
-  public:
-    NullModel(corsika::units::si::LengthType maxStepLength =
-                  corsika::units::si::meter * std::numeric_limits<double>::infinity());
-
-    template <typename Particle, typename Track>
-    process::EProcessReturn DoContinuous(Particle&, Track&) const;
-
-    template <typename Particle, typename Track>
-    corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const;
-  };
-
-} // namespace corsika::process::null_model

Processes/NullModel/testNullModel.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 <catch2/catch.hpp>
-
-#include <corsika/process/null_model/NullModel.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-using namespace corsika::units::si;
-using namespace corsika::process::null_model;
-using namespace corsika;
-
-TEST_CASE("NullModel", "[processes]") {
-
-  auto const& dummyCS =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  geometry::Point const origin(dummyCS, {0_m, 0_m, 0_m});
-  geometry::Vector<units::si::SpeedType::dimension_type> v(dummyCS, 0_m / second,
-                                                           0_m / second, 1_m / second);
-  geometry::Line line(origin, v);
-  geometry::Trajectory<geometry::Line> track(line, 10_s);
-
-  setup::Stack stack;
-  setup::Stack::ParticleType particle = stack.AddParticle(
-      std::tuple<particles::Code, units::si::HEPEnergyType,
-                 corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-          particles::Code::Electron, 100_GeV,
-          corsika::stack::MomentumVector(dummyCS, {0_GeV, 0_GeV, -1_GeV}),
-          geometry::Point(dummyCS, {0_m, 0_m, 10_km}), 0_ns});
-  SECTION("interface") {
-
-    NullModel model(10_m);
-
-    [[maybe_unused]] const process::EProcessReturn ret =
-        model.DoContinuous(particle, track);
-    LengthType const length = model.MaxStepLength(particle, track);
-
-    CHECK((length / 10_m) == Approx(1));
-  }
-}

Processes/ObservationPlane/CMakeLists.txt

-set (
-  MODEL_HEADERS
-  ObservationPlane.h
-  )
-
-set (
-  MODEL_SOURCES
-  ObservationPlane.cc
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/observation_plane
-  )
-
-add_library (ProcessObservationPlane STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessObservationPlane ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-#target dependencies on other libraries(also the header onlys)
-target_link_libraries (
-  ProcessObservationPlane
-  CORSIKAgeometry
-  CORSIKAprocesssequence
-  )
-
-target_include_directories (
-  ProcessObservationPlane 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
-
-#-- -- -- -- -- -- -- -- -- --
-#code unit testing
-CORSIKA_ADD_TEST(testObservationPlane)
-target_link_libraries (
-  testObservationPlane
-  ProcessObservationPlane
-  CORSIKAstackinterface
-  CORSIKAthirdparty # for catch2
-)

Processes/ObservationPlane/ObservationPlane.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/process/observation_plane/ObservationPlane.h>
-
-#include <fstream>
-#include <iostream>
-
-using namespace corsika::process::observation_plane;
-using namespace corsika::units::si;
-
-ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
-                                   std::string const& filename, bool deleteOnHit)
-    : plane_(obsPlane)
-    , outputStream_(filename)
-    , deleteOnHit_(deleteOnHit)
-    , energy_ground_(0_GeV)
-    , count_ground_(0) {
-  outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
-}
-
-corsika::process::EProcessReturn ObservationPlane::DoContinuous(
-    setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
-  TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetV0().dot(plane_.GetNormal());
-
-  if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
-
-  if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
-    return process::EProcessReturn::eOk;
-  }
-
-  const auto energy = particle.GetEnergy();
-  outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
-                << energy / 1_eV << ' '
-                << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
-                << std::endl;
-
-  if (deleteOnHit_) {
-    count_ground_++;
-    energy_ground_ += energy;
-    return process::EProcessReturn::eParticleAbsorbed;
-  } else {
-    return process::EProcessReturn::eOk;
-  }
-}
-
-LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
-                                           setup::Trajectory const& trajectory) {
-  TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetV0().dot(plane_.GetNormal());
-
-  if (timeOfIntersection < TimeType::zero()) {
-    return std::numeric_limits<double>::infinity() * 1_m;
-  }
-
-  auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
-  return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
-}
-
-void ObservationPlane::ShowResults() const {
-  std::cout << " ******************************" << std::endl
-            << " ObservationPlane: " << std::endl;
-  std::cout << " energy in ground (GeV)     :  " << energy_ground_ / 1_GeV << std::endl
-            << " no. of particles in ground :  " << count_ground_ << std::endl
-            << " ******************************" << std::endl;
-}
-
-void ObservationPlane::Reset() {
-  energy_ground_ = 0_GeV;
-  count_ground_ = 0;
-}