diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index a19c8bd89b0d7cf32884ed5589d08e11b299e720..e7944b8e63ef5b87ccb83101163c033f9e3c693e 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -135,8 +135,8 @@ int main(int argc, char** argv) {
auto const observationHeight = 0_km + builder.getEarthRadius();
auto const injectionHeight = 112.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
- sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
- si::detail::static_pow<2>(injectionHeight));
+ sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
+ units::static_pow<2>(injectionHeight));
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore +
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 55f316e1b5bfeb73a85aed6c2833447158d03754..b135396567b812dad56e1b6dd515736b1cd1e8f6 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -51,7 +51,7 @@ TEST_CASE("FlatExponential") {
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
LengthType const lambda = 3_m;
- auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+ auto const rho0 = 1_g / units::static_pow<3>(1_cm);
FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
protonComposition);
auto const tEnd = 5_s;
@@ -106,7 +106,7 @@ TEST_CASE("SlidingPlanarExponential") {
std::vector<float>{1.f});
LengthType const lambda = 3_m;
- auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+ auto const rho0 = 1_g / units::static_pow<3>(1_cm);
auto const tEnd = 5_s;
SlidingPlanarExponential<IMediumModel> const medium(gOrigin, rho0, lambda,
@@ -146,7 +146,7 @@ struct Exponential {
template <int N>
auto Derivative(Point const& p, Vector<dimensionless_d> const& v) const {
return v.GetComponents()[0] * (*this)(p) /
- corsika::units::si::detail::static_pow<N>(1_m);
+ corsika::units::static_pow<N>(1_m);
}
auto FirstDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 06ca5374fbfb5a8e334c1634b0c697f75819d4cd..9e9f26bde19df00c75d3f99593b471c5c705acb8 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -30,6 +30,19 @@ namespace phys::units {
*
*/
+namespace corsika::units {
+ template <int N, typename T>
+ auto constexpr static_pow([[maybe_unused]] T x) {
+ if constexpr (N == 0) {
+ return 1;
+ } else if constexpr (N > 0) {
+ return x * static_pow<N - 1, T>(x);
+ } else {
+ return 1 / static_pow<-N, T>(x);
+ }
+ }
+}
+
namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
@@ -67,19 +80,6 @@ namespace corsika::units::si {
using InverseGrammageType =
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
- namespace detail {
- template <int N, typename T>
- auto constexpr static_pow([[maybe_unused]] T x) {
- if constexpr (N == 0) {
- return 1;
- } else if constexpr (N > 0) {
- return x * static_pow<N - 1, T>(x);
- } else {
- return 1 / static_pow<-N, T>(x);
- }
- }
- } // namespace detail
-
template <typename DimFrom, typename DimTo>
auto constexpr ConversionFactorHEPToSI() {
static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 &&
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 9f033d243f66c41deb044a912d08dbc3d2d5d003..bae1ce1d7123d82c9878bc4a4bbd319f33ce9912 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -114,7 +114,7 @@ TEST_CASE("PhysicalUnits", "[Units]") {
}
SECTION("static_pow") {
- using namespace corsika::units::si::detail;
+ using namespace corsika::units;
double x = 235.7913;
REQUIRE(1 == static_pow<0, double>(x));
REQUIRE(x == static_pow<1, double>(x));
diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
index 18df797fe78d62d216df5c5e91de0b1dc0b4c400..031ea746f6651d41865d06e9cd37dc577f64b9a9 100644
--- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
@@ -57,8 +57,8 @@ TEST_CASE("CONEXSourceCut") {
auto const injectionHeight = 112.75_km + conex::earthRadius;
auto const t =
-observationHeight * cos(thetaRad) +
- sqrt(-units::si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
- units::si::detail::static_pow<2>(injectionHeight));
+ sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
+ units::static_pow<2>(injectionHeight));
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore +
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index c8f3eed3dd74a39d1de39e40be9112c5dab5f54e..dbece1ac127c2f3f3d0970971ce7fde6aba7aee5 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -50,7 +50,7 @@ units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
for (size_t i = 0; i < fractions.size(); ++i) {
auto const targetMass = particles::GetMass(components[i]);
- auto const s = detail::static_pow<2>(projectileEnergy + targetMass) -
+ auto const s = units::static_pow<2>(projectileEnergy + targetMass) -
projectileMomentumSquaredNorm;
avgCrossSection += CrossSection(s) * fractions[i];
}
@@ -92,7 +92,7 @@ process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupView& vie
for (size_t i = 0; i < components.size(); ++i) {
auto const targetMass = particles::GetMass(components[i]);
- auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) -
+ auto const s = units::static_pow<2>(projectileEnergy + targetMass) -
projectileMomentumSquaredNorm;
cross_section_of_components[i] = CrossSection(s);
}
@@ -119,7 +119,7 @@ process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupView& vie
auto const eTargetCoM = targetCoM.GetTimeLikeComponent();
auto const sqrtS = eProjectileCoM + eTargetCoM;
- auto const s = units::si::detail::static_pow<2>(sqrtS);
+ auto const s = units::static_pow<2>(sqrtS);
auto const B = this->B(s);
std::cout << B << std::endl;
@@ -158,7 +158,7 @@ process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupView& vie
p.SetMomentum(projectileScatteredLab.GetSpaceLikeComponents());
p.SetEnergy(
sqrt(projectileScatteredLab.GetSpaceLikeComponents().squaredNorm() +
- units::si::detail::static_pow<2>(particles::GetMass(
+ units::static_pow<2>(particles::GetMass(
p.GetPID())))); // Don't use energy from boost. It can be smaller than
// the momentum due to limited numerical accuracy.
@@ -185,7 +185,7 @@ units::si::CrossSectionType HadronicElasticInteraction::CrossSection(
// according to Schuler & Sjöstrand, PRD 49, 2257 (1994)
// (we ignore rho because rho^2 is just ~2 %)
auto const sigmaElastic =
- units::si::detail::static_pow<2>(sigmaTotal) /
+ units::static_pow<2>(sigmaTotal) /
(16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s)));
std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl;
diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc
index c63579e8819920a298dd39ddc90ffa8973f31454..08b67182b1f881a48ce291bc37ca1c36359a08e5 100644
--- a/Processes/InteractionCounter/testInteractionCounter.cc
+++ b/Processes/InteractionCounter/testInteractionCounter.cc
@@ -60,7 +60,7 @@ auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
- sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
+ sqrt(units::static_pow<2>(mN * vA) + pLab.squaredNorm());
auto particle =
stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
@@ -82,7 +82,7 @@ auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
- sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
+ sqrt(units::static_pow<2>(particles::GetMass(vProjectileType)) +
pLab.squaredNorm());
auto particle = stack->AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 84e04ba21079451cd2e8e2641eaaa641e508ded2..d8e9e14fc4573aa64bc99b3f47e4dfdf5e9f0c76 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -87,9 +87,7 @@ using namespace corsika::units;
template <typename TStackView>
auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
-
for (auto const& p : view) { sum += p.GetMomentum(); }
-
return sum;
}
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index cec9beb3ee2608a4d0a0314ba857e8ff47cd926b..0a795e1ee1c5809707edeb68505a671413b53303 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -126,8 +126,8 @@ CrossSectionType UrQMD::GetCrossSection(particles::Code projectileCode,
!IsNucleus(targetCode)) { // both particles are "special"
auto const mProj = particles::GetMass(projectileCode);
auto const mTar = particles::GetMass(targetCode);
- double sqrtS = sqrt(units::si::detail::static_pow<2>(mProj) +
- units::si::detail::static_pow<2>(mTar) + 2 * labEnergy * mTar) *
+ double sqrtS = sqrt(units::static_pow<2>(mProj) +
+ units::static_pow<2>(mTar) + 2 * labEnergy * mTar) *
(1 / 1_GeV);
// we must set some UrQMD globals first...
@@ -179,7 +179,7 @@ CrossSectionType UrQMD::GetCrossSection(particles::Code projectileCode,
int const At = IsNucleus(targetCode) ? particles::GetNucleusA(targetCode) : 1;
double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
- return 10_mb * M_PI * units::si::detail::static_pow<2>(maxImpact);
+ return 10_mb * M_PI * units::static_pow<2>(maxImpact);
// is a constant cross-section really reasonable?
}
}
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 71cdd8068864a3c4eda760696374f9b61e193634..5cb6be54ff7b6a0842a1729a2fa7fad2b8ae4a4f 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -86,7 +86,7 @@ auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
- sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
+ sqrt(units::static_pow<2>(mN * vA) + pLab.squaredNorm());
auto particle =
stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
@@ -107,7 +107,7 @@ auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
- sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
+ sqrt(units::static_pow<2>(particles::GetMass(vProjectileType)) +
pLab.squaredNorm());
auto particle = stack->AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
diff --git a/Stack/History/testHistoryView.cc b/Stack/History/testHistoryView.cc
index bba09b71c1dc9f08073f76dd215973132cc8c1fe..9079c815b3b625ce21b799c767e7a8b315b50ec8 100644
--- a/Stack/History/testHistoryView.cc
+++ b/Stack/History/testHistoryView.cc
@@ -75,29 +75,30 @@ TEST_CASE("HistoryStackExtension", "[stack]") {
corsika::history::EventPtr evt = p0.GetEvent();
CHECK(evt == nullptr);
- // add secondaries, 1st generation
- history::HistorySecondaryView<TestStackView> hview0(p0);
-
- auto const ev0 = p0.GetEvent();
- CHECK(ev0 == nullptr);
- // CHECK(ev0->projectile(s.begin()) == p0);
-
- // add 5 secondaries
- for (int i = 0; i < 5; ++i) {
- auto sec = hview0.AddSecondary(
+ SECTION("interface test, view") {
+
+ // add secondaries, 1st generation
+ history::HistorySecondaryView<TestStackView> hview0(p0);
+
+ auto const ev0 = p0.GetEvent();
+ CHECK(ev0 == nullptr);
+
+ // add 5 secondaries
+ for (int i = 0; i < 5; ++i) {
+ auto sec = hview0.AddSecondary(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
- CHECK(sec.GetParentEventIndex() == i);
- CHECK(sec.GetEvent() != nullptr);
- CHECK(sec.GetEvent()->parentEvent() == nullptr);
- }
+ CHECK(sec.GetParentEventIndex() == i);
+ CHECK(sec.GetEvent() != nullptr);
+ CHECK(sec.GetEvent()->parentEvent() == nullptr);
+ }
- // read 1st genertion particle particle
- auto p1 = stack.GetNextParticle();
+ // read 1st genertion particle particle
+ auto p1 = stack.GetNextParticle();
history::HistorySecondaryView<TestStackView> hview1(p1);
@@ -176,4 +177,58 @@ TEST_CASE("HistoryStackExtension", "[stack]") {
CHECK(test_proj1.GetEvent()->parentEvent()->parentEvent() == nullptr);
}
*/
+ }
+
+ SECTION("also test projectile access") {
+
+ // add secondaries, 1st generation
+ history::HistorySecondaryView<TestStackView> hview0(p0);
+ auto proj0 = hview0.GetProjectile();
+ auto const ev0 = p0.GetEvent();
+ CHECK(ev0 == nullptr);
+
+ // add 5 secondaries
+ for (int i = 0; i < 5; ++i) {
+ auto sec = proj0.AddSecondary(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
+
+ CHECK(sec.GetParentEventIndex() == i);
+ CHECK(sec.GetEvent() != nullptr);
+ CHECK(sec.GetEvent()->parentEvent() == nullptr);
+ }
+ CHECK(stack.getEntries()==6);
+
+
+ // read 1st genertion particle particle
+ auto p1 = stack.GetNextParticle();
+
+ history::HistorySecondaryView<TestStackView> hview1(p1);
+ auto proj1 = hview1.GetProjectile();
+ auto const ev1 = p1.GetEvent();
+
+ // add second generation of secondaries
+ // add 10 secondaries
+ for (int i = 0; i < 10; ++i) {
+ auto sec = proj1.AddSecondary(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Electron, 1.5_GeV,
+ corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
+
+ CHECK(sec.GetParentEventIndex() == i);
+ CHECK(sec.GetEvent()->parentEvent() == ev1);
+ CHECK(sec.GetEvent()->parentEvent()->secondaries().size() == 10);
+ CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev0);
+
+ CHECK((stack.begin() + sec.GetEvent()->projectileIndex()).GetEvent() ==
+ sec.GetEvent()->parentEvent());
+ }
+ CHECK(stack.getEntries()==16);
+
+ }
}