diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 3236170171fa4af6512bd3a8240089f6d15d272c..f3ebb57d541bf585312f47083d92b745bf82fc36 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -7,7 +7,7 @@ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of * the license. */ - + #include <corsika/cascade/Cascade.h> #include <corsika/process/ProcessSequence.h> #include <corsika/process/stack_inspector/StackInspector.h> @@ -31,6 +31,9 @@ #include <corsika/random/RNGManager.h> +#include <boost/type_index.hpp> +using boost::typeindex::type_id_with_cvr; + #include <iostream> #include <limits> #include <typeinfo> @@ -223,7 +226,10 @@ int main() { corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat"); // assemble all processes into an ordered process list - const auto sequence = /*p0 <<*/ sibyll << decay << cut << trackWriter; + auto sequence = p0 << sibyll << decay << cut << trackWriter; + + // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() + // << "\n"; // setup particle stack, and add primary particle setup::Stack stack; @@ -232,13 +238,14 @@ int main() { { auto particle = stack.NewParticle(); particle.SetPID(Code::Proton); - hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV); + hep::MomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass()); auto plab = stack::super_stupid::MomentumVector(rootCS, 0_GeV, 0_GeV, -P0); particle.SetEnergy(E0); particle.SetMomentum(plab); particle.SetTime(0_ns); Point p(rootCS, 0_m, 0_m, 0_m); particle.SetPosition(p); + cout << particle.GetEnergy() / 1_GeV << endl; } // define air shower object, run simulation diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc index e2836ff1c1b89b5afb67905b4c9814040ef19912..0dfed040c85468501f02ce0cd07ac5a44dddcedb 100644 --- a/Documentation/Examples/staticsequence_example.cc +++ b/Documentation/Examples/staticsequence_example.cc @@ -86,7 +86,7 @@ void modular() { Process3 m3; Process4 m4(0.9); - const auto sequence = m1 << m2 << m3 << m4; + auto sequence = m1 << m2 << m3 << m4; DummyData p; DummyStack s; diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 081325f96038dc8c733f4599603e0ee97a22d04f..e1b742624010bdf549008c82669ff637de54443b 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -114,7 +114,7 @@ TEST_CASE("Cascade", "[Cascade]") { stack_inspector::StackInspector<setup::Stack> p0(true); ProcessSplit p1; - const auto sequence = p0 << p1; + auto sequence = p0 << p1; setup::Stack stack; corsika::cascade::Cascade EAS(env, tracking, sequence, stack); diff --git a/Framework/Geometry/BaseVector.h b/Framework/Geometry/BaseVector.h index c16e83bb87c2d2677e87b2b11b06bd1aa958be22..7e3d8ef6112546940311b808fe3ca798b0988a79 100644 --- a/Framework/Geometry/BaseVector.h +++ b/Framework/Geometry/BaseVector.h @@ -31,10 +31,8 @@ namespace corsika::geometry { BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) : qVector(pQVector) , cs(&pCS) {} - - auto const& GetCoordinateSystem() const { - return *cs; - } + + auto const& GetCoordinateSystem() const { return *cs; } }; } // namespace corsika::geometry diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h index d338119fdc02f5b0f61d0b1f022e4a0d529c5269..bcaa45fba325367048a2a7b93c6ddb0752428300 100644 --- a/Framework/Geometry/Trajectory.h +++ b/Framework/Geometry/Trajectory.h @@ -43,9 +43,9 @@ namespace corsika::geometry { assert(t >= 0 * corsika::units::si::second); return T::ArcLength(0, t); } - + void LimitEndTo(corsika::units::si::LengthType limit) { - fTimeLength = T::TimeFromArclength(limit); + fTimeLength = T::TimeFromArclength(limit); } }; diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h index aa82c0b5d43f3c1c066150ca77d1901b5857918e..03e04c9bdc22f4bf9deba36c6a1d2964e5e60798 100644 --- a/Framework/ProcessSequence/BaseProcess.h +++ b/Framework/ProcessSequence/BaseProcess.h @@ -25,12 +25,21 @@ namespace corsika::process { */ - template <typename derived> + template <typename Derived> struct BaseProcess { - derived& GetRef() { return static_cast<derived&>(*this); } - const derived& GetRef() const { return static_cast<const derived&>(*this); } + private: + BaseProcess() {} + friend Derived; + + public: + Derived& GetRef() { return static_cast<Derived&>(*this); } + const Derived& GetRef() const { return static_cast<const Derived&>(*this); } }; + // overwrite the default trait class, to mark BaseProcess<T> as useful process + template <class T> + std::true_type is_process_impl(const BaseProcess<T>* impl); + } // namespace corsika::process #endif diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h index dd9572ac94cbb8b24352df1c77a22227679ea1c2..19090b04e6867a62c47c6644cac2368e168484ef 100644 --- a/Framework/ProcessSequence/ContinuousProcess.h +++ b/Framework/ProcessSequence/ContinuousProcess.h @@ -41,6 +41,10 @@ namespace corsika::process { corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const; }; + // overwrite the default trait class, to mark BaseProcess<T> as useful process + template <class T> + std::true_type is_process_impl(const ContinuousProcess<T>* impl); + } // namespace corsika::process #endif diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h index 0175cee60f8a4b8f8e0e4fd033cf07d64c239209..4cd83595e4a74e0a3f584e76434cae246767009b 100644 --- a/Framework/ProcessSequence/DecayProcess.h +++ b/Framework/ProcessSequence/DecayProcess.h @@ -35,17 +35,21 @@ namespace corsika::process { /// here starts the interface-definition part // -> enforce derived to implement DoDecay... template <typename Particle, typename Stack> - EProcessReturn DoDecay(Particle&, Stack&) const; + EProcessReturn DoDecay(Particle&, Stack&); template <typename Particle> - corsika::units::si::TimeType GetLifetime(Particle& p) const; + corsika::units::si::TimeType GetLifetime(Particle& p); template <typename Particle> - corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) const { + corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) { return 1. / GetRef().GetLifetime(p); } }; + // overwrite the default trait class, to mark DecayProcess<T> as useful process + template <class T> + std::true_type is_process_impl(const DecayProcess<T>* impl); + } // namespace corsika::process #endif diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h index a56edffa5f3def9dc866abead57142ec26db7890..ff0cabd08a311dc5574c44cedddf8205d8545b19 100644 --- a/Framework/ProcessSequence/InteractionProcess.h +++ b/Framework/ProcessSequence/InteractionProcess.h @@ -36,18 +36,22 @@ namespace corsika::process { /// here starts the interface-definition part // -> enforce derived to implement DoInteraction... template <typename P, typename S> - inline EProcessReturn DoInteraction(P&, S&) const; + inline EProcessReturn DoInteraction(P&, S&); template <typename Particle, typename Track> - corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t) const; + corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t); template <typename Particle, typename Track> corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p, - Track& t) const { + Track& t) { return 1. / GetRef().GetInteractionLength(p, t); } }; + // overwrite the default trait class, to mark BaseProcess<T> as useful process + template <class T> + std::true_type is_process_impl(const InteractionProcess<T>* impl); + } // namespace corsika::process #endif diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index a75ad6adc3d9aeae6159028f0d905683ff8d88d0..3f37626c692b55d77a3a2339d41c6aaada7a5567 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -21,6 +21,7 @@ #include <cmath> #include <limits> +#include <type_traits> namespace corsika::process { @@ -35,20 +36,38 @@ namespace corsika::process { https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern */ + // define a marker (trait class) to tag any class that qualifies as "Process" for the + // "ProcessSequence" + std::false_type is_process_impl(...); + template <class T> + using is_process = decltype(is_process_impl(std::declval<T*>())); + // this is a marker to track which BaseProcess is also a ProcessSequence template <typename T> struct is_process_sequence { static const bool value = false; }; + /** + T1 and T2 are both references if possible (lvalue), otherwise + (rvalue) they are just classes. This allows us to handle both, + rvalue as well as lvalue Processes in the ProcessSequence. + */ template <typename T1, typename T2> class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > { + using T1type = typename std::decay<T1>::type; + using T2type = typename std::decay<T2>::type; + public: - const T1& A; - const T2& B; + T1 A; // this is a reference, if possible + T2 B; // this is a reference, if possible + + // ProcessSequence(ProcessSequence<T1,T2>&& v) : A(v.A), B(v.B) {} + // ProcessSequence<T1,T2>& operator=(ProcessSequence<T1,T2>&& v) { A=v.A; B=v.B; + // return *this; } - ProcessSequence(const T1& in_A, const T2& in_B) + ProcessSequence(T1 in_A, T2 in_B) : A(in_A) , B(in_B) {} @@ -56,13 +75,13 @@ namespace corsika::process { // void Hello() const { detail::CallHello<T1,T2>::Call(A, B); } template <typename Particle, typename Track, typename Stack> - EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const { + EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) { EProcessReturn ret = EProcessReturn::eOk; - if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value || + if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value || is_process_sequence<T1>::value) { ret |= A.DoContinuous(p, t, s); } - if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value || + if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value || is_process_sequence<T2>::value) { ret |= B.DoContinuous(p, t, s); } @@ -70,17 +89,17 @@ namespace corsika::process { } template <typename Particle, typename Track> - corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const { + corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) { corsika::units::si::LengthType 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<T1>, T1>::value || + if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value || is_process_sequence<T1>::value) { corsika::units::si::LengthType const len = A.MaxStepLength(p, track); max_length = std::min(max_length, len); } - if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value || + if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value || is_process_sequence<T2>::value) { corsika::units::si::LengthType const len = B.MaxStepLength(p, track); max_length = std::min(max_length, len); @@ -89,29 +108,28 @@ namespace corsika::process { } template <typename Particle, typename Track> - corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p, - Track& t) const { + corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p, Track& t) { return 1. / GetInverseInteractionLength(p, t); } template <typename Particle, typename Track> - corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength( - Particle& p, Track& t) const { + corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(Particle& p, + Track& t) { return GetInverseInteractionLength(p, t); } template <typename Particle, typename Track> corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p, - Track& t) const { + Track& t) { using namespace corsika::units::si; InverseGrammageType tot = 0 * meter * meter / gram; - if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value || + if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value || is_process_sequence<T1>::value) { tot += A.GetInverseInteractionLength(p, t); } - if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value || + if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value || is_process_sequence<T2>::value) { tot += B.GetInverseInteractionLength(p, t); } @@ -121,15 +139,15 @@ namespace corsika::process { template <typename Particle, typename Stack> EProcessReturn SelectInteraction( Particle& p, Stack& s, corsika::units::si::InverseGrammageType lambda_select, - corsika::units::si::InverseGrammageType& lambda_inv_count) const { + corsika::units::si::InverseGrammageType& lambda_inv_count) { - if constexpr (is_process_sequence<T1>::value) { + if constexpr (is_process_sequence<T1type>::value) { // if A is a process sequence --> check inside const EProcessReturn ret = A.SelectInteraction(p, s, 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<T1>, T1>::value) { + } else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_count += A.GetInverseInteractionLength(p, s); // check if we should execute THIS process and then EXIT @@ -145,7 +163,7 @@ namespace corsika::process { B.SelectInteraction(p, s, 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<T2>, T2>::value) { + } else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_count += B.GetInverseInteractionLength(p, s); // check if we should execute THIS process and then EXIT @@ -158,26 +176,26 @@ namespace corsika::process { } template <typename Particle> - corsika::units::si::TimeType GetTotalLifetime(Particle& p) const { + corsika::units::si::TimeType GetTotalLifetime(Particle& p) { return 1. / GetInverseLifetime(p); } template <typename Particle> - corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) const { + corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) { return GetInverseLifetime(p); } template <typename Particle> - corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) const { + corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) { using namespace corsika::units::si; corsika::units::si::InverseTimeType tot = 0 / second; - if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value || + if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value || is_process_sequence<T1>::value) { tot += A.GetInverseLifetime(p); } - if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value || + if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value || is_process_sequence<T2>::value) { tot += B.GetInverseLifetime(p); } @@ -186,15 +204,15 @@ namespace corsika::process { // select decay process template <typename Particle, typename Stack> - EProcessReturn SelectDecay( - Particle& p, Stack& s, corsika::units::si::InverseTimeType decay_select, - corsika::units::si::InverseTimeType& decay_inv_count) const { + EProcessReturn SelectDecay(Particle& p, Stack& s, + corsika::units::si::InverseTimeType decay_select, + corsika::units::si::InverseTimeType& decay_inv_count) { if constexpr (is_process_sequence<T1>::value) { // if A is a process sequence --> check inside const EProcessReturn ret = A.SelectDecay(p, s, 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<T1>, T1>::value) { + } else if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_count += A.GetInverseLifetime(p); // check if we should execute THIS process and then EXIT @@ -210,7 +228,7 @@ namespace corsika::process { const EProcessReturn ret = B.SelectDecay(p, s, 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<T2>, T2>::value) { + } else if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_count += B.GetInverseLifetime(p); // check if we should execute THIS process and then EXIT @@ -222,10 +240,9 @@ namespace corsika::process { return EProcessReturn::eOk; } - /// TODO the const_cast is not nice, think about the constness here - void Init() const { - const_cast<T1*>(&A)->Init(); - const_cast<T2*>(&B)->Init(); + void Init() { + A.Init(); + B.Init(); } }; @@ -234,28 +251,46 @@ namespace corsika::process { /// must be allowed, this is why we define a macro to define all /// combinations here: -#define OPSEQ(C1, C2) \ - template <typename T1, typename T2> \ - inline const ProcessSequence<T1, T2> operator<<(const C1<T1>& A, const C2<T2>& B) { \ - return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); \ + template < + typename P1, typename P2, + typename std::enable_if<is_process<typename std::decay<P1>::type>::value && + is_process<typename std::decay<P2>::type>::value>::type...> + inline auto operator<<(P1&& A, P2&& B) -> ProcessSequence<P1, P2> { + return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef()); } - OPSEQ(BaseProcess, BaseProcess) - OPSEQ(BaseProcess, InteractionProcess) - OPSEQ(BaseProcess, ContinuousProcess) - OPSEQ(BaseProcess, DecayProcess) - OPSEQ(ContinuousProcess, BaseProcess) - OPSEQ(ContinuousProcess, InteractionProcess) - OPSEQ(ContinuousProcess, ContinuousProcess) - OPSEQ(ContinuousProcess, DecayProcess) - OPSEQ(InteractionProcess, BaseProcess) - OPSEQ(InteractionProcess, InteractionProcess) - OPSEQ(InteractionProcess, ContinuousProcess) - OPSEQ(InteractionProcess, DecayProcess) - OPSEQ(DecayProcess, BaseProcess) - OPSEQ(DecayProcess, InteractionProcess) - OPSEQ(DecayProcess, ContinuousProcess) - OPSEQ(DecayProcess, DecayProcess) + /* #define OPSEQ(C1, C2) \ */ + /* template < \ */ + /* typename P1, typename P2, \ */ + /* typename std::enable_if<is_process<typename std::decay<P1>::type>::value && + * \ */ + /* is_process<typename + * std::decay<P2>::type>::value>::type...> \ */ + /* inline auto operator+(P1&& A, P2&& B)->ProcessSequence<P1, P2> { \ */ + /* return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef()); \ */ + /* } */ + + /* /\*template <typename T1, typename T2> \ */ + /* inline ProcessSequence<T1, T2> operator%(C1<T1>& A, C2<T2>& B) { \ */ + /* return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); \ */ + /* }*\/ */ + + /* OPSEQ(BaseProcess, BaseProcess) */ + /* OPSEQ(BaseProcess, InteractionProcess) */ + /* OPSEQ(BaseProcess, ContinuousProcess) */ + /* OPSEQ(BaseProcess, DecayProcess) */ + /* OPSEQ(ContinuousProcess, BaseProcess) */ + /* OPSEQ(ContinuousProcess, InteractionProcess) */ + /* OPSEQ(ContinuousProcess, ContinuousProcess) */ + /* OPSEQ(ContinuousProcess, DecayProcess) */ + /* OPSEQ(InteractionProcess, BaseProcess) */ + /* OPSEQ(InteractionProcess, InteractionProcess) */ + /* OPSEQ(InteractionProcess, ContinuousProcess) */ + /* OPSEQ(InteractionProcess, DecayProcess) */ + /* OPSEQ(DecayProcess, BaseProcess) */ + /* OPSEQ(DecayProcess, InteractionProcess) */ + /* OPSEQ(DecayProcess, ContinuousProcess) */ + /* OPSEQ(DecayProcess, DecayProcess) */ /// marker to identify objectas ProcessSequence template <typename A, typename B> diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc index 1723a6893cf10eabb0b55bb54dbd7374cc477089..f46b917b537f452efa7897403beb85aa635bfb03 100644 --- a/Framework/ProcessSequence/testProcessSequence.cc +++ b/Framework/ProcessSequence/testProcessSequence.cc @@ -188,7 +188,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Process3 m3(2); Process4 m4(3); - const auto sequence = m1 << m2 << m3 << m4; + auto sequence = m1 << m2 << m3 << m4; globalCount = 0; sequence.Init(); @@ -208,7 +208,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { DummyStack s; DummyTrajectory t; - const auto sequence2 = cp1 << m2 << m3; + auto sequence2 = cp1 << m2 << m3; GrammageType const tot = sequence2.GetTotalInteractionLength(s, t); InverseGrammageType const tot_inv = sequence2.GetTotalInverseInteractionLength(s, t); cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl; @@ -222,7 +222,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { DummyStack s; - const auto sequence2 = cp1 << m2 << m3 << d3; + auto sequence2 = cp1 << m2 << m3 << d3; TimeType const tot = sequence2.GetTotalLifetime(s); InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(s); cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl; @@ -235,7 +235,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") { Process2 m2(1); Process3 m3(2); - const auto sequence2 = cp1 << m2 << m3 << cp2; + auto sequence2 = cp1 << m2 << m3 << cp2; DummyData p; DummyStack s; diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc index 814d63e1f2f44463e2861f08d4b6ae5d3730af8e..faac8c6d112e26448502910cf516eeaee949a674 100644 --- a/Framework/Utilities/testCOMBoost.cc +++ b/Framework/Utilities/testCOMBoost.cc @@ -55,7 +55,7 @@ TEST_CASE("boosts") { // boost projecticle auto const [eProjectileCoM, pProjectileCoM] = boost.toCoM(eProjectileLab, pProjectileLab); - + // boost target auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab); diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h index 93516c141ad13888e47b57d72d247c897c2ebdbf..03763ddd0a8263af71f4049ffd964d7aa30faf55 100644 --- a/Processes/Sibyll/Decay.h +++ b/Processes/Sibyll/Decay.h @@ -25,11 +25,11 @@ namespace corsika::process { setTrackedParticlesStable(); } - void setTrackedParticlesStable() const { + void setTrackedParticlesStable() { /* - Sibyll is hadronic generator - only hadrons decay - */ + Sibyll is hadronic generator + only hadrons decay + */ // set particles unstable setHadronsUnstable(); // make tracked particles stable @@ -45,17 +45,17 @@ namespace corsika::process { } } - void setUnstable(const corsika::particles::Code pCode) const { + void setUnstable(const corsika::particles::Code pCode) { int s_id = process::sibyll::ConvertToSibyllRaw(pCode); s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]); } - void setStable(const corsika::particles::Code pCode) const { + void setStable(const corsika::particles::Code pCode) { int s_id = process::sibyll::ConvertToSibyllRaw(pCode); s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]); } - void setAllStable() const { + void setAllStable() { // name? also makes EM particles stable using std::cout; @@ -75,7 +75,7 @@ namespace corsika::process { } } - void setHadronsUnstable() const { + void setHadronsUnstable() { using std::cout; using std::endl; @@ -118,7 +118,7 @@ namespace corsika::process { } template <typename Particle> - corsika::units::si::TimeType GetLifetime(Particle& p) const { + corsika::units::si::TimeType GetLifetime(Particle& p) { using std::cout; using std::endl; using namespace corsika::units::si; @@ -142,7 +142,7 @@ namespace corsika::process { } template <typename Particle, typename Stack> - void DoDecay(Particle& p, Stack& s) const { + void DoDecay(Particle& p, Stack& s) { using corsika::geometry::Point; using namespace corsika::units::si; fCount++; diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h index fe8cfecef38eeef1ddec703e7c400e505b1f70f5..48222bfe364d0cb92f5bfb1d76087fc4d13b5400 100644 --- a/Processes/Sibyll/Interaction.h +++ b/Processes/Sibyll/Interaction.h @@ -35,7 +35,7 @@ namespace corsika::process::sibyll { } template <typename Particle, typename Track> - corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const { + corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) { using namespace corsika::units; using namespace corsika::units::hep; @@ -115,7 +115,7 @@ namespace corsika::process::sibyll { } template <typename Particle, typename Stack> - corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const { + corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) { using namespace corsika::units; using namespace corsika::units::hep; @@ -137,7 +137,7 @@ namespace corsika::process::sibyll { // FOR NOW: hard coded z-axis for corsika frame QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m}; QuantityVector<length_d> const xAxis{1_m, 0_m, 0_m}; - auto pt = [](MomentumVector p) { + [[maybe_unused]] auto pt = [](MomentumVector p) { return sqrt(p.GetComponents()[0] * p.GetComponents()[0] + p.GetComponents()[1] * p.GetComponents()[1]); }; diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h index d78273822c700edaaae9e2c6bbefdfbfabb2a9e8..0be72a7692f9871d2a54ea7107f1eea560f609c7 100644 --- a/Processes/Sibyll/SibStack.h +++ b/Processes/Sibyll/SibStack.h @@ -80,9 +80,8 @@ namespace corsika::process::sibyll { corsika::units::hep::EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); } - bool HasDecayed() const - { - return abs(GetStackData().GetId(GetIndex()))>100 ? true : false; + bool HasDecayed() const { + return abs(GetStackData().GetId(GetIndex())) > 100 ? true : false; } void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); } corsika::process::sibyll::SibyllCode GetPID() const { diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc index abb29da21f41637e268fb62d6874f5dbfd1c55b3..39b59505eaba16853a06807d42ed91a56dcb9750 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -37,7 +37,7 @@ StackInspector<Stack>::~StackInspector() {} template <typename Stack> process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&, - Stack& s) const { + Stack& s) { if (!fReport) return process::EProcessReturn::eOk; [[maybe_unused]] int i = 0; @@ -60,8 +60,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr } template <typename Stack> -corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength( - Particle&, setup::Trajectory&) const { +corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength(Particle&, + setup::Trajectory&) { return std::numeric_limits<double>::infinity() * meter; } diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h index 07b987b46f8a03a2683b2d160ac86e24ac251560..6a14b4e6eb9578e8c283637c3afbaee71f97256b 100644 --- a/Processes/StackInspector/StackInspector.h +++ b/Processes/StackInspector/StackInspector.h @@ -31,14 +31,14 @@ namespace corsika::process { ~StackInspector(); void Init(); - EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const; + EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s); //~ template <typename Particle> corsika::units::si::LengthType MaxStepLength(Particle&, - corsika::setup::Trajectory&) const; + corsika::setup::Trajectory&); private: bool fReport; - mutable int fCountStep = 0; + int fCountStep = 0; }; } // namespace stack_inspector diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h index 6b554feba48fb65aadd8dcc31f09ddfea5d7288c..31966195751d0a1dfaee6d61a60421a0f80d4676 100644 --- a/Processes/TrackingLine/TrackingLine.h +++ b/Processes/TrackingLine/TrackingLine.h @@ -84,9 +84,13 @@ namespace corsika::process { p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c; auto const currentPosition = p.GetPosition(); - std::cout << "TrackingLine pid: " << p.GetPID() << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl; + std::cout << "TrackingLine pid: " << p.GetPID() + << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl; std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates() << std::endl; + std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl; + std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV + << " GeV " << std::endl; std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl; // to do: include effect of magnetic field