diff --git a/Environment/IMagneticFieldModel.h b/Environment/IMagneticFieldModel.h index a8c13b4e129d24916d75bd1f5f4f66121fc20be5..d27eda38ca03d1c0c6a35c6318e7dda7bb32da46 100644 --- a/Environment/IMagneticFieldModel.h +++ b/Environment/IMagneticFieldModel.h @@ -40,7 +40,7 @@ namespace corsika::environment { /** * A virtual default destructor. */ - virtual ~IMagneticFieldModel() = default; + virtual ~IMagneticFieldModel() = default; // LCOV_EXCL_LINE }; // END: class MagneticField diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h index c26650a055f9075f3d034820b4081f2dca11f6ae..da6c9ca47a41ad14d58b1cc991a473ee970520bd 100644 --- a/Environment/IMediumModel.h +++ b/Environment/IMediumModel.h @@ -18,7 +18,7 @@ namespace corsika::environment { class IMediumModel { public: - virtual ~IMediumModel() = default; + virtual ~IMediumModel() = default; // LCOV_EXCL_LINE virtual corsika::units::si::MassDensityType GetMassDensity( corsika::geometry::Point const&) const = 0; diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt index bb8cbb326152902709a34337fd18124f6fdd257d..1d1e7e60da4e25e6558683e5ca0183e09e183fe4 100644 --- a/Framework/Cascade/CMakeLists.txt +++ b/Framework/Cascade/CMakeLists.txt @@ -15,6 +15,19 @@ add_library (CORSIKAcascade INTERFACE) CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS}) +target_link_libraries( + CORSIKAcascade + INTERFACE + CORSIKArandom + CORSIKAstackinterface + CORSIKAparticles + CORSIKAgeometry + CORSIKAenvironment + CORSIKAprocesssequence + CORSIKAunits + CORSIKAlogging + ) + # include directive for upstream code target_include_directories ( CORSIKAcascade @@ -34,19 +47,9 @@ install ( CORSIKA_ADD_TEST(testCascade) target_link_libraries ( testCascade - # CORSIKAutls - CORSIKArandom - ProcessSibyll CORSIKAcascade ProcessStackInspector ProcessTrackingLine ProcessNullModel - CORSIKAstackinterface - CORSIKAprocesses - CORSIKAparticles - CORSIKAgeometry - CORSIKAenvironment - CORSIKAprocesssequence - CORSIKAunits CORSIKAtesting ) diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 389ec0557e8bcc6da87fd57adafcf608d48a2a8e..e7f0cec36185615f18d70ddf2b94113e5968b10a 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -9,6 +9,7 @@ #pragma once #include <corsika/environment/Environment.h> +#include <corsika/logging/Logging.h> #include <corsika/process/ProcessReturn.h> #include <corsika/random/ExponentialDistribution.h> #include <corsika/random/RNGManager.h> @@ -175,7 +176,7 @@ namespace corsika::cascade { currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step, next_interact); - // determine the maximum geometric step length + // determine the maximum geometric step length from continuous processes LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step); std::cout << "distance_max=" << distance_max << std::endl; @@ -250,7 +251,9 @@ namespace corsika::cascade { // make sure particle actually did decay if it should have done so if (secondaries.GetSize() == 1 && projectile.GetPID() == secondaries.GetNextParticle().GetPID()) - throw std::runtime_error("Cascade::Step: Particle decays into itself!"); + throw std::runtime_error( + fmt::format("Cascade: {} decayed into itself!", + particles::GetName(projectile.GetPID()))); } fProcessSequence.DoSecondaries(secondaries); @@ -263,20 +266,31 @@ namespace corsika::cascade { } else { // step-length limitation within volume std::cout << "step-length limitation" << std::endl; - fProcessSequence.DoSecondaries(secondaries); + // no extra physics happens here. just proceed to next step. } [[maybe_unused]] auto const assertion = [&] { auto const* numericalNodeAfterStep = fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition()); + C8LOG_TRACE(fmt::format( + "Geometry check: numericalNodeAfterStep={} currentLogicalNode={}", + fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode))); return numericalNodeAfterStep == currentLogicalNode; }; assert(assertion()); // numerical and logical nodes don't match - } else { // boundary crossing, step is limited by volume boundary + + } else { // boundary crossing, step is limited by volume boundary + std::cout << "boundary crossing! next node = " << nextVol << std::endl; vParticle.SetNode(nextVol); - // DoBoundary may delete the particle (or not) + /* + DoBoundary may delete the particle (or not) + + small caveat: any changes to vParticle, or even the production + of new secondaries is currently not passed to ParticleCut, + thus, particles outside the desired phase space may be produced. + */ fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol); } } diff --git a/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h index eeb4da07633674ed049b33778839bef486633316..8249d1472efd1d80187f7d3ae8b0f82935750881 100644 --- a/Framework/ProcessSequence/BaseProcess.h +++ b/Framework/ProcessSequence/BaseProcess.h @@ -13,6 +13,8 @@ namespace corsika::process { + class TDerived; // fwd decl + /** \class BaseProcess @@ -22,19 +24,17 @@ namespace corsika::process { */ - template <typename Derived> - struct BaseProcess { - private: - BaseProcess() {} - friend Derived; + template <typename TDerived> + class BaseProcess { + protected: + friend TDerived; - public: - Derived& GetRef() { return static_cast<Derived&>(*this); } - const Derived& GetRef() const { return static_cast<const Derived&>(*this); } - }; + BaseProcess() = default; // protected constructor will allow only + // derived classes to be created, not + // BaseProcess itself - // 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); + TDerived& GetRef() { return static_cast<TDerived&>(*this); } + const TDerived& GetRef() const { return static_cast<const TDerived&>(*this); } + }; } // namespace corsika::process diff --git a/Framework/ProcessSequence/BoundaryCrossingProcess.h b/Framework/ProcessSequence/BoundaryCrossingProcess.h index c0793e2b6b126ed5c2cac1b6765ce01221143312..24c2d93a8a5df5a8174e9e38cea5bca69a04f201 100644 --- a/Framework/ProcessSequence/BoundaryCrossingProcess.h +++ b/Framework/ProcessSequence/BoundaryCrossingProcess.h @@ -8,23 +8,21 @@ #pragma once -#include <corsika/environment/Environment.h> +#include <corsika/process/BaseProcess.h> #include <corsika/process/ProcessReturn.h> namespace corsika::process { + template <typename TDerived> - struct BoundaryCrossingProcess { - auto& GetRef() { return static_cast<TDerived&>(*this); } - auto const& GetRef() const { return static_cast<const TDerived&>(*this); } + struct BoundaryCrossingProcess : public BaseProcess<TDerived> { /** * This method is called when a particle crosses the boundary between the nodes * \p from and \p to. */ - template <typename Particle, typename VTNType> - EProcessReturn DoBoundaryCrossing(Particle&, VTNType const& from, VTNType const& to); + template <typename TParticle, typename TVolumeNode> + EProcessReturn DoBoundaryCrossing(TParticle&, TVolumeNode const& from, + TVolumeNode const& to); }; - template <class T> - std::true_type is_process_impl(BoundaryCrossingProcess<T> const* impl); } // namespace corsika::process diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h index b8d2f52211da1c6f6d56e04a398258278b40c758..9171eeeebaf6cb0a700bd45b2e8c5775238f00d6 100644 --- a/Framework/ProcessSequence/ContinuousProcess.h +++ b/Framework/ProcessSequence/ContinuousProcess.h @@ -8,6 +8,7 @@ #pragma once +#include <corsika/process/BaseProcess.h> #include <corsika/process/ProcessReturn.h> // for convenience #include <corsika/units/PhysicalUnits.h> @@ -22,23 +23,17 @@ namespace corsika::process { */ - template <typename derived> - struct ContinuousProcess { - derived& GetRef() { return static_cast<derived&>(*this); } - const derived& GetRef() const { return static_cast<const derived&>(*this); } + template <typename TDerived> + struct ContinuousProcess : public BaseProcess<TDerived> { // here starts the interface part - // -> enforce derived to implement DoContinuous... - template <typename Particle, typename Track> - EProcessReturn DoContinuous(Particle&, Track const&) const; + // -> enforce TDerived to implement DoContinuous... + template <typename TParticle, typename TTrack> + EProcessReturn DoContinuous(TParticle&, TTrack const&) const; - // -> enforce derived to implement MaxStepLength... - template <typename Particle, typename Track> - units::si::LengthType MaxStepLength(Particle const& p, Track const& track) const; + // -> enforce TDerived to implement MaxStepLength... + template <typename TParticle, typename TTrack> + units::si::LengthType MaxStepLength(TParticle const& p, TTrack const& 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 diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h index fb5ae9eb9ca881bc41c076e9f852a865059ff269..779abbc94940c4804872515556f4416cf263e190 100644 --- a/Framework/ProcessSequence/DecayProcess.h +++ b/Framework/ProcessSequence/DecayProcess.h @@ -8,6 +8,7 @@ #pragma once +#include <corsika/process/BaseProcess.h> #include <corsika/process/ProcessReturn.h> // for convenience #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> @@ -23,27 +24,23 @@ namespace corsika::process { */ - template <typename derived> - struct DecayProcess { - derived& GetRef() { return static_cast<derived&>(*this); } - const derived& GetRef() const { return static_cast<const derived&>(*this); } + template <typename TDerived> + struct DecayProcess : BaseProcess<TDerived> { + + using BaseProcess<TDerived>::GetRef; /// here starts the interface-definition part - // -> enforce derived to implement DoDecay... - template <typename Particle> - EProcessReturn DoDecay(Particle&); + // -> enforce TDerived to implement DoDecay... + template <typename TParticle> + EProcessReturn DoDecay(TParticle&); - template <typename Particle> - corsika::units::si::TimeType GetLifetime(Particle& p); + template <typename TParticle> + corsika::units::si::TimeType GetLifetime(TParticle& p); - template <typename Particle> - corsika::units::si::InverseTimeType GetInverseLifetime(Particle& vP) { + template <typename TParticle> + corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& vP) { return 1. / GetRef().GetLifetime(vP); } }; - // 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 diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h index ec3a086297caf09cdfddb29ea49acbb1346f36e7..204d6db476f77a96841a9acb77622b1e3f74658b 100644 --- a/Framework/ProcessSequence/InteractionProcess.h +++ b/Framework/ProcessSequence/InteractionProcess.h @@ -8,6 +8,7 @@ #pragma once +#include <corsika/process/BaseProcess.h> #include <corsika/process/ProcessReturn.h> // for convenience #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> @@ -23,16 +24,15 @@ namespace corsika::process { */ - template <typename derived> - struct InteractionProcess { + template <typename TDerived> + struct InteractionProcess : public BaseProcess<TDerived> { - derived& GetRef() { return static_cast<derived&>(*this); } - const derived& GetRef() const { return static_cast<const derived&>(*this); } + using BaseProcess<TDerived>::GetRef; /// here starts the interface-definition part - // -> enforce derived to implement DoInteraction... - template <typename Particle> - EProcessReturn DoInteraction(Particle&); + // -> enforce TDerived to implement DoInteraction... + template <typename TParticle> + EProcessReturn DoInteraction(TParticle&); template <typename TParticle> corsika::units::si::GrammageType GetInteractionLength(TParticle& p); @@ -43,8 +43,4 @@ namespace corsika::process { } }; - // 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 diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h index bef2c12d70858a7626ec4384468f5e10feed0ca3..ae1e45b951d7a2ff144c6d683bc7b88ff93eb831 100644 --- a/Framework/ProcessSequence/ProcessReturn.h +++ b/Framework/ProcessSequence/ProcessReturn.h @@ -27,12 +27,20 @@ namespace corsika::process { return static_cast<EProcessReturn>(static_cast<int>(a) | static_cast<int>(b)); } - inline EProcessReturn& operator|=(EProcessReturn& a, EProcessReturn b) { + inline EProcessReturn& operator|=(EProcessReturn& a, const EProcessReturn b) { return a = a | b; } - inline bool operator==(EProcessReturn a, EProcessReturn b) { + inline EProcessReturn operator&(const EProcessReturn a, const EProcessReturn b) { + return static_cast<EProcessReturn>(static_cast<int>(a) & static_cast<int>(b)); + } + + inline bool operator==(const EProcessReturn a, const EProcessReturn b) { return (static_cast<int>(a) & static_cast<int>(b)) != 0; } + inline bool isAbsorbed(const EProcessReturn a) { + return static_cast<int>(a & EProcessReturn::eParticleAbsorbed); + } + } // namespace corsika::process diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index c54ea09f5f85b244f7021fffb5fe19e5c6619bd2..7555f1800a720d81855edecbd4e3f83728cd444a 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -35,33 +35,29 @@ 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> + template <typename TClass> struct is_process_sequence : std::false_type {}; - template <typename T> - bool constexpr is_process_sequence_v = is_process_sequence<T>::value; + template <typename TClass> + bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value; + // we also need a marker to identify SwitchProcess namespace switch_process { - template <typename A, typename B> + template <typename TProcess1, typename TProcess2> class SwitchProcess; // fwd-decl. } // to detect SwitchProcesses inside the ProcessSequence - template <typename T> + template <typename TClass> struct is_switch_process : std::false_type {}; - template <typename T> - bool constexpr is_switch_process_v = is_switch_process<T>::value; + template <typename TClass> + bool constexpr is_switch_process_v = is_switch_process<TClass>::value; - template <typename A, typename B> - struct is_switch_process<switch_process::SwitchProcess<A, B>> : std::true_type {}; + template <typename Process1, typename Process2> + struct is_switch_process<switch_process::SwitchProcess<Process1, Process2>> + : std::true_type {}; /** T1 and T2 are both references if possible (lvalue), otherwise @@ -88,9 +84,6 @@ namespace corsika::process { : A(in_A) , B(in_B) {} - // example for a trait-based call: - // void Hello() const { detail::CallHello<T1,T2>::Call(A, B); } - template <typename Particle, typename VTNType> EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from, VTNType const& to) { @@ -113,10 +106,10 @@ namespace corsika::process { EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) { EProcessReturn ret = EProcessReturn::eOk; if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) { - ret |= A.DoContinuous(vP, vT); + if (!isAbsorbed(ret)) { ret |= A.DoContinuous(vP, vT); } } if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) { - ret |= B.DoContinuous(vP, vT); + if (!isAbsorbed(ret)) { ret |= B.DoContinuous(vP, vT); } } return ret; } @@ -317,20 +310,28 @@ namespace corsika::process { } }; - /// the << operator assembles many BaseProcess, ContinuousProcess, and - /// Interaction/DecayProcess objects into a ProcessSequence, all combinatorics - /// must be allowed, this is why we define a macro to define all - /// combinations here: - - 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&& vA, P2&& vB) -> ProcessSequence<P1, P2> { - return ProcessSequence<P1, P2>(vA.GetRef(), vB.GetRef()); + // the << operator assembles many BaseProcess, ContinuousProcess, and + // Interaction/DecayProcess objects into a ProcessSequence, all combinatorics + // must be allowed, this is why we define a macro to define all + // combinations here: + + // enable the << operator to construct ProcessSequence from two + // Processes, only if poth Processes derive from BaseProcesses + + template <typename TProcess1, typename TProcess2> + inline typename std::enable_if< + std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>, + typename std::decay<TProcess1>::type>::value && + std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>, + typename std::decay<TProcess2>::type>::value, + ProcessSequence<TProcess1, TProcess2>>::type + operator<<(TProcess1&& vA, TProcess2&& vB) { + return ProcessSequence<TProcess1, TProcess2>(vA, vB); } /// marker to identify objectas ProcessSequence - template <typename A, typename B> - struct is_process_sequence<corsika::process::ProcessSequence<A, B>> : std::true_type {}; + template <typename TProcess1, typename TProcess2> + struct is_process_sequence<corsika::process::ProcessSequence<TProcess1, TProcess2>> + : std::true_type {}; + } // namespace corsika::process diff --git a/Framework/ProcessSequence/SecondariesProcess.h b/Framework/ProcessSequence/SecondariesProcess.h index 3a8f35c74fc6718dd03d81d8edcf840be1d9b316..62133059fba81acb084da8dbb8bd7997121e0a6b 100644 --- a/Framework/ProcessSequence/SecondariesProcess.h +++ b/Framework/ProcessSequence/SecondariesProcess.h @@ -8,6 +8,7 @@ #pragma once +#include <corsika/process/BaseProcess.h> #include <corsika/process/ProcessReturn.h> // for convenience #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> @@ -23,20 +24,13 @@ namespace corsika::process { */ - template <typename derived> - struct SecondariesProcess { - - derived& GetRef() { return static_cast<derived&>(*this); } - const derived& GetRef() const { return static_cast<const derived&>(*this); } + template <typename TDerived> + struct SecondariesProcess : public BaseProcess<TDerived> { /// here starts the interface-definition part - // -> enforce derived to implement DoSecondaries... + // -> enforce TDerived to implement DoSecondaries... template <typename TSecondaries> inline EProcessReturn DoSecondaries(TSecondaries&); }; - // overwrite the default trait class, to mark BaseProcess<T> as useful process - template <class T> - std::true_type is_process_impl(const SecondariesProcess<T>* impl); - } // namespace corsika::process diff --git a/Framework/ProcessSequence/StackProcess.h b/Framework/ProcessSequence/StackProcess.h index 05d3583963f8e0d882290c4e84924728023a3ee4..d54c7944d12b03ee09160c3aec0665bbcc54d469 100644 --- a/Framework/ProcessSequence/StackProcess.h +++ b/Framework/ProcessSequence/StackProcess.h @@ -8,6 +8,7 @@ #pragma once +#include <corsika/process/BaseProcess.h> #include <corsika/process/ProcessReturn.h> // for convenience #include <corsika/setup/SetupTrajectory.h> #include <corsika/units/PhysicalUnits.h> @@ -23,19 +24,16 @@ namespace corsika::process { */ - template <typename derived> - class StackProcess { + template <typename TDerived> + class StackProcess : public BaseProcess<TDerived> { public: StackProcess() = delete; StackProcess(const unsigned int nStep) : fNStep(nStep) {} - derived& GetRef() { return static_cast<derived&>(*this); } - const derived& GetRef() const { return static_cast<const derived&>(*this); } - /// here starts the interface-definition part - // -> enforce derived to implement DoStack... + // -> enforce TDerived to implement DoStack... template <typename TStack> inline EProcessReturn DoStack(TStack&); diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc index 9907d3ef8ef7755031f6c2bd8b6473bcf94ac17f..2f58e96aa5593ddd0ceec6d723d1cbc404ec5d09 100644 --- a/Processes/ObservationPlane/ObservationPlane.cc +++ b/Processes/ObservationPlane/ObservationPlane.cc @@ -38,11 +38,8 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous( << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m << std::endl; - if (deleteOnHit_) { - return process::EProcessReturn::eParticleAbsorbed; - } else { - return process::EProcessReturn::eOk; - } + if (deleteOnHit_) { return process::EProcessReturn::eParticleAbsorbed; } + return process::EProcessReturn::eOk; } LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&, diff --git a/Processes/ParticleCut/CMakeLists.txt b/Processes/ParticleCut/CMakeLists.txt index eed2c18f29e715c8794310f420ed58cc0dff03c1..db1ee34c6db71a4013d19f3d35b2d9a5799ba24f 100644 --- a/Processes/ParticleCut/CMakeLists.txt +++ b/Processes/ParticleCut/CMakeLists.txt @@ -21,7 +21,6 @@ set_target_properties ( PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION 1 -# PUBLIC_HEADER "${MODEL_HEADERS}" ) # target dependencies on other libraries (also the header onlys) @@ -31,6 +30,7 @@ target_link_libraries ( CORSIKAparticles CORSIKAprocesssequence CORSIKAsetup + CORSIKAlogging ) target_include_directories ( @@ -44,7 +44,6 @@ install ( TARGETS ProcessParticleCut LIBRARY DESTINATION lib ARCHIVE DESTINATION lib -# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE} ) # -------------------- diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc index defa0c1014dfaa3afd0e52b459d2bafa963a36c4..bd5a5fb61ad5a681bb6940e5808100d7526bdefc 100644 --- a/Processes/ParticleCut/ParticleCut.cc +++ b/Processes/ParticleCut/ParticleCut.cc @@ -6,6 +6,7 @@ * the license. */ +#include <corsika/logging/Logging.h> #include <corsika/process/particle_cut/ParticleCut.h> using namespace std; @@ -57,39 +58,57 @@ namespace corsika::process { } } + template <typename TParticle> + bool ParticleCut::checkCutParticle(const TParticle& particle) { + + const Code pid = particle.GetPID(); + HEPEnergyType energy = particle.GetEnergy(); + C8LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid, + energy / 1_GeV, + (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV)); + if (ParticleIsEmParticle(pid)) { + C8LOG_DEBUG("removing em. particle..."); + fEmEnergy += energy; + fEmCount += 1; + return true; + } else if (ParticleIsInvisible(pid)) { + C8LOG_DEBUG("removing inv. particle..."); + fInvEnergy += energy; + fInvCount += 1; + return true; + } else if (ParticleIsBelowEnergyCut(particle)) { + C8LOG_DEBUG("removing low en. particle..."); + fEnergy += energy; + return true; + } else if (particle.GetTime() > 10_ms) { + C8LOG_DEBUG("removing OLD particle..."); + fEnergy += energy; + return true; + } + return false; // this particle will not be removed/cut + } + EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& vS) { - auto p = vS.begin(); - while (p != vS.end()) { - const Code pid = p.GetPID(); - HEPEnergyType energy = p.GetEnergy(); - cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy - << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" - << endl; - if (ParticleIsEmParticle(pid)) { - cout << "removing em. particle..." << endl; - fEmEnergy += energy; - fEmCount += 1; - p.Delete(); - } else if (ParticleIsInvisible(pid)) { - cout << "removing inv. particle..." << endl; - fInvEnergy += energy; - fInvCount += 1; - p.Delete(); - } else if (ParticleIsBelowEnergyCut(p)) { - cout << "removing low en. particle..." << endl; - fEnergy += energy; - p.Delete(); - } else if (p.GetTime() > 10_ms) { - cout << "removing OLD particle..." << endl; - fEnergy += energy; - p.Delete(); + auto particle = vS.begin(); + while (particle != vS.end()) { + if (checkCutParticle(particle)) { + particle.Delete(); } else { - ++p; // next entry in SecondaryView + ++particle; // next entry in SecondaryView } } return EProcessReturn::eOk; } + process::EProcessReturn ParticleCut::DoContinuous(Particle& particle, Track const&) { + C8LOG_TRACE("ParticleCut::DoContinuous"); + if (checkCutParticle(particle)) { + C8LOG_TRACE("removing during continuous"); + return process::EProcessReturn::eParticleAbsorbed; + } + return process::EProcessReturn::eOk; + } + ParticleCut::ParticleCut(const units::si::HEPEnergyType vCut) : fECut(vCut) { @@ -101,14 +120,16 @@ namespace corsika::process { } void ParticleCut::ShowResults() { - cout << " ******************************" << endl - << " ParticleCut: " << endl - << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl - << " no. of em. particles injected: " << fEmCount << endl - << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl - << " no. of inv. particles injected: " << fInvCount << endl - << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl - << " ******************************" << endl; + C8LOG_INFO(fmt::format( + " ******************************\n" + " ParticleCut: \n" + " energy in em. component (GeV): {}\n" + " no. of em. particles injected: {}\n" + " energy in inv. component (GeV): {}\n" + " no. of inv. particles injected: {}\n" + " energy below particle cut (GeV): {}\n" + " ******************************", + fEmEnergy / 1_GeV, fEmCount, fInvEnergy / 1_GeV, fInvCount, fEnergy / 1_GeV)); } } // namespace particle_cut } // namespace corsika::process diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h index 53571551c445852648e1e86a4a713ee22bd7f64e..4c0f6371ac023798cb739dfce1b2f94e37a3fce6 100644 --- a/Processes/ParticleCut/ParticleCut.h +++ b/Processes/ParticleCut/ParticleCut.h @@ -9,13 +9,18 @@ #pragma once #include <corsika/particles/ParticleProperties.h> +#include <corsika/process/ContinuousProcess.h> #include <corsika/process/SecondariesProcess.h> #include <corsika/setup/SetupStack.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::process { namespace particle_cut { - class ParticleCut : public process::SecondariesProcess<ParticleCut> { + class ParticleCut : public process::SecondariesProcess<ParticleCut>, + public corsika::process::ContinuousProcess<ParticleCut> { + + using Particle = corsika::setup::Stack::ParticleType; + using Track = corsika::setup::Trajectory; units::si::HEPEnergyType const fECut; @@ -28,13 +33,14 @@ namespace corsika::process { public: ParticleCut(const units::si::HEPEnergyType vCut); - bool ParticleIsInvisible(particles::Code) const; EProcessReturn DoSecondaries(corsika::setup::StackView&); - template <typename TParticle> - bool ParticleIsBelowEnergyCut(TParticle const&) const; + EProcessReturn DoContinuous(Particle& vParticle, Track const& vTrajectory); - bool ParticleIsEmParticle(particles::Code) const; + corsika::units::si::LengthType MaxStepLength( + corsika::setup::Stack::ParticleType const&, corsika::setup::Trajectory const&) { + return units::si::meter * std::numeric_limits<double>::infinity(); + } void ShowResults(); @@ -43,6 +49,16 @@ namespace corsika::process { units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; } unsigned int GetNumberEmParticles() const { return fEmCount; } unsigned int GetNumberInvParticles() const { return fInvCount; } + + protected: + template <typename TParticle> + bool checkCutParticle(const TParticle& p); + + template <typename TParticle> + bool ParticleIsBelowEnergyCut(TParticle const&) const; + + bool ParticleIsEmParticle(particles::Code) const; + bool ParticleIsInvisible(particles::Code) const; }; } // namespace particle_cut } // namespace corsika::process diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc index 8c4506120a98297582ee5d1a09154b084450fa87..a3d3edbf39eb27f4aaaf66c6255f6cb879572a9a 100644 --- a/Processes/TrackWriter/TrackWriter.cc +++ b/Processes/TrackWriter/TrackWriter.cc @@ -53,9 +53,4 @@ namespace corsika::process::track_writer { return process::EProcessReturn::eOk; } - template <> - units::si::LengthType TrackWriter::MaxStepLength(Particle&, Track&) { - return units::si::meter * std::numeric_limits<double>::infinity(); - } - } // namespace corsika::process::track_writer diff --git a/Processes/TrackWriter/TrackWriter.h b/Processes/TrackWriter/TrackWriter.h index f253c6d42fcd4af6893ae02df84361d26a445abc..b62e09feac758cf1a30ef1966d4e02379313076e 100644 --- a/Processes/TrackWriter/TrackWriter.h +++ b/Processes/TrackWriter/TrackWriter.h @@ -25,7 +25,9 @@ namespace corsika::process::track_writer { corsika::process::EProcessReturn DoContinuous(Particle&, Track&); template <typename Particle, typename Track> - corsika::units::si::LengthType MaxStepLength(Particle&, Track&); + corsika::units::si::LengthType MaxStepLength(Particle&, Track&) { + return units::si::meter * std::numeric_limits<double>::infinity(); + } private: std::string const fFilename; diff --git a/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h b/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h new file mode 100644 index 0000000000000000000000000000000000000000..d33997b8eecddca3d53be9fd4f1f7b1b95ad4835 --- /dev/null +++ b/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h @@ -0,0 +1,99 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/logging/Logging.h> +#include <corsika/stack/Stack.h> + +#include <tuple> +#include <utility> +#include <vector> + +namespace corsika::stack::node { + + /** + * @class GeometryDataInterface + * + * corresponding defintion of a stack-readout object, the iteractor + * dereference operator will deliver access to these function + // defintion of a stack-readout object, the iteractor dereference + // operator will deliver access to these function + */ + template <typename T, typename TEnvType> + class GeometryDataInterface : public T { + + protected: + using T::GetStack; + using T::GetStackData; + + public: + using T::GetIndex; + using BaseNodeType = typename TEnvType::BaseNodeType; + + public: + // default version for particle-creation from input data + void SetParticleData(const std::tuple<BaseNodeType const*> v) { + SetNode(std::get<0>(v)); + } + void SetParticleData(GeometryDataInterface& parent, + const std::tuple<BaseNodeType const*>) { + SetNode(parent.GetNode()); // copy Node from parent particle! + } + void SetParticleData() { SetNode(nullptr); } + void SetParticleData(GeometryDataInterface& parent) { + SetNode(parent.GetNode()); // copy Node from parent particle! + } + + std::string as_string() const { return fmt::format("node={}", fmt::ptr(GetNode())); } + + void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); } + BaseNodeType const* GetNode() const { return GetStackData().GetNode(GetIndex()); } + }; + + // definition of stack-data object to store geometry information + template <typename TEnvType> + + /** + * @class GeometryData + * + * definition of stack-data object to store geometry information + */ + class GeometryData { + + public: + using BaseNodeType = typename TEnvType::BaseNodeType; + + // these functions are needed for the Stack interface + void Clear() { fNode.clear(); } + unsigned int GetSize() const { return fNode.size(); } + unsigned int GetCapacity() const { return fNode.size(); } + void Copy(const int i1, const int i2) { fNode[i2] = fNode[i1]; } + void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); } + + // custom data access function + void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; } + BaseNodeType const* GetNode(const int i) const { return fNode[i]; } + + // these functions are also needed by the Stack interface + void IncrementSize() { fNode.push_back(nullptr); } + void DecrementSize() { + if (fNode.size() > 0) { fNode.pop_back(); } + } + + // custom private data section + private: + std::vector<const BaseNodeType*> fNode; + }; + + template <typename T, typename TEnv> + struct MakeGeometryDataInterface { + typedef GeometryDataInterface<T, TEnv> type; + }; + +} // namespace corsika::stack::node diff --git a/Stack/History/HistoryObservationPlane.cc b/Stack/History/HistoryObservationPlane.cc new file mode 100644 index 0000000000000000000000000000000000000000..d9374eee7e00ee8177cd3027b86a9fac7bc6b9cd --- /dev/null +++ b/Stack/History/HistoryObservationPlane.cc @@ -0,0 +1,82 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <corsika/logging/Logging.h> +#include <corsika/history/HistoryObservationPlane.hpp> + +#include <boost/histogram/ostream.hpp> + +#include <fstream> +#include <iostream> + +using namespace corsika::units::si; +using namespace corsika::history; +using namespace corsika; + +HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack, + geometry::Plane const& obsPlane, + bool deleteOnHit) + : stack_{stack} + , plane_{obsPlane} + , deleteOnHit_{deleteOnHit} {} + +corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous( + setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) { + TimeType const timeOfIntersection = + (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) / + trajectory.GetV0().dot(plane_.GetNormal()); + + if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; } + + if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) { + return process::EProcessReturn::eOk; + } + + C8LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}", + particle.GetPID())); + + auto const pid = particle.GetPID(); + if (particles::IsMuon(pid)) { fillHistoryHistogram(particle); } + + if (deleteOnHit_) { + return process::EProcessReturn::eParticleAbsorbed; + } else { + return process::EProcessReturn::eOk; + } +} + +LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&, + setup::Trajectory const& trajectory) { + TimeType const timeOfIntersection = + (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) / + trajectory.GetV0().dot(plane_.GetNormal()); + + if (timeOfIntersection < TimeType::zero()) { + return std::numeric_limits<double>::infinity() * 1_m; + } + + auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection); + return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001; +} + +void HistoryObservationPlane::fillHistoryHistogram( + setup::Stack::ParticleType const& muon) { + // double const muonEnergy = muon.GetEnergy() / 1_eV; + + // auto parent = stack_.begin() + muon.GetEvent()->projectileIndex(); + Event* event = muon.GetEvent().get(); + + int intCounter = 0; + while (event) { + event = event->parentEvent().get(); + intCounter++; + } + histogram_(intCounter); +} + +void HistoryObservationPlane::print() { std::cout << histogram_ << std::endl; } diff --git a/Stack/History/HistorySecondaryProducer.hpp b/Stack/History/HistorySecondaryProducer.hpp new file mode 100644 index 0000000000000000000000000000000000000000..2d3d560a32e9a242ee32488476f988e939c54d27 --- /dev/null +++ b/Stack/History/HistorySecondaryProducer.hpp @@ -0,0 +1,71 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/stack/SecondaryView.h> +#include <corsika/history/Event.hpp> + +#include <corsika/logging/Logging.h> + +#include <memory> +#include <type_traits> +#include <utility> + +namespace corsika::history { + + template <class T1, template <class> class T2> + class HistorySecondaryProducer { + + using TView = corsika::stack::SecondaryView<T1, T2, HistorySecondaryProducer>; + + public: + EventPtr event_; + + public: + HistorySecondaryProducer() : + event_{std::make_shared<Event>()} { + C8LOG_TRACE("HistorySecondaryProducer::HistorySecondaryProducer"); + } + + /** + * Method is called after a new SecondaryView has been + * created. Extra logic can be introduced here. + * + * The input Particle is a reference object into the original + * parent stack! It is not a reference into the SecondaryView + * itself. + */ + template <typename Particle> + void new_view(Particle& p) { + C8LOG_TRACE("HistorySecondaryProducer::new_view"); + event_->setProjectileIndex(p.GetIndex()); + event_->setParentEvent(p.GetEvent()); + } + + /** + * Method is called after a new Secondary has been created on the + * SecondaryView. Extra logic can be introduced here. + * + * The input Particle is the new secondary that was produced and + * is of course a reference into the SecondaryView itself. + */ + template <typename Particle> + auto new_secondary(Particle& sec) { + C8LOG_TRACE("HistorySecondaryProducer::new_secondary(sec)"); + + // store particles at production time in Event here + auto const sec_index = event_->addSecondary( + sec.GetEnergy(), sec.GetMomentum(), sec.GetPID()); + sec.SetParentEventIndex(sec_index); + sec.SetEvent(event_); + } + + }; + +} // namespace corsika::history