diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 75f29e4ed51bc70d7d31d4917ef4c186baee4bef..f9ec0beb0e40df08a16595bee6be8dc0048cff94 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -153,7 +153,7 @@ namespace corsika::cascade { // determine combined total interaction length (inverse) InverseGrammageType const total_inv_lambda = - fProcessSequence.GetTotalInverseInteractionLength(vParticle, step); + fProcessSequence.GetTotalInverseInteractionLength(vParticle); // sample random exponential step length in grammage corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda); @@ -242,13 +242,14 @@ namespace corsika::cascade { if (min_distance == distance_interact) { std::cout << "collide" << std::endl; - InverseGrammageType const actual_inv_length = - fProcessSequence.GetTotalInverseInteractionLength(vParticle, step); + InverseGrammageType const current_inv_length = + fProcessSequence.GetTotalInverseInteractionLength(vParticle); - random::UniformRealDistribution<InverseGrammageType> uniDist(actual_inv_length); + random::UniformRealDistribution<InverseGrammageType> uniDist( + current_inv_length); const auto sample_process = uniDist(fRNG); InverseGrammageType inv_lambda_count = 0. * meter * meter / gram; - fProcessSequence.SelectInteraction(vParticle, projectile, step, sample_process, + fProcessSequence.SelectInteraction(vParticle, projectile, sample_process, inv_lambda_count); } else if (min_distance == distance_decay) { std::cout << "decay" << std::endl; diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 89fd65ef42150149a9c0ffbd7183fed865ba659d..f29d1af9fba4c01c6452fc2dc026df7ff99f2c42 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -69,8 +69,8 @@ public: ProcessSplit(GrammageType const X0) : fX0(X0) {} - template <typename Particle, typename Track> - corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const { + template <typename Particle> + corsika::units::si::GrammageType GetInteractionLength(Particle const&) const { return fX0; } diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h index 27ad9c222d5ce644d6b344a469edaf8e717fcf87..a9b4b9fca306e865a85ac94ae6a9878ce7f0057d 100644 --- a/Framework/ProcessSequence/InteractionProcess.h +++ b/Framework/ProcessSequence/InteractionProcess.h @@ -38,13 +38,12 @@ namespace corsika::process { template <typename Particle> EProcessReturn DoInteraction(Particle&); - template <typename Particle, typename Track> - corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t); + template <typename TParticle> + corsika::units::si::GrammageType GetInteractionLength(TParticle& p); - template <typename Particle, typename Track> - corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p, - Track& t) { - return 1. / GetRef().GetInteractionLength(p, t); + template <typename TParticle> + corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) { + return 1. / GetRef().GetInteractionLength(p); } }; diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index 3bd759c04a97e55360c87473cc4457159bd967fa..55eb43ba795c35d04654a5256232221fa4c900b5 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -47,9 +47,19 @@ namespace corsika::process { // this is a marker to track which BaseProcess is also a ProcessSequence template <typename T> - struct is_process_sequence { - static const bool value = false; - }; + struct is_process_sequence : std::false_type {}; + + template <typename T> + bool constexpr is_process_sequence_v = is_process_sequence<T>::value; + + namespace switch_process { + template <typename A, typename B> + class SwitchProcess; // fwd-decl. + } + + template <typename A, typename B> + struct is_process_sequence<switch_process::SwitchProcess<A, B>> + : std::true_type {}; /** T1 and T2 are both references if possible (lvalue), otherwise @@ -57,11 +67,14 @@ namespace corsika::process { rvalue as well as lvalue Processes in the ProcessSequence. */ template <typename T1, typename T2> - class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > { + class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2>> { using T1type = typename std::decay<T1>::type; using T2type = typename std::decay<T2>::type; + static bool constexpr t1ProcSeq = is_process_sequence_v<T1type>; + static bool constexpr t2ProcSeq = is_process_sequence_v<T2type>; + public: T1 A; // this is a reference, if possible T2 B; // this is a reference, if possible @@ -78,13 +91,13 @@ namespace corsika::process { VTNType const& to) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of<BoundaryCrossingProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T1type>, T1type> || + t1ProcSeq) { ret |= A.DoBoundaryCrossing(p, from, to); } - if constexpr (std::is_base_of<BoundaryCrossingProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T2type>, T2type> || + t2ProcSeq) { ret |= B.DoBoundaryCrossing(p, from, to); } @@ -94,12 +107,10 @@ namespace corsika::process { template <typename TParticle, typename TTrack> EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) { ret |= A.DoContinuous(vP, vT); } - if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) { ret |= B.DoContinuous(vP, vT); } return ret; @@ -108,12 +119,10 @@ namespace corsika::process { template <typename TSecondaries> EProcessReturn DoSecondaries(TSecondaries& vS) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of<SecondariesProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<SecondariesProcess<T1type>, T1type> || t1ProcSeq) { ret |= A.DoSecondaries(vS); } - if constexpr (std::is_base_of<SecondariesProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<SecondariesProcess<T2type>, T2type> || t2ProcSeq) { ret |= B.DoSecondaries(vS); } return ret; @@ -121,12 +130,10 @@ namespace corsika::process { bool CheckStep() { bool ret = false; - if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) { ret |= A.CheckStep(); } - if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) { ret |= B.CheckStep(); } return ret; @@ -135,12 +142,10 @@ namespace corsika::process { template <typename TStack> EProcessReturn DoStack(TStack& vS) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of<StackProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) { if (A.CheckStep()) { ret |= A.DoStack(vS); } } - if constexpr (std::is_base_of<StackProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) { if (B.CheckStep()) { ret |= B.DoStack(vS); } } return ret; @@ -152,64 +157,58 @@ namespace corsika::process { max_length = // if no other process in the sequence implements it std::numeric_limits<double>::infinity() * corsika::units::si::meter; - if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) { corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack); max_length = std::min(max_length, len); } - if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) { corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack); max_length = std::min(max_length, len); } return max_length; } - template <typename TParticle, typename TTrack> - corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP, - TTrack& vT) { - return 1. / GetInverseInteractionLength(vP, vT); + template <typename TParticle> + corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) { + return 1. / GetInverseInteractionLength(vP); } - template <typename TParticle, typename TTrack> + template <typename TParticle> corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength( - TParticle& vP, TTrack& vT) { - return GetInverseInteractionLength(vP, vT); + TParticle& vP) { + return GetInverseInteractionLength(vP); } - template <typename TParticle, typename TTrack> - corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP, - TTrack& vT) { + template <typename TParticle> + corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) { using namespace corsika::units::si; InverseGrammageType tot = 0 * meter * meter / gram; - if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { - tot += A.GetInverseInteractionLength(vP, vT); + if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> || t1ProcSeq) { + tot += A.GetInverseInteractionLength(vP); } - if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { - tot += B.GetInverseInteractionLength(vP, vT); + if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> || t2ProcSeq) { + tot += B.GetInverseInteractionLength(vP); } return tot; } - template <typename TParticle, typename TSecondaries, typename TTrack> + template <typename TParticle, typename TSecondaries> EProcessReturn SelectInteraction( - TParticle& vP, TSecondaries& vS, TTrack& vT, + TParticle& vP, TSecondaries& vS, [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select, corsika::units::si::InverseGrammageType& lambda_inv_count) { - if constexpr (is_process_sequence<T1type>::value) { + if constexpr (t1ProcSeq) { // if A is a process sequence --> check inside const EProcessReturn ret = - A.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count); + A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count); // if A did succeed, stop routine if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) { + } else if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type>) { // if this is not a ContinuousProcess --> evaluate probability - lambda_inv_count += A.GetInverseInteractionLength(vP, vT); + lambda_inv_count += A.GetInverseInteractionLength(vP); // check if we should execute THIS process and then EXIT if (lambda_select < lambda_inv_count) { A.DoInteraction(vS); @@ -217,15 +216,15 @@ namespace corsika::process { } } // end branch A - if constexpr (is_process_sequence<T2>::value) { + if constexpr (t2ProcSeq) { // if A is a process sequence --> check inside const EProcessReturn ret = - B.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count); + B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count); // if A did succeed, stop routine if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) { + } else if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type>) { // if this is not a ContinuousProcess --> evaluate probability - lambda_inv_count += B.GetInverseInteractionLength(vP, vT); + lambda_inv_count += B.GetInverseInteractionLength(vP); // check if we should execute THIS process and then EXIT if (lambda_select < lambda_inv_count) { B.DoInteraction(vS); @@ -251,12 +250,10 @@ namespace corsika::process { corsika::units::si::InverseTimeType tot = 0 / second; - if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value || - is_process_sequence<T1>::value) { + if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type> || t1ProcSeq) { tot += A.GetInverseLifetime(p); } - if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value || - is_process_sequence<T2>::value) { + if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type> || t2ProcSeq) { tot += B.GetInverseLifetime(p); } return tot; @@ -268,12 +265,12 @@ namespace corsika::process { TParticle& vP, TSecondaries& vS, [[maybe_unused]] corsika::units::si::InverseTimeType decay_select, corsika::units::si::InverseTimeType& decay_inv_count) { - if constexpr (is_process_sequence<T1>::value) { + if constexpr (t1ProcSeq) { // if A is a process sequence --> check inside const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count); // if A did succeed, stop routine if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value) { + } else if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type>) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_count += A.GetInverseLifetime(vP); // check if we should execute THIS process and then EXIT @@ -284,12 +281,12 @@ namespace corsika::process { } } // end branch A - if constexpr (is_process_sequence<T2>::value) { + if constexpr (t2ProcSeq) { // if A is a process sequence --> check inside const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count); // if A did succeed, stop routine if (ret != EProcessReturn::eOk) { return ret; } - } else if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value) { + } else if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type>) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_count += B.GetInverseLifetime(vP); // check if we should execute THIS process and then EXIT @@ -322,9 +319,7 @@ namespace corsika::process { /// marker to identify objectas ProcessSequence template <typename A, typename B> - struct is_process_sequence<corsika::process::ProcessSequence<A, B> > { - static const bool value = true; - }; + struct is_process_sequence<corsika::process::ProcessSequence<A, B>> : std::true_type {}; } // namespace corsika::process diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc index aa78ee89c4c012eb2df3d322b6b1ddf6fa44ba1e..e8bfe7c1c835762e64ccb3afed48962775dc24ae 100644 --- a/Framework/ProcessSequence/testProcessSequence.cc +++ b/Framework/ProcessSequence/testProcessSequence.cc @@ -100,8 +100,8 @@ public: cout << "Process2::DoInteraction" << endl; return EProcessReturn::eOk; } - template <typename Particle, typename Track> - GrammageType GetInteractionLength(Particle&, Track&) const { + template <typename Particle> + GrammageType GetInteractionLength(Particle&) const { cout << "Process2::GetInteractionLength" << endl; return 3_g / (1_cm * 1_cm); } @@ -123,8 +123,8 @@ public: cout << "Process3::DoInteraction" << endl; return EProcessReturn::eOk; } - template <typename Particle, typename Track> - GrammageType GetInteractionLength(Particle&, Track&) const { + template <typename Particle> + GrammageType GetInteractionLength(Particle&) const { cout << "Process3::GetInteractionLength" << endl; return 1_g / (1_cm * 1_cm); } @@ -220,12 +220,11 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Process3 m3(2); DummyData particle; - DummyTrajectory track; auto sequence2 = cp1 << m2 << m3; - GrammageType const tot = sequence2.GetTotalInteractionLength(particle, track); + GrammageType const tot = sequence2.GetTotalInteractionLength(particle); InverseGrammageType const tot_inv = - sequence2.GetTotalInverseInteractionLength(particle, track); + sequence2.GetTotalInverseInteractionLength(particle); cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl; } diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index 2835a7043cc80101ef22d1301af981b5735de6bb..6afe1cc51acdc689b97a90dec334d727b61c7f60 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -10,6 +10,7 @@ add_subdirectory (TrackingLine) add_subdirectory (HadronicElasticModel) add_subdirectory (EnergyLoss) add_subdirectory (UrQMD) +add_subdirectory (SwitchProcess) #add_custom_target(CORSIKAprocesses) add_library (CORSIKAprocesses INTERFACE) diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc index 33f162e41531bf4008f91b21bb98e8afc13a7abb..9e05c521320787b1f8a42ca3c154cfa708285769 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.cc +++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc @@ -18,14 +18,12 @@ #include <corsika/utl/COMBoost.h> #include <corsika/setup/SetupStack.h> -#include <corsika/setup/SetupTrajectory.h> #include <iomanip> #include <iostream> using namespace corsika::setup; -using Particle = corsika::setup::Stack::ParticleType; -using Track = corsika::setup::Trajectory; +using SetupParticle = corsika::setup::Stack::ParticleType; namespace corsika::process::HadronicElasticModel { @@ -38,7 +36,7 @@ namespace corsika::process::HadronicElasticModel { template <> units::si::GrammageType HadronicElasticInteraction::GetInteractionLength( - Particle const& p, Track&) { + SetupParticle const& p) { using namespace units::si; if (p.GetPID() == particles::Code::Proton) { auto const* currentNode = p.GetNode(); @@ -79,7 +77,7 @@ namespace corsika::process::HadronicElasticModel { } template <> - process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p) { + process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupParticle& p) { if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; } using namespace units::si; diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h index 4180c04c514d00089fd044ab39027b9128a853ff..b3fd15fd2cdc6feac967182cff342b60f30e954d 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.h +++ b/Processes/HadronicElasticModel/HadronicElasticModel.h @@ -56,8 +56,8 @@ namespace corsika::process::HadronicElasticModel { units::si::CrossSectionType y = 0.05608 * units::si::barn); void Init(); - template <typename Particle, typename Track> - corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&); + template <typename Particle> + corsika::units::si::GrammageType GetInteractionLength(Particle const& p); template <typename Particle> corsika::process::EProcessReturn DoInteraction(Particle&); diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc index 6cdf54973a2d553a79f6b6ebcd97b4d5ecf6cf24..2a0cce386edc6dfc5a14f0f6fb84fd94aa5a85a7 100644 --- a/Processes/Pythia/Interaction.cc +++ b/Processes/Pythia/Interaction.cc @@ -15,7 +15,6 @@ #include <corsika/environment/NuclearComposition.h> #include <corsika/geometry/FourVector.h> #include <corsika/setup/SetupStack.h> -#include <corsika/setup/SetupTrajectory.h> #include <corsika/utl/COMBoost.h> #include <tuple> @@ -28,7 +27,6 @@ using namespace corsika; using namespace corsika::setup; using Projectile = corsika::setup::StackView::ParticleType; using Particle = corsika::setup::Stack::ParticleType; -using Track = Trajectory; namespace corsika::process::pythia { @@ -168,7 +166,7 @@ namespace corsika::process::pythia { } template <> - units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) { + units::si::GrammageType Interaction::GetInteractionLength(Particle& p) { using namespace units; using namespace units::si; diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h index ac27825acf0561d0d9e7e8e8b66961bdfc674e2f..3594358e0b39f528991183b4c414e8edaa3b5828 100644 --- a/Processes/Pythia/Interaction.h +++ b/Processes/Pythia/Interaction.h @@ -51,16 +51,16 @@ namespace corsika::process::pythia { const corsika::particles::Code TargetId, const corsika::units::si::HEPEnergyType CoMenergy); - template <typename TParticle, typename TTrack> - corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&); + template <typename TParticle> + corsika::units::si::GrammageType GetInteractionLength(TParticle&); /** In this function PYTHIA is called to produce one event. The event is copied (and boosted) into the shower lab frame. */ - template <typename TProjctile> - corsika::process::EProcessReturn DoInteraction(TProjctile&); + template <typename TProjectile> + corsika::process::EProcessReturn DoInteraction(TProjectile&); private: corsika::random::RNG& fRNG = diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc index 3fb2df13ecd1c24e674997a3c3aa468164f21d06..f6f09b84c17f80f8b579d1a11eba843e6ddbc964 100644 --- a/Processes/Pythia/testPythia.cc +++ b/Processes/Pythia/testPythia.cc @@ -118,14 +118,6 @@ TEST_CASE("pythia process") { auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it universe.AddChild(std::move(theMedium)); - geometry::Point const origin(cs, {0_m, 0_m, 0_m}); - geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second, - 1_m / second); - geometry::Line line(origin, v); - geometry::Trajectory<geometry::Line> track(line, 10_s); - - random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); - SECTION("pythia decay") { setup::Stack stack; @@ -173,7 +165,6 @@ TEST_CASE("pythia process") { process::pythia::Interaction model; model.Init(); [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); - [[maybe_unused]] const GrammageType length = - model.GetInteractionLength(particle, track); + [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); } } diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index b15d9e8859e86d6426f6b78b7834fe2211f191f1..c089ec083c9c7d9ec1d91136cb0c7796dc958c57 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -118,8 +118,8 @@ namespace corsika::process::sibyll { } template <> - units::si::GrammageType Interaction::GetInteractionLength(SetupParticle& vP, - Track&) const { + units::si::GrammageType Interaction::GetInteractionLength( + SetupParticle const& vP) const { using namespace units; using namespace units::si; diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h index e7d7a0902116c2fceb9240cf9266d04a2ea2fd72..643bd0d51bd4f387aafa3870f2dddfb0ed1fe9d5 100644 --- a/Processes/Sibyll/Interaction.h +++ b/Processes/Sibyll/Interaction.h @@ -52,8 +52,8 @@ namespace corsika::process::sibyll { GetCrossSection(const corsika::particles::Code, const corsika::particles::Code, const corsika::units::si::HEPEnergyType) const; - template <typename TParticle, typename TTrack> - corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&) const; + template <typename TParticle> + corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const; /** In this function SIBYLL is called to produce one event. The diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index 31abd9da350e1c80c24868eb2313e860b51bd821..907e95b271bb44cfb4d2982f0f1f40f88cbb0eca 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -217,7 +217,7 @@ namespace corsika::process::sibyll { template <> template <> units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength( - Particle& vP, Track&) { + Particle& vP) { using namespace units; using namespace units::si; diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h index 6bee9ba7a7964e74112c2bf60e5cbd4b1c23c2fa..261f931321ed98119ab02f4e3cace952308fdbb7 100644 --- a/Processes/Sibyll/NuclearInteraction.h +++ b/Processes/Sibyll/NuclearInteraction.h @@ -53,11 +53,11 @@ namespace corsika::process::sibyll { std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType> GetCrossSection(Particle& p, const corsika::particles::Code TargetId); - template <typename Particle, typename Track> - corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&); + template <typename Particle> + corsika::units::si::GrammageType GetInteractionLength(Particle&); - template <typename Projectle> - corsika::process::EProcessReturn DoInteraction(Projectle&); + template <typename Projectile> + corsika::process::EProcessReturn DoInteraction(Projectile&); private: TEnvironment const& fEnvironment; diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index d5f5d9bc8e39ac775de3fb9d541a18b596941bad..0c42c90abb86a758c45f5ad4a79ebe87d3f1d3fa 100644 --- a/Processes/Sibyll/testSibyll.cc +++ b/Processes/Sibyll/testSibyll.cc @@ -98,12 +98,6 @@ TEST_CASE("SibyllInterface", "[processes]") { const geometry::CoordinateSystem& cs = env.GetCoordinateSystem(); - geometry::Point const origin(cs, {0_m, 0_m, 0_m}); - geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second, - 1_m / second); - geometry::Line line(origin, v); - geometry::Trajectory<geometry::Line> track(line, 10_s); - random::RNGManager::GetInstance().RegisterRandomStream("s_rndm"); SECTION("InteractionInterface") { @@ -126,8 +120,7 @@ TEST_CASE("SibyllInterface", "[processes]") { model.Init(); [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); - [[maybe_unused]] const GrammageType length = - model.GetInteractionLength(particle, track); + [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); } SECTION("NuclearInteractionInterface") { @@ -153,8 +146,7 @@ TEST_CASE("SibyllInterface", "[processes]") { model.Init(); [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); - [[maybe_unused]] const GrammageType length = - model.GetInteractionLength(particle, track); + [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); } SECTION("DecayInterface") { diff --git a/Processes/SwitchProcess/CMakeLists.txt b/Processes/SwitchProcess/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..feac9f360ced5356d90ecd87cc7d4ca51939e0f6 --- /dev/null +++ b/Processes/SwitchProcess/CMakeLists.txt @@ -0,0 +1,50 @@ +set ( + MODEL_HEADERS + SwitchProcess.h + ) + +set ( + MODEL_NAMESPACE + corsika/process/switch_process + ) + +add_library (ProcessSwitch INTERFACE) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSwitch ${MODEL_NAMESPACE} ${MODEL_HEADERS}) + +set_target_properties ( + ProcessStackInspector + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 +# PUBLIC_HEADER "${MODEL_HEADERS}" + ) + +# target dependencies on other libraries (also the header onlys) +target_link_libraries ( + ProcessSwitch + INTERFACE + CORSIKAunits + CORSIKAprocesssequence + ) + +target_include_directories ( + ProcessSwitch + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE}) + +# -------------------- +# code unit testing +add_executable (testSwitchProcess testSwitchProcess.cc) + +target_link_libraries ( + testSwitchProcess + ProcessSwitch + CORSIKAstackinterface + CORSIKAthirdparty # for catch2 +) + +CORSIKA_ADD_TEST(testSwitchProcess) diff --git a/Processes/SwitchProcess/SwitchProcess.h b/Processes/SwitchProcess/SwitchProcess.h new file mode 100644 index 0000000000000000000000000000000000000000..31e59599bc6dceb68edf4ad030d284eb37c69bce --- /dev/null +++ b/Processes/SwitchProcess/SwitchProcess.h @@ -0,0 +1,98 @@ +#ifndef _corsika_SwitchProcess_h +#define _corsika_SwitchProcess_h + +#include <corsika/process/InteractionProcess.h> +#include <corsika/process/ProcessSequence.h> +#include <corsika/setup/SetupStack.h> +#include <corsika/units/PhysicalUnits.h> + +namespace corsika::process::switch_process { + + /** + * This process provides an energy-based switch between two interaction processes P1 and + * P1. For energies below the threshold, P1 is invoked, otherwise P2. Both can be either + * single interaction processes or multiple ones combined in a ProcessSequence. A + * SwitchProcess itself will always be regarded as a ProcessSequence rather than an + * InteractionProcess when assembled into a greater ProcessSequence. + */ + + template <class TLowEProcess, class THighEProcess> + class SwitchProcess : public BaseProcess<SwitchProcess<TLowEProcess, THighEProcess>> { + TLowEProcess& fLowEProcess; + THighEProcess& fHighEProcess; + units::si::HEPEnergyType const fThresholdEnergy; + + public: + SwitchProcess(TLowEProcess& vLowEProcess, THighEProcess& vHighEProcess, + units::si::HEPEnergyType vThresholdEnergy) + : fLowEProcess(vLowEProcess) + , fHighEProcess(vHighEProcess) + , fThresholdEnergy(vThresholdEnergy) {} + + void Init() { + fLowEProcess.Init(); + fHighEProcess.Init(); + } + + template <typename TParticle> + corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) { + return 1 / GetInteractionLength(p); + } + + template <typename TParticle> + units::si::GrammageType GetInteractionLength(TParticle& vParticle) { + if (vParticle.GetEnergy() < fThresholdEnergy) { + if constexpr (is_process_sequence_v<TLowEProcess>) { + return fLowEProcess.GetTotalInteractionLength(vParticle); + } else { + return fLowEProcess.GetInteractionLength(vParticle); + } + } else { + if constexpr (is_process_sequence_v<THighEProcess>) { + return fHighEProcess.GetTotalInteractionLength(vParticle); + } else { + return fHighEProcess.GetInteractionLength(vParticle); + } + } + } + + // required to conform to ProcessSequence interface. We cannot just + // implement DoInteraction() because we want to call SelectInteraction + // in case a member process is a ProcessSequence. + template <typename TParticle, typename TSecondaries> + EProcessReturn SelectInteraction( + TParticle& vP, TSecondaries& vS, + [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select, + corsika::units::si::InverseGrammageType& lambda_inv_count) { + if (vP.GetEnergy() < fThresholdEnergy) { + if constexpr (is_process_sequence_v<TLowEProcess>) { + return fLowEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count); + } else { + lambda_inv_count += fLowEProcess.GetInverseInteractionLength(vP); + // check if we should execute THIS process and then EXIT + if (lambda_select < lambda_inv_count) { + fLowEProcess.DoInteraction(vS); + return EProcessReturn::eInteracted; + } else { + return EProcessReturn::eOk; + } + } + } else { + if constexpr (is_process_sequence_v<THighEProcess>) { + return fHighEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count); + } else { + lambda_inv_count += fHighEProcess.GetInverseInteractionLength(vP); + // check if we should execute THIS process and then EXIT + if (lambda_select < lambda_inv_count) { + fHighEProcess.DoInteraction(vS); + return EProcessReturn::eInteracted; + } else { + return EProcessReturn::eOk; + } + } + } + } + }; +} // namespace corsika::process::switch_process + +#endif diff --git a/Processes/SwitchProcess/testSwitchProcess.cc b/Processes/SwitchProcess/testSwitchProcess.cc new file mode 100644 index 0000000000000000000000000000000000000000..4e902e8857367bfbb4d1ab7643941cdcce27c04c --- /dev/null +++ b/Processes/SwitchProcess/testSwitchProcess.cc @@ -0,0 +1,256 @@ +#include <corsika/process/switch_process/SwitchProcess.h> +#include <corsika/stack/SecondaryView.h> +#include <corsika/stack/Stack.h> +#include <corsika/units/PhysicalUnits.h> + +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file +#include <catch2/catch.hpp> + +#include <algorithm> +#include <random> + +using namespace corsika; +using namespace corsika::process; +using namespace corsika::units::si; + +class TestStackData { + +public: + // these functions are needed for the Stack interface + void Init() {} + void Clear() { fData.clear(); } + unsigned int GetSize() const { return fData.size(); } + unsigned int GetCapacity() const { return fData.size(); } + void Copy(int i1, int i2) { fData[i2] = fData[i1]; } + void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); } + + // custom data access function + void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; } + HEPEnergyType GetData(const unsigned int i) const { return fData[i]; } + + // these functions are also needed by the Stack interface + void IncrementSize() { fData.resize(fData.size() + 1); } + void DecrementSize() { + if (fData.size() > 0) { fData.pop_back(); } + } + + // custom private data section +private: + std::vector<HEPEnergyType> fData; +}; + +/** + * From static_cast of a StackIterator over entries in the + * TestStackData class you get and object of type + * TestParticleInterface defined here + * + * It provides Get/Set methods to read and write data to the "Data" + * storage of TestStackData obtained via + * "StackIteratorInterface::GetStackData()", given the index of the + * iterator "StackIteratorInterface::GetIndex()" + * + */ +template <typename StackIteratorInterface> +class TestParticleInterface + : public corsika::stack::ParticleBase<StackIteratorInterface> { + +public: + using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData; + using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex; + + /* + The SetParticleData methods are called for creating new entries + on the stack. You can specifiy various parametric versions to + perform this task: + */ + + // default version for particle-creation from input data + void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); } + void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/, + std::tuple<HEPEnergyType> v) { + SetEnergy(std::get<0>(v)); + } + + // here are the fundamental methods for access to TestStackData data + void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); } + HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); } +}; + +using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>; + +// see issue 161 +#if defined(__clang__) +using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>; +#elif defined(__GNUC__) || defined(__GNUG__) +using StackTestView = corsika::stack::MakeView<SimpleStack>::type; +#endif + +auto constexpr kgMSq = 1_kg / (1_m * 1_m); + +template <int N> +struct DummyProcess : InteractionProcess<DummyProcess<N>> { + void Init() {} + + template <typename TParticle> + corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const { + return N * kgMSq; + } + + template <typename TSecondaries> + corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) { + // to figure out which process was selected in the end, we produce N + // secondaries for DummyProcess<N> + + for (int i = 0; i < N; ++i) { + vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N}); + } + + return EProcessReturn::eOk; + } +}; + +using DummyLowEnergyProcess = DummyProcess<1>; +using DummyHighEnergyProcess = DummyProcess<2>; +using DummyAdditionalProcess = DummyProcess<3>; + +TEST_CASE("SwitchProcess from InteractionProcess") { + DummyLowEnergyProcess low; + DummyHighEnergyProcess high; + DummyAdditionalProcess proc; + + switch_process::SwitchProcess switchProcess(low, high, 1_TeV); + auto seq = switchProcess << proc; + + SimpleStack stack; + + SECTION("low energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV}); + auto p = stack.GetNextParticle(); + + // low energy process returns 1 kg/m² + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1)); + REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4)); + } + + // low energy process creates 1 secondary + //~ SECTION("SelectInteraction") { + //~ typename SimpleStack::ParticleType theParticle = + //~ stack.GetNextParticle(); // as in corsika::Cascade + //~ StackTestView view(theParticle); + //~ auto projectile = view.GetProjectile(); + + //~ InverseGrammageType invLambda = 0 / kgMSq; + //~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda); + + //~ REQUIRE(view.GetSize() == 1); + //~ } + } + + SECTION("high energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV}); + auto p = stack.GetNextParticle(); + + // high energy process returns 2 kg/m² + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2)); + REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(6. / 5)); + } + + // high energy process creates 2 secondaries + SECTION("SelectInteraction") { + typename SimpleStack::ParticleType theParticle = + stack.GetNextParticle(); // as in corsika::Cascade + StackTestView view(theParticle); + auto projectile = view.GetProjectile(); + + InverseGrammageType invLambda = 0 / kgMSq; + switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda); + + REQUIRE(view.GetSize() == 2); + } + } +} + +TEST_CASE("SwitchProcess from ProcessSequence") { + DummyProcess<1> innerA; + DummyProcess<2> innerB; + DummyProcess<3> outer; + DummyProcess<4> additional; + + auto seq = innerA << innerB; + + switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV); + auto completeSeq = switchProcess << additional; + + SimpleStack stack; + + SECTION("low energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV}); + auto p = stack.GetNextParticle(); + + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3)); + REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(4. / 7)); + } + + SECTION("SelectInteraction") { + std::vector<int> numberOfSecondaries; + + for (int i = 0; i < 1000; ++i) { + typename SimpleStack::ParticleType theParticle = + stack.GetNextParticle(); // as in corsika::Cascade + StackTestView view(theParticle); + auto projectile = view.GetProjectile(); + + double r = i / 1000.; + InverseGrammageType invLambda = r * 7./4 / kgMSq; + + InverseGrammageType accumulator = 0 / kgMSq; + completeSeq.SelectInteraction(p, projectile, invLambda, accumulator); + + numberOfSecondaries.push_back(view.GetSize()); + } + + auto const mean = + std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) / + numberOfSecondaries.size(); + REQUIRE(mean == Approx(12./7.).margin(0.01)); + } + } + + SECTION("high energy") { + stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV}); + auto p = stack.GetNextParticle(); + + SECTION("interaction length") { + REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3)); + REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(12. / 7.)); + } + + SECTION("SelectInteraction") { + std::vector<int> numberOfSecondaries; + + for (int i = 0; i < 1000; ++i) { + typename SimpleStack::ParticleType theParticle = + stack.GetNextParticle(); // as in corsika::Cascade + StackTestView view(theParticle); + auto projectile = view.GetProjectile(); + + double r = i / 1000.; + InverseGrammageType invLambda = r * 7. / 12. / kgMSq; + + InverseGrammageType accumulator = 0 / kgMSq; + completeSeq.SelectInteraction(p, projectile, invLambda, accumulator); + + numberOfSecondaries.push_back(view.GetSize()); + } + + auto const mean = + std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) / + numberOfSecondaries.size(); + REQUIRE(mean == Approx(24. / 7.).margin(0.01)); + } + } +} diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc index df4456bfb965f398720e33aafd286ed56bbdb307..98f7a7385727d30836e7590bfb6bd54f8de4f542 100644 --- a/Processes/UrQMD/UrQMD.cc +++ b/Processes/UrQMD/UrQMD.cc @@ -2,9 +2,6 @@ #include <corsika/geometry/Vector.h> #include <corsika/particles/ParticleProperties.h> #include <corsika/process/urqmd/UrQMD.h> -#include <corsika/random/RNGManager.h> -#include <corsika/setup/SetupStack.h> -#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> #include <algorithm> @@ -20,7 +17,6 @@ UrQMD::UrQMD() { iniurqmd_(); } using SetupStack = corsika::setup::Stack; using SetupParticle = corsika::setup::Stack::StackIterator; using SetupProjectile = corsika::setup::StackView::StackIterator; -using SetupTrack = corsika::setup::Trajectory; template <typename TParticle> // need template here, as this is called both with // SetupParticle as well as SetupProjectile @@ -82,7 +78,7 @@ bool UrQMD::CanInteract(particles::Code vCode) const { vCode) != std::cend(validProjectileCodes); } -GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const { +GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle) const { if (!CanInteract(vParticle.GetPID())) { // we could do the canInteract check in GetCrossSection, too but if // we do it here we have the advantage of avoiding the loop diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h index 11c8a0eb9a14268f1729a2fbf23508bcf024bcb4..a77e6e843b38055bc96c470c79f2902fcfa58ce5 100644 --- a/Processes/UrQMD/UrQMD.h +++ b/Processes/UrQMD/UrQMD.h @@ -5,7 +5,6 @@ #include <corsika/process/InteractionProcess.h> #include <corsika/random/RNGManager.h> #include <corsika/setup/SetupStack.h> -#include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> #include <array> @@ -17,7 +16,7 @@ namespace corsika::process::UrQMD { UrQMD(); void Init() {} corsika::units::si::GrammageType GetInteractionLength( - corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const; + corsika::setup::Stack::StackIterator&) const; template <typename TParticle> corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,