diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 05863e5431c44100e84c7fa2b764a4811adb3b6b..651cc981a2a0da9bf31a5a8f7e79eaeeea5a4ce1 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -80,20 +80,30 @@ int main(int argc, char** argv) {
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer<MEnv>(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addLinearLayer<environment::UniformMediumType<
+ environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
+ 1e9_cm, 112.8_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
@@ -174,16 +184,16 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
- em_continuous.ShowResults();
+ em_continuous.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
- cut.GetEmEnergy() + em_continuous.GetEnergyLost() +
+ cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
- em_continuous.Reset();
+ em_continuous.reset();
auto const hists = proposalCounted.GetHistogram();
hists.saveLab("inthist_lab_emShower.npz");
diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index bad8cee91271f546456ee61b2af398ee352b964d..51b34dc46501897c52659503a7242e4b3dd7c96d 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -72,9 +72,6 @@ void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().SeedAll(seed);
}
-template <typename T>
-using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
-
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
@@ -104,20 +101,30 @@ int main(int argc, char** argv) {
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer<MEnv>(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addLinearLayer<environment::UniformMediumType<
+ environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
+ 1e9_cm, 112.8_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 2d2a1b7cb031dc574337dddc654b90a904802cdb..ea3ff330e64446fac7aa88c2ca78088713ac34db 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -75,9 +75,6 @@ void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().SeedAll(seed);
}
-template <typename T>
-using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
-
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
@@ -107,20 +104,30 @@ int main(int argc, char** argv) {
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer<MEnv>(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addLinearLayer<environment::UniformMediumType<
+ environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
+ 1e9_cm, 112.8_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
@@ -257,16 +264,16 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
- em_continuous.ShowResults();
+ em_continuous.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
- cut.GetEmEnergy() + em_continuous.GetEnergyLost() +
+ cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
- em_continuous.Reset();
+ em_continuous.reset();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
diff --git a/Environment/LayeredSphericalAtmosphereBuilder.h b/Environment/LayeredSphericalAtmosphereBuilder.h
index 01f4855afe4cd97eb70d780fc0da8b5476847ba6..12db29715a3ae254d8de0e4b0bb5b5f248164a80 100644
--- a/Environment/LayeredSphericalAtmosphereBuilder.h
+++ b/Environment/LayeredSphericalAtmosphereBuilder.h
@@ -62,28 +62,13 @@ namespace corsika::environment {
composition_ = std::make_unique<NuclearComposition>(composition);
}
- void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
- units::si::LengthType upperBoundary) {
- auto const radius = earthRadius_ + upperBoundary;
- checkRadius(radius);
- previousRadius_ = radius;
-
- auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
- std::make_unique<geometry::Sphere>(center_, radius));
-
- auto const rho0 = b / c;
- std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
-
- node->template SetModelProperties<
- environment::SlidingPlanarExponential<TMediumInterface>>(
- center_, rho0, -c, *composition_, earthRadius_);
-
- layers_.push(std::move(node));
- }
-
- template <template <typename> typename MExtraModels, typename... TArgs>
+ template <
+ typename TMediumModel = environment::SlidingPlanarExponential<TMediumInterface>,
+ typename... TArgs>
void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
units::si::LengthType upperBoundary, TArgs&&... args) {
+ using namespace units::si;
+
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
@@ -94,16 +79,14 @@ namespace corsika::environment {
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
- node->template SetModelProperties<
- MExtraModels<SlidingPlanarExponential<TMediumInterface>>>(
- args..., center_, rho0, -c, *composition_, earthRadius_);
+ node->template SetModelProperties<TMediumModel>(args..., center_, rho0, -c,
+ *composition_, earthRadius_);
layers_.push(std::move(node));
}
- int size() const { return layers_.size(); }
-
- template <template <typename> typename MExtraModels, typename... TArgs>
+ template <typename TMediumModel = environment::HomogeneousMedium<TMediumInterface>,
+ typename... TArgs>
void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary,
TArgs&&... args) {
using namespace units::si;
@@ -120,32 +103,12 @@ namespace corsika::environment {
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<geometry::Sphere>(center_, radius));
- node->template SetModelProperties<
- MExtraModels<HomogeneousMedium<TMediumInterface>>>(args..., rho0,
- *composition_);
+ node->template SetModelProperties<TMediumModel>(args..., rho0, *composition_);
layers_.push(std::move(node));
}
- void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary) {
- using namespace units::si;
-
- auto const radius = earthRadius_ + upperBoundary;
- checkRadius(radius);
- previousRadius_ = radius;
-
- units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
- auto const rho0 = b / c;
-
- std::cout << "rho0 = " << rho0;
-
- auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
- std::make_unique<geometry::Sphere>(center_, radius));
- node->template SetModelProperties<HomogeneousMedium<TMediumInterface>>(
- rho0, *composition_);
-
- layers_.push(std::move(node));
- }
+ int size() const { return layers_.size(); }
void assemble(Environment<TMediumInterface>& env) {
auto& universe = env.GetUniverse();
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 5480aec9426406fa152f4c9d57cd45b0f97a65c8..6c0062c711aa567f2014b8e16c3f72334716829d 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -298,7 +298,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
// setup our interface types
- using IModelInterface = IMagneticFieldModel<IMediumModel>;
+ using ModelInterface = IMagneticFieldModel<IMediumModel>;
// the composition we use for the homogenous medium
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
@@ -308,14 +308,15 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
Vector B0(gCS, 0_T, 0_T, 1_T);
Vector B1(gCS, 1_T, 1_T, 0_T);
- // LayeredSphericalAtmosphereBuilder<IMediumModel> builder(gOrigin);
- LayeredSphericalAtmosphereBuilder<IModelInterface> builder(gOrigin);
+ LayeredSphericalAtmosphereBuilder<ModelInterface> builder(gOrigin);
builder.setNuclearComposition(
{{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
- builder.addLinearLayer<UniformMagneticField>(1_km, 10_km, B0);
- builder.addExponentialLayer<UniformMagneticField>(1222.6562_g / (1_cm * 1_cm),
- 994186.38_cm, 20_km, B1);
+ builder.addLinearLayer<UniformMagneticField<HomogeneousMedium<ModelInterface>>>(
+ 1_km, 10_km, B0);
+ builder.addExponentialLayer<
+ UniformMagneticField<SlidingPlanarExponential<ModelInterface>>>(
+ 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km, B1);
CHECK(builder.size() == 2);
diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
index 30831eed61e3e1f2e4168903a953366cf6dda674..6699e419c44edc920f261c605b601205f52e4444 100644
--- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
@@ -35,8 +35,8 @@ using namespace corsika::environment;
using namespace corsika::geometry;
using namespace corsika::units::si;
-template <typename T>
-using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
+// template <typename T>
+// using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
TEST_CASE("CONEXSourceCut") {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
@@ -56,20 +56,30 @@ TEST_CASE("CONEXSourceCut") {
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer<MEnv>(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
- environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
- builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir,
- geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<
+ environment::UniformMediumType<environment::UniformMagneticField<
+ SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
+ 540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addLinearLayer<environment::UniformMediumType<
+ environment::UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>(
+ 1e9_cm, 112.8_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
builder.assemble(env);
@@ -95,8 +105,7 @@ TEST_CASE("CONEXSourceCut") {
[[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
corsika::process::conex_source_cut::CONEXSourceCut conex(
- center, showerAxis, t, injectionHeight, E0,
- particles::GetPDG(particles::Code::Proton));
+ center, showerAxis, t, injectionHeight, E0, particles::Code::Proton);
HEPEnergyType const Eem{1_PeV};
auto const momentum = showerAxis.GetDirection() * Eem;
@@ -111,7 +120,8 @@ TEST_CASE("CONEXSourceCut") {
std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " "
<< emPosition.GetCoordinates(rootCS) << std::endl;
- conex.addParticle(0, Eem, 0_eV, emPosition, momentum.normalized(), 0_s);
+ conex.addParticle(particles::Code::Proton, Eem, 0_eV, emPosition, momentum.normalized(),
+ 0_s);
conex.SolveCE();
}