Resolve "FPE due to zero track length in energy loss for very short lived particles"

Documentation/Examples/cascade_example.cc

@@ -124,7 +124,7 @@ int main() {
     cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
     stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
                                  corsika::stack::MomentumVector, geometry::Point,
-                                 units::si::TimeType, unsigned short, unsigned short>{
+		      units::si::TimeType, unsigned short, unsigned short>{
         beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ});
   }
 

Processes/EnergyLoss/EnergyLoss.cc

@@ -181,7 +181,7 @@ namespace corsika::process::energy_loss {
     }
     p.SetEnergy(Enew);
     MomentumUpdate(p, Enew);
-    fEnergyLossTot += dE;
+    EnergyLossTot_ += dE;
     FillProfile(p, t, dE);
     return status;
   }
@@ -215,13 +215,13 @@ namespace corsika::process::energy_loss {
 
     using namespace corsika::geometry;
 
-    auto const toStart = vTrack.GetPosition(0) - fInjectionPoint;
-    auto const toEnd = vTrack.GetPosition(1) - fInjectionPoint;
+    auto const toStart = vTrack.GetPosition(0) - InjectionPoint_;
+    auto const toEnd = vTrack.GetPosition(1) - InjectionPoint_;
 
-    auto const v1 = (toStart * 1_Hz).dot(fShowerAxisDirection);
-    auto const v2 = (toEnd * 1_Hz).dot(fShowerAxisDirection);
-    geometry::Line const lineToStartBin(fInjectionPoint, fShowerAxisDirection * v1);
-    geometry::Line const lineToEndBin(fInjectionPoint, fShowerAxisDirection * v2);
+    auto const v1 = (toStart * 1_Hz).dot(ShowerAxisDirection_);
+    auto const v2 = (toEnd * 1_Hz).dot(ShowerAxisDirection_);
+    geometry::Line const lineToStartBin(InjectionPoint_, ShowerAxisDirection_ * v1);
+    geometry::Line const lineToEndBin(InjectionPoint_, ShowerAxisDirection_ * v2);
 
     SetupTrack const trajToStartBin(lineToStartBin, 1_s);
     SetupTrack const trajToEndBin(lineToEndBin, 1_s);
@@ -233,8 +233,8 @@ namespace corsika::process::energy_loss {
         vP.GetNode()->GetModelProperties().IntegratedGrammage(trajToEndBin,
                                                               trajToEndBin.GetLength());
 
-    const int binStart = grammageStart / fdX;
-    const int binEnd = grammageEnd / fdX;
+    const int binStart = grammageStart / dX_;
+    const int binEnd = grammageEnd / dX_;
 
     std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
               << grammageEnd << std::endl;
@@ -242,21 +242,24 @@ namespace corsika::process::energy_loss {
     auto energyCount = HEPEnergyType::zero();
 
     auto fill = [&](int bin, GrammageType weight) {
-      auto const increment = -dE * weight / (grammageEnd - grammageStart);
-      fProfile[bin] += increment;
-      energyCount += increment;
-
-      std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
-                << std::endl;
+      const auto dX = grammageEnd - grammageStart;
+      if (dX > dX_threshold_) {
+        auto const increment = -dE * weight / (grammageEnd - grammageStart);
+        Profile_[bin] += increment;
+        energyCount += increment;
+
+        std::cout << "filling bin " << bin << " with weight " << weight << ": "
+                  << increment << std::endl;
+      }
     };
 
     // fill longitudinal profile
-    fill(binStart, (1 + binStart) * fdX - grammageStart);
-    fill(binEnd, grammageEnd - binEnd * fdX);
+    fill(binStart, (1 + binStart) * dX_ - grammageStart);
+    fill(binEnd, grammageEnd - binEnd * dX_);
 
-    if (binStart == binEnd) { fill(binStart, -fdX); }
+    if (binStart == binEnd) { fill(binStart, -dX_); }
 
-    for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, fdX); }
+    for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
 
     std::cout << "total energy added to histogram: " << energyCount << std::endl;
   }
@@ -264,8 +267,8 @@ namespace corsika::process::energy_loss {
   void EnergyLoss::PrintProfile() const {
     std::ofstream file("EnergyLossProfile.dat");
     cout << "# EnergyLoss PrintProfile  X-bin [g/cm2]  dE/dX [GeV/g/cm2]  " << endl;
-    double const deltaX = fdX / 1_g * square(1_cm);
-    for (auto v : fProfile) {
+    double const deltaX = dX_ / 1_g * square(1_cm);
+    for (auto v : Profile_) {
       file << v.first * deltaX << " " << v.second / (deltaX * 1_GeV) << endl;
     }
   }

Processes/EnergyLoss/EnergyLoss.h

@@ -34,8 +34,8 @@ namespace corsika::process::energy_loss {
     template <typename TDim>
     EnergyLoss(geometry::Point const& injectionPoint,
                geometry::Vector<TDim> const& direction)
-        : fInjectionPoint(injectionPoint)
-        , fShowerAxisDirection(direction.normalized()) {}
+        : InjectionPoint_(injectionPoint)
+        , ShowerAxisDirection_(direction.normalized()) {}
 
     EnergyLoss(setup::Trajectory const& trajectory)
         : EnergyLoss(trajectory.GetPosition(0), trajectory.GetV0()){};
@@ -45,7 +45,7 @@ namespace corsika::process::energy_loss {
                                          setup::Trajectory const&);
     units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
                                         setup::Trajectory const&) const;
-    units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
+    units::si::HEPEnergyType GetTotal() const { return EnergyLossTot_; }
     void PrintProfile() const;
     static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const&,
                                                const units::si::GrammageType);
@@ -60,16 +60,20 @@ namespace corsika::process::energy_loss {
     // void FillProfileAbsorbed(setup::Stack::ParticleType const&, setup::Trajectory
     // const&);
 
-    units::si::HEPEnergyType fEnergyLossTot = units::si::HEPEnergyType::zero();
-    units::si::GrammageType const fdX = std::invoke([]() {
+    units::si::HEPEnergyType EnergyLossTot_ = units::si::HEPEnergyType::zero();
+    units::si::GrammageType const dX_ = std::invoke([]() {
       using namespace units::si;
       return 10_g / square(1_cm);
     });                                               // profile binning
-    std::map<int, units::si::HEPEnergyType> fProfile; // longitudinal profile
-    geometry::Point const fInjectionPoint;
-    geometry::Vector<units::si::dimensionless_d> const fShowerAxisDirection;
+    std::map<int, units::si::HEPEnergyType> Profile_; // longitudinal profile
+    geometry::Point const InjectionPoint_;
+    geometry::Vector<units::si::dimensionless_d> const ShowerAxisDirection_;
   };
 
+  const units::si::GrammageType dX_threshold_ = std::invoke([]() {
+    using namespace units::si;
+    return 0.0001_g / square(1_cm);
+  });
 } // namespace corsika::process::energy_loss
 
 #endif

Processes/StackInspector/StackInspector.cc

@@ -32,9 +32,9 @@ template <typename TStack>
 StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
                                        const HEPEnergyType vE0)
     : StackProcess<StackInspector<TStack>>(vNStep)
-    , fReportStack(vReportStack)
-    , fE0(vE0)
-    , fStartTime(std::chrono::system_clock::now()) {}
+    , ReportStack_(vReportStack)
+    , E0_(vE0)
+    , StartTime_(std::chrono::system_clock::now()) {}
 
 template <typename TStack>
 StackInspector<TStack>::~StackInspector() {}
@@ -47,7 +47,7 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
   for (const auto& iterP : vS) {
     HEPEnergyType E = iterP.GetEnergy();
     Etot += E;
-    if (fReportStack) {
+    if (ReportStack_) {
       geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
                                                .GetRootCoordinateSystem(); // for printout
       auto pos = iterP.GetPosition().GetCoordinates(rootCS);
@@ -60,12 +60,15 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
   }
 
   auto const now = std::chrono::system_clock::now();
-  const std::chrono::duration<double> elapsed_seconds = now - fStartTime;
+  const std::chrono::duration<double> elapsed_seconds = now - StartTime_;
   std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
-  double const progress = (fE0 - Etot) / fE0;
+  auto const dE = E0_ - Etot;
+  if (dE < dE_threshold_) return process::EProcessReturn::eOk;
+  double const progress = dE / E0_;
+
   double const eta_seconds = elapsed_seconds.count() / progress;
   std::time_t const eta_time = std::chrono::system_clock::to_time_t(
-      fStartTime + std::chrono::seconds((int)eta_seconds));
+      StartTime_ + std::chrono::seconds((int)eta_seconds));
 
   cout << "StackInspector: "
        << " time=" << std::put_time(std::localtime(&now_time), "%T")
@@ -79,8 +82,8 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
 
 template <typename TStack>
 void StackInspector<TStack>::Init() {
-  fReportStack = false;
-  fStartTime = std::chrono::system_clock::now();
+  ReportStack_ = false;
+  StartTime_ = std::chrono::system_clock::now();
 }
 
 #include <corsika/cascade/testCascade.h>

Processes/StackInspector/StackInspector.h

@@ -39,12 +39,16 @@ namespace corsika::process {
       /**
        * To set a new E0, for example when a new shower event is started
        */
-      void SetE0(const corsika::units::si::HEPEnergyType vE0) { fE0 = vE0; }
+      void SetE0(const corsika::units::si::HEPEnergyType vE0) { E0_ = vE0; }
 
     private:
-      bool fReportStack;
-      corsika::units::si::HEPEnergyType fE0;
-      decltype(std::chrono::system_clock::now()) fStartTime;
+      bool ReportStack_;
+      corsika::units::si::HEPEnergyType E0_;
+      const corsika::units::si::HEPEnergyType dE_threshold_ = std::invoke([]() {
+        using namespace units::si;
+        return 1_eV;
+      });
+      decltype(std::chrono::system_clock::now()) StartTime_;
     };
 
   } // namespace stack_inspector