Merge branch 'master' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika

Framework/Cascade/Cascade.h

@@ -71,8 +71,9 @@ namespace corsika::cascade {
           fProcessSequence.GetTotalInverseInteractionLength(particle, step);
 
       // sample random exponential step length in grammage
-      std::exponential_distribution expDist((1_m * 1_m / 1_g) / total_inv_lambda);
-      GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
+      auto constexpr grammageConversion = 1_g / (1_m * 1_m);
+      std::exponential_distribution expDist(1 / (grammageConversion * total_inv_lambda));
+      GrammageType const next_interact = grammageConversion * expDist(fRNG);
 
       std::cout << "total_inv_lambda=" << total_inv_lambda
                 << ", next_interact=" << next_interact << std::endl;

Framework/ProcessSequence/ContinuousProcess.h

@@ -13,6 +13,7 @@
 #define _include_corsika_continuousprocess_h_
 
 #include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::process {
 
@@ -32,8 +33,12 @@ namespace corsika::process {
 
     // here starts the interface part
     // -> enforce derived to implement DoContinuous...
-    template <typename P, typename T, typename S>
-    inline EProcessReturn DoContinuous(P&, T&, S&) const;
+    template <typename Particle, typename Track, typename Stack>
+    EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
+    
+    // -> enforce derived to implement MaxStepLength...
+    template <typename Particle, typename Track>
+    corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;
   };
 
 } // namespace corsika::process

Framework/ProcessSequence/DecayProcess.h

@@ -29,14 +29,13 @@ namespace corsika::process {
 
   template <typename derived>
   struct DecayProcess {
-
     derived& GetRef() { return static_cast<derived&>(*this); }
     const derived& GetRef() const { return static_cast<const derived&>(*this); }
 
     /// here starts the interface-definition part
     // -> enforce derived to implement DoDecay...
     template <typename Particle, typename Stack>
-    inline EProcessReturn DoDecay(Particle&, Stack&) const;
+    EProcessReturn DoDecay(Particle&, Stack&) const;
 
     template <typename Particle>
     corsika::units::si::TimeType GetLifetime(Particle& p) const;

Framework/ProcessSequence/ProcessSequence.h

@@ -71,8 +71,9 @@ namespace corsika::process {
 
     template <typename Particle, typename Track>
     corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const {
-      corsika::units::si::LengthType max_length =
+      corsika::units::si::LengthType max_length = // if no other process in the sequence implements it
           std::numeric_limits<double>::infinity() * corsika::units::si::meter;
+          
       if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
                     is_process_sequence<T1>::value) {
         corsika::units::si::LengthType const len = A.MaxStepLength(p, track);

Processes/NullModel/NullModel.cc

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -11,37 +10,4 @@
 
 #include <corsika/process/null_model/NullModel.h>
 
-#include <corsika/logging/Logger.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <iostream>
-#include <limits>
-using namespace std;
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::process::null_model;
-
-template <typename Stack>
-NullModel<Stack>::NullModel() {}
-
-template <typename Stack>
-NullModel<Stack>::~NullModel() {}
-
-template <typename Stack>
-process::EProcessReturn NullModel<Stack>::DoContinuous(Particle&, setup::Trajectory&,
-                                                       Stack&) const {
-  return EProcessReturn::eOk;
-}
-
-template <typename Stack>
-double NullModel<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
-  return std::numeric_limits<double>::infinity();
-}
-
-template <typename Stack>
-void NullModel<Stack>::Init() {}
-
-#include <corsika/setup/SetupStack.h>
-
-template class process::null_model::NullModel<setup::Stack>;
+void corsika::process::null_model::NullModel::Init() {}

Processes/NullModel/NullModel.h

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -15,26 +14,29 @@
 #include <corsika/process/ContinuousProcess.h>
 #include <corsika/setup/SetupTrajectory.h>
 
-namespace corsika::process {
-
-  namespace null_model {
+namespace corsika::process::null_model {
 
-    template <typename Stack>
-    class NullModel {
+  class NullModel : public corsika::process::ContinuousProcess<NullModel> {
+    corsika::units::si::LengthType fMaxStepLength{
+        corsika::units::si::meter * std::numeric_limits<double>::infinity()};
 
-      typedef typename Stack::ParticleType Particle;
+  public:
+    NullModel(corsika::units::si::LengthType maxStepLength)
+        : fMaxStepLength(maxStepLength) {}
 
-    public:
-      NullModel();
-      ~NullModel();
+    void Init();
 
-      void Init();
-      EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
-      double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
-    };
+    template <typename Particle, typename Track, typename Stack>
+    process::EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
+      return EProcessReturn::eOk;
+    }
 
-  } // namespace null_model
+    template <typename Particle, typename Track>
+    corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const {
+      return fMaxStepLength;
+    }
+  };
 
-} // namespace corsika::process
+} // namespace corsika::process::null_model
 
 #endif

Processes/NullModel/testNullModel.cc

@@ -43,11 +43,13 @@ TEST_CASE("NullModel", "[processes]") {
 
   SECTION("interface") {
 
-    NullModel<setup::Stack> model;
+    NullModel model(10_m);
 
     model.Init();
     [[maybe_unused]] const process::EProcessReturn ret =
         model.DoContinuous(particle, track, stack);
-    [[maybe_unused]] const double length = model.MaxStepLength(particle, track);
+    LengthType const length = model.MaxStepLength(particle, track);
+
+    CHECK((length / 10_m) == Approx(1));
   }
 }

Processes/Sibyll/Decay.h

@@ -86,9 +86,8 @@ namespace corsika::process {
         // i.e. corsika::particles::ListOfParticles()
         std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
         // make ALL particles unstable, then set EM stable
-        for (auto& p : corsika2sibyll) {
+        for (int sibCode : corsika2sibyll) {
           // std::cout << (int)p << std::endl;
-          const int sibCode = (int)p;
           // skip unknown and antiparticles
           if (sibCode < 1) continue;
           // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
@@ -128,19 +127,14 @@ namespace corsika::process {
 
         const corsika::units::si::TimeType t0 =
             corsika::particles::GetLifetime(p.GetPID());
+        auto const lifetime = gamma * t0;
+
         cout << "Decay: code: " << p.GetPID() << endl;
         cout << "Decay: MinStep: t0: " << t0 << endl;
-        cout << "Decay: MinStep: energy: " << E / 1_GeV << endl;
+        cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
         cout << "Decay: MinStep: gamma: " << gamma << endl;
-        // cout << "Decay: MinStep: density: " << density << endl;
-        // return as column density
-        // const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
-        // cout << "Decay: MinStep: x0: " << x0 << endl;
-        corsika::units::si::TimeType const lifetime = gamma * t0;
         cout << "Decay: MinStep: tau: " << lifetime << endl;
-        // int a = 1;
-        // const double x = -x0 * log(s_rndm_(a));
-        // cout << "Decay: next decay: " << x << endl;
+
         return lifetime;
       }
 

Processes/Sibyll/ParticleConversion.h

@@ -34,9 +34,8 @@ namespace corsika::process::sibyll {
   }
 
   SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
-    //~ assert(handledBySibyll(pCode));
     return static_cast<SibyllCode>(
-        corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
+      corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
   }
 
   corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {
@@ -44,8 +43,7 @@ namespace corsika::process::sibyll {
   }
 
   int constexpr ConvertToSibyllRaw(corsika::particles::Code pCode) {
-    return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>(
-        corsika::process::sibyll::ConvertToSibyll(pCode));
+    return static_cast<int>(ConvertToSibyll(pCode));
   }
 
   int constexpr GetSibyllXSCode(corsika::particles::Code pCode) {

Processes/Sibyll/sibyll2.3c.cc

 #include <corsika/process/sibyll/sibyll2.3c.h>
 
 #include <corsika/random/RNGManager.h>
+#include <random>
 
 double s_rndm_(int&) {
-  static corsika::random::RNG& rmng =
+  static corsika::random::RNG& rng =
       corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-  ;
-  return rmng() / (double)rmng.max();
+  
+  std::uniform_real_distribution<double> dist;
+  return dist(rng);
 }

Stack/SuperStupidStack/SuperStupidStack.h

@@ -45,15 +45,19 @@ namespace corsika::stack {
       void SetPID(const corsika::particles::Code id) {
         GetStackData().SetPID(GetIndex(), id);
       }
+
       void SetEnergy(const corsika::units::hep::EnergyType& e) {
         GetStackData().SetEnergy(GetIndex(), e);
       }
+
       void SetMomentum(const MomentumVector& v) {
         GetStackData().SetMomentum(GetIndex(), v);
       }
+
       void SetPosition(const corsika::geometry::Point& v) {
         GetStackData().SetPosition(GetIndex(), v);
       }
+
       void SetTime(const corsika::units::si::TimeType& v) {
         GetStackData().SetTime(GetIndex(), v);
       }
@@ -61,15 +65,19 @@ namespace corsika::stack {
       corsika::particles::Code GetPID() const {
         return GetStackData().GetPID(GetIndex());
       }
+
       corsika::units::hep::EnergyType GetEnergy() const {
         return GetStackData().GetEnergy(GetIndex());
       }
+
       MomentumVector GetMomentum() const {
         return GetStackData().GetMomentum(GetIndex());
       }
+
       corsika::geometry::Point GetPosition() const {
         return GetStackData().GetPosition(GetIndex());
       }
+
       corsika::units::si::TimeType GetTime() const {
         return GetStackData().GetTime(GetIndex());
       }
@@ -99,22 +107,30 @@ namespace corsika::stack {
       }
 
       int GetSize() const { return fDataPID.size(); }
+
       int GetCapacity() const { return fDataPID.size(); }
 
       void SetPID(const int i, const corsika::particles::Code id) { fDataPID[i] = id; }
+
       void SetEnergy(const int i, const corsika::units::hep::EnergyType e) {
         fDataE[i] = e;
       }
       void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
+
       void SetPosition(const int i, const corsika::geometry::Point& v) {
         fPosition[i] = v;
       }
+      
       void SetTime(const int i, const corsika::units::si::TimeType& v) { fTime[i] = v; }
 
       corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }
+
       corsika::units::hep::EnergyType GetEnergy(const int i) const { return fDataE[i]; }
+
       MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
+
       corsika::geometry::Point GetPosition(const int i) const { return fPosition[i]; }
+
       corsika::units::si::TimeType GetTime(const int i) const { return fTime[i]; }
 
       /**
@@ -156,6 +172,7 @@ namespace corsika::stack {
                             0 * corsika::units::si::meter}));
         fTime.push_back(0 * corsika::units::si::second);
       }
+
       void DecrementSize() {
         if (fDataE.size() > 0) {
           fDataPID.pop_back();