diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index 2b86077118bfbcc94ef39bbef411e098d7d055a1..adbbca31ce43a2a4cfcca7808f9930d6cbbe1b27 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 e0451e9474155c7d2fbe47b21bb1bcd6238fc988..c8fa31e6edf2c0457e5fe6617897656bac5bb5f4 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 371c4c688d0d86229a185a147b526790c23321e2..f3cd9c19aa48d5eb02c7276d10d07752f2a2cfbc 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 8898bb391cf6dab0d5c1e51aa2fcea816a64925f..594a95845c9a2e8d8640aeda4c7847ed591d9420 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 165f0b489a15abab98a8b6d8b9ca549d55d6d81a..4edfa022ebb7a128ec94b5fea172ea596e62e718 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 8fd61fcde6eec20269e0a489a1ef2e6adc588137..c6f2dbe6111e7141f2880e320c83a4bdead4a030 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 2a297005f8a726865f19a16531eefbfbdb991144..096f9b181cb9596da1b93a844792a3b07cae9ece 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 387267e7842a98dd26b4b09da1d65bdfd3ba6670..e93eb64f63835c18f50dafc9b2ce8964309e3985 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 ddbf5e9bc3634f68195558cf8ecb609cb6e2ddf4..32d3025940e78be21d15d58f5593f32d4a8f514e 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 f1b1dfcf5eee747f2007d207788ad2f686285766..9f69e73c40c61c7ad5461ac77120ec4a182e4ce2 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 eef78592967cd6a6cfe33eaf9526a5577e4de377..5ee19f4bfd6ceb74400c5a939a84e3dfe7ca00f4 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);