First very limited experiment of the electro magnetic process with proposal.

The proposal library has to be installed manually and corresponding cmake variables have to be exported.

Documentation/Examples/CMakeLists.txt

@@ -14,72 +14,26 @@ 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
-  CORSIKAprocesses
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  )
-
-CORSIKA_ADD_EXAMPLE (cascade_example)
-target_link_libraries (cascade_example
-  SuperStupidStack
-  CORSIKAunits
-  CORSIKAlogging
-  CORSIKArandom
-  ProcessSibyll
-  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
+    ProcessProposal
     CORSIKAcascade
-    ProcessEnergyLoss
     ProcessTrackWriter
-    ProcessStackInspector
-    ProcessTrackingLine
     ProcessParticleCut
-    ProcessHadronicElasticModel
-    ProcessStackInspector
+    ProcessTrackingLine
     CORSIKAprocesses
-    CORSIKAcascade
     CORSIKAparticles
     CORSIKAgeometry
     CORSIKAenvironment
     CORSIKAprocesssequence
     )
 
-  CORSIKA_ADD_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
-  target_link_libraries (vertical_EAS
+CORSIKA_ADD_EXAMPLE (cascade_example)
+target_link_libraries (cascade_example
     SuperStupidStack
     CORSIKAunits
     CORSIKAlogging
@@ -91,14 +45,14 @@ if (Pythia8_FOUND)
     CORSIKAcascade
     ProcessCONEXSourceCut
     ProcessEnergyLoss
-    ProcessObservationPlane
-    ProcessInteractionCounter
     ProcessTrackWriter
+    ProcessStackInspector
     ProcessTrackingLine
+    ProcessProposal
     ProcessParticleCut
     ProcessOnShellCheck
+    ProcessHadronicElasticModel
     ProcessStackInspector
-    ProcessLongitudinalProfile
     CORSIKAprocesses
     CORSIKAcascade
     CORSIKAparticles
@@ -106,21 +60,96 @@ if (Pythia8_FOUND)
     CORSIKAenvironment
     CORSIKAprocesssequence
     )
+
+if (Pythia8_FOUND)
+    CORSIKA_ADD_EXAMPLE (cascade_proton_example)
+    target_link_libraries (cascade_proton_example
+        SuperStupidStack
+        CORSIKAunits
+        CORSIKAlogging
+        CORSIKArandom
+        ProcessSibyll
+        ProcessPythia
+        CORSIKAcascade
+        ProcessEnergyLoss
+        ProcessTrackWriter
+        ProcessStackInspector
+        ProcessTrackingLine
+        ProcessParticleCut
+        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
+        ProcessSibyll
+        ProcessPythia
+        ProcessUrQMD
+        ProcessSwitch
+        CORSIKAcascade
+        ProcessEnergyLoss
+        ProcessObservationPlane
+        ProcessInteractionCounter
+        ProcessTrackWriter
+        ProcessTrackingLine
+        ProcessParticleCut
+        ProcessStackInspector
+        ProcessLongitudinalProfile
+        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
-  )
+    SuperStupidStack
+    CORSIKAunits
+    ProcessEnergyLoss
+    CORSIKAparticles
+    CORSIKAgeometry
+    CORSIKAenvironment
+    )
 
 CORSIKA_ADD_EXAMPLE (staticsequence_example)
 target_link_libraries (staticsequence_example
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAgeometry
-  CORSIKAlogging)
+    CORSIKAprocesssequence
+    CORSIKAunits
+    CORSIKAgeometry
+    CORSIKAlogging)
+
+CORSIKA_ADD_EXAMPLE (proposal_example RUN_OPTIONS 100.)
+target_link_libraries (proposal_example
+    SuperStupidStack
+    CORSIKAunits
+    CORSIKAlogging
+    CORSIKArandom
+    CORSIKAcascade
+    ProcessObservationPlane
+    ProcessInteractionCounter
+    ProcessTrackWriter
+    ProcessProposal
+    ProcessTrackingLine
+    ProcessParticleCut
+    ProcessStackInspector
+    ProcessLongitudinalProfile
+    CORSIKAprocesses
+    CORSIKAcascade
+    CORSIKAparticles
+    CORSIKAgeometry
+    CORSIKAenvironment
+    CORSIKAprocesssequence
+    )

Documentation/Examples/cascade_example.cc

@@ -23,6 +23,8 @@
 
 #include <corsika/geometry/Sphere.h>
 
+#include <corsika/process/proposal/Interaction.h>
+
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
@@ -135,31 +137,34 @@ int main() {
 
   random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
   random::RNGManager::GetInstance().RegisterRandomStream("pythia");
+  random::RNGManager::GetInstance().RegisterRandomStream("proposal");
   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);
+  process::proposal::Interaction proposal(env, cut);
 
   process::track_writer::TrackWriter trackWriter("tracks.dat");
   process::energy_loss::EnergyLoss eLoss{showerAxis};
 
   // assemble all processes into an ordered process list
-  auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
-                               << trackWriter;
+  auto sequence = stackInspect << sibyll << sibyllNuc << proposal << decay
+      /* << eLoss */
+      << cut << trackWriter;
 
   // define air shower object, run simulation
   cascade::Cascade EAS(env, tracking, sequence, stack);
 
   EAS.Run();
 
-  eLoss.PrintProfile(); // print longitudinal profile
+  /* 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;
+  /* cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl */
+  /*      << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl; */
 }

Documentation/Examples/proposal_example.cc

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/cascade/Cascade.h>
+#include <corsika/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/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 <corsika/process/track_writer/TrackWriter.h>
+#include <corsika/process/proposal/Interaction.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");
+  // add PROPOSAL here (?)
+  random::RNGManager::GetInstance().SeedAll();
+}
+
+int main(int argc, char** argv) {
+  if (argc != 2) {
+    std::cerr << "usage: proposal_example <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.earthRadius;
+  auto const injectionHeight = 112.75_km + builder.earthRadius;
+  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{...};
+  PROPOSAL::InterpolationDef::path_to_tables = "~/.local/share/PROPOSAL/tables/";
+  PROPOSAL::InterpolationDef::path_to_tables_readonly = "~/.local/share/PROPOSAL/tables/";
+
+  process::particle_cut::ParticleCut cut(10_GeV);
+  process::proposal::Interaction proposal(env, cut);
+  process::interaction_counter::InteractionCounter proposalCounted(proposal);
+  process::track_writer::TrackWriter trackWriter("tracks.dat");
+
+  // energy cut; n.b. ParticleCut needs to be modified not to discard EM particles!
+  /* process::particle_cut::ParticleCut cut{60_GeV}; */
+
+  // 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 << longprof << proposal << cut << observationLevel << trackWriter;
+  // define air shower object, run simulation
+  tracking_line::TrackingLine tracking;
+  cascade::Cascade EAS(env, tracking, sequence, stack);
+  EAS.Init();
+
+  // to fix the point of first interaction, uncomment the following two lines:
+  //  EAS.SetNodes();
+  //  EAS.forceInteraction();
+
+  EAS.Run();
+
+  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;
+
+  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;
+}

Processes/CMakeLists.txt

 # general
 add_subdirectory (NullModel)
 # tracking
-add_subdirectory (TrackingLine) 
+add_subdirectory (TrackingLine)
 # hadron interaction models
 add_subdirectory (Sibyll)
 add_subdirectory (QGSJetII)
@@ -13,6 +13,8 @@ if (CONEX_FOUND)
 endif (CONEX_FOUND)
 add_subdirectory (HadronicElasticModel)
 add_subdirectory (UrQMD)
+add_subdirectory (SwitchProcess)
+add_subdirectory (Proposal)
 
 # continuous physics
 add_subdirectory (EnergyLoss)
@@ -34,6 +36,7 @@ add_subdirectory (SwitchProcess)
 add_library (CORSIKAprocesses INTERFACE)
 add_dependencies(CORSIKAprocesses ProcessNullModel)
 add_dependencies(CORSIKAprocesses ProcessSibyll)
+add_dependencies(CORSIKAprocesses ProcessProposal)
 if (Pythia8_FOUND)
   add_dependencies(CORSIKAprocesses ProcessPythia)
 endif (Pythia8_FOUND)

Processes/ParticleCut/ParticleCut.cc

@@ -66,12 +66,13 @@ namespace corsika::process {
       C8LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid,
                               energy / 1_GeV,
                               (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV));
-      if (ParticleIsEmParticle(pid)) {
-        C8LOG_DEBUG("removing em. particle...");
-        fEmEnergy += energy;
-        fEmCount += 1;
-        return true;
-      } else if (ParticleIsInvisible(pid)) {
+      /* if (ParticleIsEmParticle(pid)) { */
+      /*   C8LOG_DEBUG("removing em. particle..."); */
+      /*   fEmEnergy += energy; */
+      /*   fEmCount += 1; */
+      /*   return true; */
+      /* } else */
+      if (ParticleIsInvisible(pid)) {
         C8LOG_DEBUG("removing inv. particle...");
         fInvEnergy += energy;
         fInvCount += 1;

Processes/ParticleCut/ParticleCut.h

@@ -31,7 +31,8 @@ namespace corsika::process {
       unsigned int fInvCount = 0;
 
     public:
-      ParticleCut(const units::si::HEPEnergyType vCut);
+      ParticleCut(const units::si::HEPEnergyType eCut)
+          : fECut(eCut) {}
 
       EProcessReturn DoSecondaries(corsika::setup::StackView&);
 
@@ -44,6 +45,7 @@ namespace corsika::process {
 
       void ShowResults();
 
+      units::si::HEPEnergyType GetECut() const { return fECut; }
       units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
       units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
       units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }

Processes/Proposal/CMakeLists.txt

+# find_package(PROPOSAL REQUIRED ${LIB_INCLUDE})
+
+
+FILE (GLOB MODEL_SOURCES *.cc)
+FILE (GLOB MODEL_HEADERS *.h)
+SET (MODEL_NAMESPACE corsika/process/proposal)
+
+ADD_LIBRARY (ProcessProposal STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessProposal ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+SET_TARGET_PROPERTIES ( ProcessProposal PROPERTIES VERSION ${PROJECT_VERSION}
+    SOVERSION 1
+    #  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+TARGET_LINK_LIBRARIES (
+  ProcessProposal
+  CORSIKAparticles
+  CORSIKAutilities
+  CORSIKAunits
+  CORSIKAthirdparty
+  CORSIKAgeometry
+  CORSIKAenvironment
+  ${PROPOSAL_LIBRARY}
+  )
+
+TARGET_INCLUDE_DIRECTORIES (
+  ProcessProposal
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+TARGET_INCLUDE_DIRECTORIES (
+  ProcessProposal
+  SYSTEM
+  PUBLIC ${PROPOSAL_INCLUDE_DIR}
+  INTERFACE ${PROPOSAL_INCLUDE_DIR}
+  )
+
+INSTALL (
+  TARGETS ProcessProposal
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )

Processes/Proposal/ContinuousProcess.cc

+#include <PROPOSAL/PROPOSAL.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/proposal/ContinuousProcess.h>
+#include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+
+using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupTrack = corsika::setup::Trajectory;
+
+namespace corsika::process::proposal {
+  using namespace corsika::setup;
+  using namespace corsika::environment;
+  using namespace corsika::units::si;
+
+  unordered_map<particles::Code, PROPOSAL::ParticleDef> ContinuousProcess::particles{
+      {particles::Code::Gamma, PROPOSAL::GammaDef()},
+      {particles::Code::Electron, PROPOSAL::EMinusDef()},
+      {particles::Code::Positron, PROPOSAL::EPlusDef()},
+      {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
+      {particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
+      {particles::Code::TauPlus, PROPOSAL::TauPlusDef()},
+      {particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
+  };
+
+  bool ContinuousProcess::CanInteract(particles::Code pcode) const noexcept {
+    auto search = particles.find(pcode);
+    if (search != particles.end()) return true;
+    return false;
+  }
+
+  template <>
+  ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
+                                       CORSIKA_ParticleCut const& _cut)
+      : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
+                                                           false)) {
+    auto all_compositions = std::vector<const NuclearComposition*>();
+    _env.GetUniverse()->walk([&](auto& vtn) {
+      if (vtn.HasModelProperties())
+        all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+    });
+    for (auto& ncarg : all_compositions) {
+      auto comp_vec = std::vector<PROPOSAL::Components::Component>();
+      auto frac_iter = ncarg->GetFractions().cbegin();
+      for (auto& pcode : ncarg->GetComponents()) {
+        comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
+                              *frac_iter);
+        ++frac_iter;
+      }
+      media[ncarg] = PROPOSAL::Medium(
+          "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
+          PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
+          PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
+    }
+  }
+
+  HEPEnergyType ContinuousProcess::TotalEnergyLoss(SetupParticle const& vP,
+                                                   GrammageType const vDX) {
+    auto calc_ptr = GetCalculator(vP);
+    auto upper_energy = calc_ptr->second->UpperLimitTrackIntegral(
+        vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm);
+    std::cout << "upper_energy: " << upper_energy << "MeV" << std::endl;
+    return upper_energy * 1_MeV;
+  }
+
+   void ContinuousProcess::Init() {}
+  template <>
+  EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) {
+    if (vP.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
+    std::cout << "DoContinuous..." << std::endl;
+    GrammageType const dX =
+        vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
+    HEPEnergyType dE = TotalEnergyLoss(vP, dX);
+    auto E = vP.GetEnergy();
+    const auto Ekin = E - vP.GetMass();
+    auto Enew = E + dE;
+    auto status = process::EProcessReturn::eOk;
+    if (-dE > Ekin) {
+      dE = -Ekin;
+      Enew = vP.GetMass();
+      status = process::EProcessReturn::eParticleAbsorbed;
+    }
+    vP.SetEnergy(Enew);
+    auto pnew = vP.GetMomentum();
+    vP.SetMomentum(pnew * Enew / pnew.GetNorm());
+    return status;
+  }
+
+  template <>
+  units::si::LengthType ContinuousProcess::MaxStepLength(SetupParticle const& vP,
+                                                         SetupTrack const& vT) {
+    auto constexpr dX = 1_g / square(1_cm);
+    auto const dE = -TotalEnergyLoss(vP, dX); // dE > 0
+    auto const maxLoss = 0.01 * vP.GetEnergy();
+    auto const maxGrammage = maxLoss / dE * dX;
+
+    return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT,
+                                                                    maxGrammage) *
+           1.0001; // to make sure particle gets absorbed when DoContinuous() is called
+  }
+
+} // namespace corsika::process::proposal

Processes/Proposal/ContinuousProcess.h

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _corsika_process_proposal_interaction_h_
+#define _corsika_process_proposal_interaction_h_
+
+#include <PROPOSAL/PROPOSAL.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <unordered_map>
+#include "PROPOSAL/PROPOSAL.h"
+
+using std::unordered_map;
+
+using namespace corsika::environment;
+
+using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
+
+namespace corsika::process::proposal {
+
+  class ContinuousProcess
+      : public corsika::process::ContinuousProcess<ContinuousProcess> {
+  private:
+    shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
+
+    static unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+    unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
+
+    corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream("proposal");
+
+    bool CanInteract(particles::Code pcode) const noexcept;
+
+    struct interaction_hash {
+        size_t operator()(const std::pair<const NuclearComposition*, particles::Code>& p) const
+        {
+            auto hash1 = std::hash<const NuclearComposition*>{}(p.first);
+            auto hash2 = std::hash<particles::Code>{}(p.second);
+            return hash1 ^ hash2;
+        }
+    };
+
+    unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc;
+
+    auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
+        auto medium = media.at(&comp);
+        if (corsika_code == particles::Code::Gamma) {
+            auto cross = GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+                    calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::Electron) {
+            auto cross = GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::Positron) {
+            auto cross = GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::MuMinus) {
+            auto cross = GetStdCrossSections(PROPOSAL::MuMinusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::MuPlus) {
+            auto cross = GetStdCrossSections(PROPOSAL::MuPlusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::TauMinus) {
+            auto cross = GetStdCrossSections(PROPOSAL::TauMinusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::TauPlus) {
+            auto cross = GetStdCrossSections(PROPOSAL::TauPlusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        throw std::runtime_error("PROPOSAL could not find corresponding builder");
+    }
+    
+    template <typename Particle>
+    auto GetCalculator(Particle& vP) {
+      auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+      auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID()));
+      if (calc_it != calc.end()) return calc_it;
+      return BuildCalculator(vP.GetPID(), comp);
+    }
+
+    units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
+                                             const units::si::GrammageType);
+
+  public:
+    template <typename TEnvironment>
+    ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+
+    void Init();
+
+    template <typename Particle, typename Track>
+    EProcessReturn DoContinuous(Particle&, Track const&) ;
+
+    template <typename Particle, typename Track>
+    units::si::LengthType MaxStepLength(Particle const& p, Track const& track) ;
+  };
+} // namespace corsika::process::proposal
+
+#endif

Processes/Proposal/Interaction.cc

+
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+#include <limits>
+#include <memory>
+#include <random>
+#include <tuple>
+
+using Component_PROPOSAL = PROPOSAL::Components::Component;
+
+namespace corsika::process::proposal {
+  using namespace corsika::setup;
+  using namespace corsika::environment;
+  using namespace corsika::units::si;
+
+  std::unordered_map<particles::Code, PROPOSAL::ParticleDef> Interaction::particles{
+      {particles::Code::Gamma, PROPOSAL::GammaDef()},
+      {particles::Code::Electron, PROPOSAL::EMinusDef()},
+      {particles::Code::Positron, PROPOSAL::EPlusDef()},
+      {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
+      {particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
+      {particles::Code::TauPlus, PROPOSAL::TauPlusDef()},
+      {particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
+  };
+
+  bool Interaction::CanInteract(particles::Code pcode) const noexcept {
+    if (particles.find(pcode) != particles.end()) return true;
+    return false;
+  }
+
+  template <>
+  Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut)
+      : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
+                                                           false))
+      , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
+    auto all_compositions = std::vector<const NuclearComposition*>();
+    _env.GetUniverse()->walk([&](auto& vtn) {
+      if (vtn.HasModelProperties()) {
+        all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+      }
+    });
+    for (auto& ncarg : all_compositions) {
+      auto comp_vec = std::vector<Component_PROPOSAL>();
+      auto frac_iter = ncarg->GetFractions().cbegin();
+      for (auto& pcode : ncarg->GetComponents()) {
+        comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
+                              *frac_iter);
+        ++frac_iter;
+      }
+      media[ncarg] = PROPOSAL::Medium(
+          "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
+          PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
+          PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
+    }
+  }
+
+  template <>
+  corsika::process::EProcessReturn Interaction::DoInteraction(
+      setup::StackView::StackIterator& vP) {
+    if (CanInteract(vP.GetPID())) {
+      auto calc = GetCalculator(vP); // [CrossSections]
+      std::uniform_real_distribution<double> distr(0., 1.);
+      auto [type, comp_ptr, v] =
+          get<INTERACTION>(calc->second)
+              ->TypeInteraction(vP.GetEnergy() / 1_MeV, distr(fRNG));
+      auto rnd =
+          vector<double>(get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type));
+      for (auto& it : rnd) it = distr(fRNG);
+      auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
+                                      vP.GetPosition().GetY() / 1_cm,
+                                      vP.GetPosition().GetZ() / 1_cm);
+      auto d = vP.GetDirection().GetComponents();
+      auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
+                                          d.GetZ().magnitude());
+      auto loss = make_tuple(static_cast<int>(type), point, direction,
+                             v * vP.GetEnergy() / 1_MeV, 0.);
+      auto sec = get<SECONDARIES>(calc->second)
+                     ->CalculateSecondaries(vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
+      for (auto& s : sec) {
+        auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
+        auto vec = corsika::geometry::QuantityVector(
+            get<PROPOSAL::Loss::DIRECTION>(s).GetX() * E,
+            get<PROPOSAL::Loss::DIRECTION>(s).GetY() * E,
+            get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * E);
+        auto p = corsika::stack::MomentumVector(
+            corsika::geometry::RootCoordinateSystem::GetInstance()
+                .GetRootCoordinateSystem(),
+            vec);
+        auto sec_code = corsika::particles::ConvertFromPDG(
+            static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+        vP.AddSecondary(make_tuple(sec_code, E, p, vP.GetPosition(), vP.GetTime()));
+      }
+    }
+    return process::EProcessReturn::eOk;
+  }
+
+  template <>
+  corsika::units::si::GrammageType Interaction::GetInteractionLength(
+      setup::Stack::StackIterator const& vP) {
+    if (CanInteract(vP.GetPID())) {
+      auto calc = GetCalculator(vP);
+      return get<INTERACTION>(calc->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
+             (1_cm * 1_cm);
+    }
+    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+  }
+} // namespace corsika::process::proposal

Processes/Proposal/Interaction.h

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _corsika_process_proposal_interaction_h_
+#define _corsika_process_proposalythia_interaction_h_
+
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <unordered_map>
+#include "PROPOSAL/PROPOSAL.h"
+
+using namespace corsika::environment;
+
+using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
+using std::make_pair;
+using std::make_tuple;
+
+namespace corsika::process::proposal {
+
+  class Interaction : public corsika::process::InteractionProcess<Interaction> {
+
+    shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
+
+    static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+    std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
+
+    corsika::random::RNG& fRNG;
+
+    bool CanInteract(particles::Code pcode) const noexcept;
+
+    using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
+                               unique_ptr<PROPOSAL::Interaction>>;
+    using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
+
+    struct interaction_hash {
+      size_t operator()(const calc_key_t& p) const {
+        return std::hash<const NuclearComposition*>{}(p.first) ^
+               std::hash<particles::Code>{}(p.second);
+      }
+    };
+
+    std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
+
+    template <typename Particle>
+    auto BuildCalculator(particles::Code corsika_code, Particle p_def,
+                         NuclearComposition const& comp) {
+      auto cross = GetStdCrossSections(p_def, media.at(&comp), cut, true);
+      auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
+      auto [insert_it, success] = calculators.insert(
+          {make_pair(&comp, corsika_code),
+           make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
+                      PROPOSAL::make_interaction(cross, true))});
+      return insert_it;
+    }
+
+    auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
+      if (corsika_code == particles::Code::Gamma)
+        return BuildCalculator(particles::Code::Gamma, PROPOSAL::GammaDef(), comp);
+      if (corsika_code == particles::Code::Electron)
+        return BuildCalculator(particles::Code::Electron, PROPOSAL::EMinusDef(), comp);
+      if (corsika_code == particles::Code::Positron)
+        return BuildCalculator(particles::Code::Positron, PROPOSAL::EPlusDef(), comp);
+      if (corsika_code == particles::Code::MuMinus)
+        return BuildCalculator(particles::Code::MuMinus, PROPOSAL::MuMinusDef(), comp);
+      if (corsika_code == particles::Code::MuPlus)
+        return BuildCalculator(particles::Code::MuPlus, PROPOSAL::MuPlusDef(), comp);
+      if (corsika_code == particles::Code::TauMinus)
+        return BuildCalculator(particles::Code::TauMinus, PROPOSAL::TauMinusDef(), comp);
+      if (corsika_code == particles::Code::TauPlus)
+        return BuildCalculator(particles::Code::TauPlus, PROPOSAL::TauPlusDef(), comp);
+      throw std::runtime_error("PROPOSAL could not find corresponding builder");
+    }
+
+    enum { SECONDARIES, INTERACTION };
+    template <typename Particle>
+    auto GetCalculator(Particle& vP) {
+      auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+      auto calc_it = calculators.find(make_pair(&comp, vP.GetPID()));
+      if (calc_it != calculators.end()) return calc_it;
+      return BuildCalculator(vP.GetPID(), comp);
+    }
+
+  public:
+    template <typename TEnvironment>
+    Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+
+    void Init() {};
+
+    template <typename Particle>
+    corsika::process::EProcessReturn DoInteraction(Particle&);
+
+    template <typename TParticle>
+    corsika::units::si::GrammageType GetInteractionLength(TParticle const& p);
+  }; // namespace corsika::process::proposal
+} // namespace corsika::process::proposal
+#endif