From 3dfcaef7841f81b97b9f51b89032219181af6910 Mon Sep 17 00:00:00 2001 From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu> Date: Thu, 21 Mar 2019 21:43:27 +0100 Subject: [PATCH] begin making Environment templated --- Documentation/Examples/cascade_example.cc | 64 +++++++++++-------- Environment/CMakeLists.txt | 1 + Environment/Environment.h | 18 +++--- Environment/NameModel.h | 35 ++++++++++ Environment/VolumeTreeNode.h | 12 ++-- Framework/Cascade/CMakeLists.txt | 1 + Framework/Cascade/Cascade.h | 7 +- Framework/Cascade/testCascade.cc | 16 ++--- Framework/Cascade/testCascade.h | 19 ++++++ Framework/Geometry/Trajectory.h | 2 +- .../ProcessSequence/BoundaryCrossingProcess.h | 8 +-- Framework/ProcessSequence/ProcessSequence.h | 6 +- Framework/StackInterface/ParticleBase.h | 2 +- .../HadronicElasticModel.cc | 6 +- .../HadronicElasticModel.h | 4 +- Processes/Sibyll/Interaction.cc | 16 ++--- Processes/Sibyll/Interaction.h | 7 +- Processes/Sibyll/NuclearInteraction.cc | 13 ++-- Processes/Sibyll/NuclearInteraction.h | 8 +-- Processes/Sibyll/testSibyll.cc | 10 +-- Processes/StackInspector/CMakeLists.txt | 1 + Processes/StackInspector/StackInspector.cc | 4 +- Processes/TrackingLine/TrackingLine.cc | 19 +----- Processes/TrackingLine/TrackingLine.h | 33 +++++----- .../TrackingLine/testTrackingLineStack.h | 2 +- Setup/SetupEnvironment.h | 5 +- Setup/SetupStack.h | 32 ++++++---- 27 files changed, 200 insertions(+), 151 deletions(-) create mode 100644 Environment/NameModel.h create mode 100644 Framework/Cascade/testCascade.h diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index cdec9d648..9ff2e73cd 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -16,8 +16,10 @@ #include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/setup/SetupStack.h> +#include <corsika/setup/SetupEnvironment.h> #include <corsika/setup/SetupTrajectory.h> +#include <corsika/environment/NameModel.h> #include <corsika/environment/Environment.h> #include <corsika/environment/HomogeneousMedium.h> #include <corsika/environment/NuclearComposition.h> @@ -214,16 +216,25 @@ public: HEPEnergyType GetEmEnergy() const { return fEmEnergy; } }; -struct MyBoundaryCrossingProcess - : public BoundaryCrossingProcess<MyBoundaryCrossingProcess> { - template <typename Particle> - EProcessReturn DoBoundaryCrossing(Particle&, environment::BaseNodeType const& from, - environment::BaseNodeType const& to) { - std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl; - return EProcessReturn::eOk; - } - - void Init() {} +struct MyBoundaryCrossingProcess : public BoundaryCrossingProcess<MyBoundaryCrossingProcess> { + //~ environment::BaseNodeType const& fA, fB; + + MyBoundaryCrossingProcess() {} + + //~ MyBoundaryCrossingProcess(environment::BaseNodeType const& a, environment::BaseNodeType const& b) : fA(a), fB(b) {} + + template <typename Particle> + EProcessReturn DoBoundaryCrossing(Particle& p, environment::BaseNodeType const& from, environment::BaseNodeType const& to) { + std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl; + + //~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) { + p.Delete(); + //~ } + + return EProcessReturn::eOk; + } + + void Init() {} }; // @@ -244,24 +255,25 @@ int main() { // fraction of oxygen const float fox = 0.20946; - using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>; - outerMedium->SetModelProperties<MyHomogeneousModel>( + using MyHomogeneousModel = environment::HomogeneousMedium<setup::IEnvironmentModel>; + outerMedium->SetModelProperties<setup::IEnvironmentModel>("outer", 1_kg / (1_m * 1_m * 1_m), environment::NuclearComposition( std::vector<particles::Code>{particles::Code::Nitrogen, particles::Code::Oxygen}, std::vector<float>{(float)1. - fox, fox})); - + auto innerMedium = environment::Environment::CreateNode<Sphere>( - Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 10_m); + Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, + 2000_m); - innerMedium->SetModelProperties<MyHomogeneousModel>( + innerMedium->SetModelProperties<setup::IEnvironmentModel>("inner", 1_kg / (1_m * 1_m * 1_m), environment::NuclearComposition( std::vector<particles::Code>{particles::Code::Nitrogen, particles::Code::Oxygen}, std::vector<float>{(float)1. - fox, fox})); - + outerMedium->AddChild(std::move(innerMedium)); universe.AddChild(std::move(outerMedium)); @@ -279,26 +291,26 @@ int main() { ProcessCut cut(20_GeV); random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel"); - process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env); + process::HadronicElasticModel::HadronicElasticInteraction + hadronicElastic(env); process::TrackWriter::TrackWriter trackWriter("tracks.dat"); // assemble all processes into an ordered process list // auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter; - auto sequence = p0 << sibyll << sibyllNuc << decay << cut << MyBoundaryCrossingProcess() - << trackWriter; + auto sequence = p0 /*<< sibyll << sibyllNuc << decay << cut*/ << hadronicElastic << MyBoundaryCrossingProcess() << trackWriter; // setup particle stack, and add primary particle setup::Stack stack; stack.Clear(); - const Code beamCode = Code::Nucleus; + const Code beamCode = Code::Proton; const int nuclA = 56; const int nuclZ = int(nuclA / 2.15 + 0.7); const HEPMassType mass = particles::Proton::GetMass() * nuclZ + (nuclA - nuclZ) * particles::Neutron::GetMass(); - const HEPEnergyType E0 = nuclA * 100_GeV; - double theta = 27.234; - double phi = 0.; + const HEPEnergyType E0 = nuclA * 100_TeV; + double theta = 0; + double phi = 0; { auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { @@ -316,10 +328,12 @@ int main() { cout << "input angles: theta=" << theta << " phi=" << phi << endl; cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl; Point pos(rootCS, 0_m, 0_m, 0_m); + for (int l = 0; l < 100; ++l) { stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, geometry::Point, - units::si::TimeType, unsigned short, unsigned short>{ - beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ}); + units::si::TimeType>{ + beamCode, E0, plab, pos, 0_ns}); + } } // define air shower object, run simulation diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt index ef1d24035..b0e4a43f9 100644 --- a/Environment/CMakeLists.txt +++ b/Environment/CMakeLists.txt @@ -9,6 +9,7 @@ set ( LinearApproximationIntegrator.h DensityFunction.h Environment.h + NameModel.h ) set ( diff --git a/Environment/Environment.h b/Environment/Environment.h index 97931aeff..338da588a 100644 --- a/Environment/Environment.h +++ b/Environment/Environment.h @@ -9,20 +9,17 @@ * the license. */ -#ifndef _include_Environment_h -#define _include_Environment_h +#ifndef _include_environment_Environment_h +#define _include_environment_Environment_h #include <corsika/environment/IMediumModel.h> #include <corsika/environment/VolumeTreeNode.h> #include <corsika/geometry/Point.h> #include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Sphere.h> -#include <corsika/setup/SetupEnvironment.h> #include <limits> namespace corsika::environment { - using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>; - struct Universe : public corsika::geometry::Sphere { Universe(corsika::geometry::CoordinateSystem const& pCS) : corsika::geometry::Sphere( @@ -34,16 +31,18 @@ namespace corsika::environment { bool Contains(corsika::geometry::Point const&) const override { return true; } }; - // template <typename IEnvironmentModel> + template <typename IEnvironmentModel> class Environment { public: + using BaseNodeType = VolumeTreeNode<IEnvironmentModel>; + Environment() : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance() .GetRootCoordinateSystem()} , fUniverse(std::make_unique<BaseNodeType>( std::make_unique<Universe>(fCoordinateSystem))) {} - using IEnvironmentModel = corsika::setup::IEnvironmentModel; + // using IEnvironmentModel = corsika::setup::IEnvironmentModel; auto& GetUniverse() { return fUniverse; } auto const& GetUniverse() const { return fUniverse; } @@ -63,9 +62,12 @@ namespace corsika::environment { private: corsika::geometry::CoordinateSystem const& fCoordinateSystem; - BaseNodeType::VTNUPtr fUniverse; + typename BaseNodeType::VTNUPtr fUniverse; }; + //using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>; + //using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>; + } // namespace corsika::environment #endif diff --git a/Environment/NameModel.h b/Environment/NameModel.h new file mode 100644 index 000000000..29321e5ce --- /dev/null +++ b/Environment/NameModel.h @@ -0,0 +1,35 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * 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. + */ + +#ifndef _include_NameModel_h +#define _include_NameModel_h + +#include <string> +#include <utility> + +namespace corsika::environment { + + template <typename T> + struct NameModel : public T { + + template <typename... Args> + NameModel(std::string const& name, Args&&... args) : T(std::forward<Args>(args)...), fName(name) {} + + std::string const& GetName() const { + return fName; + } + + private: + std::string fName; + }; + +} // namespace corsika::environment + +#endif diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h index e125b243d..b2ad36203 100644 --- a/Environment/VolumeTreeNode.h +++ b/Environment/VolumeTreeNode.h @@ -21,9 +21,10 @@ namespace corsika::environment { class Empty {}; //<! intended for usage as default template argument - template <typename IModelProperties = Empty> + template <typename TModelProperties = Empty> class VolumeTreeNode { public: + using IModelProperties = TModelProperties; using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>; using IMPSharedPtr = std::shared_ptr<IModelProperties>; using VolUPtr = std::unique_ptr<corsika::geometry::Volume>; @@ -31,6 +32,7 @@ namespace corsika::environment { VolumeTreeNode(VolUPtr pVolume = nullptr) : fGeoVolume(std::move(pVolume)) {} + //! convenience function equivalent to Volume::Contains bool Contains(corsika::geometry::Point const& p) const { return fGeoVolume->Contains(p); } @@ -45,7 +47,7 @@ namespace corsika::environment { } /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given - * \class Point \arg p, or nullptr iff \arg p is not contained in this volume. + * \class Point \p p, or nullptr iff \p p is not contained in this volume. */ VolumeTreeNode<IModelProperties> const* GetContainingNode( corsika::geometry::Point const& p) const { @@ -89,12 +91,12 @@ namespace corsika::environment { auto const& GetModelProperties() const { return *fModelProperties; } - template <typename TModelProperties, typename... Args> + template <typename ModelProperties, typename... Args> auto SetModelProperties(Args&&... args) { - static_assert(std::is_base_of_v<IModelProperties, TModelProperties>, + static_assert(std::is_base_of_v<IModelProperties, ModelProperties>, "unusable type provided"); - fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...); + fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...); return fModelProperties; } diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt index 646692248..3e294b214 100644 --- a/Framework/Cascade/CMakeLists.txt +++ b/Framework/Cascade/CMakeLists.txt @@ -9,6 +9,7 @@ set ( set ( CORSIKAcascade_HEADERS Cascade.h + testCascade.h ) add_library (CORSIKAcascade INTERFACE) diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index e4a79f7ca..4a82d38be 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -22,6 +22,7 @@ #include <cmath> #include <iostream> +#include <type_traits> /** * The cascade namespace assembles all objects needed to simulate full particles cascades. @@ -53,6 +54,8 @@ namespace corsika::cascade { template <typename Tracking, typename ProcessList, typename Stack> class Cascade { using Particle = typename Stack::ParticleType; + using VolumeTreeNode = std::remove_pointer_t<decltype(((Particle*) nullptr)->GetNode())>; + using MediumInterface = typename VolumeTreeNode::IModelProperties; // we only want fully configured objects Cascade() = delete; @@ -62,7 +65,7 @@ namespace corsika::cascade { * Cascade class cannot be default constructed, but needs a valid * list of physics processes for configuration at construct time. */ - Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl, + Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr, ProcessList& pl, Stack& stack) : fEnvironment(env) , fTracking(tr) @@ -239,7 +242,7 @@ namespace corsika::cascade { } private: - corsika::environment::Environment const& fEnvironment; + corsika::environment::Environment<MediumInterface> const& fEnvironment; Tracking& fTracking; ProcessList& fProcessSequence; Stack& fStack; diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 1e8a7a4de..7d473eecc 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -11,7 +11,7 @@ #include <limits> -#include <corsika/environment/Environment.h> +#include <corsika/cascade/testCascade.h> #include <corsika/cascade/Cascade.h> @@ -25,7 +25,6 @@ #include <corsika/geometry/RootCoordinateSystem.h> #include <corsika/geometry/Vector.h> -#include <corsika/environment/Environment.h> #include <corsika/environment/HomogeneousMedium.h> #include <corsika/environment/NuclearComposition.h> @@ -46,11 +45,12 @@ using namespace corsika::geometry; #include <iostream> using namespace std; -environment::Environment MakeDummyEnv() { - environment::Environment env; // dummy environment + +auto MakeDummyEnv() { + TestEnvironmentType env; // dummy environment auto& universe = *(env.GetUniverse()); - auto theMedium = environment::Environment::CreateNode<Sphere>( + auto theMedium = TestEnvironmentType::CreateNode<Sphere>( Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); @@ -113,14 +113,14 @@ TEST_CASE("Cascade", "[Cascade]") { rmng.RegisterRandomStream("cascade"); auto env = MakeDummyEnv(); - tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env); + tracking_line::TrackingLine tracking; - stack_inspector::StackInspector<setup::Stack> p0(true); + stack_inspector::StackInspector<TestCascadeStack> p0(true); const HEPEnergyType Ecrit = 85_MeV; ProcessSplit p1(Ecrit); auto sequence = p0 << p1; - setup::Stack stack; + TestCascadeStack stack; cascade::Cascade EAS(env, tracking, sequence, stack); CoordinateSystem const& rootCS = diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h new file mode 100644 index 000000000..d8c553a0f --- /dev/null +++ b/Framework/Cascade/testCascade.h @@ -0,0 +1,19 @@ +#ifndef _Framework_Cascade_testCascade_h +#define _Framework_Cascade_testCascade_h + +#include <corsika/environment/Environment.h> +#include <corsika/setup/SetupStack.h> + +using TestEnvironmentType = corsika::environment::Environment<corsika::environment::IMediumModel>; + +template <typename T> +using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>; + +// combine particle data stack with geometry information for tracking +template <typename StackIter> +using StackWithGeometryInterface = + corsika::stack::CombinedParticleInterface<corsika::setup::detail::ParticleDataStack::PIType, + SetupGeometryDataInterface, StackIter>; +using TestCascadeStack = corsika::stack::CombinedStack<typename corsika::setup::detail::ParticleDataStack::StackImpl, GeometryData<TestEnvironmentType>, StackWithGeometryInterface>; + +#endif diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h index 320630610..dd47208fc 100644 --- a/Framework/Geometry/Trajectory.h +++ b/Framework/Geometry/Trajectory.h @@ -1,4 +1,3 @@ - /* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * @@ -13,6 +12,7 @@ #define _include_TRAJECTORY_H #include <corsika/geometry/Point.h> +#include <corsika/geometry/Line.h> #include <corsika/units/PhysicalUnits.h> namespace corsika::geometry { diff --git a/Framework/ProcessSequence/BoundaryCrossingProcess.h b/Framework/ProcessSequence/BoundaryCrossingProcess.h index 6157aeeed..18dcd85dd 100644 --- a/Framework/ProcessSequence/BoundaryCrossingProcess.h +++ b/Framework/ProcessSequence/BoundaryCrossingProcess.h @@ -22,11 +22,11 @@ namespace corsika::process { /** * This method is called when a particle crosses the boundary between the nodes - * \arg from and \arg to. + * \p from and \p to. */ - template <typename Particle> - EProcessReturn DoBoundaryCrossing(Particle&, environment::BaseNodeType const& from, - environment::BaseNodeType const& to); + template <typename Particle, typename VTNType> + EProcessReturn DoBoundaryCrossing(Particle&, VTNType const& from, + VTNType const& to); }; template <class T> diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index f055b826c..c697d09f4 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -71,9 +71,9 @@ namespace corsika::process { // example for a trait-based call: // void Hello() const { detail::CallHello<T1,T2>::Call(A, B); } - template <typename Particle> - EProcessReturn DoBoundaryCrossing(Particle& p, environment::BaseNodeType const& from, - environment::BaseNodeType const& to) { + template <typename Particle, typename VTNType> + EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from, + VTNType const& to) { EProcessReturn ret = EProcessReturn::eOk; if constexpr (std::is_base_of<BoundaryCrossingProcess<T1type>, T1type>::value || diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h index d145c1e88..b5704e762 100644 --- a/Framework/StackInterface/ParticleBase.h +++ b/Framework/StackInterface/ParticleBase.h @@ -90,7 +90,7 @@ namespace corsika::stack { // protected: /** - @name Access to underlying stack dfata, these are service + @name Access to underlying stack fData, these are service function for user classes. User code can only rely on GetIndex and GetStackData to retrieve data @{ diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc index f68e31410..681d8fc0d 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.cc +++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc @@ -35,16 +35,14 @@ namespace corsika::process::HadronicElasticModel { environment::Environment const& env, units::si::CrossSectionType x, units::si::CrossSectionType y) : fX(x) - , fY(y) - , fEnvironment(env) {} + , fY(y) {} template <> units::si::GrammageType HadronicElasticInteraction::GetInteractionLength( Particle const& p, Track&) { using namespace units::si; if (p.GetPID() == particles::Code::Proton) { - auto const* currentNode = - fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition()); + auto const* currentNode = p.GetNode(); auto const& mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h index f03963872..9760569ee 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.h +++ b/Processes/HadronicElasticModel/HadronicElasticModel.h @@ -50,14 +50,12 @@ namespace corsika::process::HadronicElasticModel { corsika::random::RNG& fRNG = corsika::random::RNGManager::GetInstance().GetRandomStream( "HadronicElasticModel"); - corsika::environment::Environment const& fEnvironment; inveV2 B(eV2 s) const; corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const; public: - HadronicElasticInteraction(corsika::environment::Environment const&, - // x & y values taken from DL for pp collisions + HadronicElasticInteraction(// x & y values taken from DL for pp collisions units::si::CrossSectionType x = 0.0217 * units::si::barn, units::si::CrossSectionType y = 0.05608 * units::si::barn); void Init(); diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index 4c5f60448..683730554 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -34,8 +34,7 @@ using Track = Trajectory; namespace corsika::process::sibyll { - Interaction::Interaction(environment::Environment const& env) - : fEnvironment(env) {} + Interaction::Interaction() {} Interaction::~Interaction() { cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl; @@ -127,8 +126,7 @@ namespace corsika::process::sibyll { ideally as full particle object so that the four momenta and the boosts can be defined.. */ - const auto currentNode = - fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition()); + auto const* currentNode = p.GetNode(); const auto mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); // determine average interaction length @@ -254,8 +252,8 @@ namespace corsika::process::sibyll { HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm()); // sample target mass number - const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig); - const auto& mediumComposition = + auto const* currentNode = p.GetNode(); + auto const& mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); // get cross sections for target materials /* @@ -268,13 +266,9 @@ namespace corsika::process::sibyll { for (size_t i = 0; i < compVec.size(); ++i) { auto const targetId = compVec[i]; - const auto [sigProd, sigEla, nNuc] = + [[maybe_unsused]] const auto [sigProd, sigEla, nNuc] = GetCrossSection(corsikaBeamId, targetId, Ecm); cross_section_of_components[i] = sigProd; - [[maybe_unused]] int ideleteme = - nNuc; // to avoid not used warning in array binding - [[maybe_unused]] auto sigElaCopy = - sigEla; // to avoid not used warning in array binding } const auto targetCode = mediumComposition.SampleTarget( diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h index 40463f665..4157738ff 100644 --- a/Processes/Sibyll/Interaction.h +++ b/Processes/Sibyll/Interaction.h @@ -18,10 +18,6 @@ #include <corsika/units/PhysicalUnits.h> #include <tuple> -namespace corsika::environment { - class Environment; -} - namespace corsika::process::sibyll { class Interaction : public corsika::process::InteractionProcess<Interaction> { @@ -31,7 +27,7 @@ namespace corsika::process::sibyll { bool fInitialized = false; public: - Interaction(corsika::environment::Environment const& env); + Interaction(); ~Interaction(); void Init(); @@ -60,7 +56,6 @@ namespace corsika::process::sibyll { corsika::process::EProcessReturn DoInteraction(Particle&, Stack&); private: - corsika::environment::Environment const& fEnvironment; corsika::random::RNG& fRNG = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm"); }; diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index 6a3736323..56a0c67f7 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -34,10 +34,8 @@ using Track = Trajectory; namespace corsika::process::sibyll { - NuclearInteraction::NuclearInteraction(environment::Environment const& env, - process::sibyll::Interaction& hadint) - : fEnvironment(env) - , fHadronicInteraction(hadint) {} + NuclearInteraction::NuclearInteraction(process::sibyll::Interaction& hadint) + : fHadronicInteraction(hadint) {} NuclearInteraction::~NuclearInteraction() { cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl; @@ -175,9 +173,8 @@ namespace corsika::process::sibyll { ideally as full particle object so that the four momenta and the boosts can be defined.. */ - const auto currentNode = - fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition()); - const auto mediumComposition = + auto const* const currentNode = p.GetNode(); + auto const& mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); // determine average interaction length // weighted sum @@ -332,7 +329,7 @@ namespace corsika::process::sibyll { // // proton stand-in for nucleon const auto beamId = particles::Proton::GetCode(); - const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig); + auto const* const currentNode = p.GetNode(); const auto& mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); cout << "get nucleon-nucleus cross sections for target materials.." << endl; diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h index cf1f65933..920904eda 100644 --- a/Processes/Sibyll/NuclearInteraction.h +++ b/Processes/Sibyll/NuclearInteraction.h @@ -16,10 +16,6 @@ #include <corsika/process/InteractionProcess.h> #include <corsika/random/RNGManager.h> -namespace corsika::environment { - class Environment; -} - namespace corsika::process::sibyll { class Interaction; // fwd-decl @@ -36,8 +32,7 @@ namespace corsika::process::sibyll { int fNucCount = 0; public: - NuclearInteraction(corsika::environment::Environment const& env, - corsika::process::sibyll::Interaction& hadint); + NuclearInteraction(corsika::process::sibyll::Interaction& hadint); ~NuclearInteraction(); void Init(); @@ -52,7 +47,6 @@ namespace corsika::process::sibyll { corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s); private: - corsika::environment::Environment const& fEnvironment; corsika::process::sibyll::Interaction& fHadronicInteraction; corsika::random::RNG& fRNG = corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm"); diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index 93a2840e0..6d286624e 100644 --- a/Processes/Sibyll/testSibyll.cc +++ b/Processes/Sibyll/testSibyll.cc @@ -85,10 +85,10 @@ using namespace corsika::units; TEST_CASE("SibyllInterface", "[processes]") { // setup environment, geometry - environment::Environment env; + environment::Environment<environment::IMediumModel> env; auto& universe = *(env.GetUniverse()); - auto theMedium = environment::Environment::CreateNode<geometry::Sphere>( + auto theMedium = environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>( geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); @@ -123,7 +123,7 @@ TEST_CASE("SibyllInterface", "[processes]") { corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ particles::Code::Proton, E0, plab, pos, 0_ns}); - Interaction model(env); + Interaction model; model.Init(); [[maybe_unused]] const process::EProcessReturn ret = @@ -147,8 +147,8 @@ TEST_CASE("SibyllInterface", "[processes]") { units::si::TimeType, unsigned short, unsigned short>{ particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2}); - Interaction hmodel(env); - NuclearInteraction model(env, hmodel); + Interaction hmodel; + NuclearInteraction model(hmodel); model.Init(); [[maybe_unused]] const process::EProcessReturn ret = diff --git a/Processes/StackInspector/CMakeLists.txt b/Processes/StackInspector/CMakeLists.txt index 96b502f60..f00f44c2a 100644 --- a/Processes/StackInspector/CMakeLists.txt +++ b/Processes/StackInspector/CMakeLists.txt @@ -28,6 +28,7 @@ set_target_properties ( # target dependencies on other libraries (also the header onlys) target_link_libraries ( ProcessStackInspector + CORSIKAcascade CORSIKAunits CORSIKAgeometry CORSIKAsetup diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc index fb749c0b8..8552da191 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -17,6 +17,7 @@ #include <corsika/logging/Logger.h> #include <corsika/setup/SetupTrajectory.h> +#include <corsika/cascade/testCascade.h> #include <iostream> #include <limits> @@ -51,7 +52,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr auto pos = iterP.GetPosition().GetCoordinates(rootCS); cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30) << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, " - << " pos=" << pos; + << " pos=" << pos << " node = " << iterP.GetNode(); // if (iterP.GetPID()==Code::Nucleus) // cout << " nuc_ref=" << iterP.GetNucleusRef(); cout << endl; @@ -76,3 +77,4 @@ void StackInspector<Stack>::Init() { #include <corsika/setup/SetupStack.h> template class process::stack_inspector::StackInspector<setup::Stack>; +template class process::stack_inspector::StackInspector<TestCascadeStack>; diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc index dae65e44b..4806e7110 100644 --- a/Processes/TrackingLine/TrackingLine.cc +++ b/Processes/TrackingLine/TrackingLine.cc @@ -10,7 +10,6 @@ */ #include <corsika/process/tracking_line/TrackingLine.h> - #include <corsika/environment/Environment.h> #include <corsika/geometry/Point.h> #include <corsika/geometry/QuantityVector.h> @@ -26,9 +25,8 @@ using namespace corsika; namespace corsika::process::tracking_line { - template <class Stack, class Trajectory> std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>> - TrackingLine<Stack, Trajectory>::TimeOfIntersection(corsika::geometry::Line const& line, + TimeOfIntersection(corsika::geometry::Line const& line, geometry::Sphere const& sphere) { auto const delta = line.GetR0() - sphere.GetCenter(); auto const v = line.GetV0(); @@ -52,20 +50,5 @@ namespace corsika::process::tracking_line { return {}; } } - - template <class Stack, class Trajectory> - TrackingLine<Stack, Trajectory>::TrackingLine( - corsika::environment::Environment const& pEnv) - : fEnvironment(pEnv) {} - } // namespace corsika::process::tracking_line -#include <corsika/setup/SetupStack.h> -#include <corsika/setup/SetupTrajectory.h> -using namespace corsika::setup; -template class corsika::process::tracking_line::TrackingLine<setup::Stack, - setup::Trajectory>; - -#include "testTrackingLineStack.h" -template class corsika::process::tracking_line::TrackingLine<DummyStack, - setup::Trajectory>; diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h index a6d749d82..55c865295 100644 --- a/Processes/TrackingLine/TrackingLine.h +++ b/Processes/TrackingLine/TrackingLine.h @@ -12,15 +12,20 @@ #ifndef _include_corsika_processes_TrackingLine_h_ #define _include_corsika_processes_TrackingLine_h_ -#include <corsika/environment/Environment.h> -#include <corsika/environment/VolumeTreeNode.h> +//~ #include <corsika/environment/Environment.h> +//~ #include <corsika/environment/VolumeTreeNode.h> #include <corsika/units/PhysicalUnits.h> +#include <corsika/geometry/Vector.h> +#include <corsika/geometry/Trajectory.h> +#include <corsika/geometry/Line.h> +#include <corsika/geometry/Sphere.h> #include <optional> #include <type_traits> #include <utility> namespace corsika::environment { - class Environment; + template <typename IEnvironmentModel> class Environment; + template <typename IEnvironmentModel> class VolumeTreeNode; } namespace corsika::geometry { class Line; @@ -30,24 +35,18 @@ namespace corsika::geometry { namespace corsika::process { namespace tracking_line { - - template <typename Stack, typename Trajectory> - class TrackingLine { // - - using Particle = typename Stack::StackIterator; - - corsika::environment::Environment const& fEnvironment; - - public: + std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>> TimeOfIntersection(corsika::geometry::Line const& line, geometry::Sphere const& sphere); - TrackingLine(corsika::environment::Environment const& pEnv); + class TrackingLine { + + public: + TrackingLine() {}; + template <typename Particle> // was Stack previously, and argument was Stack::StackIterator auto GetTrack(Particle const& p) { - using std::cout; - using std::endl; using namespace corsika::units::si; using namespace corsika::geometry; geometry::Vector<SpeedType::dimension_type> const velocity = @@ -75,7 +74,7 @@ namespace corsika::process { auto const& children = currentLogicalVolumeNode->GetChildNodes(); auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes(); - std::vector<std::pair<TimeType, environment::BaseNodeType const*>> intersections; + std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections; auto addIfIntersects = [&](auto const& vtn, auto const& nextNode) { static_assert(std::is_same_v<decltype(vtn), decltype(nextNode)>); @@ -133,7 +132,7 @@ namespace corsika::process { std::cout << " t-intersect: " << min << std::endl; - return std::make_tuple(Trajectory(line, min), velocity.norm() * min, + return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min), velocity.norm() * min, minIter->second); } }; diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h index 819934d48..45975d032 100644 --- a/Processes/TrackingLine/testTrackingLineStack.h +++ b/Processes/TrackingLine/testTrackingLineStack.h @@ -24,7 +24,7 @@ struct DummyParticle { corsika::units::si::HEPEnergyType fEnergy; corsika::geometry::Vector<MOMENTUM> fMomentum; corsika::geometry::Point fPosition; - corsika::environment::BaseNodeType const* fNodePtr; + corsika::environment::VolumeTreeNode<> const* fNodePtr; DummyParticle(corsika::units::si::HEPEnergyType pEnergy, corsika::geometry::Vector<MOMENTUM> pMomentum, diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h index 08780b73b..f55cda350 100644 --- a/Setup/SetupEnvironment.h +++ b/Setup/SetupEnvironment.h @@ -13,9 +13,12 @@ #define _include_corsika_setup_environment_h_ #include <corsika/environment/IMediumModel.h> +#include <corsika/environment/NameModel.h> +#include <corsika/environment/Environment.h> namespace corsika::setup { - using IEnvironmentModel = corsika::environment::IMediumModel; + using IEnvironmentModel = environment::NameModel<environment::IMediumModel>; + using SetupEnvironment = environment::Environment<IEnvironmentModel>; } #endif diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h index 9f42d76c9..d8dbac906 100644 --- a/Setup/SetupStack.h +++ b/Setup/SetupStack.h @@ -21,14 +21,19 @@ // extension with geometry information for tracking #include <corsika/environment/Environment.h> #include <corsika/stack/CombinedStack.h> +#include <corsika/setup/SetupEnvironment.h> #include <tuple> +#include <utility> #include <vector> // definition of stack-data object to store geometry information +template <typename TEnvType> class GeometryData { public: + using BaseNodeType = typename TEnvType::BaseNodeType; + // these functions are needed for the Stack interface void Init() {} void Clear() { fNode.clear(); } @@ -38,8 +43,8 @@ public: void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); } // custom data access function - void SetNode(const int i, const corsika::environment::BaseNodeType* v) { fNode[i] = v; } - const corsika::environment::BaseNodeType* GetNode(const int i) const { + void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; } + auto const* GetNode(const int i) const { return fNode[i]; } @@ -51,35 +56,35 @@ public: // custom private data section private: - std::vector<const corsika::environment::BaseNodeType*> fNode; + std::vector<const BaseNodeType*> fNode; }; // defintion of a stack-readout object, the iteractor dereference // operator will deliver access to these function -template <typename T> +template <typename T, typename TEnvType> class GeometryDataInterface : public T { public: using T::GetIndex; using T::GetStackData; using T::SetParticleData; + using BaseNodeType = typename TEnvType::BaseNodeType; // default version for particle-creation from input data - void SetParticleData(const std::tuple<const corsika::environment::BaseNodeType*> v) { + void SetParticleData(const std::tuple<BaseNodeType const*> v) { SetNode(std::get<0>(v)); } - void SetParticleData(GeometryDataInterface<T>& parent, - const std::tuple<const corsika::environment::BaseNodeType*>) { + void SetParticleData(GeometryDataInterface& parent, const std::tuple<BaseNodeType const*>) { SetNode(parent.GetNode()); // copy Node from parent particle! } void SetParticleData() { SetNode(nullptr); } - void SetParticleData(GeometryDataInterface<T>& parent) { + void SetParticleData(GeometryDataInterface& parent) { SetNode(parent.GetNode()); // copy Node from parent particle! } - void SetNode(const corsika::environment::BaseNodeType* v) { + void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); } - const corsika::environment::BaseNodeType* GetNode() const { + auto const* GetNode() const { return GetStackData().GetNode(GetIndex()); } }; @@ -101,14 +106,17 @@ namespace corsika::setup { using ParticleDataStack = corsika::stack::nuclear_extension::NuclearStackExtension< corsika::stack::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>; + template <typename T> + using SetupGeometryDataInterface = GeometryDataInterface<T, setup::SetupEnvironment>; + // combine particle data stack with geometry information for tracking template <typename StackIter> using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<ParticleDataStack::PIType, - GeometryDataInterface, StackIter>; + SetupGeometryDataInterface, StackIter>; using StackWithGeometry = - corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl, GeometryData, + corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl, GeometryData<setup::SetupEnvironment>, StackWithGeometryInterface>; } // namespace detail -- GitLab