Merge branch '109-add-shower-visualization' into 'master'

Resolve "add shower visualization"

Closes #109

See merge request AirShowerPhysics/corsika!40

Documentation/Examples/CMakeLists.txt

@@ -28,6 +28,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogg
    ProcessSibyll
   CORSIKAcascade
   ProcessStackInspector
+  ProcessTrackWriter
   CORSIKAprocesses
   CORSIKAparticles
   CORSIKAgeometry

Documentation/Examples/cascade_example.cc

@@ -25,6 +25,8 @@
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
 
+#include <corsika/process/track_writer/TrackWriter.h>
+
 #include <corsika/units/PhysicalUnits.h>
 
 #include <corsika/random/RNGManager.h>
@@ -192,7 +194,6 @@ public:
 // The example main program for a particle cascade
 //
 int main() {
-
   // initialize random number sequence(s)
   corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
 
@@ -222,9 +223,10 @@ int main() {
   corsika::process::sibyll::Interaction sibyll;
   corsika::process::sibyll::Decay decay;
   ProcessCut cut(8_GeV);
+  corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
 
   // assemble all processes into an ordered process list
-  auto sequence = p0 << sibyll << decay << cut;
+  auto sequence = p0 << sibyll << decay << cut << trackWriter;
 
   // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
   // << "\n";
@@ -232,7 +234,7 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const hep::EnergyType E0 = 1_TeV;
+  const hep::EnergyType E0 = 10_TeV;
   {
     auto particle = stack.NewParticle();
     particle.SetPID(Code::Proton);
@@ -241,7 +243,7 @@ int main() {
     particle.SetEnergy(E0);
     particle.SetMomentum(plab);
     particle.SetTime(0_ns);
-    Point p(rootCS, 0_m, 0_m, 10_km);
+    Point p(rootCS, 0_m, 0_m, 0_m);
     particle.SetPosition(p);
     cout << particle.GetEnergy() / 1_GeV << endl;
   }

Framework/Cascade/Cascade.h

@@ -57,13 +57,7 @@ namespace corsika::cascade {
     }
 
     void Step(Particle& particle) {
-
       using namespace corsika::units::si;
-      using std::cout;
-      using std::endl;
-      using std::log;
-
-      cout << particle.GetEnergy() / 1_GeV << endl;
 
       // determine geometric tracking
       corsika::setup::Trajectory step = fTracking.GetTrack(particle);
@@ -80,11 +74,15 @@ namespace corsika::cascade {
                 << ", next_interact=" << next_interact << std::endl;
 
       // convert next_step from grammage to length
+      auto const* currentNode =
+          fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+
+      if (currentNode == &*fEnvironment.GetUniverse()) {
+        throw std::runtime_error("particle entered void universe");
+      }
+
       LengthType const distance_interact =
-          fEnvironment.GetUniverse()
-              ->GetContainingNode(particle.GetPosition())
-              ->GetModelProperties()
-              .ArclengthFromGrammage(step, next_interact);
+          currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
 
       // determine the maximum geometric step length
       LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
@@ -116,6 +114,8 @@ namespace corsika::cascade {
       particle.SetPosition(step.PositionFromArclength(min_distance));
       // .... also update time, momentum, direction, ...
 
+      step.LimitEndTo(min_distance);
+
       // apply all continuous processes on particle + track
       corsika::process::EProcessReturn status =
           fProcessSequence.DoContinuous(particle, step, fStack);
@@ -127,11 +127,11 @@ namespace corsika::cascade {
         return;
       }
 
-      std::cout << "sth. happening before geometric limit ?"
+      std::cout << "sth. happening before geometric limit ? "
                 << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
 
       if (min_distance < distance_max) { // interaction to happen within geometric limit
-        // check weather decay or interaction limits this step
+        // check whether decay or interaction limits this step
 
         if (min_distance == distance_interact) {
           std::cout << "collide" << std::endl;

Framework/Geometry/Trajectory.h

@@ -43,6 +43,10 @@ namespace corsika::geometry {
       assert(t >= 0 * corsika::units::si::second);
       return T::ArcLength(0, t);
     }
+
+    void LimitEndTo(corsika::units::si::LengthType limit) {
+      fTimeLength = T::TimeFromArclength(limit);
+    }
   };
 
 } // namespace corsika::geometry

Processes/CMakeLists.txt

@@ -2,6 +2,7 @@
 add_subdirectory (NullModel)
 add_subdirectory (Sibyll)
 add_subdirectory (StackInspector)
+add_subdirectory (TrackWriter)
 add_subdirectory (TrackingLine)
 
 #add_custom_target(CORSIKAprocesses)

Processes/Sibyll/Interaction.h

@@ -75,8 +75,9 @@ namespace corsika::process::sibyll {
       // FOR NOW: assume target is oxygen
       const int kTarget = corsika::particles::Oxygen::GetNucleusA();
 
-      hep::EnergyType Etot =
-          p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+      const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+                                                corsika::particles::Neutron::GetMass());
+      hep::EnergyType Etot = p.GetEnergy() + nucleon_mass;
       MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
       // FOR NOW: assume target is at rest
       MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
@@ -143,6 +144,18 @@ namespace corsika::process::sibyll {
 
         Point pOrig = p.GetPosition();
         TimeType tOrig = p.GetTime();
+        // sibyll CS has z along particle momentum
+        // FOR NOW: hard coded z-axis for corsika frame
+        QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
+        QuantityVector<length_d> const xAxis{1_m, 0_m, 0_m};
+        [[maybe_unused]] auto pt = [](MomentumVector p) {
+          return sqrt(p.GetComponents()[0] * p.GetComponents()[0] +
+                      p.GetComponents()[1] * p.GetComponents()[1]);
+        };
+        double theta = acos(p.GetMomentum().GetComponents()[2] / p.GetMomentum().norm());
+        cout << "ProcessSibyll: polar angle between sibyllCS and rootCS: " << theta
+             << endl;
+        CoordinateSystem sibyllCS = rootCS.rotate(xAxis, theta);
 
         /*
            the target should be defined by the Environment,
@@ -157,11 +170,16 @@ namespace corsika::process::sibyll {
             GetNucleusA(); // env.GetTargetParticle().GetPID();
 
         // FOR NOW: target is always at rest
-        const EnergyType Etarget = 0_GeV + corsika::particles::Proton::GetMass();
+        const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+                                                  corsika::particles::Neutron::GetMass());
+        const EnergyType Etarget = 0_GeV + nucleon_mass;
         const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
         cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
         cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
              << endl;
+        cout << "beam momentum in sibyll frame (GeV/c): "
+             << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV << endl;
+
         cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
         cout << "time: " << tOrig << endl;
 
@@ -220,23 +238,31 @@ namespace corsika::process::sibyll {
 
           // add particles from sibyll to stack
           // link to sibyll stack
+          // here we need to pass projectile momentum and energy to define the local
+          // sibyll frame and the boosts to the lab. frame
           SibStack ss;
 
-          // SibStack does not know about momentum yet so we need counter to access
           // momentum array in Sibyll
-          int i = -1;
           MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
           EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
           for (auto& psib : ss) {
-            ++i;
+
             // skip particles that have decayed in Sibyll
-            if (abs(s_plist_.llist[i]) > 100) continue;
+            if (psib.HasDecayed()) 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();
+            // FOR NOW: fill vector in sibCS and then rotate into rootCS
+            // can be done in SibStack by passing sibCS
+            // get momentum vector in sibyllCS
+            const auto pSibyllComponentsSibCS = psib.GetMomentum().GetComponents();
+            // temporary vector in sibyllCS
+            auto SibVector = MomentumVector(sibyllCS, pSibyllComponentsSibCS);
+            // rotatate to rootCS
+            const auto pSibyllComponents = SibVector.GetComponents(rootCS);
+            // boost to lab. frame
             hep::EnergyType en_lab = 0. * 1_GeV;
             hep::MomentumType p_lab_components[3];
             en_lab = psib.GetEnergy() * gamma;
@@ -255,7 +281,6 @@ namespace corsika::process::sibyll {
             auto pnew = s.NewParticle();
             pnew.SetEnergy(en_lab);
             pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
-
             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));

Processes/Sibyll/SibStack.h

@@ -31,7 +31,7 @@ namespace corsika::process::sibyll {
     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; }
@@ -91,6 +91,9 @@ namespace corsika::process::sibyll {
     corsika::units::hep::EnergyType GetEnergy() const {
       return GetStackData().GetEnergy(GetIndex());
     }
+    bool HasDecayed() const {
+      return abs(GetStackData().GetId(GetIndex())) > 100 ? true : false;
+    }
     void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
     corsika::process::sibyll::SibyllCode GetPID() const {
       return static_cast<corsika::process::sibyll::SibyllCode>(

Processes/TrackWriter/CMakeLists.txt

+
+set (
+  MODEL_SOURCES
+  TrackWriter.cc
+  )
+
+set (
+  MODEL_HEADERS
+  TrackWriter.h
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/track_writer
+  )
+
+add_library (ProcessTrackWriter STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackWriter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+  ProcessTrackWriter
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+#  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessTrackWriter
+  CORSIKAunits
+  CORSIKAparticles
+  CORSIKAgeometry
+  CORSIKAsetup
+  )
+
+target_include_directories (
+  ProcessTrackWriter 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS ProcessTrackWriter
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )
+
+
+# --------------------
+# code unit testing
+#add_executable (testNullModel testNullModel.cc)
+
+#target_link_libraries (
+#  testNullModel  ProcessNullModel
+#  CORSIKAsetup
+#  CORSIKAgeometry
+#  CORSIKAunits
+#  CORSIKAthirdparty # for catch2
+#  )
+
+#add_test (NAME testNullModel COMMAND testNullModel)
+

Processes/TrackWriter/TrackWriter.cc

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/track_writer/TrackWriter.h>
+#include <string>
+
+void corsika::process::TrackWriter::TrackWriter::Init() {
+  using namespace std::string_literals;
+
+  fFile.open(fFilename);
+  fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
+        << '\n';
+}

Processes/TrackWriter/TrackWriter.h

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _Processes_TrackWriter_h_
+#define _Processes_TrackWriter_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <fstream>
+#include <limits>
+#include <string>
+
+namespace corsika::process::TrackWriter {
+
+  class TrackWriter : public corsika::process::ContinuousProcess<TrackWriter> {
+
+  public:
+    TrackWriter(std::string const& filename)
+        : fFilename(filename) {}
+
+    void Init();
+
+    template <typename Particle, typename Track, typename Stack>
+    corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&) {
+      using namespace corsika::units::si;
+      auto const start = t.GetPosition(0).GetCoordinates();
+      auto const delta = t.GetPosition(1).GetCoordinates() - start;
+      auto const& name = corsika::particles::GetName(p.GetPID());
+
+      fFile << name << "    " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
+            << start[1] / 1_m << ' ' << start[2] / 1_m << "   " << delta[0] / 1_m << ' '
+            << delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
+
+      return corsika::process::EProcessReturn::eOk;
+    }
+
+    template <typename Particle, typename Track>
+    corsika::units::si::LengthType MaxStepLength(Particle&, Track&) {
+      return corsika::units::si::meter * std::numeric_limits<double>::infinity();
+    }
+
+  private:
+    std::string const fFilename;
+    std::ofstream fFile;
+  };
+
+} // namespace corsika::process::TrackWriter
+
+#endif

Processes/TrackingLine/TrackingLine.h

@@ -84,7 +84,8 @@ namespace corsika::process {
             p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
 
         auto const currentPosition = p.GetPosition();
-
+        std::cout << "TrackingLine pid: " << p.GetPID()
+                  << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
         std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
                   << std::endl;
         std::cout << "TrackingLine   E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;

Processes/TrackingLine/testTrackingLine.cc

@@ -9,6 +9,7 @@
  */
 
 #include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
 
 #include <corsika/process/tracking_line/TrackingLine.h>
 
@@ -47,6 +48,7 @@ struct DummyParticle {
   auto GetEnergy() const { return fEnergy; }
   auto GetMomentum() const { return fMomentum; }
   auto GetPosition() const { return fPosition; }
+  auto GetPID() const { return corsika::particles::Code::Unknown; }
 };
 
 struct DummyStack {

Tools/plot_tracks.sh

+#!/bin/sh
+
+# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+#
+# See file AUTHORS for a list of contributors.
+#
+# This software is distributed under the terms of the GNU General Public
+# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+# the license.
+
+# with this script you can plot an animation of output of TrackWriter
+
+track_dat=$1
+if [ -z "$track_dat" ]; then
+  echo "usage: $0 <track.dat> [output.gif]" >&2
+  exit 1
+fi
+
+output=$2
+if [ -z "$output" ]; then
+  output="$track_dat.gif"
+fi
+
+cat <<EOF | gnuplot
+set term gif animate size 600,600
+set output "$output"
+
+set xlabel "x / m"
+set ylabel "y / m"
+set zlabel "z / m"
+set title "CORSIKA 8 preliminary"
+
+do for [t=0:360:1] {
+	set view 60, t
+        splot "$track_dat" u 3:4:5:6:7:8 w vectors nohead t ""
+}
+EOF
+
+exit $?