Now stores interpolation tables in $CORSIKA_DATA/PROPOSAL.

Remove proposal from cascade example, improve error message.
Improve output, and shower diagnostics: hunting all the energy Produce REAL error message.

Documentation/Examples/cascade_example.cc

@@ -23,8 +23,6 @@
 
 #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>
@@ -137,35 +135,32 @@ 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, true, true);
-  // 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 /* << proposal*/
-                               << decay
-                               /* << eLoss */
-                               << cut << trackWriter;
+  auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
+                               << trackWriter;
 
   // define air shower object, run simulation
   cascade::Cascade EAS(env, tracking, sequence, stack);
 
   EAS.Run();
 
-  /* eLoss.PrintProfile(); // print longitudinal profile */
+  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;
+  cut.Reset();  
 }

Documentation/Examples/em_shower.cc

@@ -130,9 +130,6 @@ int main(int argc, char** argv) {
   // 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, false, true);
   process::proposal::Interaction proposal(env, cut);
   process::proposal::ContinuousProcess em_continuous(env, cut);
@@ -147,7 +144,7 @@ int main(int argc, char** argv) {
   process::observation_plane::ObservationPlane observationLevel(obsPlane,
                                                                 "particles.dat");
 
-  auto sequence = proposalCounted << cut << em_continuous << longprof << observationLevel
+  auto sequence = proposalCounted << em_continuous << longprof << cut << observationLevel
                                   << trackWriter;
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
@@ -160,11 +157,16 @@ int main(int argc, char** argv) {
   EAS.Run();
 
   cut.ShowResults();
+  em_continuous.ShowResults();
+  observationLevel.ShowResults();
   const HEPEnergyType Efinal =
-      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround();
   cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-
+  observationLevel.Reset();
+  cut.Reset();
+  em_continuous.Reset();
+  
   auto const hists = proposalCounted.GetHistogram();
   hists.saveLab("inthist_lab.txt");
   hists.saveCMS("inthist_cms.txt");

Documentation/Examples/vertical_EAS.cc

@@ -184,9 +184,6 @@ int main(int argc, char** argv) {
   decaySibyll.PrintDecayConfig();
 
   // 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{60_GeV, false, true};
   process::proposal::Interaction proposal(env, cut);
   process::proposal::ContinuousProcess em_continuous(env, cut);
@@ -210,8 +207,8 @@ int main(int argc, char** argv) {
                                                        55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
 
-  auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << cut << em_continuous
-				<< longprof << observationLevel;
+  auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << em_continuous
+				<< cut << longprof << observationLevel;
 
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
@@ -225,10 +222,15 @@ int main(int argc, char** argv) {
 
 
   cut.ShowResults();
+  em_continuous.ShowResults();
+  observationLevel.ShowResults();
   const HEPEnergyType Efinal =
-      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround();
   cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+  observationLevel.Reset();
+  cut.Reset();
+  em_continuous.Reset();
 
   auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
     urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();

Processes/ObservationPlane/ObservationPlane.cc

@@ -9,6 +9,7 @@
 #include <corsika/process/observation_plane/ObservationPlane.h>
 
 #include <fstream>
+#include <iostream>
 
 using namespace corsika::process::observation_plane;
 using namespace corsika::units::si;
@@ -17,7 +18,9 @@ ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
                                    std::string const& filename, bool deleteOnHit)
     : plane_(obsPlane)
     , outputStream_(filename)
-    , deleteOnHit_(deleteOnHit) {
+    , deleteOnHit_(deleteOnHit)
+    , energy_ground_(0_GeV)
+    , count_ground_(0) {
   outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
 }
 
@@ -33,13 +36,19 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
     return process::EProcessReturn::eOk;
   }
 
+  const auto energy = particle.GetEnergy();
   outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
-                << particle.GetEnergy() * (1 / 1_eV) << ' '
+                << energy / 1_eV << ' '
                 << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
                 << std::endl;
 
-  if (deleteOnHit_) { return process::EProcessReturn::eParticleAbsorbed; }
-  return process::EProcessReturn::eOk;
+  if (deleteOnHit_) {
+    count_ground_++;
+    energy_ground_ += energy;
+    return process::EProcessReturn::eParticleAbsorbed;
+  } else {
+    return process::EProcessReturn::eOk;
+  }
 }
 
 LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
@@ -55,3 +64,16 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
   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;
+}

Processes/ObservationPlane/ObservationPlane.h

@@ -36,9 +36,16 @@ namespace corsika::process::observation_plane {
         corsika::setup::Stack::ParticleType const&,
         corsika::setup::Trajectory const& vTrajectory);
 
+    void ShowResults() const;
+    void Reset();
+    corsika::units::si::HEPEnergyType GetEnergyGround() const { return energy_ground_; }
+    
   private:
     geometry::Plane const plane_;
     std::ofstream outputStream_;
     bool const deleteOnHit_;
+
+    units::si::HEPEnergyType energy_ground_ = 0 * units::si::electronvolt;
+    unsigned int count_ground_ = 0;
   };
 } // namespace corsika::process::observation_plane

Processes/ParticleCut/ParticleCut.cc

@@ -131,5 +131,14 @@ namespace corsika::process {
           " ******************************",
           fEmEnergy / 1_GeV, uiEmCount, fInvEnergy / 1_GeV, uiInvCount, fEnergy / 1_GeV));
     }
+
+    void ParticleCut::Reset() {
+      fEmEnergy = 0_GeV;
+      uiEmCount = 0;
+      fInvEnergy = 0_GeV;
+      uiInvCount = 0;
+      fEnergy = 0_GeV;
+    }
+
   } // namespace particle_cut
 } // namespace corsika::process

Processes/ParticleCut/ParticleCut.h

@@ -52,7 +52,8 @@ namespace corsika::process {
 
       bool ParticleIsEmParticle(particles::Code) const;
 
-      void ShowResults();
+      void ShowResults() const;
+      void Reset();
 
       units::si::HEPEnergyType GetECut() const { return fECut; }
       units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }

Processes/Proposal/ContinuousProcess.cc

@@ -17,6 +17,8 @@
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/COMBoost.h>
 
+#include <iostream>
+
 namespace corsika::process::proposal {
 
   using namespace corsika::units::si;
@@ -97,10 +99,11 @@ namespace corsika::process::proposal {
     auto final_energy = get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
                             vP.GetEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
                         1_MeV;
-
+    auto dE = vP.GetEnergy() - final_energy;
+    energy_lost_ += dE;
+    
     // if the particle has a charge take multiple scattering into account
-    if (vP.GetChargeNumber() != 0)
-      Scatter(vP, vP.GetEnergy() - final_energy, dX);
+    if (vP.GetChargeNumber() != 0) Scatter(vP, dE, dX);
 
     // Update the energy and absorbe the particle if it's below the energy
     // threshold, because it will no longer propagated.
@@ -131,4 +134,14 @@ namespace corsika::process::proposal {
     return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
   }
 
+  void ContinuousProcess::ShowResults() const {
+    std::cout << " ******************************" << std::endl
+              << " PROCESS::ContinuousProcess: " << std::endl;
+    std::cout << " energy lost dE (GeV)      :  " << energy_lost_ / 1_GeV << std::endl;
+  }
+
+  void ContinuousProcess::Reset() {
+    energy_lost_ = 0_GeV;
+  }
+
 } // namespace corsika::process::proposal

Processes/Proposal/ContinuousProcess.h

@@ -20,8 +20,6 @@
 
 namespace corsika::process::proposal {
 
-  using namespace corsika::units::si;
-
   //!
   //! Electro-magnetic and gamma continous losses produced by proposal. It makes
   //! use of interpolation tables which are runtime intensive calculation, but can be
@@ -37,6 +35,8 @@ namespace corsika::process::proposal {
     std::unordered_map<calc_key_t, calc_t, hash>
         calc; //!< Stores the displacement and scattering calculators.
 
+    units::si::HEPEnergyType energy_lost_ = 0 * units::si::electronvolt;
+    
     //!
     //! Build the displacement and scattering calculators and add it to calc.
     //!
@@ -74,5 +74,9 @@ namespace corsika::process::proposal {
     //!
     template <typename Particle, typename Track>
     corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&);
+
+    void ShowResults() const;
+    void Reset();
+    corsika::units::si::HEPEnergyType GetEnergyLost() const { return energy_lost_; }
   };
 } // namespace corsika::process::proposal

Processes/Proposal/ProposalProcessBase.cc

@@ -14,6 +14,7 @@
 #include <corsika/setup/SetupTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/COMBoost.h>
+#include <cstdlib>
 #include <limits>
 #include <memory>
 #include <random>
@@ -50,6 +51,18 @@ namespace corsika::process::proposal {
     PROPOSAL::InterpolationDef::order_of_interpolation = 2;
     PROPOSAL::InterpolationDef::nodes_cross_section = 100;
     PROPOSAL::InterpolationDef::nodes_propagate = 1000;
+
+    //! If corsika data exist store interpolation tables to the corresponding
+    //! path, otherwise interpolation tables would only stored in main memory if
+    //! no explicit intrpolation def is specified.
+    if (auto data_path = std::getenv("CORSIKA_DATA")) {
+      PROPOSAL::InterpolationDef::path_to_tables = std::string(data_path) + "/PROPOSAL";
+    } else {
+      throw std::runtime_error(
+          "It is not recommended to run PROPOSAL without its tables in "
+          "$CORSIKA_DATA/PROPOSAL. This would be extremely slow. Please provide the "
+          "table directory. ");
+    }
   }
 
   size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {