Merge branch '378-weird-longitudinalprofile' into 'master'

Resolve "Weird LongitudinalProfile"

Closes #378

See merge request AirShowerPhysics/corsika!321

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,

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};

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;

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_);

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;

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;
     }

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;

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 {

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

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,

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;
             }

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();
 

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; }

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

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)

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 {

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);