Merge branch 'improvements' into 'master'

Some improvements here and there

See merge request !116

Framework/Geometry/CMakeLists.txt

@@ -10,6 +10,7 @@ set (
   Point.h
   Line.h
   Sphere.h
+  Plane.h
   Volume.h
   CoordinateSystem.h
   RootCoordinateSystem.h

Framework/Geometry/Plane.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 _include_Framework_Geometry_Plane_h_
+#define _include_Framework_Geometry_Plane_h_
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::geometry {
+  class Plane {
+
+    using DimLessVec = Vector<corsika::units::si::dimensionless_d>;
+
+    Point const fCenter;
+    DimLessVec const fNormal;
+
+  public:
+    Plane(Point const& vCenter, DimLessVec const& vNormal)
+        : fCenter(vCenter)
+        , fNormal(vNormal) {}
+    bool IsAbove(Point const& vP) const {
+      return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
+    }
+
+    Point const& GetCenter() const { return fCenter; }
+    DimLessVec const& GetNormal() const { return fNormal; }
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/ProcessSequence/ContinuousProcess.h

@@ -34,11 +34,11 @@ namespace corsika::process {
     // here starts the interface part
     // -> enforce derived to implement DoContinuous...
     template <typename Particle, typename Track>
-    EProcessReturn DoContinuous(Particle&, Track&) const;
+    EProcessReturn DoContinuous(Particle&, Track const&) const;
 
     // -> enforce derived to implement MaxStepLength...
     template <typename Particle, typename Track>
-    corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;
+    units::si::LengthType MaxStepLength(Particle const& p, Track const& track) const;
   };
 
   // overwrite the default trait class, to mark BaseProcess<T> as useful process

Framework/Units/PhysicalConstants.h

@@ -50,6 +50,9 @@ namespace corsika::units::constants {
   constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 
   auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
+  
+  // molar gas constant
+  auto constexpr R = Rep(8.314'459'8) * joule / (mole * kelvin);
 
   // etc.
 

Processes/EnergyLoss/EnergyLoss.cc

@@ -24,9 +24,8 @@ using namespace std;
 
 using namespace corsika;
 using namespace corsika::units::si;
-using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
+using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupTrack = corsika::setup::Trajectory;
 
 namespace corsika::process::energy_loss {
 
@@ -53,7 +52,7 @@ namespace corsika::process::energy_loss {
    * For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
    *
    */
-  HEPEnergyType EnergyLoss::BetheBloch(Particle& p, GrammageType const dX) {
+  HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const dX) {
 
     // all these are material constants and have to come through Environment
     // right now: values for nitrogen_D
@@ -70,20 +69,20 @@ namespace corsika::process::energy_loss {
     // end of material constants
 
     // this is the Bethe-Bloch coefficiet 4pi N_A r_e^2 m_e c^2
-    auto const K = 0.307075_MeV / 1_mol * square(1_cm);
+    auto constexpr K = 0.307075_MeV / 1_mol * square(1_cm);
     HEPEnergyType const E = p.GetEnergy();
     HEPMassType const m = p.GetMass();
     double const gamma = E / m;
-    double const Z = p.GetChargeNumber();
-    double const Z2 = pow(Z, 2);
-    HEPMassType const me = particles::Electron::GetMass();
+    int const Z = p.GetChargeNumber();
+    int const Z2 = Z * Z;
+    HEPMassType constexpr me = particles::Electron::GetMass();
     auto const m2 = m * m;
-    auto const me2 = me * me;
-    double const gamma2 = pow(gamma, 2);
+    auto constexpr me2 = me * me;
+    double const gamma2 = gamma * gamma;
 
     double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2    (1-1/gamma)*(1+1/gamma);
                                                 // (gamma_2-1)/gamma_2 = (1-1/gamma2);
-    double const c2 = 1;                        // HEP convention here c=c2=1
+    double constexpr c2 = 1;                    // HEP convention here c=c2=1
     cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
     [[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
     HEPMassType const Wmax =
@@ -148,7 +147,8 @@ namespace corsika::process::energy_loss {
            (0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
   }
 
-  process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t) {
+  process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p,
+                                                   SetupTrack const& t) {
     if (p.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
     GrammageType const dX =
         p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
@@ -170,12 +170,25 @@ namespace corsika::process::energy_loss {
     p.SetEnergy(Enew);
     MomentumUpdate(p, Enew);
     fEnergyLossTot += dE;
-    GetXbin(p, dE);
+    GetXbin(p, t, dE);
     return status;
   }
 
-  units::si::LengthType EnergyLoss::MaxStepLength(Particle&, Track&) {
-    return units::si::meter * std::numeric_limits<double>::infinity();
+  units::si::LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
+                                                  SetupTrack const& vTrack) const {
+    if (vParticle.GetChargeNumber() == 0) {
+      return units::si::meter * std::numeric_limits<double>::infinity();
+    }
+
+    auto constexpr dX = 1_g / square(1_cm);
+    auto const dE = -BetheBloch(vParticle, dX); // dE > 0
+    //~ auto const Ekin = vParticle.GetEnergy() - vParticle.GetMass();
+    auto const maxLoss = 0.01 * vParticle.GetEnergy();
+    auto const maxGrammage = maxLoss / dE * dX;
+
+    return vParticle.GetNode()->GetModelProperties().ArclengthFromGrammage(vTrack,
+                                                                           maxGrammage) *
+           1.0001; // to make sure particle gets absorbed when DoContinuous() is called
   }
 
   void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& vP,
@@ -187,17 +200,17 @@ namespace corsika::process::energy_loss {
 
 #include <corsika/geometry/CoordinateSystem.h>
 
-  int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& vP,
+  int EnergyLoss::GetXbin(SetupParticle const& vP, SetupTrack const& vTrack,
                           const HEPEnergyType dE) {
 
     using namespace corsika::geometry;
 
     CoordinateSystem const& rootCS =
         RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-    Point pos1(rootCS, 0_m, 0_m, 0_m);
-    Point pos2(rootCS, 0_m, 0_m, vP.GetPosition().GetCoordinates()[2]);
-    Vector delta = (pos2 - pos1) / 1_s;
-    Trajectory t(Line(pos1, delta), 1_s);
+    Point const pos1(rootCS, 0_m, 0_m, 0_m);
+    Point const pos2(rootCS, 0_m, 0_m, vTrack.GetPosition(0).GetCoordinates()[2]);
+    auto const delta = (pos2 - pos1) / 1_s;
+    Trajectory const t(Line(pos1, delta), 1_s);
 
     GrammageType const grammage =
         vP.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());

Processes/EnergyLoss/EnergyLoss.h

@@ -25,34 +25,30 @@ namespace corsika::process::energy_loss {
   class EnergyLoss : public corsika::process::ContinuousProcess<EnergyLoss> {
 
     using MeVgcm2 = decltype(1e6 * units::si::electronvolt / units::si::gram *
-                             corsika::units::si::square(1e-2 * units::si::meter));
+                             units::si::square(1e-2 * units::si::meter));
 
-    void MomentumUpdate(corsika::setup::Stack::ParticleType&,
-                        corsika::units::si::HEPEnergyType Enew);
+    void MomentumUpdate(setup::Stack::ParticleType&, units::si::HEPEnergyType Enew);
 
   public:
     EnergyLoss();
     void Init() {}
+    process::EProcessReturn DoContinuous(setup::Stack::ParticleType&,
+                                         setup::Trajectory const&);
+    units::si::LengthType MaxStepLength(setup::Stack::ParticleType const&,
+                                        setup::Trajectory const&) const;
 
-    corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
-                                                  corsika::setup::Trajectory&);
-    corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&,
-                                                 corsika::setup::Trajectory&);
-
-    corsika::units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
+    units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
     void PrintProfile() const;
 
   private:
-    corsika::units::si::HEPEnergyType BetheBloch(
-        corsika::setup::Stack::ParticleType& p,
-        const corsika::units::si::GrammageType dX);
+    static units::si::HEPEnergyType BetheBloch(setup::Stack::ParticleType const& p,
+                                               const units::si::GrammageType dX);
 
-    int GetXbin(corsika::setup::Stack::ParticleType& p,
-                const corsika::units::si::HEPEnergyType dE);
+    int GetXbin(setup::Stack::ParticleType const&, setup::Trajectory const&, units::si::HEPEnergyType);
 
-    corsika::units::si::HEPEnergyType fEnergyLossTot;
-    corsika::units::si::GrammageType fdX; // profile binning
-    std::map<int, double> fProfile;       // longitudinal profile
+    units::si::HEPEnergyType fEnergyLossTot;
+    units::si::GrammageType fdX;    // profile binning
+    std::map<int, double> fProfile; // longitudinal profile
   };
 
 } // namespace corsika::process::energy_loss

Processes/Pythia/Decay.cc

@@ -148,9 +148,6 @@ namespace corsika::process::pythia {
     // set particle stable
     Decay::SetStable(vP.GetPID());
 
-    // remove original particle from corsika stack
-    vP.Delete();
-    //    if (fCount>10) throw std::runtime_error("stop here");
   }
 
 } // namespace corsika::process::pythia

Processes/Pythia/Interaction.cc

@@ -411,8 +411,6 @@ namespace corsika::process::pythia {
              << "Elab_final=" << Elab_final / 1_GeV
              << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
       }
-      // delete current particle
-      vP.Delete();
     }
     return process::EProcessReturn::eOk;
   }

Processes/Pythia/Random.cc

@@ -12,9 +12,6 @@
 
 namespace corsika::process::pythia {
 
-  double Random::flat() {
-    std::uniform_real_distribution<double> dist;
-    return dist(fRNG);
-  }
+  double Random::flat() { return fDist(fRNG); }
 
 } // namespace corsika::process::pythia

Processes/Pythia/Random.h

@@ -23,6 +23,7 @@ namespace corsika::process {
       double flat();
 
     private:
+      std::uniform_real_distribution<double> fDist;
       corsika::random::RNG& fRNG =
           corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
     };

Processes/Sibyll/Interaction.cc

@@ -366,8 +366,6 @@ namespace corsika::process::sibyll {
              << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
       }
     }
-    // delete current particle
-    vP.Delete();
     return process::EProcessReturn::eOk;
   }
 

Processes/StackInspector/StackInspector.cc

@@ -37,7 +37,7 @@ template <typename TStack>
 StackInspector<TStack>::~StackInspector() {}
 
 template <typename TStack>
-process::EProcessReturn StackInspector<TStack>::DoStack(TStack& vS) {
+process::EProcessReturn StackInspector<TStack>::DoStack(TStack const& vS) {
   if (!fReport) return process::EProcessReturn::eOk;
   [[maybe_unused]] int i = 0;
   HEPEnergyType Etot = 0_GeV;

Processes/StackInspector/StackInspector.h

@@ -30,7 +30,7 @@ namespace corsika::process {
       ~StackInspector();
 
       void Init();
-      EProcessReturn DoStack(TStack&);
+      EProcessReturn DoStack(TStack const&);
 
     private:
       bool fReport;

Processes/TrackWriter/TrackWriter.cc

@@ -37,9 +37,9 @@ namespace corsika::process::track_writer {
     using namespace units::si;
     auto const start = vT.GetPosition(0).GetCoordinates();
     auto const delta = vT.GetPosition(1).GetCoordinates() - start;
-    auto const& name = particles::GetName(vP.GetPID());
+    auto const pdg = static_cast<int>(particles::GetPDG(vP.GetPID()));
 
-    fFile << name << "    " << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
+    fFile << pdg << ' ' << vP.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';
 

Processes/TrackingLine/TrackingLine.cc

@@ -16,31 +16,27 @@
 #include <corsika/geometry/Vector.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
 
-#include <algorithm>
-#include <iostream>
+#include <limits>
 #include <stdexcept>
 #include <utility>
 
-using namespace corsika;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
 
 namespace corsika::process::tracking_line {
 
-  std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
-  TimeOfIntersection(corsika::geometry::Line const& line,
-                     geometry::Sphere const& sphere) {
+  std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const& line,
+                                                                  Sphere const& sphere) {
     auto const delta = line.GetR0() - sphere.GetCenter();
     auto const v = line.GetV0();
-    auto const vSqNorm = v.squaredNorm();
+    auto const vSqNorm =
+        v.squaredNorm(); // todo: get rid of this by having V0 normalized always
     auto const R = sphere.GetRadius();
 
     auto const vDotDelta = v.dot(delta);
     auto const discriminant =
         vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
 
-    //~ std::cout << "discriminant: " << discriminant << std::endl;
-    //~ std::cout << "alpha: " << alpha << std::endl;
-    //~ std::cout << "beta: " << beta << std::endl;
-
     if (discriminant.magnitude() > 0) {
       auto const sqDisc = sqrt(discriminant);
       auto const invDenom = 1 / vSqNorm;
@@ -50,4 +46,17 @@ namespace corsika::process::tracking_line {
       return {};
     }
   }
+
+  TimeType TimeOfIntersection(Line const& vLine, Plane const& vPlane) {
+    auto const delta = vPlane.GetCenter() - vLine.GetR0();
+    auto const v = vLine.GetV0();
+    auto const n = vPlane.GetNormal();
+    auto const c = n.dot(v);
+
+    if (c.magnitude() == 0) {
+      return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
+    } else {
+      return n.dot(delta) / c;
+    }
+  }
 } // namespace corsika::process::tracking_line

Processes/TrackingLine/TrackingLine.h

@@ -13,6 +13,7 @@
 #define _include_corsika_processes_TrackingLine_h_
 
 #include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/geometry/Vector.h>
@@ -33,8 +34,10 @@ namespace corsika::process {
   namespace tracking_line {
 
     std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
-    TimeOfIntersection(corsika::geometry::Line const& line,
-                       geometry::Sphere const& sphere);
+    TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
+
+    corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
+                                                    geometry::Plane const&);
 
     class TrackingLine {
 

README.md

@@ -60,7 +60,7 @@ CORSIKA 8 is tested regularly at least on gcc7.3.0 and clang-6.0.0.
 Additional software prerequisites: eigen3, boost, cmake, g++, git. 
 On a bare Ubuntu 18.04, just add:
 ```
-sudo apt-get install libeigen3-dev libboost-dev cmake g++ git
+sudo apt install libeigen3-dev cmake g++ git
 ```
 
 Follow these steps to download and install CORSIKA 8 milestone2