Merge branch 'sibyll' into 'master'

Sibyll

Closes #72

See merge request !28

CMakeLists.txt

@@ -7,6 +7,7 @@ project (
   LANGUAGES CXX
   )
 
+# as long as there still are modules using it:
 enable_language (Fortran)
 
 # ignore many irrelevant Up-to-date messages during install

Framework/Cascade/SibStack.h

@@ -7,9 +7,12 @@
 #include <corsika/cascade/sibyll2.3c.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 #include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
 
 using namespace std;
 using namespace corsika::stack;
+using namespace corsika::units;
+using namespace corsika::geometry;
 
 class SibStackData {
 
@@ -19,13 +22,30 @@ public:
   void Clear() { s_plist_.np = 0; }
 
   int GetSize() const { return s_plist_.np; }
+#warning check actual capacity of sibyll stack
   int GetCapacity() const { return 8000; }
 
   void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
-  void SetEnergy(const int i, const double v) { s_plist_.p[3][i] = v; }
+  void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
+  void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
+    auto tmp = v.GetComponents();
+    for (int idx = 0; idx < 3; ++idx)
+      s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;
+  }
 
   int GetId(const int i) const { return s_plist_.llist[i]; }
-  double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
+
+  EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+  super_stupid::MomentumVector GetMomentum(const int i) const {
+    CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+    corsika::geometry::QuantityVector<momentum_d> components{
+        s_plist_.p[0][i] * 1_GeV / si::constants::c,
+        s_plist_.p[1][i] * 1_GeV / si::constants::c,
+        s_plist_.p[2][i] * 1_GeV / si::constants::c};
+    super_stupid::MomentumVector v1(rootCS, components);
+    return v1;
+  }
 
   void Copy(const int i1, const int i2) {
     s_plist_.llist[i1] = s_plist_.llist[i2];
@@ -45,13 +65,19 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
   using ParticleBase<StackIteratorInterface>::GetIndex;
 
 public:
-  void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
-  double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+  void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
+  EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
   void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
   corsika::process::sibyll::SibyllCode GetPID() const {
     return static_cast<corsika::process::sibyll::SibyllCode>(
         GetStackData().GetId(GetIndex()));
   }
+  super_stupid::MomentumVector GetMomentum() const {
+    return GetStackData().GetMomentum(GetIndex());
+  }
+  void SetMomentum(const super_stupid::MomentumVector& v) {
+    GetStackData().SetMomentum(GetIndex(), v);
+  }
 };
 
 typedef Stack<SibStackData, ParticleInterface> SibStack;

Framework/Geometry/BaseTrajectory.h

@@ -46,7 +46,7 @@ namespace corsika::geometry {
         corsika::units::si::LengthType) const = 0;
 
     virtual LengthType ArcLength(corsika::units::si::TimeType t1,
-                                   corsika::units::si::TimeType t2) const = 0;
+                                 corsika::units::si::TimeType t2) const = 0;
 
     virtual corsika::units::si::TimeType GetDuration(
         corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {

Framework/Geometry/Helix.h

@@ -58,8 +58,8 @@ namespace corsika::geometry {
 
     auto GetRadius() const { return radius; }
 
-    corsika::units::si::LengthType ArcLength(
-      corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
+    corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
+                                             corsika::units::si::TimeType t2) const {
       return (vPar + vPerp).norm() * (t2 - t1);
     }
 

Framework/Geometry/Trajectory.h

@@ -12,8 +12,8 @@
 #ifndef _include_TRAJECTORY_H
 #define _include_TRAJECTORY_H
 
-#include <corsika/units/PhysicalUnits.h>
 #include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
 
 using corsika::units::si::LengthType;
 using corsika::units::si::TimeType;

Processes/Sibyll/CMakeLists.txt

@@ -20,6 +20,7 @@ set (
 set (
   MODEL_HEADERS
   ParticleConversion.h
+  ProcessDecay.h
   ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
   )
 

Processes/Sibyll/ProcessDecay.h

+#ifndef _include_ProcessDecay_h_
+#define _include_ProcessDecay_h_
+
+#include <corsika/process/ContinuousProcess.h>
+
+#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+// using namespace corsika::particles;
+
+namespace corsika::process {
+
+  namespace sibyll {
+
+    void setHadronsUnstable() {
+      // name? also makes EM particles stable
+
+      // loop over all particles in sibyll
+      // should be changed to loop over human readable list
+      // 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) {
+        // 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(
+        // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
+        s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
+        // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
+        // std::endl;
+      }
+      // set Leptons and Proton and Neutron stable
+      // use stack to loop over particles
+      setup::Stack ds;
+      ds.NewParticle().SetPID(corsika::particles::Code::Proton);
+      ds.NewParticle().SetPID(corsika::particles::Code::Neutron);
+      ds.NewParticle().SetPID(corsika::particles::Code::Electron);
+      ds.NewParticle().SetPID(corsika::particles::Code::Positron);
+      ds.NewParticle().SetPID(corsika::particles::Code::NuE);
+      ds.NewParticle().SetPID(corsika::particles::Code::NuEBar);
+      ds.NewParticle().SetPID(corsika::particles::Code::MuMinus);
+      ds.NewParticle().SetPID(corsika::particles::Code::MuPlus);
+      ds.NewParticle().SetPID(corsika::particles::Code::NuMu);
+      ds.NewParticle().SetPID(corsika::particles::Code::NuMuBar);
+
+      for (auto& p : ds) {
+        int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+        // set particle stable by setting table value negative
+        //	cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
+        //     << " stable in Sibyll .." << endl;
+        s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+        p.Delete();
+      }
+    }
+
+    class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+    public:
+      ProcessDecay() {}
+      void Init() {
+        // setHadronsUnstable();
+      }
+
+      void setAllStable() {
+        // name? also makes EM particles stable
+
+        // loop over all particles in sibyll
+        // should be changed to loop over human readable list
+        // i.e. corsika::particles::ListOfParticles()
+        for (auto& p : 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(static_cast<SibyllCode>(sibCode)) << " stable"
+                    << std::endl;
+          s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
+          std::cout << "decay table value: " << s_csydec_.idb[sibCode - 1] << std::endl;
+        }
+      }
+
+      friend void setHadronsUnstable();
+
+      template <typename Particle>
+      double MinStepLength(Particle& p, setup::Trajectory&) const {
+        corsika::units::hep::EnergyType E = p.GetEnergy();
+        corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
+        // env.GetDensity();
+
+        const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
+
+        const double gamma = E / m;
+
+        TimeType t0 = GetLifetime(p.GetPID());
+        cout << "ProcessDecay: code: " << (p.GetPID()) << endl;
+        cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
+        cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
+        cout << "ProcessDecay: MinStep: density: " << density << endl;
+        // return as column density
+        const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
+        cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
+        int a = 1;
+        const double x = -x0 * log(s_rndm_(a));
+        cout << "ProcessDecay: next decay: " << x << endl;
+        return x;
+      }
+
+      template <typename Particle, typename Stack>
+      void DoDiscrete(Particle& p, Stack& s) const {
+        SibStack ss;
+        ss.Clear();
+        // copy particle to sibyll stack
+        auto pin = ss.NewParticle();
+        pin.SetPID(process::sibyll::ConvertToSibyllRaw(p.GetPID()));
+        pin.SetEnergy(p.GetEnergy());
+        pin.SetMomentum(p.GetMomentum());
+        // remove original particle from corsika stack
+        p.Delete();
+        // set all particles/hadrons unstable
+        setHadronsUnstable();
+        // call sibyll decay
+        std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
+        decsib_();
+        // print output
+        int print_unit = 6;
+        sib_list_(print_unit);
+        // copy particles from sibyll stack to corsika
+        int i = -1;
+        for (auto& psib : ss) {
+          ++i;
+          // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
+          if (abs(s_plist_.llist[i]) > 100) continue;
+          // add to corsika stack
+          // cout << "decay product: " << process::sibyll::ConvertFromSibyll(
+          // psib.GetPID() ) << endl;
+          auto pnew = s.NewParticle();
+          pnew.SetEnergy(psib.GetEnergy());
+          pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+          pnew.SetMomentum(psib.GetMomentum());
+        }
+        // empty sibyll stack
+        ss.Clear();
+      }
+
+      template <typename Particle, typename Stack>
+      EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+        return EProcessReturn::eOk;
+      }
+    };
+  } // namespace sibyll
+} // namespace corsika::process
+
+#endif

Processes/Sibyll/code_generator.py

@@ -159,8 +159,6 @@ if __name__ == "__main__":
     pythia_db = load_pythiadb(sys.argv[1])
     read_sibyll_codes(sys.argv[2], pythia_db)
     
-    print (str(pythia_db))
-    
     with open("Generated.inc", "w") as f:
         print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
         print(generate_sibyll_enum(pythia_db), file=f)

Processes/StackInspector/StackInspector.cc

@@ -10,9 +10,9 @@
  */
 
 #include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/particles/ParticleProperties.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/units/PhysicalUnits.h>
-#include <corsika/particles/ParticleProperties.h>
 
 #include <corsika/logging/Logger.h>
 
@@ -28,7 +28,8 @@ using namespace corsika::process::stack_inspector;
 
 template <typename Stack>
 StackInspector<Stack>::StackInspector(const bool aReport)
-    : fReport(aReport) {}
+    : fReport(aReport)
+    , fCountStep(0) {}
 
 template <typename Stack>
 StackInspector<Stack>::~StackInspector() {}
@@ -36,7 +37,6 @@ StackInspector<Stack>::~StackInspector() {}
 template <typename Stack>
 process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
                                                             Stack& s) const {
-  static int countStep = 0;
   if (!fReport) return EProcessReturn::eOk;
   [[maybe_unused]] int i = 0;
   EnergyType Etot = 0_GeV;
@@ -51,8 +51,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
          << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
          << " pos=" << pos << endl;
   }
-  countStep++;
-  cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize()
+  fCountStep++;
+  cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
        << " Estack=" << Etot / 1_GeV << " GeV" << endl;
   return EProcessReturn::eOk;
 }
@@ -63,7 +63,9 @@ double StackInspector<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const
 }
 
 template <typename Stack>
-void StackInspector<Stack>::Init() {}
+void StackInspector<Stack>::Init() {
+  fCountStep = 0;
+}
 
 #include <corsika/setup/SetupStack.h>
 

Processes/StackInspector/StackInspector.h

@@ -40,6 +40,7 @@ namespace corsika::process {
 
     private:
       bool fReport;
+      mutable int fCountStep = 0;
     };
 
   } // namespace stack_inspector

Processes/TrackingLine/testTrackingLine.cc

@@ -13,8 +13,8 @@
 #include <corsika/process/tracking_line/TrackingLine.h>
 
 #include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
 #include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
 
 #include <corsika/setup/SetupTrajectory.h>
 using corsika::setup::Trajectory;
@@ -33,66 +33,80 @@ using namespace std;
 using namespace corsika::units::si;
 
 struct DummyParticle {
-    EnergyType fEnergy;
-    Vector<momentum_d> fMomentum;
-    Point fPosition;
-    
-    DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition) : fEnergy(pEnergy), fMomentum(pMomentum), fPosition(pPosition) {}
-    
-    auto GetEnergy() const { return fEnergy; }
-    auto GetMomentum() const { return fMomentum; }
-    auto GetPosition() const { return fPosition; }
+  EnergyType fEnergy;
+  Vector<momentum_d> fMomentum;
+  Point fPosition;
+
+  DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition)
+      : fEnergy(pEnergy)
+      , fMomentum(pMomentum)
+      , fPosition(pPosition) {}
+
+  auto GetEnergy() const { return fEnergy; }
+  auto GetMomentum() const { return fMomentum; }
+  auto GetPosition() const { return fPosition; }
 };
 
 struct DummyStack {
-    using ParticleType = DummyParticle;
+  using ParticleType = DummyParticle;
 };
 
 TEST_CASE("TrackingLine") {
   corsika::environment::Environment env; // dummy environment
   auto const& cs = env.GetCoordinateSystem();
-  
+
   tracking_line::TrackingLine<DummyStack> tracking(env);
 
   SECTION("intersection with sphere") {
     Point const origin(cs, {0_m, 0_m, 0_m});
     Point const center(cs, {0_m, 0_m, 10_m});
     Sphere const sphere(center, 1_m);
-    Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+    Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+                                                            0_m / second, 1_m / second);
     Line line(origin, v);
-    
+
     geometry::Trajectory<Line> traj(line, 12345_s);
-    
-    auto const opt = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+
+    auto const opt =
+        tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
     REQUIRE(opt.has_value());
-    
+
     auto [t1, t2] = opt.value();
     REQUIRE(t1 / 9_s == Approx(1));
     REQUIRE(t2 / 11_s == Approx(1));
-    
-    auto const optNoIntersection = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+
+    auto const optNoIntersection =
+        tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
     REQUIRE_FALSE(optNoIntersection.has_value());
   }
-  
+
   SECTION("maximally possible propagation") {
-      auto& universe = *(env.GetUniverse());
-      
-      //~ std::cout << env.GetUniverse().get() << std::endl;
-      
-      DummyParticle p(1_J, Vector<momentum_d>(cs, 0*kilogram*meter/second, 0*kilogram*meter/second, 1*kilogram*meter/second), Point(cs, 0_m, 0_m,0_m));
-      
-      auto const radius = 20_m;
-      
-      auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
-      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);      
-      universe.AddChild(std::move(theMedium));
-      
-      Point const origin(cs, {0_m, 0_m, 0_m});
-      Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
-      Line line(origin, v);
-      
-      auto const traj = tracking.GetTrack(p);
-      
-      REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius)).GetComponents(cs).norm().magnitude() == Approx(0).margin(1e-4));
+    auto& universe = *(env.GetUniverse());
+
+    //~ std::cout << env.GetUniverse().get() << std::endl;
+
+    DummyParticle p(
+        1_J,
+        Vector<momentum_d>(cs, 0 * kilogram * meter / second,
+                           0 * kilogram * meter / second, 1 * kilogram * meter / second),
+        Point(cs, 0_m, 0_m, 0_m));
+
+    auto const radius = 20_m;
+
+    auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+        Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+    universe.AddChild(std::move(theMedium));
+
+    Point const origin(cs, {0_m, 0_m, 0_m});
+    Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+                                                            0_m / second, 1_m / second);
+    Line line(origin, v);
+
+    auto const traj = tracking.GetTrack(p);
+
+    REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
+                .GetComponents(cs)
+                .norm()
+                .magnitude() == Approx(0).margin(1e-4));
   }
 }

Setup/SetupEnvironment.h

@@ -15,7 +15,7 @@
 #include <corsika/environment/IMediumModel.h>
 
 namespace corsika::setup {
-    using IEnvironmentModel = corsika::environment::IMediumModel;
+  using IEnvironmentModel = corsika::environment::IMediumModel;
 }
 
 #endif

Stack/SuperStupidStack/SuperStupidStack.h

@@ -137,7 +137,7 @@ namespace corsika::stack {
       void IncrementSize() {
         fDataPID.push_back(Code::Unknown);
         fDataE.push_back(0 * joule);
-#warning this here makes no sense: see issue #48
+        //#TODO this here makes no sense: see issue #48
         geometry::CoordinateSystem& dummyCS =
             geometry::RootCoordinateSystem::GetInstance().GetRootCS();
         fMomentum.push_back(MomentumVector(