diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index c12c79b01d9663c8caf08a8068ec0dc13b529e3c..510d8e3922f2c7c9a1bbd0578bea34f74641e500 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -10,6 +10,7 @@ set (
Point.h
Line.h
Sphere.h
+ Plane.h
Volume.h
CoordinateSystem.h
RootCoordinateSystem.h
diff --git a/Framework/Geometry/Plane.h b/Framework/Geometry/Plane.h
new file mode 100644
index 0000000000000000000000000000000000000000..2c3f103f52bb9278ffcca63766512ceda13d32c2
--- /dev/null
+++ b/Framework/Geometry/Plane.h
@@ -0,0 +1,41 @@
+
+/*
+ * (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
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 58d378aa3eb5bdad5f8029555ea8825afa450007..27f6c8fae4cb52e037f7ac128dd5fde161dd5566 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/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
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 75dfc51cedd83bd266f07f971aef48ca717b7d81..78f9ecd35a5872ccc05ddebab3dd5f6873a65d1f 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/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.
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 990d6172501dc67f0df4958cf651d587003f1923..aaa6dfd6146deb0f34cfbdd79f7e78d76aafa4be 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/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());
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index deaff4faf0837f6cac44ce2da975637393035d49..711e135f716099fc3204a9be14239239a97a3016 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/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
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 222c68c8dd57d7fd57c453b62f53777e85ad062a..670d24196423606bbce79329ce95272630412633 100644
--- a/Processes/Pythia/Decay.cc
+++ b/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
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index ae8104350086f5cfda94b174305d0f0366faa2ce..6cdf54973a2d553a79f6b6ebcd97b4d5ecf6cf24 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/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;
}
diff --git a/Processes/Pythia/Random.cc b/Processes/Pythia/Random.cc
index 8fcbcf4ee7274ba36252f7969306073f3cffddd5..0d08403ca1d890a062e2632124079c6a3cb220f6 100644
--- a/Processes/Pythia/Random.cc
+++ b/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
diff --git a/Processes/Pythia/Random.h b/Processes/Pythia/Random.h
index 276fc532ca53ac84b708e8ad12cb16fbaed01c9a..cd35fc4eeafcbdf4ec507cf54ab966271a8183f0 100644
--- a/Processes/Pythia/Random.h
+++ b/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");
};
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index d674f9a3a775970bf60de08c0cf33b55fc73a17d..b15d9e8859e86d6426f6b78b7834fe2211f191f1 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/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;
}
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 814763d5cbfcbc0bd4201db70f220dc4e8ea123f..e56503378ed5abfc15578c474ed0afd4d54d6412 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/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;
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 57a929dd2362b302d43d05175878a257c5c8a832..62f36eeebae774f74d26953697707e698c75da6e 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -30,7 +30,7 @@ namespace corsika::process {
~StackInspector();
void Init();
- EProcessReturn DoStack(TStack&);
+ EProcessReturn DoStack(TStack const&);
private:
bool fReport;
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index e422bce8af8f054819cd03c174eaf594e85d6c76..d3247e87c8ad0e66620afcad2b2dc18d3e1c6049 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/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';
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 319030f91d5347d47af58cddbdea607684a16677..149556bf245b9fb418fb3db03f267e64b688fb5f 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/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
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index a801a2ee2baa206424211151c6a7fb9a616d0b70..df02264fe6f4c5640fe40c799255e87ce26a958e 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/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 {
diff --git a/README.md b/README.md
index 67dd505f0d4dafc3f5ee73e055df7e4754429c99..d5cf8f7c2764daf2fe25f5c861502318a7ea4026 100644
--- a/README.md
+++ b/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