Merge branch '97-cascade-example-debugging-2' into 'master'

Resolve "cascade example debugging"

Closes #97

See merge request AirShowerPhysics/corsika!36

#.gitlab-ci.yml#

+image: ubuntu:bionic
+
+variables:
+  GIT_SSL_NO_VERIFY: "1"
+
+before_script:
+  - apt-get update --yes
+  - apt-get install --yes cmake libboost-dev libeigen3-dev python3 gfortran
+
+build:
+  stage: build
+  tags:
+    - run1
+  script:
+    - mkdir build
+    - cd build
+    - cmake .. 
+    - cmake --build .
+    - ctest -V
+
+# code_quality:
+#   image: docker:stable
+#   variables:
+#     DOCKER_DRIVER: overlay2
+#   allow_failure: true
+#   services:
+#     - docker:stable-dind
+#   script:
+#     - export SP_VERSION=$(echo "$CI_SERVER_VERSION" | sed 's/^\([0-9]*\)\.\([0-9]*\).*/\1-\2-stable/')
+#     - docker run
+#         --env SOURCE_CODE="$PWD"
+#         --volume "$PWD":/code
+#         --volume /var/run/docker.sock:/var/run/docker.sock
+#         "registry.gitlab.com/gitlab-org/security-products/codequality:$SP_VERSION" /code
+#   artifacts:
+#     reports:
+#       codequality: gl-code-quality-report.json

CMakeLists.txt

@@ -28,7 +28,7 @@ if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
 endif()
  
 if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
-  message(STATUS "Setting build type to '${default_bluild_type}' as no other was specified.")
+  message(STATUS "Setting build type to '${default_build_type}' as no other was specified.")
   set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
       STRING "Choose the type of build." FORCE)
   # Set the possible values of build type for cmake-gui

Documentation/Examples/cascade_example.cc

@@ -119,8 +119,10 @@ public:
 
   template <typename Particle>
   LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
+    cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
+    cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
     const Code pid = p.GetPID();
-    if (isEmParticle(pid) || isInvisible(pid)) {
+    if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
       cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
       return 0_m;
     } else {
@@ -134,40 +136,26 @@ public:
   EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
     cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
     const Code pid = p.GetPID();
+    EProcessReturn ret = EProcessReturn::eOk;
     if (isEmParticle(pid)) {
       cout << "removing em. particle..." << endl;
       fEmEnergy += p.GetEnergy();
       fEmCount += 1;
       p.Delete();
+      ret = EProcessReturn::eParticleAbsorbed;
     } else if (isInvisible(pid)) {
       cout << "removing inv. particle..." << endl;
       fInvEnergy += p.GetEnergy();
       fInvCount += 1;
       p.Delete();
+      ret = EProcessReturn::eParticleAbsorbed;
     } else if (isBelowEnergyCut(p)) {
       cout << "removing low en. particle..." << endl;
       fEnergy += p.GetEnergy();
       p.Delete();
+      ret = EProcessReturn::eParticleAbsorbed;
     }
-    // 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;
+    return ret;
   }
 
   void Init() {
@@ -182,10 +170,11 @@ public:
   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
+         << " 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
+         << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
          << " ******************************" << endl;
   }
 

Processes/Sibyll/Decay.h

@@ -42,13 +42,23 @@ namespace corsika::process {
         }
       }
 
+      void setUnstable(const corsika::particles::Code pCode) const {
+        int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+        s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
+      }
+
+      void setStable(const corsika::particles::Code pCode) const {
+        int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+        s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+      }
+
       void setAllStable() const {
+        // name? also makes EM particles stable
 
-        // using namespace corsika::io;
         using std::cout;
         using std::endl;
 
-        // name? also makes EM particles stable
+        std::cout << "Decay: setting all particles stable.." << std::endl;
 
         // loop over all particles in sibyll
         // should be changed to loop over human readable list
@@ -58,11 +68,7 @@ namespace corsika::process {
           const int sibCode = (int)p;
           // skip unknown and antiparticles
           if (sibCode < 1) continue;
-          const corsika::particles::Code pid =
-              ConvertFromSibyll(static_cast<SibyllCode>(sibCode));
-          std::cout << "Sibyll: Decay: setting " << pid << " 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;
         }
       }
 
@@ -124,6 +130,7 @@ namespace corsika::process {
             corsika::particles::GetLifetime(p.GetPID());
         cout << "Decay: code: " << p.GetPID() << endl;
         cout << "Decay: MinStep: t0: " << t0 << endl;
+        cout << "Decay: MinStep: energy: " << E / 1_GeV << endl;
         cout << "Decay: MinStep: gamma: " << gamma << endl;
         // cout << "Decay: MinStep: density: " << density << endl;
         // return as column density
@@ -139,20 +146,29 @@ namespace corsika::process {
 
       template <typename Particle, typename Stack>
       void DoDecay(Particle& p, Stack& s) const {
+        using corsika::geometry::Point;
+        using namespace corsika::units::si;
         SibStack ss;
         ss.Clear();
         // copy particle to sibyll stack
         auto pin = ss.NewParticle();
-        pin.SetPID(process::sibyll::ConvertToSibyllRaw(p.GetPID()));
+        const corsika::particles::Code pCode = p.GetPID();
+        pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
         pin.SetEnergy(p.GetEnergy());
         pin.SetMomentum(p.GetMomentum());
+        // remember position
+        Point decayPoint = p.GetPosition();
+        TimeType t0 = p.GetTime();
         // remove original particle from corsika stack
         p.Delete();
         // set all particles/hadrons unstable
-        setHadronsUnstable();
+        // setHadronsUnstable();
+        setUnstable(pCode);
         // call sibyll decay
         std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
         decsib_();
+        // reset to stable
+        setStable(pCode);
         // print output
         int print_unit = 6;
         sib_list_(print_unit);
@@ -169,6 +185,8 @@ namespace corsika::process {
           pnew.SetEnergy(psib.GetEnergy());
           pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
           pnew.SetMomentum(psib.GetMomentum());
+          pnew.SetPosition(decayPoint);
+          pnew.SetTime(t0);
         }
         // empty sibyll stack
         ss.Clear();

Processes/Sibyll/Interaction.h

@@ -28,12 +28,6 @@ namespace corsika::process::sibyll {
       corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
       rmng.RegisterRandomStream("s_rndm");
 
-      // test random number generator
-      std::cout << "Interaction: "
-                << " test sequence of random numbers." << std::endl;
-      int a = 0;
-      for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
-
       // initialize Sibyll
       sibyll_ini_();
     }
@@ -56,9 +50,8 @@ namespace corsika::process::sibyll {
 
       // beam particles for sibyll : 1, 2, 3 for p, pi, k
       // read from cross section code table
-      int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
-
-      bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
+      const int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
+      const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
 
       /*
          the target should be defined by the Environment,
@@ -67,7 +60,7 @@ namespace corsika::process::sibyll {
        */
       // target nuclei: A < 18
       // FOR NOW: assume target is oxygen
-      int kTarget = 16;
+      const int kTarget = corsika::particles::Oxygen::GetNucleusA();
 
       hep::EnergyType Etot =
           p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
@@ -77,8 +70,8 @@ namespace corsika::process::sibyll {
       Ptot += p.GetMomentum();
       Ptot += pTarget;
       // calculate cm. energy
-      hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
-      double Ecm = sqs / 1_GeV;
+      const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
+      const double Ecm = sqs / 1_GeV;
 
       std::cout << "Interaction: "
                 << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
@@ -149,8 +142,8 @@ namespace corsika::process::sibyll {
         CoordinateSystem& rootCS =
             RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
-        QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
-        Point pOrig(rootCS, coordinates);
+        Point pOrig = p.GetPosition();
+        TimeType tOrig = p.GetTime();
 
         /*
            the target should be defined by the Environment,
@@ -160,8 +153,9 @@ namespace corsika::process::sibyll {
            here we need: GetTargetMassNumber() or GetTargetPID()??
                          GetTargetMomentum() (zero in EAS)
         */
-        // FOR NOW: set target to proton
-        int kTarget = 1; // env.GetTargetParticle().GetPID();
+        // FOR NOW: set target to oxygen
+        const int kTarget = corsika::particles::Oxygen::
+            GetNucleusA(); // env.GetTargetParticle().GetPID();
 
         cout << "defining target momentum.." << endl;
         // FOR NOW: target is always at rest
@@ -171,6 +165,8 @@ namespace corsika::process::sibyll {
              << pTarget.GetComponents() / 1_GeV * constants::c << endl;
         cout << "beam momentum (GeV/c): "
              << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+        cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
+        cout << "time: " << tOrig << endl;
 
         // get energy of particle from stack
         /*
@@ -184,7 +180,7 @@ namespace corsika::process::sibyll {
         EnergyType E = p.GetEnergy();
         EnergyType Etot = E + Etarget;
         // total momentum
-        MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+        MomentumVector Ptot = p.GetMomentum();
         // invariant mass, i.e. cm. energy
         EnergyType Ecm =
             sqrt(Etot * Etot - Ptot.squaredNorm()); // sqrt( 2. * E * 0.93827_GeV );
@@ -208,8 +204,8 @@ namespace corsika::process::sibyll {
                   << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
         if (E < 8.5_GeV || Ecm < 10_GeV) {
           std::cout << "Interaction: "
-                    << " DoInteraction: dropping particle.." << std::endl;
-          p.Delete();
+                    << " DoInteraction: should have dropped particle.." << std::endl;
+          // p.Delete(); delete later... different process
         } else {
           // Sibyll does not know about units..
           double sqs = Ecm / 1_GeV;
@@ -265,6 +261,8 @@ namespace corsika::process::sibyll {
             corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
                 p_lab_components[0], p_lab_components[1], p_lab_components[2]};
             pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
+            pnew.SetPosition(pOrig);
+            pnew.SetTime(tOrig);
             Ptot_final += pnew.GetMomentum();
           }
           // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV

Processes/Sibyll/sibyll2.3c.h

@@ -59,7 +59,7 @@ extern struct {
 // extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
 
 // sibyll main subroutine
-void sibyll_(int&, int&, double&);
+void sibyll_(const int&, const int&, const double&);
 
 // subroutine to initiate sibyll
 void sibyll_ini_();
@@ -79,8 +79,9 @@ void decsib_();
 // interaction length
 // double fpni_(double&, int&);
 
-void sib_sigma_hnuc_(int&, int&, double&, double&, double&);
-void sib_sigma_hp_(int&, double&, double&, double&, double&, double*, double&, double&);
+void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&);
+void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
+                   double&);
 
 double s_rndm_(int&);
 

Stack/SuperStupidStack/SuperStupidStack.h

@@ -74,12 +74,9 @@ namespace corsika::stack {
         return GetStackData().GetTime(GetIndex());
       }
 
-      //#warning this does not really work, nor makes sense:
       corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>
       GetDirection() const {
-        auto P = GetMomentum();
-        return P / P.norm() * 1e10 *
-               (corsika::units::si::meter / corsika::units::si::second);
+        return GetMomentum() / GetEnergy() * corsika::units::constants::c;
       }
     };