improve output, and shower diagnostics: hunting all the energy

Documentation/Examples/cascade_example.cc

@@ -162,4 +162,5 @@ int main() {
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
   cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
        << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+  cut.Reset();  
 }

Documentation/Examples/em_shower.cc

@@ -144,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;
@@ -157,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

@@ -207,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;
@@ -222,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;
 }
 
@@ -32,13 +35,16 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
   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())) << ' '
-                << particle.GetEnergy() * (1 / 1_eV) << ' '
+                << 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;
@@ -58,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

@@ -94,15 +94,26 @@ namespace corsika::process {
       energy_ = 0_GeV;
     }
 
-    void ParticleCut::ShowResults() {
+    void ParticleCut::ShowResults() const {
       cout << " ******************************" << endl
-           << " ParticleCut: " << endl
-           << " energy in em.  component (GeV):  " << emEnergy_ / 1_GeV << endl
-           << " no. of em.  particles injected:  " << emCount_ << endl
-           << " energy in inv. component (GeV):  " << invEnergy_ / 1_GeV << endl
-           << " no. of inv. particles injected:  " << invCount_ << endl
-           << " energy below particle cut (GeV): " << energy_ / 1_GeV << endl
+           << " ParticleCut: " << endl;
+      if (discardEm_) 
+	cout << " energy in em.  component (GeV)     :  " << emEnergy_ / 1_GeV << endl
+	     << " no. of em.  particles removed      :  " << emCount_ << endl;
+      if (discardInv_)
+	cout << " energy in inv. component (GeV)     :  " << invEnergy_ / 1_GeV << endl
+	     << " no. of inv. particles removed      :  " << invCount_ << endl;      
+      cout << " energy below energy threshold (GeV): " << energy_ / 1_GeV << endl
            << " ******************************" << endl;
     }
+
+    void ParticleCut::Reset() {
+      emEnergy_ = 0_GeV;
+      emCount_ = 0;
+      invEnergy_ = 0_GeV;
+      invCount_ = 0;
+      energy_ = 0_GeV;
+    }
+    
   } // namespace particle_cut
 } // namespace corsika::process

Processes/ParticleCut/ParticleCut.h

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

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