diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 5f218b8780c063abfbbcbcd8339641ff5846f3f9..28306c590bb7da30d19d70f7ddfa19917ce525a0 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -121,11 +121,22 @@ namespace corsika {
ContinuousProcessStepLength const continuousMaxStep =
sequence_.getMaxStepLength(vParticle, step);
LengthType const continuous_max_dist = continuousMaxStep;
- ContinuousProcessIndex const limitingId = continuousMaxStep;
// take minimum of geometry, interaction, decay for next step
- auto const min_distance =
- std::min({distance_interact, distance_decay, continuous_max_dist, geomMaxLength});
+ LengthType const min_discrete = std::min(distance_interact, distance_decay);
+ LengthType const min_non_continuous = std::min(min_discrete, geomMaxLength);
+ LengthType const min_distance = std::min(min_non_continuous, continuous_max_dist);
+
+ // inform ContinuousProcesses (if applicable) that it is responsible for step-limit
+ // this would become simpler if we follow the idea of Max to enumerate ALL types of
+ // processes. Then non-contunuous are included and no further logic is needed to
+ // distinguish between continuous and non-continuous limit.
+ ContinuousProcessIndex limitingId;
+ bool const isContinuous = continuous_max_dist < min_non_continuous;
+ if (isContinuous) {
+ limitingId =
+ continuousMaxStep; // the current step IS limited by a known continuous process
+ }
CORSIKA_LOG_DEBUG(
"transport particle by : {} m "
@@ -136,10 +147,12 @@ namespace corsika {
min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
distance_interact / 1_m, continuous_max_dist / 1_m);
- // here the particle is actually moved along the trajectory to new position:
+ // move particle along the trajectory to new position
+ // also update momentum/direction/time
step.setLength(min_distance);
vParticle.setPosition(step.getPosition(1));
- // assumption: tracking does not change absolute momentum:
+ // assumption: tracking does not change absolute momentum (continuous physics can and
+ // will):
vParticle.setMomentum(step.getDirection(1) * vParticle.getMomentum().getNorm());
vParticle.setTime(vParticle.getTime() + step.getDuration());
@@ -148,54 +161,52 @@ namespace corsika {
ProcessReturn::ParticleAbsorbed) {
CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.getPID(), vParticle.getEnergy() / 1_GeV);
- if (!vParticle.isErased()) vParticle.erase();
+ if (vParticle.isErased()) {
+ CORSIKA_LOG_WARN(
+ "Particle marked as Absorbed in doContinuous, but prematurely erased. This "
+ "may be bug. Check.");
+ } else {
+ vParticle.erase();
+ }
return;
}
+ if (isContinuous) {
+ return; // there is nothing further, step is finished
+ }
- CORSIKA_LOG_DEBUG("sth. happening before geometric limit ? {}",
+ CORSIKA_LOG_DEBUG("discrete process before geometric limit ? {}",
((min_distance < geomMaxLength) ? "yes" : "no"));
- if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
+ if (geomMaxLength < min_discrete) {
+ // geometric / tracking limit
- // check whether decay or interaction limits this step the
- // outcome of decay or interaction MAY be a) new particles in
- // secondaries, b) the projectile particle deleted (or
- // changed)
+ if (nextVol != currentLogicalNode) {
+ // boundary crossing, step is limited by volume boundary
- TStackView secondaries(vParticle);
+ CORSIKA_LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
- if (min_distance < continuous_max_dist) {
+ if (nextVol == environment_.getUniverse().get()) {
+ CORSIKA_LOG_DEBUG(
+ "particle left physics world, is now in unknown space -> delete");
+ vParticle.erase();
+ }
+ vParticle.setNode(nextVol);
/*
- Create SecondaryView object on Stack. The data container
- remains untouched and identical, and 'projectil' is identical
- to 'vParticle' above this line. However,
- projectil.AddSecondaries populate the SecondaryView, which can
- then be used afterwards for further processing. Thus: it is
- important to use projectle/view (and not vParticle) for Interaction,
- and Decay!
- */
+ doBoundary may delete the particle (or not)
- [[maybe_unused]] auto projectile = secondaries.getProjectile();
-
- if (distance_interact < distance_decay) {
- interaction(secondaries);
- } else {
- decay(secondaries);
- // make sure particle actually did decay if it should have done so
- if (secondaries.getSize() == 1 &&
- projectile.getPID() == secondaries.getNextParticle().getPID())
- throw std::runtime_error(fmt::format("Particle {} decays into itself!",
- get_name(projectile.getPID())));
- }
+ caveat: any changes to vParticle, or even the production
+ of new secondaries is currently not passed to ParticleCut,
+ thus, particles outside the desired phase space may be produced.
- sequence_.doSecondaries(secondaries);
- vParticle.erase();
- } else { // step-length limitation within volume
+ \todo: this must be fixed.
+ */
- CORSIKA_LOG_DEBUG("step-length limitation");
- // no further physics happens here. just proceed to next step.
+ sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ return; // step finished
}
+ CORSIKA_LOG_DEBUG("step limit reached (e.g. deflection). nothing further happens.");
+
[[maybe_unused]] auto const assertion = [&] {
auto const* numericalNodeAfterStep =
environment_.getUniverse()->getContainingNode(vParticle.getPosition());
@@ -208,31 +219,43 @@ namespace corsika {
assert(assertion()); // numerical and logical nodes should match, since
// we did not cross any volume boundary
- } else { // boundary crossing, step is limited by volume boundary
-
- if (nextVol != currentLogicalNode) {
-
- CORSIKA_LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
-
- if (nextVol == environment_.getUniverse().get()) {
- CORSIKA_LOG_DEBUG(
- "particle left physics world, is now in unknown space -> delete");
- vParticle.erase();
- }
- vParticle.setNode(nextVol);
- /*
- doBoundary may delete the particle (or not)
+ // step length limit
+ return;
+ }
- caveat: any changes to vParticle, or even the production
- of new secondaries is currently not passed to ParticleCut,
- thus, particles outside the desired phase space may be produced.
+ // interaction or decay to happen in this step
+ // the outcome of decay or interaction MAY be a) new particles in
+ // secondaries, b) the projectile particle deleted (or
+ // changed)
- \todo: this must be fixed.
- */
+ TStackView secondaries(vParticle);
- sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
- }
+ /*
+ Create SecondaryView object on Stack. The data container
+ remains untouched and identical, and 'projectil' is identical
+ to 'vParticle' above this line. However,
+ projectil.AddSecondaries populate the SecondaryView, which can
+ then be used afterwards for further processing. Thus: it is
+ important to use projectle/view (and not vParticle) for Interaction,
+ and Decay!
+ */
+
+ [[maybe_unused]] auto projectile = secondaries.getProjectile();
+
+ if (distance_interact < distance_decay) {
+ interaction(secondaries);
+ } else {
+ decay(secondaries);
+ // make sure particle actually did decay if it should have done so
+ // \todo this should go to a validation code and not be included here
+ if (secondaries.getSize() == 1 &&
+ projectile.getPID() == secondaries.getNextParticle().getPID())
+ throw std::runtime_error(fmt::format("Particle {} decays into itself!",
+ get_name(projectile.getPID())));
}
+
+ sequence_.doSecondaries(secondaries);
+ vParticle.erase();
}
template <typename TTracking, typename TProcessList, typename TStack,
diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 7d0aa253b838adf8a06d5f0cdfbb71022d896695..47504f4ea5f0dd02672418e260eb61dd6971c5c3 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -53,15 +53,15 @@ namespace corsika {
return particle::detail::isHadron[static_cast<CodeIntType>(p)];
}
- inline bool constexpr is_em(Code c) {
+ inline bool constexpr is_em(Code const c) {
return c == Code::Electron || c == Code::Positron || c == Code::Gamma;
}
- inline bool constexpr is_muon(Code c) {
+ inline bool constexpr is_muon(Code const c) {
return c == Code::MuPlus || c == Code::MuMinus;
}
- inline bool constexpr is_neutrino(Code c) {
+ inline bool constexpr is_neutrino(Code const c) {
return c == Code::NuE || c == Code::NuMu || c == Code::NuTau || c == Code::NuEBar ||
c == Code::NuMuBar || c == Code::NuTauBar;
}
@@ -88,7 +88,7 @@ namespace corsika {
return stream << get_name(code);
}
- inline Code convert_from_PDG(PDGCode p) {
+ inline Code convert_from_PDG(PDGCode const p) {
static_assert(particle::detail::conversionArray.size() % 2 == 1);
// this will fail, for the strange case where the maxPDG is negative...
int constexpr maxPDG{(particle::detail::conversionArray.size() - 1) >> 1};
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 7a0d6409feadf182a1e8402c581d118009a75b2f..b93331d05079962f5e710ac8de952935eb289517 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -67,12 +67,10 @@ namespace corsika {
}
if constexpr (t2ProcSeq) {
- if (!isAbsorbed(ret)) { ret |= B_.doContinuous(particle, vT, limitId); }
+ ret |= B_.doContinuous(particle, vT, limitId);
} else if constexpr (is_continuous_process_v<process2_type>) {
- if (!isAbsorbed(ret)) {
- ret |= B_.doContinuous(particle, vT,
- limitId == ContinuousProcessIndex(IndexProcess2));
- }
+ ret |=
+ B_.doContinuous(particle, vT, limitId == ContinuousProcessIndex(IndexProcess2));
}
return ret;
diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl
index d332913b5d995c83135fdbcb2adf4d383cf13456..48deafe6e79a474f06ecf664d9b70640b8e7756b 100644
--- a/corsika/detail/modules/LongitudinalProfile.inl
+++ b/corsika/detail/modules/LongitudinalProfile.inl
@@ -35,9 +35,11 @@ namespace corsika {
GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
- CORSIKA_LOG_INFO(
- "pos1={} m, pos2={}, X={} g/cm2", vTrack.getPosition(0).getCoordinates() / 1_m,
- vTrack.getPosition(1).getCoordinates() / 1_m, grammageStart / 1_g * square(1_cm));
+ CORSIKA_LOG_TRACE("pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2",
+ vTrack.getPosition(0).getCoordinates() / 1_m,
+ vTrack.getPosition(1).getCoordinates() / 1_m,
+ grammageStart / 1_g * square(1_cm),
+ grammageEnd / 1_g * square(1_cm));
// Note: particle may go also "upward", thus, grammageEnd<grammageStart
const int binStart = std::ceil(grammageStart / dX_);
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 33cb556ec0021a6fedf2d254b2d35ff559df18bb..c4da50e940f88e3b8a3d3e9f36ac0018eeb09466 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -48,7 +48,6 @@ namespace corsika {
if (deleteOnHit_) {
count_ground_++;
energy_ground_ += energy;
- particle.erase();
return ProcessReturn::ParticleAbsorbed;
} else {
return ProcessReturn::Ok;
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index f4d75667518b743a421b429a6afb8170c71e9869..ffbbbed8ddd51cb99ffe45f92217ed426fe2ee91 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -153,7 +153,6 @@ namespace corsika {
CORSIKA_LOG_TRACE("ParticleCut::DoContinuous");
if (checkCutParticle(particle)) {
CORSIKA_LOG_TRACE("removing during continuous");
- particle.erase();
// signal to upstream code that this particle was deleted
return ProcessReturn::ParticleAbsorbed;
}
diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl
index 10b0f1582a31fcc01ada4b638f428917fc369747..caae6c59c8cb43b01cd1e9789366b51db2f9ec14 100644
--- a/corsika/detail/modules/StackInspector.inl
+++ b/corsika/detail/modules/StackInspector.inl
@@ -24,8 +24,8 @@
namespace corsika {
template <typename TStack>
- inline StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
- const HEPEnergyType vE0)
+ inline StackInspector<TStack>::StackInspector(int const vNStep, bool const vReportStack,
+ HEPEnergyType const vE0)
: StackProcess<StackInspector<TStack>>(vNStep)
, ReportStack_(vReportStack)
, E0_(vE0)
@@ -35,12 +35,12 @@ namespace corsika {
inline StackInspector<TStack>::~StackInspector() {}
template <typename TStack>
- inline void StackInspector<TStack>::doStack(const TStack& vS) {
+ inline void StackInspector<TStack>::doStack(TStack const& vS) {
[[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV;
- for (const auto& iterP : vS) {
+ for (auto const& iterP : vS) {
HEPEnergyType E = iterP.getEnergy();
Etot += E;
if (ReportStack_) {
@@ -60,7 +60,7 @@ namespace corsika {
}
auto const now = std::chrono::system_clock::now();
- const std::chrono::duration<double> elapsed_seconds = now - StartTime_;
+ std::chrono::duration<double> const elapsed_seconds = now - StartTime_;
std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
auto const dE = E0_ - Etot;
if (dE < dE_threshold_) return;
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index dff3b2bfb7dddd48309fe91a0b4ca94c0f1fad18..73fbf351707c35c894dd8cf48e0d250589909527 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -64,18 +64,18 @@ namespace corsika {
double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma);
// (gamma_2-1)/gamma_2 = (1-1/gamma2);
double constexpr c2 = 1; // HEP convention here c=c2=1
- CORSIKA_LOG_DEBUG("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
+ CORSIKA_LOG_TRACE("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
[[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
HEPMassType const Wmax =
2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
// approx, but <<1% HEPMassType const Wmax = 2*me*c2*beta2*gamma2; for HEAVY
// PARTICLES Wmax ~ 2me v2 for non-relativistic particles
- CORSIKA_LOG_DEBUG("BetheBloch Wmax={}", Wmax);
+ CORSIKA_LOG_TRACE("BetheBloch Wmax={}", Wmax);
// Sternheimer parameterization, density corrections towards high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
// MISSING
- CORSIKA_LOG_DEBUG("BetheBloch p.getMomentum().getNorm()/m{}=",
+ CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m{}=",
p.getMomentum().getNorm() / m);
double const x = log10(p.getMomentum().getNorm() / m);
double delta = 0;
@@ -86,7 +86,7 @@ namespace corsika {
} else if (x < x0) { // and IF conductor (otherwise, this is 0)
delta = delta0 * pow(100, 2 * (x - x0));
}
- CORSIKA_LOG_DEBUG("BetheBloch delta={}", delta);
+ CORSIKA_LOG_TRACE("BetheBloch delta={}", delta);
// with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
@@ -142,17 +142,27 @@ namespace corsika {
inline ProcessReturn BetheBlochPDG::doContinuous(setup::Stack::particle_type& p,
setup::Trajectory const& t,
bool const) {
+
+ // if this step was limiting the CORSIKA stepping, the particle is lost
+ /* see Issue https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/389
+ if (limitStep) {
+ fillProfile(t, p.getEnergy());
+ p.setEnergy(p.getMass());
+ return ProcessReturn::ParticleAbsorbed;
+ }
+ */
+
if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
GrammageType const dX =
p.getNode()->getModelProperties().getIntegratedGrammage(t, t.getLength());
- CORSIKA_LOG_DEBUG("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
+ CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
p.getChargeNumber(), dX / 1_g * square(1_cm));
HEPEnergyType dE = getTotalEnergyLoss(p, dX);
auto E = p.getEnergy();
[[maybe_unused]] const auto Ekin = E - p.getMass();
auto Enew = E + dE;
- CORSIKA_LOG_DEBUG("EnergyLoss dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV",
+ CORSIKA_LOG_TRACE("EnergyLoss dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV",
dE / 1_MeV, E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV);
p.setEnergy(Enew);
updateMomentum(p, Enew);
@@ -168,22 +178,17 @@ namespace corsika {
}
auto constexpr dX = 1_g / square(1_cm);
- auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX; // dE > 0
- //~ auto const Ekin = vParticle.getEnergy() - vParticle.getMass();
-
- auto const emCut = get_energy_threshold(vParticle.getPID());
- // 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 3D position where it reaches emCut is not very
- // 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;
+ auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX;
+ auto const energy = vParticle.getEnergy();
+ auto const energy_lim = std::max(
+ energy * 0.9, // either 10% relative loss max., or
+ get_energy_threshold(vParticle.getPID()) // energy thresholds globally defined for
+ // individual particles
+ *
+ 0.99 // need to go 1% below global e-cut to assure removal in ParticleCut. The
+ // 1% does not matter since at cut-time the entire energy is removed.
+ );
+ auto const maxGrammage = (vParticle.getEnergy() - energy_lim) / dEdX;
return vParticle.getNode()->getModelProperties().getArclengthFromGrammage(
vTrack, maxGrammage);
@@ -219,7 +224,7 @@ namespace corsika {
if (binEnd < 0) binEnd = 0;
if (binEnd > maxBin) binEnd = maxBin;
- CORSIKA_LOG_DEBUG("energy deposit of -dE={} between {} and {}", -dE, grammageStart,
+ CORSIKA_LOG_TRACE("energy deposit of -dE={} between {} and {}", -dE, grammageStart,
grammageEnd);
auto energyCount = HEPEnergyType::zero();
@@ -230,7 +235,7 @@ namespace corsika {
profile_[bin] += increment;
energyCount += increment;
- CORSIKA_LOG_DEBUG("filling bin {} with weight {} : {} ", bin, weight, increment);
+ CORSIKA_LOG_TRACE("filling bin {} with weight {} : {} ", bin, weight, increment);
};
// fill longitudinal profile
@@ -242,7 +247,7 @@ namespace corsika {
for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
}
- CORSIKA_LOG_DEBUG("total energy added to histogram: {} ", energyCount);
+ CORSIKA_LOG_TRACE("total energy added to histogram: {} ", energyCount);
}
inline void BetheBlochPDG::printProfile() const {
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index b4b0bf9638fbc6dd00d4d258d465cc6da4e65ca1..0e230a762aae1ca039e06570ea0bb2f87fa54a50 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -124,23 +124,20 @@ namespace corsika::proposal {
// Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
// hyper parameter which must be adjusted.
//
- auto const emCut = get_energy_threshold(
- code); //! energy thresholds globally defined for individual particles
-
- // 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);
+ auto const energy = vP.getEnergy();
+ auto const energy_lim = std::max(
+ energy * 0.9, // either 10% relative loss max., or
+ get_energy_threshold(
+ code) // energy thresholds globally defined for individual particles
+ *
+ 0.99 // need to go 1% below global e-cut to assure removal in ParticleCut. The
+ // 1% does not matter since at cut-time the entire energy is removed.
+ );
// solving the track integral for giving energy lim
auto c = getCalculator(vP, calc);
auto grammage = get<eDISPLACEMENT>(c->second)->SolveTrackIntegral(
- vP.getEnergy() / 1_MeV, energy_lim / 1_MeV) *
+ energy / 1_MeV, energy_lim / 1_MeV) *
1_g / square(1_cm);
// return it in distance aequivalent
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index d41c149024ff2d932fb259a909e07e6ac2ac2e63..d58ac21d851ff8b19be258e09e9a63a500784b1a 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -82,7 +82,7 @@ namespace corsika::pythia8 {
inline void Decay::setHandleDecay(Code const vParticleCode) {
handleAllDecays_ = false;
- CORSIKA_LOG_INFO("Pythia::Decay: set to handle decay of {} ", vParticleCode);
+ CORSIKA_LOG_DEBUG("Pythia::Decay: set to handle decay of {} ", vParticleCode);
if (Decay::canHandleDecay(vParticleCode))
handledDecays_.insert(vParticleCode);
else
@@ -106,12 +106,12 @@ namespace corsika::pythia8 {
}
inline void Decay::setUnstable(Code const pCode) {
- CORSIKA_LOG_INFO("Pythia::Decay: setting {} unstable..", pCode);
+ CORSIKA_LOG_DEBUG("Pythia::Decay: setting {} unstable..", pCode);
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
}
inline void Decay::setStable(Code const pCode) {
- CORSIKA_LOG_INFO("Pythia::Decay: setting {} stable..", pCode);
+ CORSIKA_LOG_DEBUG("Pythia::Decay: setting {} stable..", pCode);
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
@@ -122,8 +122,8 @@ namespace corsika::pythia8 {
inline bool Decay::canDecay(Code const pCode) {
bool const ans =
Pythia8::Pythia::particleData.canDecay(static_cast<int>(get_PDG(pCode)));
- CORSIKA_LOG_INFO("Pythia::Decay: checking if particle: {} can decay in PYTHIA? {} ",
- pCode, ans);
+ CORSIKA_LOG_DEBUG("Pythia::Decay: checking if particle: {} can decay in PYTHIA? {} ",
+ pCode, ans);
return ans;
}
@@ -153,13 +153,13 @@ namespace corsika::pythia8 {
TimeType const t0 = get_lifetime(pid);
auto const lifetime = gamma * t0;
- CORSIKA_LOG_INFO("Pythia::Decay: code: {}", particle.getPID());
- CORSIKA_LOG_INFO("Pythia::Decay: MinStep: t0: {}", t0);
- CORSIKA_LOG_INFO("Pythia::Decay: MinStep: energy: {} GeV", E / 1_GeV);
- CORSIKA_LOG_INFO("Pythia::Decay: momentum: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_INFO("Pythia::Decay: MinStep: gamma: {}", gamma);
- CORSIKA_LOG_INFO("Pythia::Decay: MinStep: tau: {} ", lifetime);
+ CORSIKA_LOG_TRACE("Pythia::Decay: code: {}", particle.getPID());
+ CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: t0: {}", t0);
+ CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: energy: {} GeV", E / 1_GeV);
+ CORSIKA_LOG_TRACE("Pythia::Decay: momentum: {} GeV",
+ particle.getMomentum().getComponents() / 1_GeV);
+ CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: gamma: {}", gamma);
+ CORSIKA_LOG_TRACE("Pythia::Decay: MinStep: tau: {} ", lifetime);
return lifetime;
} else
@@ -209,7 +209,7 @@ namespace corsika::pythia8 {
if (!Pythia8::Pythia::next())
throw std::runtime_error("Pythia::Decay: decay failed!");
else
- CORSIKA_LOG_INFO("Pythia::Decay: particles after decay: {} ", event.size());
+ CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
if (print_listing_) {
// list final state
@@ -227,9 +227,10 @@ namespace corsika::pythia8 {
FourVector const fourMomRest{Erest, pRest};
auto const fourMomLab = boost.fromCoM(fourMomRest);
- CORSIKA_LOG_INFO("particle: id={} momentum={} energy={} ", pyId,
- fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
- fourMomLab.getTimeLikeComponent());
+ CORSIKA_LOG_TRACE(
+ "particle: id={} momentum={} energy={} ", pyId,
+ fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
+ fourMomLab.getTimeLikeComponent());
view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(),
fourMomLab.getSpaceLikeComponents(), decayPoint,
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 415c1ada386c5840dd2c9ded27a28fdf7c3f2236..5770fe439dbc309588e990da8669c9f85ceb619c 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -176,8 +176,8 @@ namespace corsika {
CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution
- if (dist < 0.0001_m) { // security margin to assure transfer to next
- // logical volume
+ if (dist < -0.0001_m) { // security margin to assure transfer to next
+ // logical volume
if (first_entry == 0) {
d_enter = dist;
} else {
@@ -200,7 +200,7 @@ namespace corsika {
// both physical solutions (entry, exit) must be positive, and as small as
// possible
- if (dist < 0.0001_m) { // need small numerical margin, to assure transport
+ if (dist < -0.0001_m) { // need small numerical margin, to assure transport
// into next logical volume
continue;
}
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index dafc1f02437fcecc0c926b9b9b94a283997c03b6..e133d68679c27bc597050a46d84c5bba2bfed70e 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -48,9 +48,9 @@ namespace corsika {
using PDGCodeType = std::underlying_type<PDGCode>::type;
// forward declarations to be used in GeneratedParticleProperties
- int16_t constexpr get_charge_number(Code); //!< electric charge in units of e
- ElectricChargeType constexpr get_charge(Code); //!< electric charge
- HEPMassType constexpr get_mass(Code); //!< mass
+ int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e
+ ElectricChargeType constexpr get_charge(Code const); //!< electric charge
+ HEPMassType constexpr get_mass(Code const); //!< mass
HEPEnergyType constexpr get_energy_threshold(
Code const); //!< get energy threshold below which the particle is discarded, by
//!< default set to particle mass
@@ -67,24 +67,26 @@ namespace corsika {
}
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
- PDGCode constexpr get_PDG(Code);
- std::string_view constexpr get_name(Code); //!< name of the particle as string
- TimeType constexpr get_lifetime(Code); //!< lifetime
+ PDGCode constexpr get_PDG(Code const);
+ std::string_view constexpr get_name(Code const); //!< name of the particle as string
+ TimeType constexpr get_lifetime(Code const); //!< lifetime
//! true iff the particle is a hard-coded nucleus or Code::Nucleus
- bool constexpr is_nucleus(Code);
- bool constexpr is_hadron(Code); //!< true iff particle is hadron
- bool constexpr is_em(Code); //!< true iff particle is electron, positron or gamma
- bool constexpr is_muon(Code); //!< true iff particle is mu+ or mu-
- bool constexpr is_neutrino(Code); //!< true iff particle is (anti-) neutrino
- int constexpr get_nucleus_A(Code); //!< returns A for hard-coded nucleus, otherwise 0
- int constexpr get_nucleus_Z(Code); //!< returns Z for hard-coded nucleus, otherwise 0
+ bool constexpr is_nucleus(Code const);
+ bool constexpr is_hadron(Code const); //!< true iff particle is hadron
+ bool constexpr is_em(Code const); //!< true iff particle is electron, positron or gamma
+ bool constexpr is_muon(Code const); //!< true iff particle is mu+ or mu-
+ bool constexpr is_neutrino(Code const); //!< true iff particle is (anti-) neutrino
+ int constexpr get_nucleus_A(
+ Code const); //!< returns A for hard-coded nucleus, otherwise 0
+ int constexpr get_nucleus_Z(
+ Code const); //!< returns Z for hard-coded nucleus, otherwise 0
//! returns mass of (A,Z) nucleus, disregarding binding energy
HEPMassType get_nucleus_mass(unsigned int const, unsigned int const);
//! convert PDG code to CORSIKA 8 internal code
- Code convert_from_PDG(PDGCode);
+ Code convert_from_PDG(PDGCode const);
std::initializer_list<Code> constexpr get_all_particles();
diff --git a/corsika/framework/process/ContinuousProcessIndex.hpp b/corsika/framework/process/ContinuousProcessIndex.hpp
index bc2617d2bb6ed5366431d1369480ef9abf593426..9d6e84d0f39baaa9a1e0ad7dbf369aba6a967fb3 100644
--- a/corsika/framework/process/ContinuousProcessIndex.hpp
+++ b/corsika/framework/process/ContinuousProcessIndex.hpp
@@ -18,6 +18,8 @@ namespace corsika {
class ContinuousProcessIndex {
public:
+ ContinuousProcessIndex()
+ : id_(-1) {} // default
ContinuousProcessIndex(int const id)
: id_(id) {}
void setIndex(int const id) { id_ = id; }
diff --git a/corsika/setup/SetupTrajectory.hpp b/corsika/setup/SetupTrajectory.hpp
index 95a91115b0aa422c860448fa8057d70c35f95456..c40351d26f1c7bdef9d9b8ac18d3653caaf1927d 100644
--- a/corsika/setup/SetupTrajectory.hpp
+++ b/corsika/setup/SetupTrajectory.hpp
@@ -38,15 +38,15 @@ namespace corsika::setup {
The default tracking algorithm.
*/
- typedef corsika::tracking_leapfrog_curved::Tracking Tracking;
+ // typedef corsika::tracking_leapfrog_curved::Tracking Tracking;
// typedef corsika::tracking_leapfrog_straight::Tracking Tracking;
- // typedef corsika::tracking_line::Tracking Tracking;
+ typedef corsika::tracking_line::Tracking Tracking;
/**
The default trajectory.
*/
/// definition of Trajectory base class, to be used in tracking and cascades
- // typedef StraightTrajectory Trajectory;
- typedef corsika::LeapFrogTrajectory Trajectory;
+ typedef StraightTrajectory Trajectory;
+ // typedef corsika::LeapFrogTrajectory Trajectory;
} // namespace corsika::setup
diff --git a/do-clang-format.py b/do-clang-format.py
index 3bff0e619a14137fb9441b394a985875d0bbe7b2..d6ea597535e69da32de80ff31962ba998523305d 100755
--- a/do-clang-format.py
+++ b/do-clang-format.py
@@ -10,6 +10,14 @@ import os
import sys
import re
+do_progress = False
+try:
+ from progress.bar import ChargingBar
+ do_progress = True
+except ImportError as e:
+ do_progress = False
+ # no progress bar
+
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--apply', action="store_true",
help="Apply clang-format to files which need changes.")
@@ -85,17 +93,26 @@ version = subp.check_output(cmd.split() + ["--version"]).decode("utf-8")
print (version)
print ("Note: the clang-format version has an impact on the result. Make sure you are consistent with current CI. Consider \'--docker\' option.")
+
+bar = None
+if do_progress:
+ bar = ChargingBar('Processing', max=len(filelist))
+
if args.apply:
for filename in filelist:
+ if bar: bar.next()
subp.check_call(cmd.split() + ["-i", filename])
+ if bar: bar.finish()
else:
# only print files which need formatting
files_need_formatting = 0
for filename in filelist:
+ if bar: bar.next()
a = open(filename, "rb").read()
b = subp.check_output(cmd.split() + [filename])
if a != b:
files_need_formatting += 1
print(filename)
+ if bar: bar.finish()
sys.exit(1 if files_need_formatting > 0 else 0)
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 641fb95f85ce3e61e1779b21037ffc0ae6b3003c..bdad125b05e7557a4db27da54ee2eecae27d7592 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -6,8 +6,6 @@
* the license.
*/
-#define DEBUG 1
-
/* clang-format off */
// InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have
@@ -96,13 +94,15 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
- logging::set_level(logging::level::trace);
+ logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("vertical_EAS");
if (argc < 4) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl;
- std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
+ std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed] \n"
+ " if A=0, Z is interpreted as PDG code \n"
+ " if no seed is given, a random seed is chosen \n"
+ << std::endl;
return 1;
}
feenableexcept(FE_INVALID);
@@ -131,7 +131,7 @@ int main(int argc, char** argv) {
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
- builder.addLinearLayer(1e9_cm, 112.8_km+constants::EarthRadius::Mean);
+ builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
builder.assemble(env);
CORSIKA_LOG_DEBUG(
@@ -153,11 +153,20 @@ int main(int argc, char** argv) {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
- const Code beamCode = Code::Nucleus;
unsigned short const A = std::stoi(std::string(argv[1]));
- unsigned short Z = std::stoi(std::string(argv[2]));
- auto const mass = get_nucleus_mass(A, Z);
- const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
+ Code beamCode;
+ HEPEnergyType mass;
+ unsigned short Z = 0;
+ if (A > 0) {
+ beamCode = Code::Nucleus;
+ Z = std::stoi(std::string(argv[2]));
+ mass = get_nucleus_mass(A, Z);
+ } else {
+ int pdg = std::stoi(std::string(argv[2]));
+ beamCode = convert_from_PDG(PDGCode(pdg));
+ mass = get_mass(beamCode);
+ }
+ HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
@@ -187,17 +196,21 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
- if (A != 1) {
+ if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
- if (Z == 1) {
- stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
- } else if (Z == 0) {
- stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
+ if (A == 1) {
+ if (Z == 1) {
+ stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
+ } else if (Z == 0) {
+ stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
+ } else {
+ std::cerr << "illegal parameters" << std::endl;
+ return EXIT_FAILURE;
+ }
} else {
- std::cerr << "illegal parameters" << std::endl;
- return EXIT_FAILURE;
+ stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
}
}
@@ -210,8 +223,7 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
- corsika::qgsjetII::Interaction qgsjet;
-
+ // corsika::qgsjetII::Interaction qgsjet;
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
@@ -259,7 +271,7 @@ int main(int argc, char** argv) {
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
- StackInspector<setup::Stack> stackInspect(1000, false, E0);
+ StackInspector<setup::Stack> stackInspect(50000, false, E0);
// assemble all processes into an ordered process list
struct EnergySwitch {
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index ef12f58e8489341ad64adcfd0d67194c064514cb..1bc67084bc633c62539ec8a644eb42528a8e4106 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -23,7 +23,7 @@
using namespace corsika;
-TEST_CASE("ParticleCut", "[processes]") {
+TEST_CASE("ParticleCut", "processes") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
@@ -222,7 +222,10 @@ TEST_CASE("ParticleCut", "[processes]") {
for (auto proType : particleList) {
auto particle = stack.addParticle(std::make_tuple(
proType, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
- cut.doContinuous(particle, track);
+
+ if (cut.doContinuous(particle, track) == ProcessReturn::ParticleAbsorbed) {
+ particle.erase();
+ }
}
CHECK(stack.getEntries() == 9);