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

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

Documentation/Examples/cascade_example.cc

@@ -27,6 +27,8 @@
 #include <corsika/geometry/Sphere.h>
 #include <corsika/process/sibyll/ParticleConversion.h>
 
+#include <corsika/process/sibyll/ProcessDecay.h>
+
 #include <corsika/units/PhysicalUnits.h>
 
 #include <iostream>
@@ -38,6 +40,7 @@ using namespace corsika::process;
 using namespace corsika::units;
 using namespace corsika::particles;
 using namespace corsika::random;
+using namespace corsika::setup;
 using namespace corsika::geometry;
 using namespace corsika::environment;
 
@@ -45,16 +48,206 @@ using namespace std;
 using namespace corsika::units::si;
 
 static int fCount = 0;
+static EnergyType fEnergy = 0. * 1_GeV;
+
+// FOR NOW: global static variables for ParticleCut process
+// this is just wrong...
+static EnergyType fEmEnergy;
+static int fEmCount;
+
+static EnergyType fInvEnergy;
+static int fInvCount;
+
+class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
+public:
+  ProcessEMCut() {}
+  template <typename Particle>
+  bool isBelowEnergyCut(Particle& p) const {
+    // FOR NOW: center-of-mass energy hard coded
+    const EnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+    if (p.GetEnergy() < 50_GeV || Ecm < 10_GeV)
+      return true;
+    else
+      return false;
+  }
+
+  bool isEmParticle(Code pCode) const {
+    bool is_em = false;
+    // FOR NOW: switch
+    switch (pCode) {
+      case Code::Electron:
+        is_em = true;
+        break;
+      case Code::Gamma:
+        is_em = true;
+        break;
+      default:
+        break;
+    }
+    return is_em;
+  }
+
+  void defineEmParticles() const {
+    // create bool array identifying em particles
+  }
+
+  bool isInvisible(Code pCode) const {
+    bool is_inv = false;
+    // FOR NOW: switch
+    switch (pCode) {
+      case Code::NuE:
+        is_inv = true;
+        break;
+      case Code::NuEBar:
+        is_inv = true;
+        break;
+      case Code::NuMu:
+        is_inv = true;
+        break;
+      case Code::NuMuBar:
+        is_inv = true;
+        break;
+      case Code::MuPlus:
+        is_inv = true;
+        break;
+      case Code::MuMinus:
+        is_inv = true;
+        break;
+
+      default:
+        break;
+    }
+    return is_inv;
+  }
+
+  template <typename Particle>
+  double MinStepLength(Particle& p, setup::Trajectory&) const {
+    const Code pid = p.GetPID();
+    if (isEmParticle(pid) || isInvisible(pid)) {
+      cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+      return 0.;
+    } else {
+      double next_step = std::numeric_limits<double>::infinity();
+      cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+      return next_step;
+    }
+  }
+
+  template <typename Particle, typename Stack>
+  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+    // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
+    // cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
+    // cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
+    // const Code pid = p.GetPID();
+    // if( isEmParticle( pid ) ){
+    //   cout << "removing em. particle..." << endl;
+    //   fEmEnergy += p.GetEnergy();
+    //   fEmCount  += 1;
+    //   p.Delete();
+    //   return EProcessReturn::eParticleAbsorbed;
+    // }
+    // if ( isInvisible( pid ) ){
+    //   cout << "removing inv. particle..." << endl;
+    //   fInvEnergy += p.GetEnergy();
+    //   fInvCount  += 1;
+    //   p.Delete();
+    //   return EProcessReturn::eParticleAbsorbed;
+    // }
+    return EProcessReturn::eOk;
+  }
+
+  template <typename Particle, typename Stack>
+  void DoDiscrete(Particle& p, Stack&) const {
+    cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
+    const Code pid = p.GetPID();
+    if (isEmParticle(pid)) {
+      cout << "removing em. particle..." << endl;
+      fEmEnergy += p.GetEnergy();
+      fEmCount += 1;
+      p.Delete();
+    } else if (isInvisible(pid)) {
+      cout << "removing inv. particle..." << endl;
+      fInvEnergy += p.GetEnergy();
+      fInvCount += 1;
+      p.Delete();
+    } else if (isBelowEnergyCut(p)) {
+      cout << "removing low en. particle..." << endl;
+      fEnergy += p.GetEnergy();
+      fCount += 1;
+      p.Delete();
+    }
+  }
+
+  void Init() {
+    fEmEnergy = 0. * 1_GeV;
+    fEmCount = 0;
+    fInvEnergy = 0. * 1_GeV;
+    fInvCount = 0;
+    fEnergy = 0. * 1_GeV;
+    // defineEmParticles();
+  }
+
+  void ShowResults() {
+    cout << " ******************************" << endl
+         << " ParticleCut: " << endl
+         << " energy in em.  component (GeV): " << fEmEnergy / 1_GeV << endl
+         << " no. of em.  particles injected: " << fEmCount << endl
+         << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+         << " no. of inv. particles injected: " << fInvCount << endl
+         << " ******************************" << endl;
+  }
+
+  EnergyType GetInvEnergy() { return fInvEnergy; }
+
+  EnergyType GetCutEnergy() { return fEnergy; }
+
+  EnergyType GetEmEnergy() { return fEmEnergy; }
+
+private:
+};
 
 class ProcessSplit : public corsika::process::InteractionProcess<ProcessSplit> {
 public:
   ProcessSplit() {}
 
+  void setTrackedParticlesStable() const {
+    /*
+      Sibyll is hadronic generator
+      only hadrons decay
+     */
+    // set particles unstable
+    corsika::process::sibyll::setHadronsUnstable();
+    // make tracked particles stable
+    std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
+    setup::Stack ds;
+    ds.NewParticle().SetPID(Code::PiPlus);
+    ds.NewParticle().SetPID(Code::PiMinus);
+    ds.NewParticle().SetPID(Code::KPlus);
+    ds.NewParticle().SetPID(Code::KMinus);
+    ds.NewParticle().SetPID(Code::K0Long);
+    ds.NewParticle().SetPID(Code::K0Short);
+
+    for (auto& p : ds) {
+      int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+      // set particle stable by setting table value negative
+      s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+      p.Delete();
+    }
+  }
+
   template <typename Particle, typename Track>
   double GetInteractionLength(Particle& p, Track&) const {
+
+    // coordinate system, get global frame of reference
+    CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
+    const Code corsikaBeamId = p.GetPID();
+
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
-    int kBeam = 1;
+    int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
+
+    bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
 
     /*
        the target should be defined by the Environment,
@@ -63,33 +256,55 @@ public:
      */
     // target nuclei: A < 18
     // FOR NOW: assume target is oxygen
-    int kTarget = 1;
-    double beamEnergy = p.GetEnergy() / 1_GeV;
+    int kTarget = 16;
+
+    EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+    super_stupid::MomentumVector Ptot(
+        rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+    // FOR NOW: assume target is at rest
+    super_stupid::MomentumVector pTarget(
+        rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+    Ptot += p.GetMomentum();
+    Ptot += pTarget;
+    // calculate cm. energy
+    EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared);
+    double Ecm = sqs / 1_GeV;
+
     std::cout << "ProcessSplit: "
-              << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
-    double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
-    double dumdif[3];
+              << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
+              << " beam can interact:" << kBeam << endl
+              << " beam XS code:" << kBeam << endl
+              << " beam pid:" << p.GetPID() << endl
+              << " target mass number:" << kTarget << std::endl;
 
-    if (kTarget == 1)
-      sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
-    else
-      sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy);
+    double int_length = 0;
+    if (kInteraction) {
 
-    std::cout << "ProcessSplit: "
-              << "MinStep: sibyll return: " << prodCrossSection << std::endl;
-    CrossSectionType sig = prodCrossSection * 1_mbarn;
-    std::cout << "ProcessSplit: "
-              << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+      double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
+      double dumdif[3];
+
+      if (kTarget == 1)
+        sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
+      else
+        sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
+
+      std::cout << "ProcessSplit: "
+                << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+      CrossSectionType sig = prodCrossSection * 1_mbarn;
+      std::cout << "ProcessSplit: "
+                << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+      const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+      std::cout << "ProcessSplit: "
+                << "nucleon mass " << nucleon_mass << std::endl;
+      // calculate interaction length in medium
+      int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
+      // pick random step lenth
+      std::cout << "ProcessSplit: "
+                << "interaction length (g/cm2): " << int_length << std::endl;
+    } else
+      int_length = std::numeric_limits<double>::infinity();
 
-    const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
-    std::cout << "ProcessSplit: "
-              << "nucleon mass " << nucleon_mass << std::endl;
-    // calculate interaction length in medium
-    double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
-    // pick random step lenth
-    std::cout << "ProcessSplit: "
-              << "interaction length (g/cm2): " << int_length << std::endl;
-    // add exponential sampling
     /*
       what are the units of the output? slant depth or 3space length?
 
@@ -110,9 +325,38 @@ public:
 
   template <typename Particle, typename Stack>
   void DoInteraction(Particle& p, Stack& s) const {
-    cout << "DoInteraction: " << p.GetPID() << " interaction? "
+    cout << "ProcessSibyll: "
+         << "DoInteraction: " << p.GetPID() << " interaction? "
          << process::sibyll::CanInteract(p.GetPID()) << endl;
     if (process::sibyll::CanInteract(p.GetPID())) {
+      cout << "defining coordinates" << endl;
+      // coordinate system, get global frame of reference
+      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
+      QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+      Point pOrig(rootCS, coordinates);
+
+      /*
+         the target should be defined by the Environment,
+         ideally as full particle object so that the four momenta
+         and the boosts can be defined..
+
+         here we need: GetTargetMassNumber() or GetTargetPID()??
+                       GetTargetMomentum() (zero in EAS)
+      */
+      // FOR NOW: set target to proton
+      int kTarget = 1; // env.GetTargetParticle().GetPID();
+
+      cout << "defining target momentum.." << endl;
+      // FOR NOW: target is always at rest
+      const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+      const auto pTarget = super_stupid::MomentumVector(
+          rootCS, 0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c,
+          0. * 1_GeV / si::constants::c);
+      cout << "target momentum (GeV/c): "
+           << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
+      cout << "beam momentum (GeV/c): "
+           << p.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
 
       // get energy of particle from stack
       /*
@@ -121,18 +365,30 @@ public:
         (assuming proton at rest as target AND
         assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
       */
+      // total energy: E_beam + E_target
+      // in lab. frame: E_beam + m_target*c**2
       EnergyType E = p.GetEnergy();
-      EnergyType Ecm = sqrt(2. * E * 0.93827_GeV);
+      EnergyType Etot = E + Etarget;
+      // total momentum
+      super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+      // invariant mass, i.e. cm. energy
+      EnergyType Ecm = sqrt(Etot * Etot -
+                            Ptot.squaredNorm() *
+                                si::constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
+      /*
+       get transformation between Stack-frame and SibStack-frame
+       for EAS Stack-frame is lab. frame, could be different for CRMC-mode
+       the transformation should be derived from the input momenta
+     */
+      const double gamma = Etot / Ecm;
+      const auto gambet = Ptot / (Ecm / si::constants::c);
 
-      int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+      std::cout << "ProcessSplit: "
+                << " DoDiscrete: gamma:" << gamma << endl;
+      std::cout << "ProcessSplit: "
+                << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
 
-      /*
-         the target should be defined by the Environment,
-         ideally as full particle object so that the four momenta
-         and the boosts can be defined..
-       */
-      // FOR NOW: set target to proton
-      int kTarget = 1; // p.GetPID();
+      int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
 
       std::cout << "ProcessSplit: "
                 << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
@@ -148,7 +404,8 @@ public:
         // running sibyll, filling stack
         sibyll_(kBeam, kTarget, sqs);
         // running decays
-        // decsib_();
+        setTrackedParticlesStable();
+        decsib_();
         // print final state
         int print_unit = 6;
         sib_list_(print_unit);
@@ -159,49 +416,62 @@ public:
         // add particles from sibyll to stack
         // link to sibyll stack
         SibStack ss;
-        /*
-          get transformation between Stack-frame and SibStack-frame
-          for EAS Stack-frame is lab. frame, could be different for CRMC-mode
-          the transformation should be derived from the input momenta
-          in general transformation is rotation + boost
-        */
-        const EnergyType proton_mass = 0.93827_GeV;
-        const double gamma = (E + proton_mass) / (Ecm);
-        const double gambet = sqrt(E * E - proton_mass * proton_mass) / Ecm;
 
         // SibStack does not know about momentum yet so we need counter to access momentum
         // array in Sibyll
         int i = -1;
-        for (auto& p : ss) {
+        super_stupid::MomentumVector Ptot_final(
+            rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+        for (auto& psib : ss) {
           ++i;
-          // transform to lab. frame, primitve
-          const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
+          // skip particles that have decayed in Sibyll
+          if (abs(s_plist_.llist[i]) > 100) continue;
+
+          // transform energy to lab. frame, primitve
+          // compute beta_vec * p_vec
+          // arbitrary Lorentz transformation based on sibyll routines
+          const auto gammaBetaComponents = gambet.GetComponents();
+          const auto pSibyllComponents = psib.GetMomentum().GetComponents();
+          EnergyType en_lab = 0. * 1_GeV;
+          MomentumType p_lab_components[3];
+          en_lab = psib.GetEnergy() * gamma;
+          EnergyType pnorm = 0. * 1_GeV;
+          for (int j = 0; j < 3; ++j)
+            pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c) /
+                     (gamma + 1.);
+          pnorm += psib.GetEnergy();
+
+          for (int j = 0; j < 3; ++j) {
+            p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
+                                                             gammaBetaComponents[j] /
+                                                             si::constants::c;
+            en_lab -=
+                (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+          }
+
           // add to corsika stack
           auto pnew = s.NewParticle();
-          pnew.SetEnergy(en_lab * 1_GeV);
-          pnew.SetPID(process::sibyll::ConvertFromSibyll(p.GetPID()));
+          pnew.SetEnergy(en_lab);
+          pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+
+          corsika::geometry::QuantityVector<momentum_d> p_lab_c{
+              p_lab_components[0], p_lab_components[1], p_lab_components[2]};
+          pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
+          Ptot_final += pnew.GetMomentum();
         }
+        // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV *
+        // si::constants::c << endl;
       }
-    } else
-      p.Delete();
+    }
   }
 
   void Init() {
     fCount = 0;
 
-    // define reference frame? --> defines boosts between corsika stack and model stack.
-
-    // initialize random numbers for sibyll
-    // FOR NOW USE SIBYLL INTERNAL !!!
-    //    rnd_ini_();
-
     corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
 
     rmng.RegisterRandomStream("s_rndm");
 
-    // //    corsika::random::RNG srng;
-    // auto srng = rmng.GetRandomStream("s_rndm");
-
     // test random number generator
     std::cout << "ProcessSplit: "
               << " test sequence of random numbers." << std::endl;
@@ -211,29 +481,11 @@ public:
     // initialize Sibyll
     sibyll_ini_();
 
-    // set particles stable / unstable
-    // use stack to loop over particles
-    setup::Stack ds;
-    ds.NewParticle().SetPID(Code::Proton);
-    ds.NewParticle().SetPID(Code::Neutron);
-    ds.NewParticle().SetPID(Code::PiPlus);
-    ds.NewParticle().SetPID(Code::PiMinus);
-    ds.NewParticle().SetPID(Code::KPlus);
-    ds.NewParticle().SetPID(Code::KMinus);
-    ds.NewParticle().SetPID(Code::K0Long);
-    ds.NewParticle().SetPID(Code::K0Short);
-
-    for (auto& p : ds) {
-      int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
-      // set particle stable by setting table value negative
-      cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
-           << " stable in Sibyll .." << endl;
-      s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
-      p.Delete();
-    }
+    setTrackedParticlesStable();
   }
 
   int GetCount() { return fCount; }
+  EnergyType GetEnergy() { return fEnergy; }
 
 private:
 };
@@ -245,9 +497,6 @@ double s_rndm_(int&) {
   return rmng() / (double)rmng.max();
 }
 
-class A {};
-class B : public A {};
-
 int main() {
   corsika::environment::Environment env; // dummy environment
   auto& universe = *(env.GetUniverse());
@@ -266,20 +515,37 @@ int main() {
 
   universe.AddChild(std::move(theMedium));
 
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
   tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
+
   ProcessSplit p1;
-  const auto sequence = p0 + p1;
+  corsika::process::sibyll::ProcessDecay p2;
+  ProcessEMCut p3;
+  const auto sequence = /*p0 +*/ p1 + p2 + p3;
   setup::Stack stack;
 
   corsika::cascade::Cascade EAS(tracking, sequence, stack);
 
   stack.Clear();
   auto particle = stack.NewParticle();
-  EnergyType E0 = 100_TeV;
+  EnergyType E0 = 100_GeV;
+  MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV) / si::constants::c;
+  auto plab = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c,
+                                           0. * 1_GeV / si::constants::c, P0);
   particle.SetEnergy(E0);
+  particle.SetMomentum(plab);
   particle.SetPID(Code::Proton);
+  particle.SetTime(0_ns);
+  Point p(rootCS, 0_m, 0_m, 0_m);
+  particle.SetPosition(p);
   EAS.Init();
   EAS.Run();
-  cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+  cout << "Result: E0=" << E0 / 1_GeV
+       << "GeV, particles below energy threshold =" << p1.GetCount() << endl;
+  cout << "total energy below threshold (GeV): " << p1.GetEnergy() / 1_GeV << std::endl;
+  p3.ShowResults();
+  cout << "total energy (GeV): "
+       << (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
 }

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(