diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index 5ad445d0832d4cac93da9f25831c2810a8a3b241..74f731c7d324e3f81165e5d2eaf5035955833767 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -46,21 +46,6 @@ namespace corsika::process { template <class T> constexpr bool is_process_v = std::is_base_of_v<BaseProcess, T>; - namespace switch_process { - template <typename A, typename B> - class SwitchProcess; // fwd-decl. - } - - // to detect SwitchProcesses inside the ProcessSequence - template <typename T> - struct is_switch_process : std::false_type {}; - - template <typename T> - bool constexpr is_switch_process_v = is_switch_process<T>::value; - - template <typename A, typename B> - struct is_switch_process<switch_process::SwitchProcess<A, B>> : std::true_type {}; - template <class... Ts> class ProcessSequence : public std::tuple<Ts...> { using storage_t = std::tuple<Ts...>; @@ -77,10 +62,11 @@ namespace corsika::process { // void Hello() const { detail::CallHello<T1,T2>::Call(A, B); } template <template <class> class Base, class Functor> - void apply(Functor&& functor) { + void for_each(Functor&& functor) { boost::mp11::tuple_for_each(static_cast<storage_t&>(*this), [&](auto&& x) { using T = std::decay_t<decltype(x)>; - if constexpr (std::is_base_of_v<Base<T>, T>) functor(x); + if constexpr (std::is_base_of_v<Base<T>, T>) + std::forward<Functor>(functor)(std::forward<decltype(x)>(x)); }); } @@ -88,7 +74,7 @@ namespace corsika::process { EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from, VTNType const& to) { EProcessReturn ret = EProcessReturn::eOk; - apply<BoundaryCrossingProcess>( + for_each<BoundaryCrossingProcess>( [&](auto&& proc) { ret |= proc.DoBoundaryCrossing(p, from, to); }); return ret; } @@ -96,14 +82,14 @@ namespace corsika::process { template <typename TParticle, typename TTrack> EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) { EProcessReturn ret = EProcessReturn::eOk; - apply<ContinuousProcess>([&](auto&& proc) { ret |= proc.DoContinuous(vP, vT); }); + for_each<ContinuousProcess>([&](auto&& proc) { ret |= proc.DoContinuous(vP, vT); }); return ret; } template <typename TSecondaries> EProcessReturn DoSecondaries(TSecondaries& vS) { EProcessReturn ret = EProcessReturn::eOk; - apply<SecondariesProcess>([&](auto&& proc) { ret |= proc.DoSecondaries(vS); }); + for_each<SecondariesProcess>([&](auto&& proc) { ret |= proc.DoSecondaries(vS); }); return ret; } @@ -117,7 +103,7 @@ namespace corsika::process { */ bool CheckStep() { bool ret = false; - apply<StackProcess>([&](auto&& proc) { ret |= proc.CheckStep(); }); + for_each<StackProcess>([&](auto&& proc) { ret |= proc.CheckStep(); }); return ret; } @@ -127,7 +113,7 @@ namespace corsika::process { template <typename TStack> EProcessReturn DoStack(TStack& vS) { EProcessReturn ret = EProcessReturn::eOk; - apply<StackProcess>([&](auto&& proc) { + for_each<StackProcess>([&](auto&& proc) { if (proc.CheckStep()) ret |= proc.DoStack(vS); }); return ret; @@ -139,8 +125,7 @@ namespace corsika::process { // infinite if no other process in the sequence implements it LengthType max_length = std::numeric_limits<double>::infinity() * meter; - - apply<ContinuousProcess>([&](auto&& proc) { + for_each<ContinuousProcess>([&](auto&& proc) { const auto len = proc.MaxStepLength(vP, vTrack); if (len < max_length) max_length = len; }); @@ -163,52 +148,11 @@ namespace corsika::process { using namespace corsika::units::si; InverseGrammageType tot = 0 * meter * meter / gram; - - apply<InteractionProcess>( + for_each<InteractionProcess>( [&](auto&& proc) { tot += proc.GetInverseInteractionLength(vP); }); return tot; } - template <typename TParticle, typename TSecondaries> - EProcessReturn SelectInteraction( - TParticle& vP, TSecondaries& vS, - corsika::units::si::InverseGrammageType lambda_select, - corsika::units::si::InverseGrammageType& lambda_inv_count) { - boost::ignore_unused(vP, vS, lambda_select, lambda_inv_count); - // if constexpr (t1ProcSeq || t1SwitchProc) { - // // if A is a process sequence --> check inside - // const EProcessReturn ret = - // 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_v<InteractionProcess<T1type>, T1type>) { - // // if this is not a ContinuousProcess --> evaluate probability - // 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); - // return EProcessReturn::eInteracted; - // } - // } // end branch A - // - // if constexpr (t2ProcSeq || t2SwitchProc) { - // // if A is a process sequence --> check inside - // const EProcessReturn ret = - // 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_v<InteractionProcess<T2type>, T2type>) { - // // if this is not a ContinuousProcess --> evaluate probability - // 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); - // return EProcessReturn::eInteracted; - // } - // } // end branch A - return EProcessReturn::eOk; - } - template <typename TParticle> corsika::units::si::TimeType GetTotalLifetime(TParticle& p) { return 1. / GetInverseLifetime(p); @@ -223,21 +167,42 @@ namespace corsika::process { corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) { using namespace corsika::units::si; InverseTimeType tot = 0 / second; - apply<DecayProcess>([&](auto&& proc) { tot += proc.GetInverseLifetime(p); }); + for_each<DecayProcess>([&](auto&& proc) { tot += proc.GetInverseLifetime(p); }); return tot; } + // select interaction process + template <typename TParticle, typename TSecondaries> + EProcessReturn SelectInteraction( + TParticle& vP, TSecondaries& vS, + corsika::units::si::InverseGrammageType lambda_select, + corsika::units::si::InverseGrammageType& lambda_inv_count) { + EProcessReturn ret = EProcessReturn::eOk; + for_each<InteractionProcess>([&](auto&& proc) { + // check if we should execute THIS process and then EXIT + if (ret == EProcessReturn::eOk) { + // if this is not a ContinuousProcess --> evaluate probability + lambda_inv_count += proc.GetInverseInteractionLength(vP); + if (lambda_select < lambda_inv_count) { + proc.DoInteraction(vS); + ret = EProcessReturn::eInteracted; + } + } + }); + return ret; + } + // select decay process template <typename TParticle, typename TSecondaries> EProcessReturn SelectDecay(TParticle& vP, TSecondaries& vS, corsika::units::si::InverseTimeType decay_select, corsika::units::si::InverseTimeType& decay_inv_count) { EProcessReturn ret = EProcessReturn::eOk; - apply<DecayProcess>([&](auto&& proc) { + for_each<DecayProcess>([&](auto&& proc) { + // check if we should execute THIS process and then skip all others if (ret == EProcessReturn::eOk) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_count += proc.GetInverseLifetime(vP); - // check if we should execute THIS process and then EXIT if (decay_select < decay_inv_count) { proc.DoDecay(vS); ret = EProcessReturn::eDecayed; @@ -262,17 +227,6 @@ namespace corsika::process { template <class... Ts> ProcessSequence(Ts&&... ts)->ProcessSequence<Ts...>; - /// marker to identify objectas ProcessSequence - template <typename T> - struct is_process_sequence : std::false_type {}; - - template <class... Ts> - struct is_process_sequence<corsika::process::ProcessSequence<Ts...>> : std::true_type { - }; - - template <typename T> - bool constexpr is_process_sequence_v = is_process_sequence<T>::value; - } // namespace corsika::process #endif diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc index ce1e2932566da8823616f53b764ab94fd813232c..0b8c95e4a6ce6ba1cbb7411ca8350c28d01ff0f5 100644 --- a/Framework/ProcessSequence/testProcessSequence.cc +++ b/Framework/ProcessSequence/testProcessSequence.cc @@ -270,7 +270,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Stack1 s1(1); Stack1 s2(2); - ProcessSequence sequence(Stack1(1), Stack1(2)); + ProcessSequence sequence(s1, s2); DummyStack stack; @@ -281,14 +281,3 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { CHECK(s2.GetCount() == 10); } } - -/* - Note: there is a fine-grained dedicated test-suite for SwitchProcess - in Processes/SwitchProcess/testSwtichProcess - */ -TEST_CASE("SwitchProcess") { - Process1 p1(0); - Process2 p2(0); - switch_process::SwitchProcess s(p1, p2, 10_GeV); - REQUIRE(is_switch_process_v<decltype(s)>); -} diff --git a/Processes/SwitchProcess/SwitchProcess.h b/Processes/SwitchProcess/SwitchProcess.h index d5712f0209185fadc4b734ee452e20aa1e905245..8209eeed7eb40d534588c1b0167336de84d0ee7f 100644 --- a/Processes/SwitchProcess/SwitchProcess.h +++ b/Processes/SwitchProcess/SwitchProcess.h @@ -20,28 +20,40 @@ 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 distinct case when assembled into a - * (greater) ProcessSequence. + * P1. For energies below the threshold, P1 is invoked, otherwise P2. */ template <class TLowEProcess, class THighEProcess> - class SwitchProcess : public BaseProcess { - TLowEProcess& fLowEProcess; - THighEProcess& fHighEProcess; - units::si::HEPEnergyType const fThresholdEnergy; + class SwitchProcess + : public InteractionProcess<SwitchProcess<TLowEProcess, THighEProcess>> { + TLowEProcess low_process_; + THighEProcess high_process_; + units::si::HEPEnergyType const threshold_; + + template <typename TParticle, typename TSecondaries, typename Process> + static EProcessReturn impl(TParticle& vP, TSecondaries& vS, + corsika::units::si::InverseGrammageType lambda_select, + corsika::units::si::InverseGrammageType& lambda_inv_count, + Process& process) { + lambda_inv_count += process.GetInverseInteractionLength(vP); + // check if we should execute THIS process and then EXIT + if (lambda_select < lambda_inv_count) { + process.DoInteraction(vS); + return EProcessReturn::eInteracted; + } + return EProcessReturn::eOk; + }; public: - SwitchProcess(TLowEProcess& vLowEProcess, THighEProcess& vHighEProcess, + SwitchProcess(TLowEProcess vLowEProcess, THighEProcess vHighEProcess, units::si::HEPEnergyType vThresholdEnergy) - : fLowEProcess(vLowEProcess) - , fHighEProcess(vHighEProcess) - , fThresholdEnergy(vThresholdEnergy) {} + : low_process_(vLowEProcess) + , high_process_(vHighEProcess) + , threshold_(vThresholdEnergy) {} void Init() { - fLowEProcess.Init(); - fHighEProcess.Init(); + low_process_.Init(); + high_process_.Init(); } template <typename TParticle> @@ -51,19 +63,9 @@ namespace corsika::process::switch_process { 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); - } - } + return vParticle.GetEnergy() < threshold_ + ? low_process_.GetInteractionLength(vParticle) + : high_process_.GetInteractionLength(vParticle); } // required to conform to ProcessSequence interface. We cannot just @@ -72,37 +74,14 @@ namespace corsika::process::switch_process { template <typename TParticle, typename TSecondaries> EProcessReturn SelectInteraction( TParticle& vP, TSecondaries& vS, - [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select, + 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; - } - } - } + return vP.GetEnergy() < threshold_ + ? impl(vP, vS, lambda_select, lambda_inv_count, low_process_) + : impl(vP, vS, lambda_select, lambda_inv_count, high_process_); } }; + } // namespace corsika::process::switch_process #endif diff --git a/Processes/SwitchProcess/testSwitchProcess.cc b/Processes/SwitchProcess/testSwitchProcess.cc index 53257a4e7fc7cedb888600688c54d5a0776e971c..9e6b395e0296330ce18a4a767ae8198b46e07ce2 100644 --- a/Processes/SwitchProcess/testSwitchProcess.cc +++ b/Processes/SwitchProcess/testSwitchProcess.cc @@ -180,85 +180,3 @@ TEST_CASE("SwitchProcess from InteractionProcess") { } } } - -TEST_CASE("SwitchProcess from ProcessSequence") { - DummyProcess<1> innerA; - DummyProcess<2> innerB; - DummyProcess<3> outer; - DummyProcess<4> additional; - - ProcessSequence seq(innerA, innerB); - - switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV); - ProcessSequence 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)); - } - } -}