From bb15d6525bb2eb72e6a36134c5bce1ffd34dd56b Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Thu, 14 Feb 2019 10:42:43 +0100
Subject: [PATCH] first simple version of energy loss
---
Documentation/Examples/CMakeLists.txt | 8 +-
Documentation/Examples/cascade_example.cc | 19 +--
Environment/DensityFunction.h | 3 +-
Environment/InhomogeneousMedium.h | 3 +-
Environment/LinearApproximationIntegrator.h | 3 +-
Framework/Geometry/Trajectory.h | 11 +-
Framework/Geometry/Vector.h | 2 +
Framework/Geometry/testGeometry.cc | 6 +-
Framework/Particles/ParticleProperties.h | 23 ++--
Framework/Particles/pdxml_reader.py | 4 +-
Framework/Particles/testParticles.cc | 2 +-
Framework/Units/PhysicalUnits.h | 4 +-
Processes/CMakeLists.txt | 2 +
Processes/EnergyLoss/CMakeLists.txt | 64 +++++++++++
Processes/EnergyLoss/EnergyLoss.cc | 108 ++++++++++++++++++
Processes/EnergyLoss/EnergyLoss.h | 57 +++++++++
.../HadronicElasticModel.h | 4 +-
Processes/Sibyll/Interaction.cc | 4 +-
Processes/Sibyll/NuclearInteraction.cc | 3 +-
.../NuclearStackExtension.h | 19 ++-
Stack/SuperStupidStack/SuperStupidStack.h | 12 ++
21 files changed, 324 insertions(+), 37 deletions(-)
create mode 100644 Processes/EnergyLoss/CMakeLists.txt
create mode 100644 Processes/EnergyLoss/EnergyLoss.cc
create mode 100644 Processes/EnergyLoss/EnergyLoss.h
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 621cab01..f8f1152a 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -24,14 +24,16 @@ CORSIKA_ADD_TEST (stack_example)
add_executable (cascade_example cascade_example.cc)
target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
- CORSIKArandom
- ProcessSibyll
- CORSIKAcascade
+ CORSIKArandom
+# CORSIKAprocesses
+ ProcessSibyll
+ ProcessEnergyLoss
ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
CORSIKAprocesses
+ CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f9d1e95a..a8750538 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -11,6 +11,7 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
@@ -262,10 +263,12 @@ int main() {
// hadronicElastic(env);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+ process::EnergyLoss::EnergyLoss eLoss(2_MeV / 1_g * square(1_cm));
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
- auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
+ auto sequence = p0 << sibyll << sibyllNuc << decay << eLoss << cut << trackWriter;
+ // auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
@@ -276,17 +279,16 @@ int main() {
const Code beamCode = Code::Nucleus;
const int nuclA = 56;
const int nuclZ = int(nuclA / 2.15 + 0.7);
- const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
- (nuclA - nuclZ) * particles::Neutron::GetMass();
+ const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
const HEPEnergyType E0 =
nuclA *
- 100_GeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
+ 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt(Elab * Elab - m * m);
+ return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
@@ -315,6 +317,9 @@ int main() {
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
- cout << "total energy (GeV): " << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
+ << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+ eLoss.SaveSave();
}
diff --git a/Environment/DensityFunction.h b/Environment/DensityFunction.h
index 32ee30ad..92ba841b 100644
--- a/Environment/DensityFunction.h
+++ b/Environment/DensityFunction.h
@@ -1,4 +1,5 @@
-/**
+
+/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h
index feb6497c..0b1a0c1b 100644
--- a/Environment/InhomogeneousMedium.h
+++ b/Environment/InhomogeneousMedium.h
@@ -1,4 +1,5 @@
-/**
+
+/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
diff --git a/Environment/LinearApproximationIntegrator.h b/Environment/LinearApproximationIntegrator.h
index 910d2397..90e14d99 100644
--- a/Environment/LinearApproximationIntegrator.h
+++ b/Environment/LinearApproximationIntegrator.h
@@ -1,4 +1,5 @@
-/**
+
+/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 32063061..e8a1a58b 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -37,20 +37,21 @@ namespace corsika::geometry {
Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
corsika::units::si::TimeType GetDuration() const { return fTimeLength; }
+ corsika::units::si::LengthType GetLength() const { return GetDistance(fTimeLength); }
corsika::units::si::LengthType GetDistance(corsika::units::si::TimeType t) const {
- assert(t > fTimeLength);
+ assert(t <= fTimeLength);
assert(t >= 0 * corsika::units::si::second);
- return T::ArcLength(0, t);
+ return T::ArcLength(0 * corsika::units::si::second, t);
}
void LimitEndTo(corsika::units::si::LengthType limit) {
fTimeLength = T::TimeFromArclength(limit);
}
-
+
auto NormalizedDirection() const {
- static_assert(std::is_same_v<T, corsika::geometry::Line>);
- return T::GetV0().normalized();
+ static_assert(std::is_same_v<T, corsika::geometry::Line>);
+ return T::GetV0().normalized();
}
};
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 5ac1f118..737b7830 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -75,12 +75,14 @@ namespace corsika::geometry {
* think about whether squaredNorm() might be cheaper for your computation.
*/
auto norm() const { return BaseVector<dim>::qVector.norm(); }
+ auto GetNorm() const { return BaseVector<dim>::qVector.norm(); }
/*!
* returns the squared norm of the Vector. Before using this method,
* think about whether norm() might be cheaper for your computation.
*/
auto squaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); }
+ auto GetSquaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); }
/*!
* returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index d73ae99b..dae272be 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -177,8 +177,10 @@ TEST_CASE("Trajectories") {
CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
-
- CHECK((base.NormalizedDirection().GetComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}).eVector.norm() == Approx(0).margin(absMargin));
+
+ CHECK((base.NormalizedDirection().GetComponents(rootCS) -
+ QuantityVector<dimensionless_d>{1, 0, 0})
+ .eVector.norm() == Approx(0).margin(absMargin));
}
SECTION("Helix") {
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index dcbbe639..456711b4 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -44,8 +44,8 @@ namespace corsika::particles {
using PDGCodeType = std::underlying_type<PDGCode>::type;
// forward declarations to be used in GeneratedParticleProperties
- int16_t constexpr GetElectricChargeNumber(Code const);
- corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const);
+ int16_t constexpr GetChargeNumber(Code const);
+ corsika::units::si::ElectricChargeType constexpr GetCharge(Code const);
corsika::units::si::HEPMassType constexpr GetMass(Code const);
PDGCode constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
@@ -72,17 +72,18 @@ namespace corsika::particles {
}
/*!
- * returns electric charge of particle / (e/3), e.g. return 3 for a proton.
+ * returns electric charge number of particle return 1 for a proton.
*/
- int16_t constexpr GetElectricChargeNumber(Code const p) {
- return detail::electric_charges[static_cast<CodeIntType>(p)];
+ int16_t constexpr GetChargeNumber(Code const p) {
+ // electric_charges stores charges in units of (e/3), e.g. 3 for a proton
+ return detail::electric_charges[static_cast<CodeIntType>(p)] / 3;
}
/*!
* returns electric charge of particle, e.g. return 1.602e-19_C for a proton.
*/
- corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
- return GetElectricChargeNumber(p) * (corsika::units::constants::e * (1. / 3.));
+ corsika::units::si::ElectricChargeType constexpr GetCharge(Code const p) {
+ return GetChargeNumber(p) * (corsika::units::constants::e);
}
constexpr std::string const& GetName(Code const p) {
@@ -112,6 +113,14 @@ namespace corsika::particles {
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p);
Code ConvertFromPDG(PDGCode);
+
+ /**
+ * Get mass of nucleus
+ **/
+ corsika::units::si::HEPMassType constexpr GetNucleusMass(const int vA, const int vZ) {
+ return Proton::GetMass() * vZ + (vA - vZ) * Neutron::GetMass();
+ }
+
} // namespace corsika::particles
#endif
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index 946940f7..84ef0ee9 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -418,8 +418,8 @@ def gen_classes(particle_db):
string += " public:\n"
string += " static constexpr Code GetCode() { return Type; }\n"
string += " static constexpr corsika::units::si::HEPMassType GetMass() { return corsika::particles::GetMass(Type); }\n"
- string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(Type); }\n"
- string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n"
+ string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetCharge(Type); }\n"
+ string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetChargeNumber(Type); }\n"
string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n"
string += " static constexpr Code GetAntiParticle() { return AntiType; }\n"
string += " static constexpr bool IsNucleus() { return corsika::particles::IsNucleus(Type); }\n"
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 8bb73e7a..3fd99ea2 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -43,7 +43,7 @@ TEST_CASE("ParticleProperties", "[Particles]") {
SECTION("Charges") {
REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
- REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
+ REQUIRE(GetCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
}
SECTION("Names") {
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 786044ae..6bc71d70 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -45,9 +45,11 @@ namespace corsika::units::si {
using hepmomentum_d = phys::units::hepenergy_d;
using hepmass_d = phys::units::hepenergy_d;
- /// defining cross section
+ /// defining cross section as area
using sigma_d = phys::units::area_d;
+ // constexpr quantity<area_d> barn{Rep(1e-28L) * square(meter)};
+
/// add the unit-types
using LengthType = phys::units::quantity<phys::units::length_d, double>;
using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 26a44fe6..e23f660b 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -5,6 +5,7 @@ add_subdirectory (StackInspector)
add_subdirectory (TrackWriter)
add_subdirectory (TrackingLine)
add_subdirectory (HadronicElasticModel)
+add_subdirectory (EnergyLoss)
#add_custom_target(CORSIKAprocesses)
add_library (CORSIKAprocesses INTERFACE)
@@ -12,3 +13,4 @@ add_dependencies(CORSIKAprocesses ProcessNullModel)
add_dependencies(CORSIKAprocesses ProcessSibyll)
add_dependencies(CORSIKAprocesses ProcessStackInspector)
add_dependencies(CORSIKAprocesses ProcessTrackingLine)
+add_dependencies(CORSIKAprocesses ProcessEnergyLoss)
diff --git a/Processes/EnergyLoss/CMakeLists.txt b/Processes/EnergyLoss/CMakeLists.txt
new file mode 100644
index 00000000..f108ba85
--- /dev/null
+++ b/Processes/EnergyLoss/CMakeLists.txt
@@ -0,0 +1,64 @@
+
+set (
+ MODEL_SOURCES
+ EnergyLoss.cc
+ )
+
+set (
+ MODEL_HEADERS
+ EnergyLoss.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/energy_loss
+ )
+
+add_library (ProcessEnergyLoss STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessEnergyLoss ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+ ProcessEnergyLoss
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+# PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessEnergyLoss
+ CORSIKAunits
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAsetup
+ )
+
+target_include_directories (
+ ProcessEnergyLoss
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (
+ TARGETS ProcessEnergyLoss
+ 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
+# )
+# CORSIKA_ADD_TEST(testNullModel)
+
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
new file mode 100644
index 00000000..8646b4d5
--- /dev/null
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -0,0 +1,108 @@
+
+/*
+ * (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/energy_loss/EnergyLoss.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+using namespace std;
+
+using namespace corsika;
+using namespace corsika::units::si;
+using namespace corsika::setup;
+using Particle = Stack::ParticleType;
+using Track = Trajectory;
+
+namespace corsika::process::EnergyLoss {
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+
+ EnergyLoss::EnergyLoss(MeVgcm2 const vdEdX)
+ : fdEdX(vdEdX)
+ , fEnergyLossTot(0_GeV) {}
+
+ process::EProcessReturn EnergyLoss::DoContinuous(Particle& p, Track& t, Stack&) {
+ GrammageType const dX =
+ p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
+ HEPEnergyType dE = -dX * fdEdX * pow(p.GetChargeNumber(), 2);
+ auto E = p.GetEnergy();
+ const auto Ekin = E - p.GetParticleMass();
+ auto Enew = E + dE;
+ cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
+ << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2, dE=" << dE / 1_MeV << "MeV, "
+ << " E=" << E / 1_GeV << "GeV, Ekin=" << Ekin / 1_GeV
+ << ", Enew=" << Enew / 1_GeV << "GeV" << endl;
+ auto status = process::EProcessReturn::eOk;
+ if (-dE > Ekin) {
+ dE = -Ekin;
+ Enew = p.GetParticleMass();
+ status = process::EProcessReturn::eParticleAbsorbed;
+ }
+ p.SetEnergy(Enew);
+ MomentumUpdate(p, Enew);
+ fEnergyLossTot += dE;
+ GetXbin(p, dE);
+ return status;
+ }
+
+ units::si::LengthType EnergyLoss::MaxStepLength(Particle&, Track&) {
+ return units::si::meter * std::numeric_limits<double>::infinity();
+ }
+
+ void EnergyLoss::MomentumUpdate(corsika::setup::Stack::ParticleType& p,
+ corsika::units::si::HEPEnergyType Enew) {
+ HEPMomentumType Pnew = elab2plab(Enew, p.GetParticleMass());
+ auto pnew = p.GetMomentum();
+ p.SetMomentum(pnew * Pnew / pnew.GetNorm());
+ }
+
+#include <corsika/geometry/CoordinateSystem.h>
+
+ int EnergyLoss::GetXbin(corsika::setup::Stack::ParticleType& p,
+ const HEPEnergyType dE) {
+
+ using namespace corsika::geometry;
+
+ const GrammageType deltaX = 10_g / square(1_cm); // binning
+
+ CoordinateSystem const& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ Point pos1(rootCS, 0_m, 0_m, 0_m);
+ Point pos2(rootCS, 0_m, 0_m, p.GetPosition().GetCoordinates()[2]);
+ Vector delta = (pos2 - pos1) / 1_s;
+ Trajectory t(Line(pos1, delta), 1_s);
+
+ GrammageType const grammage =
+ p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
+
+ const int bin = grammage / deltaX;
+
+ if (!fSave.count(bin)) { cout << "EnergyLoss new x bin " << bin << endl; }
+ fSave[bin] += -dE / 1_GeV;
+ return bin;
+ }
+
+ void EnergyLoss::SaveSave() {
+
+ cout << "EnergyLoss Save " << endl;
+ for (auto v : fSave) { cout << v.first << " " << v.second << endl; }
+ }
+
+} // namespace corsika::process::EnergyLoss
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
new file mode 100644
index 00000000..3af10c01
--- /dev/null
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -0,0 +1,57 @@
+
+/*
+ * (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_EnergyLoss_h_
+#define _Processes_EnergyLoss_h_
+
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <map>
+
+namespace corsika::process::EnergyLoss {
+
+ 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));
+
+ void MomentumUpdate(corsika::setup::Stack::ParticleType&,
+ corsika::units::si::HEPEnergyType Enew);
+
+ public:
+ EnergyLoss(MeVgcm2 const vdEdX);
+ void Init() {}
+
+ corsika::process::EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType&,
+ corsika::setup::Trajectory&,
+ corsika::setup::Stack&);
+ corsika::units::si::LengthType MaxStepLength(corsika::setup::Stack::ParticleType&,
+ corsika::setup::Trajectory&);
+
+ corsika::units::si::HEPEnergyType GetTotal() const { return fEnergyLossTot; }
+ void SaveSave();
+
+ private:
+ int GetXbin(corsika::setup::Stack::ParticleType& p,
+ const corsika::units::si::HEPEnergyType dE);
+
+ MeVgcm2 fdEdX;
+ corsika::units::si::HEPEnergyType fEnergyLossTot;
+ std::map<int, double> fSave;
+ };
+
+} // namespace corsika::process::EnergyLoss
+
+#endif
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index f0396387..7f37905d 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -44,8 +44,8 @@ namespace corsika::process::HadronicElasticModel {
using SquaredHEPEnergyType = decltype(corsika::units::si::HEPEnergyType() *
corsika::units::si::HEPEnergyType());
- using eV2 = decltype(units::si::detail::static_pow<2>(units::si::electronvolt));
- using inveV2 = decltype(units::si::detail::static_pow<-2>(units::si::electronvolt));
+ using eV2 = decltype(units::si::square(units::si::electronvolt));
+ using inveV2 = decltype(1 / units::si::square(units::si::electronvolt));
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream(
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 5ccc1b74..3c27e4d9 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -277,8 +277,8 @@ namespace corsika::process::sibyll {
sigEla; // to avoid not used warning in array binding
}
- const auto targetCode = mediumComposition.SampleTarget(
- cross_section_of_components, fRNG);
+ const auto targetCode =
+ mediumComposition.SampleTarget(cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 6a373632..575903e0 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -356,7 +356,8 @@ namespace corsika::process::sibyll {
[[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
}
- const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG);
+ const auto targetCode =
+ mediumComposition.SampleTarget(cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
index 11923aba..d6d54523 100644
--- a/Stack/NuclearStackExtension/NuclearStackExtension.h
+++ b/Stack/NuclearStackExtension/NuclearStackExtension.h
@@ -159,7 +159,24 @@ namespace corsika::stack {
int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); }
/// @}
- // int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
+ /**
+ * Overwrite normal GetParticleMass function with nuclear version
+ */
+ corsika::units::si::HEPMassType GetParticleMass() const {
+ if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
+ corsika::particles::Code::Nucleus)
+ return corsika::particles::GetNucleusMass(GetNuclearA(), GetNuclearZ());
+ return InnerParticleInterface<StackIteratorInterface>::GetParticleMass();
+ }
+ /**
+ * Overwirte normal GetChargeNumber function with nuclear version
+ **/
+ int16_t GetChargeNumber() const {
+ if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
+ corsika::particles::Code::Nucleus)
+ return GetNuclearZ();
+ return InnerParticleInterface<StackIteratorInterface>::GetChargeNumber();
+ }
protected:
void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 365a45aa..f28b5132 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -117,10 +117,22 @@ namespace corsika::stack {
corsika::units::si::TimeType GetTime() const {
return GetStackData().GetTime(GetIndex());
}
+ /**
+ * @name derived quantities
+ *
+ * @{
+ */
corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection()
const {
return GetMomentum() / GetEnergy();
}
+ corsika::units::si::HEPMassType GetParticleMass() const {
+ return corsika::particles::GetMass(GetPID());
+ }
+ int16_t GetChargeNumber() const {
+ return corsika::particles::GetChargeNumber(GetPID());
+ }
+ ///@}
};
/**
--
GitLab