fixed problem with particles getting trapped in zero-length steps at their energy-threshold

Documentation/Examples/cascade_example.cc

@@ -142,7 +142,7 @@ int main() {
   process::particle_cut::ParticleCut cut(80_GeV, true, true);
 
   process::track_writer::TrackWriter trackWriter("tracks.dat");
-  process::energy_loss::EnergyLoss eLoss{showerAxis};
+  process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
 
   // assemble all processes into an ordered process list
   auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut

Documentation/Examples/em_shower.cc

@@ -159,8 +159,9 @@ int main(int argc, char** argv) {
   cut.ShowResults();
   em_continuous.ShowResults();
   observationLevel.ShowResults();
+  cout << "Cascade energy cut: " << EAS.GetEnergyCut()/1_GeV << " GeV" << endl;
   const HEPEnergyType Efinal =
-    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround();
+    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround() + EAS.GetEnergyCut();
   cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
   observationLevel.Reset();

Documentation/Examples/stopping_power.cc

@@ -49,7 +49,7 @@ int main() {
 
   environment::ShowerAxis showerAxis{injectionPos,
                                      Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
-  process::energy_loss::EnergyLoss eLoss{showerAxis};
+  process::energy_loss::EnergyLoss eLoss{showerAxis, 300_MeV};
 
   setup::Stack stack;
 

Documentation/Examples/vertical_EAS.cc

@@ -224,8 +224,9 @@ int main(int argc, char** argv) {
   cut.ShowResults();
   em_continuous.ShowResults();
   observationLevel.ShowResults();
+  cout << "Cascade energy cut: " << EAS.GetEnergyCut()/1_GeV << " GeV" << endl;
   const HEPEnergyType Efinal =
-    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround();
+    cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy() + em_continuous.GetEnergyLost() + observationLevel.GetEnergyGround() + EAS.GetEnergyCut();
   cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
   observationLevel.Reset();

Framework/Cascade/Cascade.h

@@ -83,8 +83,11 @@ namespace corsika::cascade {
         : fEnvironment(env)
         , fTracking(tr)
         , fProcessSequence(pl)
-        , fStack(stack) {}
+        , fStack(stack)
+        , energy_cut_(0 * corsika::units::si::electronvolt) {}
 
+    corsika::units::si::HEPEnergyType GetEnergyCut() const {  return energy_cut_; }
+    
     /**
      * set the nodes for all particles on the stack according to their numerical
      * position
@@ -209,10 +212,11 @@ namespace corsika::cascade {
 
       // apply all continuous processes on particle + track
       process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
-
+      
       if (status == process::EProcessReturn::eParticleAbsorbed) {
         std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
                   << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+	energy_cut_ += vParticle.GetEnergy();
         vParticle.Delete();
         return;
       }
@@ -316,7 +320,10 @@ namespace corsika::cascade {
     TProcessList& fProcessSequence;
     TStack& fStack;
     corsika::random::RNG& fRNG =
-        corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+    corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+    
+    corsika::units::si::HEPEnergyType energy_cut_;
+    
   }; // namespace corsika::cascade
 
 } // namespace corsika::cascade

Processes/EnergyLoss/EnergyLoss.cc

@@ -30,8 +30,10 @@ using SetupTrack = corsika::setup::Trajectory;
 
 using namespace corsika::process::energy_loss;
 
-EnergyLoss::EnergyLoss(environment::ShowerAxis const& shower_axis)
+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) {
@@ -176,9 +178,8 @@ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack co
        << " E=" << E / 1_GeV << "GeV,  Ekin=" << Ekin / 1_GeV << ", Enew=" << Enew / 1_GeV
        << "GeV" << endl;
   auto status = process::EProcessReturn::eOk;
-  if (-dE > Ekin) {
-    dE = -Ekin;
-    Enew = p.GetMass();
+  if (E<emCut_) {
+    Enew = emCut_;
     status = process::EProcessReturn::eParticleAbsorbed;
   }
   p.SetEnergy(Enew);
@@ -194,14 +195,24 @@ LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
   }
 
   auto constexpr dX = 1_g / square(1_cm);
-  auto const dE = -TotalEnergyLoss(vParticle, dX); // dE > 0
+  auto const dEdX = - TotalEnergyLoss(vParticle, dX) / dX; // dE > 0
   //~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
-  auto const maxLoss = 0.01 * vParticle.GetEnergy();
-  auto const maxGrammage = maxLoss / dE * dX;
+
+  // 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) *
-         1.0001; // to make sure particle gets absorbed when DoContinuous() is called
+                                                                         maxGrammage);
 }
 
 void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,

Processes/EnergyLoss/EnergyLoss.h

@@ -27,9 +27,9 @@ namespace corsika::process::energy_loss {
                              units::si::square(1e-2 * units::si::meter));
 
     void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
-
+    
   public:
-    EnergyLoss(environment::ShowerAxis const& showerAxis);
+    EnergyLoss(environment::ShowerAxis const& showerAxis, corsika::units::si::HEPEnergyType emCut);
 
     process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
                                          setup::Trajectory const&);
@@ -51,7 +51,10 @@ namespace corsika::process::energy_loss {
       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
   };
 

Processes/ParticleCut/ParticleCut.cc

@@ -26,9 +26,9 @@ namespace corsika::process {
       if (vP.GetPID() == particles::Code::Nucleus) {
         // calculate energy per nucleon
         auto const ElabNuc = energyLab / vP.GetNuclearA();
-        return (ElabNuc < eCut_);
+        return (ElabNuc <= eCut_);
       } else {
-        return (energyLab < eCut_);
+        return (energyLab <= eCut_);
       }
     }
 

Processes/Proposal/ContinuousProcess.cc

@@ -88,6 +88,7 @@ namespace corsika::process::proposal {
   EProcessReturn ContinuousProcess::DoContinuous(setup::Stack::ParticleType& vP,
                                                  setup::Trajectory const& vT) {
     if (!CanInteract(vP.GetPID())) return process::EProcessReturn::eOk;
+    if (vT.GetLength()==0_m) return process::EProcessReturn::eOk;
 
     // calculate passed grammage
     auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
@@ -106,9 +107,12 @@ namespace corsika::process::proposal {
 
     // Update the energy and absorbe the particle if it's below the energy
     // threshold, because it will no longer propagated.
-    vP.SetEnergy(final_energy);
-    if (vP.GetEnergy() <= emCut_)
+    if (final_energy <= emCut_) {
+      vP.SetEnergy(emCut_);
+      vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
       return process::EProcessReturn::eParticleAbsorbed;
+    }
+    vP.SetEnergy(final_energy);
     vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
     return process::EProcessReturn::eOk;
   }
@@ -121,8 +125,16 @@ namespace corsika::process::proposal {
 
     // Limit the step size of a conitnous loss. The maximal continuous loss seems to be a
     // hyper parameter which must be adjusted.
-    // in any case: never go below emCut_
-    auto energy_lim = std::max(0.9 * vP.GetEnergy(), emCut_);;
+    //
+    // 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.
+    // 
+    auto energy_lim = std::max(0.9 * vP.GetEnergy(), 0.99*emCut_);
 
     // solving the track integral for giving energy lim
     auto c = GetCalculator(vP, calc);

Processes/Proposal/Interaction.cc

@@ -91,6 +91,7 @@ namespace corsika::process::proposal {
             vec);
         auto sec_code = corsika::particles::ConvertFromPDG(
             static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+	std::cout << " proposal secondary: " << sec_code << " " << E/1_GeV << std::endl;
         vP.AddSecondary(make_tuple(sec_code, E, p, vP.GetPosition(), vP.GetTime()));
       }
     }