From c2725ac596beb4fee4ed1811845b5fd44162a3ce Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 18 Dec 2018 10:19:13 +0100
Subject: [PATCH] fixed some namespaces, use HEP energy units, changed
process-sequence operator
---
CMakeLists.txt | 2 +-
Documentation/Examples/cascade_example.cc | 9 +-
.../Examples/staticsequence_example.cc | 2 +-
Framework/Cascade/Cascade.h | 12 +-
Framework/Cascade/testCascade.cc | 6 +-
Framework/Geometry/QuantityVector.h | 18 +-
Framework/Geometry/Trajectory.h | 2 +-
Framework/Particles/ParticleProperties.cc | 9 +-
Framework/Particles/ParticleProperties.h | 17 +-
Framework/Particles/pdxml_reader.py | 2 +-
Framework/Particles/testParticles.cc | 3 +-
Framework/ProcessSequence/BaseProcess.h | 32 ----
Framework/ProcessSequence/ContinuousProcess.h | 1 -
Framework/ProcessSequence/ProcessReturn.h | 1 +
Framework/ProcessSequence/ProcessSequence.h | 113 +------------
.../ProcessSequence/testProcessSequence.cc | 8 +-
Framework/StackInterface/ParticleBase.h | 2 -
Framework/StackInterface/StackIterator.h | 3 -
Framework/Units/PhysicalConstants.h | 9 +-
Framework/Units/PhysicalUnits.h | 43 +++--
Framework/Units/testUnits.cc | 8 +-
Processes/Sibyll/Decay.h | 155 +++++++++---------
Processes/Sibyll/Interaction.h | 95 ++++++-----
Processes/Sibyll/SibStack.h | 64 +++++---
Processes/StackInspector/StackInspector.cc | 6 +-
Processes/TrackingLine/TrackingLine.h | 2 +-
Processes/TrackingLine/testTrackingLine.cc | 13 +-
Setup/SetupTrajectory.h | 7 +-
Stack/SuperStupidStack/SuperStupidStack.h | 100 ++++++-----
.../SuperStupidStack/testSuperStupidStack.cc | 3 +-
30 files changed, 339 insertions(+), 408 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 841f3a9d..567dd5ce 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -37,7 +37,7 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
endif()
# enable warnings and disallow non-standard language
-set(CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers")
+set(CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers -Wno-nonportable-include-path")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -O0 -g")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3")
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f2326e06..e221b28f 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -43,7 +43,7 @@ using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
-using namespace corsika::units::si;
+using namespace corsika::units::hep;
static EnergyType fEnergy = 0. * 1_GeV;
@@ -228,7 +228,7 @@ int main() {
corsika::process::sibyll::Interaction sibyll;
corsika::process::sibyll::Decay decay;
ProcessEMCut cut;
- const auto sequence = /*p0 +*/ sibyll + decay + cut;
+ const auto sequence = p0 << sibyll << decay << cut;
setup::Stack stack;
@@ -237,9 +237,8 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
- MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV) / si::constants::c;
- auto plab = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c,
- 0. * 1_GeV / si::constants::c, P0);
+ hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
+ auto plab = stack::super_stupid::MomentumVector(rootCS, 0. * 1_GeV, 0. * 1_GeV, P0);
particle.SetEnergy(E0);
particle.SetMomentum(plab);
particle.SetPID(Code::Proton);
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index 522c1711..e2836ff1 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/Documentation/Examples/staticsequence_example.cc
@@ -86,7 +86,7 @@ void modular() {
Process3 m3;
Process4 m4(0.9);
- const auto sequence = m1 + m2 + m3 + m4;
+ const auto sequence = m1 << m2 << m3 << m4;
DummyData p;
DummyStack s;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 0edd5050..a94de8bb 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -19,7 +19,8 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
-#include <type_traits>
+#include <cmath>
+#include <iostream>
namespace corsika::cascade {
@@ -56,7 +57,14 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
- using namespace corsika::units::si;
+
+ using std::cout;
+ using std::endl;
+ using std::log;
+
+ // get access to random number generator
+ static corsika::random::RNG& rmng =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 3d7d0741..80ebf4b2 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -113,7 +113,7 @@ TEST_CASE("Cascade", "[Cascade]") {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
- const auto sequence = p0 + p1;
+ const auto sequence = p0 << p1;
setup::Stack stack;
corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
@@ -126,8 +126,8 @@ TEST_CASE("Cascade", "[Cascade]") {
particle.SetPID(particles::Code::Electron);
particle.SetEnergy(E0);
particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
- particle.SetMomentum(corsika::stack::super_stupid::MomentumVector(
- rootCS, {0 * newton * second, 0 * newton * second, -1 * newton * second}));
+ particle.SetMomentum(
+ corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}));
particle.SetTime(0_ns);
EAS.Init();
EAS.Run();
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index c128039d..bb498f08 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -117,16 +117,16 @@ namespace corsika::geometry {
auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
};
-} // namespace corsika::geometry
+ template <typename dim>
+ auto& operator<<(std::ostream& os, corsika::geometry::QuantityVector<dim> qv) {
+ using Quantity = phys::units::quantity<dim, double>;
-template <typename dim>
-auto& operator<<(std::ostream& os, corsika::geometry::QuantityVector<dim> qv) {
- using Quantity = phys::units::quantity<dim, double>;
+ os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
+ << phys::units::to_unit_symbol<dim, double>(
+ Quantity(phys::units::detail::magnitude_tag, 1));
+ return os;
+ }
- os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
- << phys::units::to_unit_symbol<dim, double>(
- Quantity(phys::units::detail::magnitude_tag, 1));
- return os;
-}
+} // namespace corsika::geometry
#endif
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 0d142d91..5a9d1ed5 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -38,7 +38,7 @@ namespace corsika::geometry {
corsika::units::si::TimeType GetDuration() const { return fTimeLength; }
- LengthType GetDistance(corsika::units::si::TimeType t) const {
+ corsika::units::si::LengthType GetDistance(corsika::units::si::TimeType t) const {
assert(t > fTimeLength);
assert(t >= 0 * corsika::units::si::second);
return T::ArcLength(0, t);
diff --git a/Framework/Particles/ParticleProperties.cc b/Framework/Particles/ParticleProperties.cc
index 6c823056..2fff766f 100644
--- a/Framework/Particles/ParticleProperties.cc
+++ b/Framework/Particles/ParticleProperties.cc
@@ -10,11 +10,12 @@
*/
#include <corsika/particles/ParticleProperties.h>
+#include <iostream>
-namespace corsika::particles::io {
+namespace corsika::particles {
- std::ostream& operator<<(std::ostream& stream, Code const p) {
- return stream << GetName(p);
+ std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p) {
+ return stream << corsika::particles::GetName(p);
}
-} // namespace corsika::particles::io
+} // namespace corsika::particles
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index db3f10fe..44b9b55e 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -20,7 +20,7 @@
#include <array>
#include <cstdint>
-#include <iostream>
+#include <iosfwd>
#include <type_traits>
#include <corsika/units/PhysicalConstants.h>
@@ -36,8 +36,6 @@
namespace corsika::particles {
- using corsika::units::si::second;
-
enum class Code : int16_t;
using PDGCodeType = int32_t;
@@ -76,7 +74,7 @@ namespace corsika::particles {
}
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
- return GetElectricChargeNumber(p) * (corsika::units::si::constants::e / 3.);
+ return GetElectricChargeNumber(p) * (corsika::units::constants::e / 3.);
}
constexpr std::string const& GetName(Code const p) {
@@ -100,15 +98,12 @@ namespace corsika::particles {
return detail::nucleusZ[static_cast<CodeIntType const>(p)];
}
- namespace io {
-
- std::ostream& operator<<(std::ostream& stream, Code const p);
+ /**
+ * the output operator for particles
+ **/
- } // namespace io
+ std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p);
} // namespace corsika::particles
-// to inject the operator<< into the root namespace
-using namespace corsika::particles::io;
-
#endif
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index ca7701d4..b1a11ae8 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -292,7 +292,7 @@ def gen_properties(particle_db):
# particle masses table
string += "static constexpr std::array<corsika::units::hep::MassType const, size> masses = {\n"
for p in particle_db.values():
- string += " {mass:e} * (1e9 * corsika::units::si::constants::eV), // {name:s}\n".format(mass = p['mass'], name = p['name'])
+ string += " {mass:e} * (1e9 * corsika::units::constants::eV), // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
# PDG code table
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 219e8dd8..f2112605 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -16,7 +16,8 @@
// cpp file
#include <catch2/catch.hpp>
-using namespace corsika::units::si;
+using namespace corsika::units;
+using namespace corsika::units::hep;
using namespace corsika::particles;
TEST_CASE("ParticleProperties", "[Particles]") {
diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 9cc29196..aa82c0b5 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -29,39 +29,7 @@ namespace corsika::process {
struct BaseProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
- /*
- template <typename Particle, typename Stack>
- inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
-
- template <typename Particle, typename Track, typename Stack>
- inline EProcessReturn DoContinuous(Particle&, Track&, Stack&) const; // {}
- */
- /*
- template <typename Particle, typename Track>
- inline double GetInverseInteractionLength(Particle& p, Track& t) const {
- return 1./GetRef().GetInteractionLength(p, t);
- }
- */
- };
-
- /*
- template<template<typename, typename> class T, typename A, typename B>
- typename BaseProcess< T<A, B> >::is_process_sequence
- {
- static const bool value = true;
- };
- */
-
- /*
- template <typename T>
- struct is_base {
- static const bool value = false;
- };
- template <typename T>
- struct is_base<BaseProcess<T>> {
- static const bool value = true;
};
- */
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index c578f235..86779739 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -13,7 +13,6 @@
#define _include_corsika_continuousprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
-//#include <corsika/setup/SetupTrajectory.h>
namespace corsika::process {
diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index afd75d69..ae4869ea 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -26,6 +26,7 @@ namespace corsika::process {
eInteracted = 3,
eDecayed = 4,
};
+
} // namespace corsika::process
#endif
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index d58c14b5..01d8d8d6 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -15,7 +15,6 @@
#include <corsika/process/BaseProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
-//#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/units/PhysicalUnits.h>
@@ -23,112 +22,8 @@
#include <cmath>
#include <limits>
-//#include <corsika/setup/SetupTrajectory.h>
-// using corsika::setup::Trajectory;
-//#include <variant>
-
-//#include <type_traits> // still needed ?
-
namespace corsika::process {
- // namespace detail {
-
- /* /\* template<typename TT1, typename TT2, typename Type = void> *\/ */
- /* /\* struct CallHello { *\/ */
- /* /\* static void Call(const TT1&, const TT2&) { *\/ */
- /* /\* std::cout << "normal" << std::endl; *\/ */
- /* /\* } *\/ */
- /* /\* }; *\/ */
-
- /* /\* template<typename TT1, typename TT2> *\/ */
- /* /\* struct CallHello<TT1, TT2, typename
- * std::enable_if<std::is_base_of<ContinuousProcess<TT2>, TT2>::value>::type> *\/ */
- /* /\* { *\/ */
- /* /\* static void Call(const TT1&, const TT2&) { *\/ */
- /* /\* std::cout << "special" << std::endl; *\/ */
- /* /\* } *\/ */
- /* }; */
-
- /* template<typename T1, typename T2, typename Particle, typename Trajectory, typename
- * Stack> //, typename Type = void> */
- /* struct DoContinuous { */
- /* static EProcessReturn Call(const T1& A, const T2& B, Particle& p, Trajectory& t,
- * Stack& s) { */
- /* EProcessReturn ret = EProcessReturn::eOk; */
- /* if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) { */
- /* A.DoContinuous(p, t, s); */
- /* } */
- /* if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) { */
- /* B.DoContinuous(p, t, s); */
- /* } */
- /* return ret; */
- /* } */
- /* }; */
-
- /* /\* */
- /* template<typename T1, typename T2, typename Particle, typename Trajectory, typename
- * Stack> */
- /* struct DoContinuous<T1,T2,Particle,Trajectory,Stack, typename
- * std::enable_if<std::is_base_of<DiscreteProcess<T1>, T1>::value>::type> { */
- /* static EProcessReturn Call(const T1& A, const T2& B, Particle& p, Trajectory& t,
- * Stack& s) { */
- /* EProcessReturn ret = EProcessReturn::eOk; */
- /* A.DoContinuous(p, t, s); */
- /* B.DoContinuous(p, t, s); */
- /* return ret; */
- /* } */
- /* }; */
-
- /* template<typename T1, typename T2, typename Particle, typename Trajectory,
- * typename Stack> */
- /* struct DoContinuous<T1,T2,Particle,Trajectory,Stack, typename
- * std::enable_if<std::is_base_of<DiscreteProcess<T2>, T2>::value>::type> { */
- /* static EProcessReturn Call(const T1& A, const T2&, Particle& p, Trajectory& t,
- * Stack& s) { */
- /* EProcessReturn ret = EProcessReturn::eOk; */
- /* A.DoContinuous(p, t, s); */
- /* B.DoContinuous(p, t, s); */
- /* return ret; */
- /* } */
- /* }; */
- /* *\/ */
-
- /* template<typename T1, typename T2, typename Particle, typename Stack>//, typename
- * Type = void> */
- /* struct DoDiscrete { */
- /* static EProcessReturn Call(const T1& A, const T2& B, Particle& p, Stack& s) { */
- /* if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) { */
- /* A.DoDiscrete(p, s); */
- /* } */
- /* if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) { */
- /* B.DoDiscrete(p, s); */
- /* } */
- /* return EProcessReturn::eOk; */
- /* } */
- /* }; */
- /* /\* */
- /* template<typename T1, typename T2, typename Particle, typename Stack> */
- /* struct DoDiscrete<T1,T2,Particle,Stack, typename
- * std::enable_if<std::is_base_of<ContinuousProcess<T1>, T1>::value>::type> { */
- /* static EProcessReturn Call(const T1&, const T2& B, Particle& p, Stack& s) { */
- /* // A.DoDiscrete(p, s); */
- /* B.DoDiscrete(p, s); */
- /* return EProcessReturn::eOk; */
- /* } */
- /* }; */
-
- /* template<typename T1, typename T2, typename Particle, typename Stack> */
- /* struct DoDiscrete<T1,T2,Particle,Stack, typename
- * std::enable_if<std::is_base_of<ContinuousProcess<T2>, T2>::value>::type> { */
- /* static EProcessReturn Call(const T1& A, const T2&, Particle& p, Stack& s) { */
- /* A.DoDiscrete(p, s); */
- /* //B.DoDiscrete(p, s); */
- /* return EProcessReturn::eOk; */
- /* } */
- /* }; */
- /* *\/ */
- //} // end namespace detail
-
/**
\class ProcessSequence
@@ -337,10 +232,10 @@ namespace corsika::process {
/// must be allowed, this is why we define a macro to define all
/// combinations here:
-#define OPSEQ(C1, C2) \
- template <typename T1, typename T2> \
- inline const ProcessSequence<T1, T2> operator+(const C1<T1>& A, const C2<T2>& B) { \
- return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); \
+#define OPSEQ(C1, C2) \
+ template <typename T1, typename T2> \
+ inline const ProcessSequence<T1, T2> operator<<(const C1<T1>& A, const C2<T2>& B) { \
+ return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); \
}
OPSEQ(BaseProcess, BaseProcess)
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index e8a4828c..1723a689 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -188,7 +188,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process3 m3(2);
Process4 m4(3);
- const auto sequence = m1 + m2 + m3 + m4;
+ const auto sequence = m1 << m2 << m3 << m4;
globalCount = 0;
sequence.Init();
@@ -208,7 +208,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
DummyStack s;
DummyTrajectory t;
- const auto sequence2 = cp1 + m2 + m3;
+ const auto sequence2 = cp1 << m2 << m3;
GrammageType const tot = sequence2.GetTotalInteractionLength(s, t);
InverseGrammageType const tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
@@ -222,7 +222,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
DummyStack s;
- const auto sequence2 = cp1 + m2 + m3 + d3;
+ const auto sequence2 = cp1 << m2 << m3 << d3;
TimeType const tot = sequence2.GetTotalLifetime(s);
InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(s);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
@@ -235,7 +235,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process2 m2(1);
Process3 m3(2);
- const auto sequence2 = cp1 + m2 + m3 + cp2;
+ const auto sequence2 = cp1 << m2 << m3 << cp2;
DummyData p;
DummyStack s;
diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index 1d1d51ef..5ef6ca66 100644
--- a/Framework/StackInterface/ParticleBase.h
+++ b/Framework/StackInterface/ParticleBase.h
@@ -30,13 +30,11 @@ namespace corsika::stack {
template <typename StackIterator>
class ParticleBase {
- // friend class Stack<StackData, PI>; // for access to GetIterator
public:
ParticleBase() {}
private:
ParticleBase(ParticleBase&);
- // ParticleBase& operation=(ParticleBase& p);
public:
/// delete this particle on the stack. The corresponding iterator
diff --git a/Framework/StackInterface/StackIterator.h b/Framework/StackInterface/StackIterator.h
index 9e8d07e6..47424f60 100644
--- a/Framework/StackInterface/StackIterator.h
+++ b/Framework/StackInterface/StackIterator.h
@@ -14,9 +14,6 @@
#include <corsika/stack/ParticleBase.h>
-#include <iomanip>
-#include <iostream>
-
class StackData; // forward decl
namespace corsika::stack {
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 94f98ac2..9aef04ee 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -27,7 +27,7 @@
#include <phys/units/quantity.hpp>
-namespace corsika::units::si::constants {
+namespace corsika::units::constants {
using namespace phys::units;
@@ -38,12 +38,13 @@ namespace corsika::units::si::constants {
// Avogadro constant
constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole};
- // electronvolt
- constexpr quantity<energy_d> eV{Rep(1.6021766208e-19L) * joule};
// elementary charge
constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
+ // electronvolt
+ constexpr quantity<energy_d> eV{e / coulomb * joule};
+
// Planck constant
constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
@@ -59,6 +60,6 @@ namespace corsika::units::si::constants {
// etc.
-} // namespace corsika::units::si::constants
+} // namespace corsika::units::constants
#endif // PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 03029541..7565c603 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -6,20 +6,34 @@
#include <phys/units/io.hpp>
#include <phys/units/quantity.hpp>
+/*
+ It is essentially a bug of the phys/units package to define the
+ operator<< not in the same namespace as the types it is working
+ on. This breaks ADL (argument-dependent lookup). Here we "fix" this:
+ */
+namespace phys::units {
+ // using namespace phys::units::io;
+ using phys::units::io::operator<<;
+} // namespace phys::units
+
/**
* @file PhysicalUnits
*
- * Add new units and types we need
+ * Add new units and types we need. Units are compile-time. We support
+ * different system of units in parallel. Literals are used for
+ * optimal coding style.
*
- * Define _XeV literals, etc., allowing 10_GeV in the code.
*/
namespace corsika::units::hep {
using namespace phys::units;
using namespace phys::units::literals;
+ using namespace phys::units::io;
/// defining HEP energy, mass, momentum
using energy_hep_d = phys::units::energy_d;
+ constexpr phys::units::quantity<energy_hep_d> GeV{corsika::units::constants::eV};
+ // corsika::units::constants::e / phys::units::coulomb * phys::units::joule };
using MassType = phys::units::quantity<energy_hep_d, double>;
using MomentumType = phys::units::quantity<energy_hep_d, double>;
@@ -30,16 +44,20 @@ namespace corsika::units::hep {
namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
- // namespace literals = phys::units::literals;
+ using namespace phys::units::io;
+ using phys::units::io::operator<<;
/// defining momentum you suckers
/// dimensions, i.e. composition in base SI dimensions
using momentum_d = phys::units::dimensions<1, 1, -1>;
// defining the unit of momentum, so far newton-meter, maybe go to HEP?
- constexpr phys::units::quantity<momentum_d> newton_second{meter * kilogram / second};
+ constexpr phys::units::quantity<momentum_d> newton_second{
+ phys::units::meter * phys::units::kilogram / phys::units::second};
/// defining cross section
- constexpr phys::units::quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
+ using sigma_d = phys::units::dimensions<2, 0, 0>;
+ constexpr phys::units::quantity<sigma_d> barn{phys::units::Rep(1.e-28L) *
+ phys::units::meter * phys::units::meter};
/// add the unit-types
using LengthType = phys::units::quantity<phys::units::length_d, double>;
@@ -75,10 +93,16 @@ namespace phys {
namespace units {
namespace literals {
QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
- magnitude(corsika::units::si::constants::eV))
+ magnitude(corsika::units::constants::eV))
+
+ // QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::area_d,
+ // magnitude(corsika::units::si::constants::barn))
- QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::area_d,
- magnitude(corsika::units::si::constants::barn))
+ QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
+ magnitude(corsika::units::constants::barn))
+
+ QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
+ magnitude(corsika::units::constants::meter))
QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
magnitude(1_m * 1_kg / 1_s))
@@ -87,7 +111,4 @@ namespace phys {
} // namespace units
} // namespace phys
-// we want to call the operator<< without namespace... I think
-using namespace phys::units::io;
-
#endif
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index db210a77..4dc6316e 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -39,7 +39,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
[[maybe_unused]] LengthType arr1[2] = {{1_mm}, {2_cm}};
- std::array<EnergyType, 4> arr2; // empty array
+ [[maybe_unused]] std::array<EnergyType, 4> arr2; // empty array
[[maybe_unused]] std::array<EnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
@@ -97,9 +97,13 @@ TEST_CASE("PhysicalUnits", "[Units]") {
SECTION("Unit system conversion") {
const units::hep::MassType m_hep = 3_GeV;
+ std::cout << m_hep << std::endl;
+
+ const units::si::MassType m_hep2 = 3_kg;
+ std::cout << m_hep2 << std::endl;
REQUIRE(m_hep == 3_GeV); // hep::mass identical to si::energy
- auto type_check = m_hep / units::si::constants::cSquared;
+ auto type_check = m_hep / units::constants::cSquared;
REQUIRE(dynamic_cast<units::si::MassType*>(&type_check)); // hep::mass*c2 is mass unit
const units::hep::EnergyType e_hep = 4_GeV;
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 05046679..53173247 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -13,74 +13,6 @@ namespace corsika::process {
namespace sibyll {
- void setHadronsUnstable() {
- // name? also makes EM particles stable
-
- // loop over all particles in sibyll
- // should be changed to loop over human readable list
- // i.e. corsika::particles::ListOfParticles()
- std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
- // make ALL particles unstable, then set EM stable
- for (auto& p : corsika2sibyll) {
- // std::cout << (int)p << std::endl;
- const int sibCode = (int)p;
- // skip unknown and antiparticles
- if (sibCode < 1) continue;
- // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
- // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
- s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
- // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
- // std::endl;
- }
- // set Leptons and Proton and Neutron stable
- // use stack to loop over particles
- setup::Stack ds;
- ds.NewParticle().SetPID(corsika::particles::Code::Proton);
- ds.NewParticle().SetPID(corsika::particles::Code::Neutron);
- ds.NewParticle().SetPID(corsika::particles::Code::Electron);
- ds.NewParticle().SetPID(corsika::particles::Code::Positron);
- ds.NewParticle().SetPID(corsika::particles::Code::NuE);
- ds.NewParticle().SetPID(corsika::particles::Code::NuEBar);
- ds.NewParticle().SetPID(corsika::particles::Code::MuMinus);
- ds.NewParticle().SetPID(corsika::particles::Code::MuPlus);
- ds.NewParticle().SetPID(corsika::particles::Code::NuMu);
- ds.NewParticle().SetPID(corsika::particles::Code::NuMuBar);
-
- for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
- // << " stable in Sibyll .." << endl;
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- p.Delete();
- }
- }
-
- void setTrackedParticlesStable() {
- /*
- Sibyll is hadronic generator
- only hadrons decay
- */
- // set particles unstable
- setHadronsUnstable();
- // make tracked particles stable
- std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
- setup::Stack ds;
- ds.NewParticle().SetPID(particles::Code::PiPlus);
- ds.NewParticle().SetPID(particles::Code::PiMinus);
- ds.NewParticle().SetPID(particles::Code::KPlus);
- ds.NewParticle().SetPID(particles::Code::KMinus);
- ds.NewParticle().SetPID(particles::Code::K0Long);
- ds.NewParticle().SetPID(particles::Code::K0Short);
-
- for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- p.Delete();
- }
- }
-
class Decay : public corsika::process::DecayProcess<Decay> {
public:
Decay() {}
@@ -89,7 +21,32 @@ namespace corsika::process {
setTrackedParticlesStable();
}
- void setAllStable() {
+ void setTrackedParticlesStable() const {
+ /*
+ Sibyll is hadronic generator
+ only hadrons decay
+ */
+ // set particles unstable
+ setHadronsUnstable();
+ // make tracked particles stable
+ std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
+ const std::vector<corsika::particles::Code> particleList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+
+ for (auto p : particleList) {
+ // set particle stable by setting table value negative
+ const int sibid = process::sibyll::ConvertToSibyllRaw(p);
+ s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
+ }
+ }
+
+ void setAllStable() const {
+
+ // using namespace corsika::io;
+ using std::cout;
+ using std::endl;
+
// name? also makes EM particles stable
// loop over all particles in sibyll
@@ -100,27 +57,71 @@ namespace corsika::process {
const int sibCode = (int)p;
// skip unknown and antiparticles
if (sibCode < 1) continue;
- std::cout << "Sibyll: Decay: setting "
- << ConvertFromSibyll(static_cast<SibyllCode>(sibCode)) << " stable"
- << std::endl;
+ 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;
}
}
- friend void setHadronsUnstable();
+ void setHadronsUnstable() const {
+
+ using std::cout;
+ using std::endl;
+ // using namespace corsika::io;
+ using namespace corsika::units::si;
+
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. corsika::particles::ListOfParticles()
+ std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
+ // make ALL particles unstable, then set EM stable
+ for (auto& p : corsika2sibyll) {
+ // std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if (sibCode < 1) continue;
+ // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
+ // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
+ s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
+ // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
+ // std::endl;
+ }
+ // set Leptons and Proton and Neutron stable
+ // use stack to loop over particles
+ const std::vector<corsika::particles::Code> particleList = {
+ corsika::particles::Code::Proton, corsika::particles::Code::Neutron,
+ corsika::particles::Code::Electron, corsika::particles::Code::Positron,
+ corsika::particles::Code::NuE, corsika::particles::Code::NuEBar,
+ corsika::particles::Code::MuMinus, corsika::particles::Code::MuPlus,
+ corsika::particles::Code::NuMu, corsika::particles::Code::NuMuBar};
+
+ for (auto p : particleList) {
+ // set particle stable by setting table value negative
+ // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
+ // << " stable in Sibyll .." << endl;
+ const int sibid = process::sibyll::ConvertToSibyllRaw(p);
+ s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
+ }
+ }
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle& p) const {
+ using std::cout;
+ using std::endl;
+ using namespace corsika::units::si;
+
corsika::units::hep::EnergyType E = p.GetEnergy();
corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
- // const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
-
const double gamma = E / m;
- const TimeType t0 = particles::GetLifetime(p.GetPID());
- cout << "Decay: code: " << (p.GetPID()) << endl;
+ const corsika::units::si::TimeType t0 =
+ corsika::particles::GetLifetime(p.GetPID());
+ cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
// cout << "Decay: MinStep: density: " << density << endl;
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 24401ad7..80b4e449 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -3,9 +3,6 @@
#include <corsika/process/InteractionProcess.h>
-//#include <corsika/setup/SetupStack.h>
-//#include <corsika/setup/SetupTrajectory.h>
-
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
@@ -14,10 +11,6 @@
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
-using namespace corsika;
-using namespace corsika::process::sibyll;
-using namespace corsika::units::si;
-
namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
@@ -27,6 +20,11 @@ namespace corsika::process::sibyll {
~Interaction() {}
void Init() {
+
+ using corsika::random::RNGManager;
+ using std::cout;
+ using std::endl;
+
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
@@ -40,10 +38,16 @@ namespace corsika::process::sibyll {
sibyll_ini_();
}
- // void setTrackedParticlesStable();
-
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const {
+
+ using namespace corsika::units;
+ using namespace corsika::units::hep;
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ using std::cout;
+ using std::endl;
+
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -65,14 +69,15 @@ namespace corsika::process::sibyll {
// FOR NOW: assume target is oxygen
int kTarget = 16;
- EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
- super_stupid::MomentumVector Ptot(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ hep::EnergyType Etot =
+ p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+ MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// FOR NOW: assume target is at rest
- super_stupid::MomentumVector pTarget(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
- EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
+ hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
double Ecm = sqs / 1_GeV;
std::cout << "Interaction: "
@@ -94,12 +99,12 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
+ si::CrossSectionType sig = prodCrossSection * 1_mbarn;
std::cout << "Interaction: "
<< "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
- const MassType nucleon_mass =
- 0.93827_GeV / corsika::units::si::constants::cSquared;
+ const si::MassType nucleon_mass =
+ 0.93827_GeV / corsika::units::constants::cSquared;
std::cout << "Interaction: "
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
@@ -109,13 +114,33 @@ namespace corsika::process::sibyll {
<< "interaction length (g/cm2): " << int_length << std::endl;
return int_length;
- } else {
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
+
+
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+
+ /*
+ what are the units of the output? slant depth or 3space length?
+
+ */
+ //
+ // int a = 0;
+ // const double next_step = -int_length * log(s_rndm_(a));
+ // std::cout << "Interaction: "
+ // << "next step (g/cm2): " << next_step << std::endl;
+ // return next_step;
}
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
+
+ using namespace corsika::units;
+ using namespace corsika::units::hep;
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ using std::cout;
+ using std::endl;
+
cout << "ProcessSibyll: "
<< "DoInteraction: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
@@ -142,9 +167,7 @@ namespace corsika::process::sibyll {
cout << "defining target momentum.." << endl;
// FOR NOW: target is always at rest
const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
- const auto pTarget = super_stupid::MomentumVector(
- rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c,
- 0. * 1_GeV / constants::c);
+ const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
cout << "target momentum (GeV/c): "
<< pTarget.GetComponents() / 1_GeV * constants::c << endl;
cout << "beam momentum (GeV/c): "
@@ -162,18 +185,17 @@ namespace corsika::process::sibyll {
EnergyType E = p.GetEnergy();
EnergyType Etot = E + Etarget;
// total momentum
- super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ MomentumVector Ptot = p.GetMomentum(); // + pTarget;
// invariant mass, i.e. cm. energy
EnergyType Ecm =
- sqrt(Etot * Etot - Ptot.squaredNorm() *
- constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
+ sqrt(Etot * Etot - Ptot.squaredNorm()); // sqrt( 2. * E * 0.93827_GeV );
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
the transformation should be derived from the input momenta
*/
const double gamma = Etot / Ecm;
- const auto gambet = Ptot / (Ecm / constants::c);
+ const auto gambet = Ptot / Ecm;
std::cout << "Interaction: "
<< " DoDiscrete: gamma:" << gamma << endl;
@@ -195,7 +217,6 @@ namespace corsika::process::sibyll {
// running sibyll, filling stack
sibyll_(kBeam, kTarget, sqs);
// running decays
- // setTrackedParticlesStable();
decsib_();
// print final state
int print_unit = 6;
@@ -211,7 +232,7 @@ namespace corsika::process::sibyll {
// SibStack does not know about momentum yet so we need counter to access
// momentum array in Sibyll
int i = -1;
- super_stupid::MomentumVector Ptot_final(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ MomentumVector Ptot_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
for (auto& psib : ss) {
++i;
// skip particles that have decayed in Sibyll
@@ -222,20 +243,18 @@ namespace corsika::process::sibyll {
// arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents();
const auto pSibyllComponents = psib.GetMomentum().GetComponents();
- EnergyType en_lab = 0. * 1_GeV;
- MomentumType p_lab_components[3];
+ hep::EnergyType en_lab = 0. * 1_GeV;
+ hep::MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
- EnergyType pnorm = 0. * 1_GeV;
+ hep::EnergyType pnorm = 0. * 1_GeV;
for (int j = 0; j < 3; ++j)
- pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
- (gamma + 1.);
+ pnorm += (pSibyllComponents[j] * gammaBetaComponents[j]) / (gamma + 1.);
pnorm += psib.GetEnergy();
for (int j = 0; j < 3; ++j) {
- p_lab_components[j] = pSibyllComponents[j] -
- (-1) * pnorm * gammaBetaComponents[j] / constants::c;
- en_lab -=
- (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
+ p_lab_components[j] =
+ pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j];
+ en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j];
}
// add to corsika stack
@@ -243,9 +262,9 @@ namespace corsika::process::sibyll {
pnew.SetEnergy(en_lab);
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- corsika::geometry::QuantityVector<momentum_d> p_lab_c{
+ corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
- pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
+ pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
Ptot_final += pnew.GetMomentum();
}
// cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index 264cff1b..12da1bdd 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -1,17 +1,16 @@
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
+namespace corsika::process::sibyll {
-using namespace std;
-using namespace corsika::stack;
-using namespace corsika::units::si;
-using namespace corsika::geometry;
+ typedef corsika::geometry::Vector<corsika::units::hep::energy_hep_d> MomentumVector;
namespace corsika::process::sibyll {
@@ -27,25 +26,35 @@ namespace corsika::process::sibyll {
int GetCapacity() const { return 8000; }
void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
- void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
+ void SetEnergy(const int i, const corsika::units::hep::EnergyType v) {
+ using namespace corsika::units::hep;
+ s_plist_.p[3][i] = v / 1_GeV;
+ }
+ void SetMomentum(const int i, const MomentumVector& v) {
+ using namespace corsika::units;
+ using namespace corsika::units::hep;
auto tmp = v.GetComponents();
for (int idx = 0; idx < 3; ++idx)
- s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV;
}
int GetId(const int i) const { return s_plist_.llist[i]; }
- EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+ corsika::units::hep::EnergyType GetEnergy(const int i) const {
+ using namespace corsika::units::hep;
+ return s_plist_.p[3][i] * 1_GeV;
+ }
- super_stupid::MomentumVector GetMomentum(const int i) const {
+ MomentumVector GetMomentum(const int i) const {
+ using corsika::geometry::CoordinateSystem;
+ using corsika::geometry::QuantityVector;
+ using corsika::geometry::RootCoordinateSystem;
+ using namespace corsika::units::hep;
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- corsika::geometry::QuantityVector<momentum_d> components{
- s_plist_.p[0][i] * 1_GeV / constants::c,
- s_plist_.p[1][i] * 1_GeV / constants::c,
- s_plist_.p[2][i] * 1_GeV / constants::c};
- super_stupid::MomentumVector v1(rootCS, components);
+ QuantityVector<energy_hep_d> components = {
+ s_plist_.p[0][i] * 1_GeV, s_plist_.p[1][i] * 1_GeV, s_plist_.p[2][i] * 1_GeV};
+ MomentumVector v1(rootCS, components);
return v1;
}
@@ -62,28 +71,31 @@ namespace corsika::process::sibyll {
};
template <typename StackIteratorInterface>
- class ParticleInterface : public ParticleBase<StackIteratorInterface> {
- using ParticleBase<StackIteratorInterface>::GetStackData;
- using ParticleBase<StackIteratorInterface>::GetIndex;
+ class ParticleInterface : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
- EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ void SetEnergy(const corsika::units::hep::EnergyType v) {
+ GetStackData().SetEnergy(GetIndex(), v);
+ }
+ corsika::units::hep::EnergyType GetEnergy() const {
+ return GetStackData().GetEnergy(GetIndex());
+ }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const {
return static_cast<corsika::process::sibyll::SibyllCode>(
GetStackData().GetId(GetIndex()));
}
- super_stupid::MomentumVector GetMomentum() const {
- return GetStackData().GetMomentum(GetIndex());
- }
- void SetMomentum(const super_stupid::MomentumVector& v) {
+ MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
+ void SetMomentum(const MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v);
}
};
- typedef Stack<SibStackData, ParticleInterface> SibStack;
+ typedef corsika::stack::Stack<SibStackData, ParticleInterface> SibStack;
-} // namespace corsika::process::sibyll
+} // end namespace corsika::process::sibyll
#endif
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 52fe7e74..abb29da2 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -23,6 +23,7 @@
using namespace std;
using namespace corsika;
+using namespace corsika::particles;
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
@@ -37,7 +38,8 @@ StackInspector<Stack>::~StackInspector() {}
template <typename Stack>
process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
Stack& s) const {
- if (!fReport) return EProcessReturn::eOk;
+
+ if (!fReport) return process::EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
EnergyType Etot = 0_GeV;
@@ -54,7 +56,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
fCountStep++;
cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
<< " Estack=" << Etot / 1_GeV << " GeV" << endl;
- return EProcessReturn::eOk;
+ return process::EProcessReturn::eOk;
}
template <typename Stack>
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 157751ea..f92ee8ee 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -77,7 +77,7 @@ namespace corsika::process {
auto GetTrack(Particle const& p) {
using namespace corsika::units::si;
geometry::Vector<SpeedType::dimension_type> const velocity =
- p.GetMomentum() / p.GetEnergy() * corsika::units::si::constants::cSquared;
+ p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index a3df3264..737b536f 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -32,12 +32,14 @@ using namespace corsika::geometry;
using namespace std;
using namespace corsika::units::si;
+typedef corsika::units::hep::energy_hep_d MOMENTUM;
+
struct DummyParticle {
EnergyType fEnergy;
- Vector<momentum_d> fMomentum;
+ Vector<MOMENTUM> fMomentum;
Point fPosition;
- DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition)
+ DummyParticle(EnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
: fEnergy(pEnergy)
, fMomentum(pMomentum)
, fPosition(pPosition) {}
@@ -85,11 +87,8 @@ TEST_CASE("TrackingLine") {
//~ std::cout << env.GetUniverse().get() << std::endl;
- DummyParticle p(
- 1_J,
- Vector<momentum_d>(cs, 0 * kilogram * meter / second,
- 0 * kilogram * meter / second, 1 * kilogram * meter / second),
- Point(cs, 0_m, 0_m, 0_m));
+ DummyParticle p(1_J, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV),
+ Point(cs, 0_m, 0_m, 0_m));
auto const radius = 20_m;
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index acbbc373..de543484 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -18,15 +18,12 @@
#include <corsika/units/PhysicalUnits.h>
-#include <variant>
+// #include <variant>
namespace corsika::setup {
- using corsika::geometry::Helix;
- using corsika::geometry::Line;
-
/// definition of Trajectory base class, to be used in tracking and cascades
- typedef corsika::geometry::Trajectory<Line> Trajectory;
+ typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
/*
typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index d865c2d1..3a12c052 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -29,20 +29,7 @@ namespace corsika::stack {
namespace super_stupid {
- using corsika::geometry::Point;
- using corsika::geometry::Vector;
- using corsika::particles::Code;
- using corsika::units::si::energy_d;
- using corsika::units::si::EnergyType;
- using corsika::units::si::joule;
- using corsika::units::si::meter;
- using corsika::units::si::momentum_d;
- using corsika::units::si::newton_second;
- using corsika::units::si::second;
- using corsika::units::si::SpeedType;
- using corsika::units::si::TimeType;
-
- typedef Vector<momentum_d> MomentumVector;
+ typedef corsika::geometry::Vector<corsika::units::hep::energy_hep_d> MomentumVector;
/**
* Example of a particle object on the stack.
@@ -51,30 +38,48 @@ namespace corsika::stack {
template <typename StackIteratorInterface>
class ParticleInterface : public ParticleBase<StackIteratorInterface> {
- using ParticleBase<StackIteratorInterface>::GetStackData;
- using ParticleBase<StackIteratorInterface>::GetIndex;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetPID(const Code id) { GetStackData().SetPID(GetIndex(), id); }
- void SetEnergy(const EnergyType& e) { GetStackData().SetEnergy(GetIndex(), e); }
+ void SetPID(const corsika::particles::Code id) {
+ GetStackData().SetPID(GetIndex(), id);
+ }
+ void SetEnergy(const corsika::units::hep::EnergyType& e) {
+ GetStackData().SetEnergy(GetIndex(), e);
+ }
void SetMomentum(const MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v);
}
- void SetPosition(const Point& v) { GetStackData().SetPosition(GetIndex(), v); }
- void SetTime(const TimeType& v) { GetStackData().SetTime(GetIndex(), v); }
+ void SetPosition(const corsika::geometry::Point& v) {
+ GetStackData().SetPosition(GetIndex(), v);
+ }
+ void SetTime(const corsika::units::si::TimeType& v) {
+ GetStackData().SetTime(GetIndex(), v);
+ }
- Code GetPID() const { return GetStackData().GetPID(GetIndex()); }
- EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ corsika::particles::Code GetPID() const {
+ return GetStackData().GetPID(GetIndex());
+ }
+ corsika::units::hep::EnergyType GetEnergy() const {
+ return GetStackData().GetEnergy(GetIndex());
+ }
MomentumVector GetMomentum() const {
return GetStackData().GetMomentum(GetIndex());
}
- Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); }
- TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
+ corsika::geometry::Point GetPosition() const {
+ return GetStackData().GetPosition(GetIndex());
+ }
+ corsika::units::si::TimeType GetTime() const {
+ return GetStackData().GetTime(GetIndex());
+ }
-#warning this does not really work, nor makes sense:
- Vector<SpeedType::dimension_type> GetDirection() const {
+ //#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 * (units::si::meter / units::si::second);
+ return P / P.norm() * 1e10 *
+ (corsika::units::si::meter / corsika::units::si::second);
}
};
@@ -99,17 +104,21 @@ namespace corsika::stack {
int GetSize() const { return fDataPID.size(); }
int GetCapacity() const { return fDataPID.size(); }
- void SetPID(const int i, const Code id) { fDataPID[i] = id; }
- void SetEnergy(const int i, const EnergyType e) { fDataE[i] = e; }
+ void SetPID(const int i, const corsika::particles::Code id) { fDataPID[i] = id; }
+ void SetEnergy(const int i, const corsika::units::hep::EnergyType e) {
+ fDataE[i] = e;
+ }
void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
- void SetPosition(const int i, const Point& v) { fPosition[i] = v; }
- void SetTime(const int i, const TimeType& v) { fTime[i] = v; }
+ void SetPosition(const int i, const corsika::geometry::Point& v) {
+ fPosition[i] = v;
+ }
+ void SetTime(const int i, const corsika::units::si::TimeType& v) { fTime[i] = v; }
- Code GetPID(const int i) const { return fDataPID[i]; }
- EnergyType GetEnergy(const int i) const { return fDataE[i]; }
+ corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }
+ corsika::units::hep::EnergyType GetEnergy(const int i) const { return fDataE[i]; }
MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
- Point GetPosition(const int i) const { return fPosition[i]; }
- TimeType GetTime(const int i) const { return fTime[i]; }
+ corsika::geometry::Point GetPosition(const int i) const { return fPosition[i]; }
+ corsika::units::si::TimeType GetTime(const int i) const { return fTime[i]; }
/**
* Function to copy particle at location i2 in stack to i1
@@ -135,15 +144,20 @@ namespace corsika::stack {
protected:
void IncrementSize() {
+ using corsika::geometry::Point;
+ using corsika::particles::Code;
fDataPID.push_back(Code::Unknown);
- fDataE.push_back(0 * joule);
+ fDataE.push_back(0 * corsika::units::si::joule);
//#TODO this here makes no sense: see issue #48
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
fMomentum.push_back(MomentumVector(
- dummyCS, {0 * newton_second, 0 * newton_second, 0 * newton_second}));
- fPosition.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
- fTime.push_back(0 * second);
+ dummyCS, {0 * corsika::units::si::joule, 0 * corsika::units::si::joule,
+ 0 * corsika::units::si::joule}));
+ fPosition.push_back(
+ Point(dummyCS, {0 * corsika::units::si::meter, 0 * corsika::units::si::meter,
+ 0 * corsika::units::si::meter}));
+ fTime.push_back(0 * corsika::units::si::second);
}
void DecrementSize() {
if (fDataE.size() > 0) {
@@ -158,11 +172,11 @@ namespace corsika::stack {
private:
/// the actual memory to store particle data
- std::vector<Code> fDataPID;
- std::vector<EnergyType> fDataE;
+ std::vector<corsika::particles::Code> fDataPID;
+ std::vector<corsika::units::hep::EnergyType> fDataE;
std::vector<MomentumVector> fMomentum;
- std::vector<Point> fPosition;
- std::vector<TimeType> fTime;
+ std::vector<corsika::geometry::Point> fPosition;
+ std::vector<corsika::units::si::TimeType> fTime;
}; // end class SuperStupidStackImpl
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index 21f88290..092baf5f 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/Stack/SuperStupidStack/testSuperStupidStack.cc
@@ -36,8 +36,7 @@ TEST_CASE("SuperStupidStack", "[stack]") {
p.SetEnergy(1.5_GeV);
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- p.SetMomentum(MomentumVector(
- dummyCS, {1 * newton_second, 1 * newton_second, 1 * newton_second}));
+ p.SetMomentum(MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}));
p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
p.SetTime(100_s);
--
GitLab