From 2cf9818938be06b65b82c8bf5ca666d3a0c6860b Mon Sep 17 00:00:00 2001 From: ralfulrich <ralf.ulrich@kit.edu> Date: Fri, 23 Oct 2020 21:39:16 +0200 Subject: [PATCH] make_layered_spherical_atmosphere_builder --- Documentation/Examples/CMakeLists.txt | 1 - Documentation/Examples/boundary_example.cc | 16 ++++-- Documentation/Examples/cascade_example.cc | 12 ++-- Documentation/Examples/em_shower.cc | 28 ++++------ Documentation/Examples/hybrid_MC.cc | 55 ++++++++++--------- Documentation/Examples/vertical_EAS.cc | 16 ++---- .../LayeredSphericalAtmosphereBuilder.h | 38 +++++++++++-- Environment/testEnvironment.cc | 17 ++++-- Framework/ProcessSequence/ProcessSequence.h | 5 +- .../ProcessSequence/SwitchProcessSequence.h | 9 +-- .../CONEXSourceCut/testCONEXSourceCut.cc | 29 ++++------ 11 files changed, 124 insertions(+), 102 deletions(-) diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 2b8607711..adbbca31c 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -126,7 +126,6 @@ if (Pythia8_FOUND) ProcessSibyll ProcessPythia8 ProcessUrQMD - ProcessSwitch CORSIKAcascade ProcessPythia8 ProcessObservationPlane diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc index e0451e947..c8fa31e6e 100644 --- a/Documentation/Examples/boundary_example.cc +++ b/Documentation/Examples/boundary_example.cc @@ -104,12 +104,16 @@ int main() { auto world = EnvType::CreateNode<Sphere>( Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); - auto const props = - world->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>( - 1_kg / (1_m * 1_m * 1_m), - environment::NuclearComposition( - std::vector<particles::Code>{particles::Code::Proton}, - std::vector<float>{1.f})); + using MyHomogeneousModel = + environment::MediumPropertyModel<environment::UniformMagneticField< + environment::HomogeneousMedium<setup::EnvironmentInterface>>>; + + auto const props = world->SetModelProperties<MyHomogeneousModel>( + environment::Medium::AirDry1Atm, Vector(rootCS, 0_T, 0_T, 0_T), + 1_kg / (1_m * 1_m * 1_m), + environment::NuclearComposition( + std::vector<particles::Code>{particles::Code::Proton}, + std::vector<float>{1.f})); // add a "target" sphere with 5km readius at 0,0,0 auto target = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km); diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 371c4c688..f3cd9c19a 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -56,7 +56,7 @@ using namespace corsika::units::si; // int main() { - logging::SetLevel(logging::level::debug); + logging::SetLevel(logging::level::info); std::cout << "cascade_example" << std::endl; @@ -72,7 +72,7 @@ int main() { const CoordinateSystem& rootCS = env.GetCoordinateSystem(); - auto outerMedium = setup::Environment::CreateNode<Sphere>( + auto world = setup::Environment::CreateNode<Sphere>( Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); using MyHomogeneousModel = @@ -81,7 +81,7 @@ int main() { // fraction of oxygen const float fox = 0.20946; - auto const props = outerMedium->SetModelProperties<MyHomogeneousModel>( + auto const props = world->SetModelProperties<MyHomogeneousModel>( environment::Medium::AirDry1Atm, Vector(rootCS, 0_T, 0_T, 0_T), 1_kg / (1_m * 1_m * 1_m), environment::NuclearComposition( @@ -93,10 +93,8 @@ int main() { setup::Environment::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m); innerMedium->SetModelProperties(props); - - outerMedium->AddChild(std::move(innerMedium)); - - universe.AddChild(std::move(outerMedium)); + world->AddChild(std::move(innerMedium)); + universe.AddChild(std::move(world)); // setup particle stack, and add primary particle setup::Stack stack; diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc index 8898bb391..594a95845 100644 --- a/Documentation/Examples/em_shower.cc +++ b/Documentation/Examples/em_shower.cc @@ -55,7 +55,7 @@ void registerRandomStreams() { } template <typename T> -using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; +using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; int main(int argc, char** argv) { @@ -74,26 +74,20 @@ int main(int argc, char** argv) { EnvType env; const CoordinateSystem& rootCS = env.GetCoordinateSystem(); Point const center{rootCS, 0_m, 0_m, 0_m}; - environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv> - builder{center}; + auto builder = environment::make_layered_spherical_atmosphere_builder< + setup::EnvironmentInterface, + MyExtraEnv>::create(center, units::constants::EarthRadius::Mean, + environment::Medium::AirDry1Atm, + geometry::Vector{rootCS, 0_T, 0_T, 1_T}); builder.setNuclearComposition( {{particles::Code::Nitrogen, particles::Code::Oxygen}, {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now - builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); + builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km); + builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km); + builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km); + builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km); + builder.addLinearLayer(1e9_cm, 112.8_km); builder.assemble(env); // setup particle stack, and add primary particle diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc index 165f0b489..4edfa022e 100644 --- a/Documentation/Examples/hybrid_MC.cc +++ b/Documentation/Examples/hybrid_MC.cc @@ -22,6 +22,7 @@ #include <corsika/geometry/Sphere.h> #include <corsika/logging/Logging.h> #include <corsika/process/ProcessSequence.h> +#include <corsika/process/SwitchProcessSequence.h> #include <corsika/process/StackProcess.h> #include <corsika/process/conex_source_cut/CONEXSourceCut.h> #include <corsika/process/energy_loss/EnergyLoss.h> @@ -33,7 +34,6 @@ #include <corsika/process/sibyll/Decay.h> #include <corsika/process/sibyll/Interaction.h> #include <corsika/process/sibyll/NuclearInteraction.h> -#include <corsika/process/switch_process/SwitchProcess.h> #include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/process/urqmd/UrQMD.h> #include <corsika/random/RNGManager.h> @@ -73,7 +73,7 @@ void registerRandomStreams(const int seed) { } template <typename T> -using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; +using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; int main(int argc, char** argv) { @@ -98,26 +98,20 @@ int main(int argc, char** argv) { EnvType env; const CoordinateSystem& rootCS = env.GetCoordinateSystem(); Point const center{rootCS, 0_m, 0_m, 0_m}; - environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv> - builder{center}; + auto builder = environment::make_layered_spherical_atmosphere_builder< + setup::EnvironmentInterface, + MyExtraEnv>::create(center, units::constants::EarthRadius::Mean, + environment::Medium::AirDry1Atm, + geometry::Vector{rootCS, 0_T, 0_T, 1_T}); builder.setNuclearComposition( {{particles::Code::Nitrogen, particles::Code::Oxygen}, {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now - builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); + builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km); + builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km); + builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km); + builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km); + builder.addLinearLayer(1e9_cm, 112.8_km); builder.assemble(env); // setup particle stack, and add primary particle @@ -228,14 +222,23 @@ int main(int argc, char** argv) { process::interaction_counter::InteractionCounter urqmdCounted{urqmd}; // assemble all processes into an ordered process list - - auto sibyllSequence = sibyllNucCounted << sibyllCounted; - process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence, - 55_GeV); - auto decaySequence = decayPythia << decaySibyll; - - auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut - << conex << longprof << observationLevel; + struct EnergySwitch { + HEPEnergyType cutE_; + EnergySwitch(HEPEnergyType cutE) + : cutE_(cutE) {} + process::SwitchResult operator()(const setup::Stack::ParticleType& p) { + if (p.GetEnergy() < cutE_) + return process::SwitchResult::First; + else + return process::SwitchResult::Second; + } + }; + auto hadronSequence = + process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted), + EnergySwitch(55_GeV)); + auto decaySequence = process::sequence(decayPythia, decaySibyll); + auto sequence = process::sequence(hadronSequence, reset_particle_mass, decaySequence, + eLoss, cut, conex, longprof, observationLevel); // define air shower object, run simulation tracking_line::TrackingLine tracking; diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 8fd61fcde..c6f2dbe61 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -76,14 +76,7 @@ void registerRandomStreams(const int seed) { } template <typename T> -using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; - -template <typename... TArgs> -auto make_builder(geometry::Point const& center, TArgs... args) { - return environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv, - TArgs...>{ - args..., center, units::constants::EarthRadius::Mean}; -} +using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; int main(int argc, char** argv) { @@ -108,8 +101,11 @@ int main(int argc, char** argv) { EnvType env; const CoordinateSystem& rootCS = env.GetCoordinateSystem(); Point const center{rootCS, 0_m, 0_m, 0_m}; - auto builder = make_builder(center, environment::Medium::AirDry1Atm, - geometry::Vector{rootCS, 0_T, 0_T, 1_T}); + auto builder = environment::make_layered_spherical_atmosphere_builder< + setup::EnvironmentInterface, + MyExtraEnv>::create(center, units::constants::EarthRadius::Mean, + environment::Medium::AirDry1Atm, + geometry::Vector{rootCS, 0_T, 0_T, 1_T}); builder.setNuclearComposition( {{particles::Code::Nitrogen, particles::Code::Oxygen}, {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now diff --git a/Environment/LayeredSphericalAtmosphereBuilder.h b/Environment/LayeredSphericalAtmosphereBuilder.h index 2a297005f..096f9b181 100644 --- a/Environment/LayeredSphericalAtmosphereBuilder.h +++ b/Environment/LayeredSphericalAtmosphereBuilder.h @@ -26,6 +26,10 @@ namespace corsika::environment { + // forward-decl + template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment> + struct make_layered_spherical_atmosphere_builder; + /** * Helper class to setup concentric spheres of layered atmosphere * with spcified density profiles (exponential, linear, ...). @@ -78,14 +82,18 @@ namespace corsika::environment { LayeredSphericalAtmosphereBuilder& operator=( const LayeredSphericalAtmosphereBuilder&) = delete; - public: - LayeredSphericalAtmosphereBuilder( - TModelArgs... args, corsika::geometry::Point center, - units::si::LengthType earthRadius = units::constants::EarthRadius::Mean) + // friend, to allow construction + template <typename, template <typename> typename> + friend struct make_layered_spherical_atmosphere_builder; + + protected: + LayeredSphericalAtmosphereBuilder(TModelArgs... args, corsika::geometry::Point center, + units::si::LengthType earthRadius) : center_(center) , earthRadius_(earthRadius) , additionalModelArgs_{args...} {} + public: void setNuclearComposition(NuclearComposition composition) { composition_ = std::make_unique<NuclearComposition>(composition); } @@ -186,4 +194,26 @@ namespace corsika::environment { }; // end class LayeredSphericalAtmosphereBuilder + /** + * \class make_layered_spherical_atmosphere_builder + * + * Helper class to create LayeredSphericalAtmosphereBuilder, the + * extra environment models have to be passed as template-template + * argument to make_layered_spherical_atmosphere_builder, the member + * function `create` does then take an unspecified number of extra + * parameters to internalize those models for all layers later + * produced. + **/ + template <typename TMediumInterface = environment::IMediumModel, + template <typename> typename MExtraEnvirnoment = detail::NoExtraModel> + struct make_layered_spherical_atmosphere_builder { + template <typename... TArgs> + static auto create(geometry::Point const& center, units::si::LengthType earthRadius, + TArgs... args) { + return environment::LayeredSphericalAtmosphereBuilder<TMediumInterface, + MExtraEnvirnoment, TArgs...>{ + std::forward<TArgs>(args)..., center, earthRadius}; + } + }; + } // namespace corsika::environment diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc index 387267e78..e93eb64f6 100644 --- a/Environment/testEnvironment.cc +++ b/Environment/testEnvironment.cc @@ -202,7 +202,11 @@ TEST_CASE("InhomogeneousMedium") { } TEST_CASE("LayeredSphericalAtmosphereBuilder") { - LayeredSphericalAtmosphereBuilder builder(gOrigin, units::constants::EarthRadius::Mean); + + LayeredSphericalAtmosphereBuilder builder = + environment::make_layered_spherical_atmosphere_builder<>::create( + gOrigin, units::constants::EarthRadius::Mean); + builder.setNuclearComposition( {{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}}); @@ -307,12 +311,15 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") { // create magnetic field vectors Vector B0(gCS, 0_T, 0_T, 1_T); - LayeredSphericalAtmosphereBuilder<ModelInterface, UniformMagneticField> builder{ - gOrigin}; + LayeredSphericalAtmosphereBuilder builder = + environment::make_layered_spherical_atmosphere_builder< + ModelInterface, + UniformMagneticField>::create(gOrigin, units::constants::EarthRadius::Mean, B0); + builder.setNuclearComposition( {{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}}); - builder.addLinearLayer(1_km, 10_km, B0); - builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km, B0); + builder.addLinearLayer(1_km, 10_km); + builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km); CHECK(builder.size() == 2); diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index ddbf5e9bc..32d302594 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -50,8 +50,11 @@ namespace corsika::process { static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>; static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>; + static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<TProcess1type>; + static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<TProcess2type>; + // make sure only BaseProcess types TProcess1/2 are passed - static_assert(std::is_base_of_v<BaseProcess<TProcess11ype>, TProcess11ype>, + static_assert(std::is_base_of_v<BaseProcess<TProcess1type>, TProcess1type>, "can only use process derived from BaseProcess in " "ProcessSequence, for Process 1"); static_assert(std::is_base_of_v<BaseProcess<TProcess2type>, TProcess2type>, diff --git a/Framework/ProcessSequence/SwitchProcessSequence.h b/Framework/ProcessSequence/SwitchProcessSequence.h index f1b1dfcf5..9f69e73c4 100644 --- a/Framework/ProcessSequence/SwitchProcessSequence.h +++ b/Framework/ProcessSequence/SwitchProcessSequence.h @@ -72,9 +72,6 @@ namespace corsika::process { "can only use process derived from BaseProcess in " "SwitchProcessSequence, for Process 2"); - // make sure TSelect is a function - static_assert(!std::is_function_v<TSelect>, "TSelect must be a function type"); - // make sure none of TProcess1/2 is a StackProcess static_assert(!std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type>, "cannot use StackProcess in SwitchProcessSequence, for Process 1"); @@ -82,14 +79,10 @@ namespace corsika::process { "cannot use StackProcess in SwitchProcessSequence, for Process 2"); // if TProcess1/2 are already ProcessSequences, make sure they do not contain - // StackProcess - // if constexpr (t1ProcSeq) { + // any StackProcess static_assert(!contains_stack_process_v<TProcess1type>, "cannot use StackProcess in SwitchProcessSequence, remove from " "ProcessSequence 1"); - //} - - // if constexpr (t2ProcSeq) static_assert(!contains_stack_process_v<TProcess2type>, "cannot use StackProcess in SwitchProcessSequence, remove from " "ProcessSequence 2"); diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc index eef785929..5ee19f4bf 100644 --- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc +++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc @@ -36,7 +36,8 @@ using namespace corsika::geometry; using namespace corsika::units::si; template <typename T> -using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>; +using MExtraEnvirnoment = + environment::MediumPropertyModel<environment::UniformMagneticField<T>>; TEST_CASE("CONEXSourceCut") { random::RNGManager::GetInstance().RegisterRandomStream("cascade"); @@ -49,27 +50,21 @@ TEST_CASE("CONEXSourceCut") { const CoordinateSystem& rootCS = env.GetCoordinateSystem(); Point const center{rootCS, 0_m, 0_m, 0_m}; - environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv> - builder{center, conex::earthRadius}; + auto builder = environment::make_layered_spherical_atmosphere_builder< + setup::EnvironmentInterface, + MExtraEnvirnoment>::create(center, conex::earthRadius, + environment::Medium::AirDry1Atm, + geometry::Vector{rootCS, 0_T, 50_mT, 0_T}); builder.setNuclearComposition( {{particles::Code::Nitrogen, particles::Code::Oxygen}, {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now - builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km, - environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); - builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm, - geometry::Vector(rootCS, 0_T, 0_T, 1_T)); + builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km); + builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km); + builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km); + builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km); + builder.addLinearLayer(1e9_cm, 112.8_km); builder.assemble(env); -- GitLab