IAP GITLAB

Skip to content
Snippets Groups Projects
Commit bd6e3453 authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

improved/simplified interface to LayeredSp...Atm..Builder

parent c1dd7926
No related branches found
No related tags found
1 merge request!265Add medium type, as new property (air, water, rock...)
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/em_shower.cc
+ 27
− 17
View file @ bd6e3453
...@@ -80,20 +80,30 @@ int main(int argc, char** argv) { ...@@ -80,20 +80,30 @@ int main(int argc, char** argv) {
{{particles::Code::Nitrogen, particles::Code::Oxygen}, {{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now {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, builder.addExponentialLayer<
environment::EMediumType::eAir, environment::UniformMediumType<environment::UniformMagneticField<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
environment::EMediumType::eAir, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); builder.addExponentialLayer<
builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::UniformMediumType<environment::UniformMagneticField<
environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
environment::EMediumType::eAir, builder.addExponentialLayer<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); environment::UniformMediumType<environment::UniformMagneticField<
builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 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); builder.assemble(env);
...@@ -174,16 +184,16 @@ int main(int argc, char** argv) { ...@@ -174,16 +184,16 @@ int main(int argc, char** argv) {
EAS.Run(); EAS.Run();
cut.ShowResults(); cut.ShowResults();
em_continuous.ShowResults(); em_continuous.showResults();
observationLevel.ShowResults(); observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.GetEnergyLost() + cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround(); observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset(); observationLevel.Reset();
cut.Reset(); cut.Reset();
em_continuous.Reset(); em_continuous.reset();
auto const hists = proposalCounted.GetHistogram(); auto const hists = proposalCounted.GetHistogram();
hists.saveLab("inthist_lab_emShower.npz"); hists.saveLab("inthist_lab_emShower.npz");
......
Documentation/Examples/hybrid_MC.cc
+ 24
− 17
View file @ bd6e3453
...@@ -72,9 +72,6 @@ void registerRandomStreams(const int seed) { ...@@ -72,9 +72,6 @@ void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().SeedAll(seed); random::RNGManager::GetInstance().SeedAll(seed);
} }
template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) { int main(int argc, char** argv) {
logging::SetLevel(logging::level::info); logging::SetLevel(logging::level::info);
...@@ -104,20 +101,30 @@ int main(int argc, char** argv) { ...@@ -104,20 +101,30 @@ int main(int argc, char** argv) {
{{particles::Code::Nitrogen, particles::Code::Oxygen}, {{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now {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, builder.addExponentialLayer<
environment::EMediumType::eAir, environment::UniformMediumType<environment::UniformMagneticField<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
environment::EMediumType::eAir, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); builder.addExponentialLayer<
builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::UniformMediumType<environment::UniformMagneticField<
environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
environment::EMediumType::eAir, builder.addExponentialLayer<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); environment::UniformMediumType<environment::UniformMagneticField<
builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 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); builder.assemble(env);
......
Documentation/Examples/vertical_EAS.cc
+ 27
− 20
View file @ bd6e3453
...@@ -75,9 +75,6 @@ void registerRandomStreams(const int seed) { ...@@ -75,9 +75,6 @@ void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().SeedAll(seed); random::RNGManager::GetInstance().SeedAll(seed);
} }
template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) { int main(int argc, char** argv) {
logging::SetLevel(logging::level::info); logging::SetLevel(logging::level::info);
...@@ -107,20 +104,30 @@ int main(int argc, char** argv) { ...@@ -107,20 +104,30 @@ int main(int argc, char** argv) {
{{particles::Code::Nitrogen, particles::Code::Oxygen}, {{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now {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, builder.addExponentialLayer<
environment::EMediumType::eAir, environment::UniformMediumType<environment::UniformMagneticField<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
environment::EMediumType::eAir, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); builder.addExponentialLayer<
builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::UniformMediumType<environment::UniformMagneticField<
environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
environment::EMediumType::eAir, builder.addExponentialLayer<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); environment::UniformMediumType<environment::UniformMagneticField<
builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 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); builder.assemble(env);
...@@ -257,16 +264,16 @@ int main(int argc, char** argv) { ...@@ -257,16 +264,16 @@ int main(int argc, char** argv) {
EAS.Run(); EAS.Run();
cut.ShowResults(); cut.ShowResults();
em_continuous.ShowResults(); em_continuous.showResults();
observationLevel.ShowResults(); observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() + const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.GetEnergyLost() + cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround(); observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset(); observationLevel.Reset();
cut.Reset(); cut.Reset();
em_continuous.Reset(); em_continuous.reset();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() + auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram(); urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
......
Environment/LayeredSphericalAtmosphereBuilder.h
+ 11
− 48
View file @ bd6e3453
...@@ -62,28 +62,13 @@ namespace corsika::environment { ...@@ -62,28 +62,13 @@ namespace corsika::environment {
composition_ = std::make_unique<NuclearComposition>(composition); composition_ = std::make_unique<NuclearComposition>(composition);
} }
void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c, template <
units::si::LengthType upperBoundary) { typename TMediumModel = environment::SlidingPlanarExponential<TMediumInterface>,
auto const radius = earthRadius_ + upperBoundary; typename... TArgs>
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>
void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c, void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
units::si::LengthType upperBoundary, TArgs&&... args) { units::si::LengthType upperBoundary, TArgs&&... args) {
using namespace units::si;
auto const radius = earthRadius_ + upperBoundary; auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius); checkRadius(radius);
previousRadius_ = radius; previousRadius_ = radius;
...@@ -94,16 +79,14 @@ namespace corsika::environment { ...@@ -94,16 +79,14 @@ namespace corsika::environment {
auto const rho0 = b / c; auto const rho0 = b / c;
std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl; std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
node->template SetModelProperties< node->template SetModelProperties<TMediumModel>(args..., center_, rho0, -c,
MExtraModels<SlidingPlanarExponential<TMediumInterface>>>( *composition_, earthRadius_);
args..., center_, rho0, -c, *composition_, earthRadius_);
layers_.push(std::move(node)); layers_.push(std::move(node));
} }
int size() const { return layers_.size(); } template <typename TMediumModel = environment::HomogeneousMedium<TMediumInterface>,
typename... TArgs>
template <template <typename> typename MExtraModels, typename... TArgs>
void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary, void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary,
TArgs&&... args) { TArgs&&... args) {
using namespace units::si; using namespace units::si;
...@@ -120,32 +103,12 @@ namespace corsika::environment { ...@@ -120,32 +103,12 @@ namespace corsika::environment {
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>( auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<geometry::Sphere>(center_, radius)); std::make_unique<geometry::Sphere>(center_, radius));
node->template SetModelProperties< node->template SetModelProperties<TMediumModel>(args..., rho0, *composition_);
MExtraModels<HomogeneousMedium<TMediumInterface>>>(args..., rho0,
*composition_);
layers_.push(std::move(node)); layers_.push(std::move(node));
} }
void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary) { int size() const { return layers_.size(); }
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));
}
void assemble(Environment<TMediumInterface>& env) { void assemble(Environment<TMediumInterface>& env) {
auto& universe = env.GetUniverse(); auto& universe = env.GetUniverse();
......
Environment/testEnvironment.cc
+ 7
− 6
View file @ bd6e3453
...@@ -298,7 +298,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { ...@@ -298,7 +298,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") { TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
// setup our interface types // setup our interface types
using IModelInterface = IMagneticFieldModel<IMediumModel>; using ModelInterface = IMagneticFieldModel<IMediumModel>;
// the composition we use for the homogenous medium // the composition we use for the homogenous medium
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
...@@ -308,14 +308,15 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") { ...@@ -308,14 +308,15 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
Vector B0(gCS, 0_T, 0_T, 1_T); Vector B0(gCS, 0_T, 0_T, 1_T);
Vector B1(gCS, 1_T, 1_T, 0_T); Vector B1(gCS, 1_T, 1_T, 0_T);
// LayeredSphericalAtmosphereBuilder<IMediumModel> builder(gOrigin); LayeredSphericalAtmosphereBuilder<ModelInterface> builder(gOrigin);
LayeredSphericalAtmosphereBuilder<IModelInterface> builder(gOrigin);
builder.setNuclearComposition( builder.setNuclearComposition(
{{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}}); {{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
builder.addLinearLayer<UniformMagneticField>(1_km, 10_km, B0); builder.addLinearLayer<UniformMagneticField<HomogeneousMedium<ModelInterface>>>(
builder.addExponentialLayer<UniformMagneticField>(1222.6562_g / (1_cm * 1_cm), 1_km, 10_km, B0);
994186.38_cm, 20_km, B1); builder.addExponentialLayer<
UniformMagneticField<SlidingPlanarExponential<ModelInterface>>>(
1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km, B1);
CHECK(builder.size() == 2); CHECK(builder.size() == 2);
......
Processes/CONEXSourceCut/testCONEXSourceCut.cc
+ 29
− 19
View file @ bd6e3453
...@@ -35,8 +35,8 @@ using namespace corsika::environment; ...@@ -35,8 +35,8 @@ using namespace corsika::environment;
using namespace corsika::geometry; using namespace corsika::geometry;
using namespace corsika::units::si; using namespace corsika::units::si;
template <typename T> // template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>; // using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
TEST_CASE("CONEXSourceCut") { TEST_CASE("CONEXSourceCut") {
random::RNGManager::GetInstance().RegisterRandomStream("cascade"); random::RNGManager::GetInstance().RegisterRandomStream("cascade");
...@@ -56,20 +56,30 @@ TEST_CASE("CONEXSourceCut") { ...@@ -56,20 +56,30 @@ TEST_CASE("CONEXSourceCut") {
{{particles::Code::Nitrogen, particles::Code::Oxygen}, {{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now {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, builder.addExponentialLayer<
environment::EMediumType::eAir, environment::UniformMediumType<environment::UniformMagneticField<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, 1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, environment::EMediumType::eAir,
environment::EMediumType::eAir, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); builder.addExponentialLayer<
builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, environment::UniformMediumType<environment::UniformMagneticField<
environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, environment::EMediumType::eAir,
builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, geometry::Vector(rootCS, 0_T, 0_T, 1_T));
environment::EMediumType::eAir, builder.addExponentialLayer<
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); environment::UniformMediumType<environment::UniformMagneticField<
builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir, SlidingPlanarExponential<setup::EnvironmentInterface>>>>(
geometry::Vector(rootCS, 0_T, 0_T, 1_T)); 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); builder.assemble(env);
...@@ -95,8 +105,7 @@ TEST_CASE("CONEXSourceCut") { ...@@ -95,8 +105,7 @@ TEST_CASE("CONEXSourceCut") {
[[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env); [[maybe_unused]] process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
corsika::process::conex_source_cut::CONEXSourceCut conex( corsika::process::conex_source_cut::CONEXSourceCut conex(
center, showerAxis, t, injectionHeight, E0, center, showerAxis, t, injectionHeight, E0, particles::Code::Proton);
particles::GetPDG(particles::Code::Proton));
HEPEnergyType const Eem{1_PeV}; HEPEnergyType const Eem{1_PeV};
auto const momentum = showerAxis.GetDirection() * Eem; auto const momentum = showerAxis.GetDirection() * Eem;
...@@ -111,7 +120,8 @@ TEST_CASE("CONEXSourceCut") { ...@@ -111,7 +120,8 @@ TEST_CASE("CONEXSourceCut") {
std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " " std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " "
<< emPosition.GetCoordinates(rootCS) << std::endl; << 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(); conex.SolveCE();
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum