diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index 9ff2e73cd36379f40bf423a3c96b2894a2afbf7f..8b1f90be659916eb1d5e014f198bca484acdb82e 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/Documentation/Examples/cascade_example.cc @@ -15,13 +15,13 @@ #include <corsika/process/stack_inspector/StackInspector.h> #include <corsika/process/tracking_line/TrackingLine.h> -#include <corsika/setup/SetupStack.h> #include <corsika/setup/SetupEnvironment.h> +#include <corsika/setup/SetupStack.h> #include <corsika/setup/SetupTrajectory.h> -#include <corsika/environment/NameModel.h> #include <corsika/environment/Environment.h> #include <corsika/environment/HomogeneousMedium.h> +#include <corsika/environment/NameModel.h> #include <corsika/environment/NuclearComposition.h> #include <corsika/geometry/Sphere.h> @@ -216,25 +216,29 @@ public: HEPEnergyType GetEmEnergy() const { return fEmEnergy; } }; -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() {} +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, + typename Particle::BaseNodeType const& from, + typename Particle::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() {} }; // @@ -246,34 +250,34 @@ int main() { random::RNGManager::GetInstance().RegisterRandomStream("cascade"); // setup environment, geometry - environment::Environment env; + using EnvType = environment::Environment<setup::IEnvironmentModel>; + EnvType env; auto& universe = *(env.GetUniverse()); - auto outerMedium = environment::Environment::CreateNode<Sphere>( + auto outerMedium = environment::Environment<EnvType>::CreateNode<Sphere>( Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); // fraction of oxygen const float fox = 0.20946; using MyHomogeneousModel = environment::HomogeneousMedium<setup::IEnvironmentModel>; - outerMedium->SetModelProperties<setup::IEnvironmentModel>("outer", - 1_kg / (1_m * 1_m * 1_m), + 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}, - 2000_m); - innerMedium->SetModelProperties<setup::IEnvironmentModel>("inner", - 1_kg / (1_m * 1_m * 1_m), + auto innerMedium = environment::Environment<EnvType>::CreateNode<Sphere>( + Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 2000_m); + + 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)); @@ -291,14 +295,16 @@ 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*/ << hadronicElastic << MyBoundaryCrossingProcess() << trackWriter; + auto sequence = + p0 /*<< sibyll << sibyllNuc << decay << cut*/ << hadronicElastic + << MyBoundaryCrossingProcess() + << trackWriter; // setup particle stack, and add primary particle setup::Stack stack; @@ -329,10 +335,9 @@ int main() { 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>{ - beamCode, E0, plab, pos, 0_ns}); + stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, + corsika::stack::MomentumVector, geometry::Point, + units::si::TimeType>{beamCode, E0, plab, pos, 0_ns}); } } diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 7d473eeccbb022056ab816ea8d67a57b384f9b3d..5ec522a670de7ea048d33397e868590ea9e6faed 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -45,7 +45,6 @@ using namespace corsika::geometry; #include <iostream> using namespace std; - auto MakeDummyEnv() { TestEnvironmentType env; // dummy environment auto& universe = *(env.GetUniverse()); diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc index d73ae99b2f21590984aae1e90c683777549b7b32..dae272be1bc052b3174ad42ae0b2a3f50777322b 100644 --- a/Framework/Geometry/testGeometry.cc +++ b/Framework/Geometry/testGeometry.cc @@ -177,8 +177,10 @@ TEST_CASE("Trajectories") { CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates()); CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3)); - - CHECK((base.NormalizedDirection().GetComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}).eVector.norm() == Approx(0).margin(absMargin)); + + CHECK((base.NormalizedDirection().GetComponents(rootCS) - + QuantityVector<dimensionless_d>{1, 0, 0}) + .eVector.norm() == Approx(0).margin(absMargin)); } SECTION("Helix") { diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc index 681d8fc0d04125afea223cf644f05f1115e10549..42a42e0abd257c2b01f69e44c5ae8dcf3d6dba47 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.cc +++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc @@ -24,16 +24,15 @@ #include <iostream> using namespace corsika::setup; -using Particle = Stack::ParticleType; -using Track = Trajectory; +using Particle = corsika::setup::Stack::ParticleType; +using Track = corsika::setup::Trajectory; namespace corsika::process::HadronicElasticModel { void HadronicElasticInteraction::Init() {} - HadronicElasticInteraction::HadronicElasticInteraction( - environment::Environment const& env, units::si::CrossSectionType x, - units::si::CrossSectionType y) + HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x, + units::si::CrossSectionType y) : fX(x) , fY(y) {} @@ -87,8 +86,7 @@ namespace corsika::process::HadronicElasticModel { using namespace units::si; using namespace units::constants; - const auto* currentNode = - fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition()); + const auto* currentNode = p.GetNode(); const auto& composition = currentNode->GetModelProperties().GetNuclearComposition(); const auto& components = composition.GetComponents(); diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h index 9760569eeaa7c08649d5adcd811fa871c3f92109..d0826a304ec39c7f1b7fac8517cc869f8b21650d 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.h +++ b/Processes/HadronicElasticModel/HadronicElasticModel.h @@ -19,10 +19,6 @@ #include <corsika/units/PhysicalConstants.h> #include <corsika/units/PhysicalUnits.h> -namespace corsika::environment { - class Environment; -} - namespace corsika::process::HadronicElasticModel { /** * A simple model for elastic hadronic interactions based on the formulas diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index 6837305541343873041887a4f56af7b055fc6d87..290c72b2ad4abe36b79a074430aa4c53dd0d01ba 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -271,8 +271,8 @@ namespace corsika::process::sibyll { cross_section_of_components[i] = sigProd; } - const auto targetCode = mediumComposition.SampleTarget( - cross_section_of_components, fRNG); + const auto targetCode = + mediumComposition.SampleTarget(cross_section_of_components, fRNG); cout << "Interaction: target selected: " << targetCode << endl; /* FOR NOW: allow nuclei with A<18 or protons only. diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index 56a0c67f76ae3e14e613e7c620128594fb67ad71..23ec7510df721c468b044e502a996d3fbc0ecc98 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -353,7 +353,8 @@ namespace corsika::process::sibyll { [[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS } - const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG); + const auto targetCode = + mediumComposition.SampleTarget(cross_section_of_components, fRNG); cout << "Interaction: target selected: " << targetCode << endl; /* FOR NOW: allow nuclei with A<18 or protons only. diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index 6d286624e520c82a37c0d184ff1639df025a1acf..318bb5c2369129ff44e780d3b6fb5c4c2446f7f5 100644 --- a/Processes/Sibyll/testSibyll.cc +++ b/Processes/Sibyll/testSibyll.cc @@ -88,9 +88,10 @@ TEST_CASE("SibyllInterface", "[processes]") { environment::Environment<environment::IMediumModel> env; auto& universe = *(env.GetUniverse()); - 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()); + 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()); using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>; theMedium->SetModelProperties<MyHomogeneousModel>( diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc index 8552da1910f4fae141451a921eeb3f8c71c780f8..a22742165bec0a6f4a9bfe3b06e33e25ca2e9c3c 100644 --- a/Processes/StackInspector/StackInspector.cc +++ b/Processes/StackInspector/StackInspector.cc @@ -16,8 +16,8 @@ #include <corsika/logging/Logger.h> -#include <corsika/setup/SetupTrajectory.h> #include <corsika/cascade/testCascade.h> +#include <corsika/setup/SetupTrajectory.h> #include <iostream> #include <limits> diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc index 4806e7110f8961a5dc5ebd3c5eb17fc52771a6d5..6b119c890aa573b4d62fea9b693472384e1bdfaf 100644 --- a/Processes/TrackingLine/TrackingLine.cc +++ b/Processes/TrackingLine/TrackingLine.cc @@ -9,12 +9,12 @@ * the license. */ -#include <corsika/process/tracking_line/TrackingLine.h> #include <corsika/environment/Environment.h> #include <corsika/geometry/Point.h> #include <corsika/geometry/QuantityVector.h> #include <corsika/geometry/Sphere.h> #include <corsika/geometry/Vector.h> +#include <corsika/process/tracking_line/TrackingLine.h> #include <algorithm> #include <iostream> @@ -27,7 +27,7 @@ namespace corsika::process::tracking_line { std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>> TimeOfIntersection(corsika::geometry::Line const& line, - geometry::Sphere const& sphere) { + geometry::Sphere const& sphere) { auto const delta = line.GetR0() - sphere.GetCenter(); auto const v = line.GetV0(); auto const vSqNorm = v.squaredNorm(); @@ -51,4 +51,3 @@ namespace corsika::process::tracking_line { } } } // namespace corsika::process::tracking_line - diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h index 55c865295117b12d14ebe4959d97e1a746a01fc0..ef482fd3d8af4208e295b30abd73c1531a4cb112 100644 --- a/Processes/TrackingLine/TrackingLine.h +++ b/Processes/TrackingLine/TrackingLine.h @@ -74,7 +74,7 @@ namespace corsika::process { auto const& children = currentLogicalVolumeNode->GetChildNodes(); auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes(); - std::vector<std::pair<TimeType, decltype(p.GetNode())>> 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)>); diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc index 19efb14f30f13c6197f51a0659b18c4b2deaee7c..26dce27d178bae5bb7c7cb4f48b0f4d2c58df060 100644 --- a/Processes/TrackingLine/testTrackingLine.cc +++ b/Processes/TrackingLine/testTrackingLine.cc @@ -1,4 +1,3 @@ - /* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * @@ -36,10 +35,10 @@ using namespace std; using namespace corsika::units::si; TEST_CASE("TrackingLine") { - corsika::environment::Environment env; // dummy environment + environment::Environment<environment::Empty> env; // dummy environment auto const& cs = env.GetCoordinateSystem(); - tracking_line::TrackingLine<DummyStack, setup::Trajectory> tracking(env); + tracking_line::TrackingLine tracking; SECTION("intersection with sphere") { Point const origin(cs, {0_m, 0_m, 0_m}); @@ -52,7 +51,7 @@ TEST_CASE("TrackingLine") { setup::Trajectory traj(line, 12345_s); auto const opt = - tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m)); + tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m)); REQUIRE(opt.has_value()); auto [t1, t2] = opt.value(); @@ -60,25 +59,29 @@ TEST_CASE("TrackingLine") { REQUIRE(t2 / 11_s == Approx(1)); auto const optNoIntersection = - tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m)); + tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m)); REQUIRE_FALSE(optNoIntersection.has_value()); } SECTION("maximally possible propagation") { auto& universe = *(env.GetUniverse()); - //~ std::cout << env.GetUniverse().get() << std::endl; - auto const radius = 20_m; - auto theMedium = corsika::environment::Environment::CreateNode<Sphere>( + auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>( Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius); universe.AddChild(std::move(theMedium)); - - Point p0(cs, 0_m, 0_m, 0_m); - auto* const node = universe.GetContainingNode(p0); - DummyParticle p(1_GeV, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV), - p0, node); + + TestTrackingLineStack stack; + stack.AddParticle( + std::tuple<particles::Code, units::si::HEPEnergyType, + corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ + particles::Code::MuPlus, + 1_GeV, + {cs, {0_GeV, 0_GeV, 1_GeV}}, + {cs, {0_m, 0_m, 0_km}}, + 0_ns}); + auto p = stack.GetNextParticle(); Point const origin(cs, {0_m, 0_m, 0_m}); Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second, diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h index 45975d032d4b901b65c13868e95a92b79f816f11..5c3de074efeaf68ce66e4b7654795abaecaf8b91 100644 --- a/Processes/TrackingLine/testTrackingLineStack.h +++ b/Processes/TrackingLine/testTrackingLineStack.h @@ -1,4 +1,3 @@ - /* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * @@ -17,34 +16,19 @@ #include <corsika/geometry/Vector.h> #include <corsika/particles/ParticleProperties.h> #include <corsika/units/PhysicalUnits.h> +#include <corsika/setup/SetupStack.h> -typedef corsika::units::si::hepmomentum_d MOMENTUM; - -struct DummyParticle { - corsika::units::si::HEPEnergyType fEnergy; - corsika::geometry::Vector<MOMENTUM> fMomentum; - corsika::geometry::Point fPosition; - corsika::environment::VolumeTreeNode<> const* fNodePtr; +using TestEnvironmentType = corsika::environment::Environment<corsika::environment::Empty>; - DummyParticle(corsika::units::si::HEPEnergyType pEnergy, - corsika::geometry::Vector<MOMENTUM> pMomentum, - corsika::geometry::Point pPosition, - corsika::environment::BaseNodeType const* pNodePtr) - : fEnergy(pEnergy) - , fMomentum(pMomentum) - , fPosition(pPosition) - , fNodePtr(pNodePtr) {} +template <typename T> +using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>; - auto GetEnergy() const { return fEnergy; } - auto GetMomentum() const { return fMomentum; } - auto GetPosition() const { return fPosition; } - auto GetPID() const { return corsika::particles::Code::Unknown; } - auto* GetNode() const { return fNodePtr; } -}; +// 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 TestTrackingLineStack = corsika::stack::CombinedStack<typename corsika::setup::detail::ParticleDataStack::StackImpl, GeometryData<TestEnvironmentType>, StackWithGeometryInterface>; -struct DummyStack { - using ParticleType = DummyParticle; - using StackIterator = DummyParticle; -}; #endif