add further doc and some minor fixes in particle cut

Documentation/Examples/vertical_EAS.cc

@@ -195,7 +195,6 @@ int main(int argc, char** argv) {
 
   process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
 
-  process::energy_loss::EnergyLoss eLoss(showerAxis);
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
 
   Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
@@ -205,10 +204,6 @@ int main(int argc, char** argv) {
   process::UrQMD::UrQMD urqmd;
   process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
 
-  /* process::conex_source_cut::CONEXSourceCut conexSource( */
-  /*     center, showerAxis, t, injectionHeight, E0, */
-  /*     particles::GetPDG(particles::Code::Proton)); */
-
   // assemble all processes into an ordered process list
 
   auto sibyllSequence = sibyllNucCounted << sibyllCounted;
@@ -216,8 +211,6 @@ int main(int argc, char** argv) {
                                                        55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
 
-  //auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource
-  //                            << longprof << eLoss << cut << observationLevel;
   auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << cut << em_continuous
 				<< longprof << observationLevel;
 
@@ -231,16 +224,12 @@ int main(int argc, char** argv) {
 
   EAS.Run();
 
-  eLoss.PrintProfile(); // print longitudinal profile
-  /* conexSource.SolveCE(); */
 
   cut.ShowResults();
   const HEPEnergyType Efinal =
       cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
   cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-  cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
-       << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
 
   auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
     urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();

Processes/ParticleCut/ParticleCut.cc

@@ -84,7 +84,7 @@ namespace corsika::process {
 
     ParticleCut::ParticleCut(const units::si::HEPEnergyType eCut, bool discardEm, bool discardInv)
         : eCut_(eCut)
-        , discardEm_(cutEm)
+        , discardEm_(discardEm)
         , discardInv_(discardInv) {
 
       emEnergy_ = 0_GeV;

Processes/ParticleCut/ParticleCut.h

@@ -18,8 +18,8 @@ namespace corsika::process {
     class ParticleCut : public process::SecondariesProcess<ParticleCut> {
 
       units::si::HEPEnergyType const eCut_;
-      bool cutEm_;
-      bool cutInv_;
+      bool discardEm_;
+      bool discardInv_;
 
       units::si::HEPEnergyType energy_ = 0 * units::si::electronvolt;
       units::si::HEPEnergyType emEnergy_ = 0 * units::si::electronvolt;

Processes/Proposal/ContinuousProcess.cc

@@ -23,10 +23,19 @@ namespace corsika::process::proposal {
 
   void ContinuousProcess::BuildCalculator(particles::Code code,
                                           environment::NuclearComposition const& comp) {
+    // search crosssection builder for given particle
     auto p_cross = cross.find(code);
     if (p_cross == cross.end())
       throw std::runtime_error("PROPOSAL could not find corresponding builder");
+
+    // interpolate the crosssection for given media and energy cut. These may
+    // take some minutes if you have to build the tables and cannot read the
+    // from disk
     auto c = p_cross->second(media.at(&comp), cut);
+
+    // Build displacement integral and scattering object and interpolate them too and
+    // saved in the calc map by a key build out of a hash of composed of the component and
+    // particle code.
     auto disp = PROPOSAL::make_displacement(c, true);
     auto scatter = PROPOSAL::make_scattering("highland", particle[code], media.at(&comp));
     calc[std::make_pair(&comp, code)] =
@@ -38,31 +47,35 @@ namespace corsika::process::proposal {
                                        particle_cut::ParticleCut& _cut)
       : ProposalProcessBase(_env, _cut) {}
 
-  template <>
-  HEPEnergyType ContinuousProcess::TotalEnergyLoss(setup::Stack::ParticleType const& vP,
-                                                   GrammageType const& vDX) {
-    auto c = GetCalculator(vP, calc);
-    return vP.GetEnergy() - get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
-                                vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm) *
-                                1_MeV;
-  }
-
   template <>
   void ContinuousProcess::Scatter(setup::Stack::ParticleType& vP,
                                   HEPEnergyType const& loss,
                                   GrammageType const& grammage) {
+    // Get or build corresponding calculators
     auto c = GetCalculator(vP, calc);
+
+    // Cast corsika vector to proposal vector
     auto d = vP.GetDirection().GetComponents();
     auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
                                         d.GetZ().magnitude());
-    auto E_f = vP.GetEnergy() - loss; // final energy
+
+    auto E_f = vP.GetEnergy() - loss;
+
+    // draw random numbers required for scattering process
     std::uniform_real_distribution<double> distr(0., 1.);
     auto rnd = array<double, 4>();
     for (auto& it : rnd) it = distr(fRNG);
+
+    // calculate deflection based on particle energy, loss
     auto [mean_dir, final_dir] = get<SCATTERING>(c->second)->Scatter(
         grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV, direction,
         rnd);
+
+    // TODO: neglect mean direction deflection because Trajectory is a const ref
     (void)mean_dir;
+
+    // update particle direction after continuous loss caused by multiple
+    // scattering
     auto vec = corsika::geometry::QuantityVector(
         final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f);
     vP.SetMomentum(corsika::stack::MomentumVector(
@@ -74,11 +87,26 @@ namespace corsika::process::proposal {
   EProcessReturn ContinuousProcess::DoContinuous(setup::Stack::ParticleType& vP,
                                                  setup::Trajectory const& vT) {
     if (!CanInteract(vP.GetPID())) return process::EProcessReturn::eOk;
+
+    // calculate passed grammage
     auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
-    auto energy_loss = TotalEnergyLoss(vP, dX);
-    if (vP.GetChargeNumber() != 0) Scatter(vP, energy_loss, dX);
-    vP.SetEnergy(vP.GetEnergy() - energy_loss);
-    if (vP.GetEnergy() <= cut.GetECut()) return process::EProcessReturn::eParticleAbsorbed;
+
+    // Get or build corresponding track integral calculator and solve the
+    // integral
+    auto c = GetCalculator(vP, calc);
+    auto final_energy = get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
+                            vP.GetEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
+                        1_MeV;
+
+    // if the particle has a charge take multiple scattering into account
+    if (vP.GetChargeNumber() != 0 && vP.GetEnergy() > final_energy)
+      Scatter(vP, vP.GetEnergy() - final_energy, dX);
+
+    // 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 (std::fabs((vP.GetEnergy() - cut.GetECut() )/ cut.GetECut()) < 0.01)
+      return process::EProcessReturn::eParticleAbsorbed;
     vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
     return process::EProcessReturn::eOk;
   }
@@ -88,12 +116,22 @@ namespace corsika::process::proposal {
       setup::Stack::ParticleType const& vP, setup::Trajectory const& vT) {
     if (!CanInteract(vP.GetPID()))
       return units::si::meter * std::numeric_limits<double>::infinity();
+
+    // Limit the step size of a conitnous loss. The maximal continuous loss seems to be a
+    // hyper parameter which must be adjusted.
     auto energy_lim = 0.9 * vP.GetEnergy();
+
+    // If energy lim is below the cut it couldn't be descriped by these process. So the
+    // energy_lim is replaced by the cut.
     if (cut.GetECut() > energy_lim) energy_lim = cut.GetECut();
+
+    // solving the track integral for giving energy lim
     auto c = GetCalculator(vP, calc);
     auto grammage = get<DISPLACEMENT>(c->second)->SolveTrackIntegral(
                         vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) *
                     1_g / square(1_cm);
+
+    // return it in distance aequivalent
     return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
   }
 

Processes/Proposal/ContinuousProcess.h

@@ -43,19 +43,14 @@ namespace corsika::process::proposal {
     void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
 
   public:
+
+    //!
     //! Produces the continuous loss calculator for leptons based on nuclear
     //! compositions and stochastic description limited by the particle cut.
     //!
     template <typename TEnvironment>
     ContinuousProcess(TEnvironment const&, particle_cut::ParticleCut&);
 
-    //!
-    //! Calculates the energy of the particle which has passed the given
-    //! grammage.
-    //!
-    template <typename Particle>
-    corsika::units::si::HEPEnergyType TotalEnergyLoss(
-        Particle const&, corsika::units::si::GrammageType const&);
 
     //!
     //! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into

Processes/Proposal/Interaction.cc

@@ -31,10 +31,20 @@ namespace corsika::process::proposal {
 
   void Interaction::BuildCalculator(particles::Code code,
                                     environment::NuclearComposition const& comp) {
+    // search crosssection builder for given particle
     auto p_cross = cross.find(code);
     if (p_cross == cross.end())
       throw std::runtime_error("PROPOSAL could not find corresponding builder");
+
+    // interpolate the crosssection for given media and energy cut. These may
+    // take some minutes if you have to build the tables and cannot read the
+    // from disk
     auto c = p_cross->second(media.at(&comp), cut);
+
+    // Look which interactions take place and build the corresponding
+    // interaction and secondarie builder. The interaction integral will
+    // interpolated too and saved in the calc map by a key build out of a hash
+    // of composed of the component and particle code.
     auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
     calc[std::make_pair(&comp, code)] = std::make_tuple(
         PROPOSAL::make_secondaries(inter_types, particle[code], media.at(&comp)),
@@ -45,11 +55,19 @@ namespace corsika::process::proposal {
   corsika::process::EProcessReturn Interaction::DoInteraction(
       setup::StackView::StackIterator& vP) {
     if (CanInteract(vP.GetPID())) {
-      auto c = GetCalculator(vP, calc); // [CrossSections]
+      // Get or build corresponding calculators
+      auto c = GetCalculator(vP, calc);
+
+      // Get the rates of the interaction types for every component.
       std::uniform_real_distribution<double> distr(0., 1.);
+
+      // sample a interaction-type, loss and component
       auto rates = get<INTERACTION>(c->second)->Rates(vP.GetEnergy() / 1_MeV);
       auto [type, comp_ptr, v] = get<INTERACTION>(c->second)->SampleLoss(
           vP.GetEnergy() / 1_MeV, rates, distr(fRNG));
+
+      // Read how much random numbers are required to calculate the secondaries.
+      // Calculate the secondaries and deploy them on the corsika stack.
       auto rnd = vector<double>(get<SECONDARIES>(c->second)->RequiredRandomNumbers(type));
       for (auto& it : rnd) it = distr(fRNG);
       auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,

Processes/Proposal/Interaction.h

@@ -23,24 +23,44 @@ namespace corsika::process::proposal {
 
   using namespace corsika::units::si;
 
+  //!
+  //! Electro-magnetic and gamma stochastic losses produced by proposal. It makes
+  //! use of interpolation tables which are runtime intensive calculation, but can be
+  //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
+  //!
   class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
 
+    enum { SECONDARIES, INTERACTION };
     using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
                                unique_ptr<PROPOSAL::Interaction>>;
 
-    std::unordered_map<calc_key_t, calculator_t, hash> calc;
+    std::unordered_map<calc_key_t, calculator_t, hash>
+        calc; //!< Stores the secondaries and interaction calculators.
 
+    //!
+    //! Build the secondaries and interaction calculators and add it to calc.
+    //!
     void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
 
-    enum { SECONDARIES, INTERACTION };
-
   public:
+    //!
+    //! Produces the stoachastic loss calculator for leptons based on nuclear
+    //! compositions and stochastic description limited by the particle cut.
+    //!
     template <typename TEnvironment>
     Interaction(TEnvironment const& env, particle_cut::ParticleCut& cut);
 
+    //!
+    //! Calculate the rates for the different targets and interactions. Sample a
+    //! pair of interaction-type, component and rate, followed by sampling a loss and
+    //! produce the corresponding secondaries and store them on the particle stack.
+    //!
     template <typename Particle>
     corsika::process::EProcessReturn DoInteraction(Particle&);
 
+    //!
+    //! Calculates the  mean free path length
+    //!
     template <typename TParticle>
     corsika::units::si::GrammageType GetInteractionLength(TParticle const& p);
   };