Cleanup

Documentation/Examples/cascade_example.cc

@@ -229,6 +229,8 @@ int main() {
   // setup processes, decays and interactions
   tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
+  
+  corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
   corsika::process::sibyll::Interaction sibyll(env);
   corsika::process::sibyll::Decay decay;
   ProcessCut cut(8_GeV);

Environment/Environment.h

@@ -20,6 +20,8 @@
 #include <limits>
 
 namespace corsika::environment {
+  using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
+    
   struct Universe : public corsika::geometry::Sphere {
     Universe(corsika::geometry::CoordinateSystem const& pCS)
         : corsika::geometry::Sphere(
@@ -37,7 +39,7 @@ namespace corsika::environment {
     Environment()
         : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
                                 .GetRootCoordinateSystem()}
-        , fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+        , fUniverse(std::make_unique<BaseNodeType>(
               std::make_unique<Universe>(fCoordinateSystem))) {}
 
     using IEnvironmentModel = corsika::setup::IEnvironmentModel;
@@ -48,19 +50,19 @@ namespace corsika::environment {
     auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
 
     // factory method for creation of VolumeTreeNodes
-    template <class VolumeType, typename... VolumeArgs>
-    static auto CreateNode(VolumeArgs&&... args) {
-      static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
+    template <class TVolumeType, typename... TVolumeArgs>
+    static auto CreateNode(TVolumeArgs&&... args) {
+      static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
                     "unusable type provided, needs to be derived from "
                     "\"corsika::geometry::Volume\"");
 
-      return std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
-          std::make_unique<VolumeType>(std::forward<VolumeArgs>(args)...));
+      return std::make_unique<BaseNodeType>(
+          std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
     }
 
   private:
     corsika::geometry::CoordinateSystem const& fCoordinateSystem;
-    VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
+    BaseNodeType::VTNUPtr fUniverse;
   };
 
 } // namespace corsika::environment

Environment/HomogeneousMedium.h

@@ -19,6 +19,8 @@
 #include <corsika/random/RNGManager.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <cassert>
+
 /**
  * a homogeneous medium
  */
@@ -42,30 +44,6 @@ namespace corsika::environment {
     }
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
-    corsika::particles::Code const& GetTarget(
-        std::vector<corsika::units::si::CrossSectionType>& sigma) const override {
-      using namespace corsika::units::si;
-      int i = -1;
-      corsika::units::si::CrossSectionType total_weighted_sigma = 0._mbarn;
-      std::vector<float> fractions;
-      for (auto w : fNuclComp.GetFractions()) {
-        i++;
-        std::cout << "HomogeneousMedium: fraction: " << w << std::endl;
-        total_weighted_sigma += w * sigma[i];
-        fractions.push_back(w * sigma[i] / 1_mbarn);
-      }
-
-      for (auto f : fractions) {
-        f = f / (total_weighted_sigma / 1_mbarn);
-        std::cout << "HomogeneousMedium: reweighted fraction: " << f << std::endl;
-      }
-      std::discrete_distribution channelDist(fractions.begin(), fractions.end());
-      static corsika::random::RNG& kRNG =
-          corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-      const int ichannel = channelDist(kRNG);
-      return fNuclComp.GetComponents()[ichannel];
-    }
-
     corsika::units::si::GrammageType IntegratedGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const&,
         corsika::units::si::LengthType pTo) const override {

Environment/IMediumModel.h

@@ -16,6 +16,7 @@
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Trajectory.h>
+#include <random>
 #include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::environment {
@@ -38,9 +39,27 @@ namespace corsika::environment {
         corsika::units::si::GrammageType) const = 0;
 
     virtual NuclearComposition const& GetNuclearComposition() const = 0;
+    
+    template <class TRNG>
+    corsika::particles::Code SampleTarget(
+        std::vector<corsika::units::si::CrossSectionType> const& sigma, TRNG& randomStream) const {
+      using namespace corsika::units::si;
+      
+      auto const& nuclComp = GetNuclearComposition();
+      auto const& fractions = nuclComp.GetFractions();
+      assert(sigma.size() == fractions.size());
+      
+      std::vector<float> weights(fractions.size());
+      
+      for (size_t i = 0; i < fractions.size(); ++i) {
+        std::cout << "HomogeneousMedium: fraction: " << fractions[i] << std::endl;
+        weights[i] = fractions[i] * sigma[i].magnitude();
+      }
 
-    virtual corsika::particles::Code const& GetTarget(
-        std::vector<corsika::units::si::CrossSectionType>&) const = 0;
+      std::discrete_distribution channelDist(weights.cbegin(), weights.cend());
+      const int iChannel = channelDist(randomStream);
+      return nuclComp.GetComponents()[iChannel];
+    }
   };
 
 } // namespace corsika::environment

Environment/NuclearComposition.h

@@ -15,6 +15,7 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <numeric>
 #include <stdexcept>
+#include <cassert>
 #include <vector>
 
 namespace corsika::environment {
@@ -28,6 +29,8 @@ namespace corsika::environment {
                        std::vector<float> pFractions)
         : fNumberFractions(pFractions)
         , fComponents(pComponents) {
+            
+      assert(pComponents.size() == pFractions.size());
       auto const sumFractions =
           std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
 

Framework/Cascade/Cascade.h

@@ -56,6 +56,7 @@ namespace corsika::cascade {
       }
     }
 
+private:
     void Step(Particle& particle) {
       using namespace corsika::units::si;
 

Processes/Sibyll/Decay.h

@@ -21,8 +21,6 @@
 
 #include <corsika/particles/ParticleProperties.h>
 
-#include <fenv.h>
-
 namespace corsika::process {
 
   namespace sibyll {
@@ -81,7 +79,7 @@ namespace corsika::process {
         // i.e. corsika::particles::ListOfParticles()
         for (auto& p : corsika2sibyll) {
           // std::cout << (int)p << std::endl;
-          const int sibCode = (int)p;
+          const int sibCode = static_cast<int>(p);
           // skip unknown and antiparticles
           if (sibCode < 1) continue;
           s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
@@ -114,7 +112,7 @@ namespace corsika::process {
         }
         // set Leptons and Proton and Neutron stable
         // use stack to loop over particles
-        const std::vector<corsika::particles::Code> particleList = {
+        constexpr corsika::particles::Code particleList[] = {
             corsika::particles::Code::Proton,   corsika::particles::Code::Neutron,
             corsika::particles::Code::Electron, corsika::particles::Code::Positron,
             corsika::particles::Code::NuE,      corsika::particles::Code::NuEBar,
@@ -131,7 +129,7 @@ namespace corsika::process {
       }
 
       template <typename Particle>
-      corsika::units::si::TimeType GetLifetime(Particle& p) {
+      corsika::units::si::TimeType GetLifetime(Particle const& p) {
         using std::cout;
         using std::endl;
         using namespace corsika::units::si;
@@ -169,11 +167,6 @@ namespace corsika::process {
         using corsika::geometry::Point;
         using namespace corsika::units::si;
 
-        // TODO: this should be done in a central, common place. Not here..
-#ifndef CORSIKA_OSX
-        feenableexcept(FE_INVALID);
-#endif
-
         fCount++;
         SibStack ss;
         ss.Clear();
@@ -189,8 +182,8 @@ namespace corsika::process {
         // with sibyll internal values
         pin.SetMass(corsika::particles::GetMass(pCode));
         // remember position
-        Point decayPoint = p.GetPosition();
-        TimeType t0 = p.GetTime();
+        Point const decayPoint = p.GetPosition();
+        TimeType const t0 = p.GetTime();
         // remove original particle from corsika stack
         p.Delete();
         // set all particles/hadrons unstable
@@ -219,11 +212,6 @@ namespace corsika::process {
         }
         // empty sibyll stack
         ss.Clear();
-
-        // TODO: this should be done in a central, common place. Not here..
-#ifndef CORSIKA_OSX
-        fedisableexcept(FE_INVALID);
-#endif
       }
     };
   } // namespace sibyll

Processes/Sibyll/Interaction.h

@@ -45,16 +45,13 @@ namespace corsika::process::sibyll {
       using std::cout;
       using std::endl;
 
-      corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-      rmng.RegisterRandomStream("s_rndm");
-
       // initialize Sibyll
       sibyll_ini_();
     }
 
     std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
-        const corsika::particles::Code& BeamId, const corsika::particles::Code& TargetId,
-        const corsika::units::si::HEPEnergyType& CoMenergy) {
+        const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
+        const corsika::units::si::HEPEnergyType CoMenergy) {
       using namespace corsika::units::si;
       double sigProd, dummy, dum1, dum2, dum3, dum4;
       double dumdif[3];
@@ -67,19 +64,18 @@ namespace corsika::process::sibyll {
               "Sibyll target outside range. Only nuclei with A<18 are allowed.");
         sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy);
         return std::make_tuple(sigProd * 1_mbarn, iTarget);
+      } else if (TargetId == corsika::particles::Proton::GetCode()) {
+        sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
+        return std::make_tuple(sigProd * 1_mbarn, 1);
       } else {
-        if (TargetId == corsika::particles::Proton::GetCode()) {
-          sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
-          return std::make_tuple(sigProd * 1_mbarn, 1);
-        } else
-          // no interaction in sibyll possible, return infinite cross section? or throw?
-          sigProd = std::numeric_limits<double>::infinity();
+        // no interaction in sibyll possible, return infinite cross section? or throw?
+        sigProd = std::numeric_limits<double>::infinity();
         return std::make_tuple(sigProd * 1_mbarn, 0);
       }
     }
 
     template <typename Particle, typename Track>
-    corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) {
+    corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
 
       using namespace corsika::units;
       using namespace corsika::units::si;
@@ -105,11 +101,13 @@ namespace corsika::process::sibyll {
 
       // total momentum and energy
       HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
-      MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-      Ptot += p.GetMomentum();
-      Ptot += pTarget;
+      MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+      pTotLab += p.GetMomentum();
+      pTotLab += pTarget;
+      auto const pTotLabNorm = pTotLab.norm();
       // calculate cm. energy
-      const HEPEnergyType ECoM = sqrt(Elab * Elab - Ptot.squaredNorm());
+      const HEPEnergyType ECoM = sqrt(
+          (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
 
       std::cout << "Interaction: LambdaInt: \n"
                 << " input energy: " << p.GetEnergy() / 1_GeV << endl
@@ -136,7 +134,7 @@ namespace corsika::process::sibyll {
         // get weights of components from environment/medium
         const auto w = mediumComposition.GetFractions();
         // loop over components in medium
-        for (auto targetId : mediumComposition.GetComponents()) {
+        for (auto const targetId : mediumComposition.GetComponents()) {
           i++;
           cout << "Interaction: get interaction length for target: " << targetId << endl;
 
@@ -155,7 +153,7 @@ namespace corsika::process::sibyll {
 
         // calculate interaction length in medium
 #warning check interaction length units
-        GrammageType int_length =
+        GrammageType const int_length =
             avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
         std::cout << "Interaction: "
                   << "interaction length (g/cm2): "
@@ -238,7 +236,7 @@ namespace corsika::process::sibyll {
 
         // sample target mass number
         const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
-        const auto mediumComposition =
+        const auto& mediumComposition =
             currentNode->GetModelProperties().GetNuclearComposition();
         // get cross sections for target materials
         /*
@@ -246,14 +244,18 @@ namespace corsika::process::sibyll {
           should be passed from GetInteractionLength if possible
          */
 #warning reading interaction cross section again, should not be necessary
-        std::vector<si::CrossSectionType> cross_section_of_components;
-        for (auto targetId : mediumComposition.GetComponents()) {
+        auto const& compVec = mediumComposition.GetComponents();
+        std::vector<si::CrossSectionType> cross_section_of_components(
+            mediumComposition.GetComponents().size());
+
+        for (size_t i = 0; i < compVec.size(); ++i) {
+          auto const targetId = compVec[i];
           const auto [sigProd, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm);
-          cross_section_of_components.push_back(sigProd);
+          cross_section_of_components[i] = sigProd;
         }
 
-        const auto targetCode =
-            currentNode->GetModelProperties().GetTarget(cross_section_of_components);
+        const auto targetCode = currentNode->GetModelProperties().SampleTarget(
+            cross_section_of_components, fRNG);
         cout << "Interaction: target selected: " << targetCode << endl;
         /*
           FOR NOW: allow nuclei with A<18 or protons only.
@@ -336,6 +338,8 @@ namespace corsika::process::sibyll {
 
   private:
     corsika::environment::Environment const& fEnvironment;
+    corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
   };
 
 } // namespace corsika::process::sibyll

Processes/Sibyll/testSibyll.cc

@@ -13,6 +13,8 @@
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 
+#include <corsika/random/RNGManager.h>
+
 #include <corsika/particles/ParticleProperties.h>
 
 #include <corsika/geometry/Point.h>
@@ -109,6 +111,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
   geometry::Line line(origin, v);
   geometry::Trajectory<geometry::Line> track(line, 10_s);
 
+  corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+
   SECTION("InteractionInterface") {
 
     setup::Stack stack;