diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index ece23e9f534a0c5318209aef8ad0850b6bd1b940..92c52b3a2eb71a92be1d2a79d54f8968d34768f0 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/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
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 9bb1daf3f0c6eb7ff4eaeb86f5cdbfa5b9cd9618..32180ffc527a3bdcc04d02cb842b767e53909e1e 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/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();
diff --git a/Documentation/Examples/stopping_power.cc b/Documentation/Examples/stopping_power.cc
index a20e308984c286ac6cb7a4015f7d3da01edc4857..6054b38ac9c0da1ec5fe34143fbebd7683cff7aa 100644
--- a/Documentation/Examples/stopping_power.cc
+++ b/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;
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index b639f6374582a9b83109de6c1fff7fcbb0861623..8b3b8fbc4d5434919818206a3b32e63a122eba5a 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/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();
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 389ec0557e8bcc6da87fd57adafcf608d48a2a8e..21fad265637112c38920d0a2e91b05090a0afd36 100644
--- a/Framework/Cascade/Cascade.h
+++ b/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
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 0922e473edeaf99bdfce1944b620913a515d6b84..62790e9f50d38951555b3adf6bf16ecddaea342b 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/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,
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index 0059acb5d15b4f7c8418263e43fb519215b6fb7d..2e34650353491c6379f53ad8e288db3b41179607 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/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
};
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
index ab3cd652a3e193d3afd00dc994ad4e5e44fb806c..05e543674631439333b18325da12ccf4b8b0e539 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/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_);
}
}
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 0337bc8b0fc7c23cba77e4ecee186078d402ddea..1b17c5efca557713c08833dc973ba7b49ba2245a 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/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);
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index aeeaae4841c9fdda82b09439e8ce5f6d0c382d10..e5af126205a88e59606f92aae043eddb263e2016 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/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()));
}
}