diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index 313dd46f8ff3622bb5ad1b4e7ef8ef73c7fcbf2b..758b430e677109de684eaa52f081c6a74bd09c3d 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -14,6 +14,7 @@ add_subdirectory (SwitchProcess) # continuous physics add_subdirectory (EnergyLoss) add_subdirectory (TrackWriter) +add_subdirectory (ObservationPlane) # stack processes add_subdirectory (StackInspector) # secondaries process diff --git a/Processes/ObservationPlane/CMakeLists.txt b/Processes/ObservationPlane/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..2707887a156b1c59ce74370fb617e434d155eeb7 --- /dev/null +++ b/Processes/ObservationPlane/CMakeLists.txt @@ -0,0 +1,44 @@ +set ( + MODEL_HEADERS + ObservationPlane.h + ) + +set ( + MODEL_SOURCES + ObservationPlane.cc + ) + +set ( + MODEL_NAMESPACE + corsika/process/observation_plane + ) + +add_library (ProcessObservationPlane STATIC ${MODEL_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessObservationPlane ${MODEL_NAMESPACE} ${MODEL_HEADERS}) + +# target dependencies on other libraries (also the header onlys) +target_link_libraries ( + ProcessObservationPlane + CORSIKAgeometry + CORSIKAprocesssequence + ) + +target_include_directories ( + ProcessObservationPlane + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE}) + +# -------------------- +# code unit testing +CORSIKA_ADD_TEST(testObservationPlane) +target_link_libraries ( + testObservationPlane + ProcessObservationPlane + CORSIKAstackinterface + CORSIKAthirdparty # for catch2 +) + diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc new file mode 100644 index 0000000000000000000000000000000000000000..69c2f1382b3ca569fe8d62901e1289338906e4b5 --- /dev/null +++ b/Processes/ObservationPlane/ObservationPlane.cc @@ -0,0 +1,50 @@ +/* + * (c) Copyright 2019 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. + */ + +#include <corsika/process/observation_plane/ObservationPlane.h> + +#include <fstream> + +using namespace corsika::process::observation_plane; +using namespace corsika::units::si; + +ObservationPlane::ObservationPlane(geometry::Plane const& vObsPlane, + std::string const& vFilename) + : fObsPlane(vObsPlane) + , fOutputStream(vFilename) { + fOutputStream << "#PDG code, energy / eV, distance to center / m" << std::endl; +} + +corsika::process::EProcessReturn ObservationPlane::DoContinuous( + setup::Stack::ParticleType const& vParticle, setup::Trajectory const& vTrajectory) { + + if (fObsPlane.IsAbove(vTrajectory.GetPosition(1.0001))) { + return process::EProcessReturn::eOk; + } + + fOutputStream << static_cast<int>(particles::GetPDG(vParticle.GetPID())) << ' ' + << vParticle.GetEnergy() * (1 / 1_eV) << ' ' + << (vTrajectory.GetPosition(1) - fObsPlane.GetCenter()).norm() / 1_m + << std::endl; + + return process::EProcessReturn::eParticleAbsorbed; +} + +LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&, + setup::Trajectory const& vTrajectory) { + if (!fObsPlane.IsAbove(vTrajectory.GetR0())) { + return std::numeric_limits<double>::infinity() * 1_m; + } + + auto const pointOfIntersection = vTrajectory.GetPosition( + (fObsPlane.GetCenter() - vTrajectory.GetR0()).dot(fObsPlane.GetNormal()) / + vTrajectory.GetV0().dot(fObsPlane.GetNormal())); + return (vTrajectory.GetR0() - pointOfIntersection).norm(); +} diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h new file mode 100644 index 0000000000000000000000000000000000000000..9ecec9f635fb83f5bc2726a16a554ffe10902a06 --- /dev/null +++ b/Processes/ObservationPlane/ObservationPlane.h @@ -0,0 +1,49 @@ +/* + * (c) Copyright 2019 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 _Processes_ObservationPlane_h_ +#define _Processes_ObservationPlane_h_ + +#include <corsika/geometry/Plane.h> +#include <corsika/process/ContinuousProcess.h> +#include <corsika/setup/SetupStack.h> +#include <corsika/setup/SetupTrajectory.h> +#include <corsika/units/PhysicalUnits.h> + +#include <fstream> + +namespace corsika::process::observation_plane { + + /** + * The ObservationPlane writes PDG codes, energies, and distances of particles to the + * central point of the plane into its output file. The particles are considered + * "absorbed" afterwards. + */ + class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> { + + public: + ObservationPlane(geometry::Plane const& vObsPlane, std::string const& vFilename); + void Init() {} + + corsika::process::EProcessReturn DoContinuous( + corsika::setup::Stack::ParticleType const& vParticle, + corsika::setup::Trajectory const& vTrajectory); + + corsika::units::si::LengthType MaxStepLength( + corsika::setup::Stack::ParticleType const&, + corsika::setup::Trajectory const& vTrajectory); + + private: + geometry::Plane const fObsPlane; + std::ofstream fOutputStream; + }; +} // namespace corsika::process::observation_plane + +#endif diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc new file mode 100644 index 0000000000000000000000000000000000000000..1df725253e8e0f125b7a11626359ccccf3471a65 --- /dev/null +++ b/Processes/ObservationPlane/testObservationPlane.cc @@ -0,0 +1,96 @@ +/* + * (c) Copyright 2019 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. + */ + +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file +#include <catch2/catch.hpp> + +#include <corsika/process/observation_plane/ObservationPlane.h> + +#include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> +#include <corsika/geometry/Vector.h> + +#include <corsika/particles/ParticleProperties.h> +#include <corsika/units/PhysicalUnits.h> + +using namespace corsika::units::si; +using namespace corsika::process::observation_plane; +using namespace corsika; +using namespace corsika::geometry; +using namespace corsika::particles; + +TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") { + + auto const& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); + + /* + Test with downward going 1_GeV neutrino, starting at 0,1_m,10m + + ObservationPlane has origin at 0,0,0 + */ + + Point const start(rootCS, {0_m, 1_m, 10_m}); + Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second, + -units::constants::c); + Line line(start, vec); + Trajectory<Line> track(line, 10_m / units::constants::c); + + // setup particle stack, and add primary particle + setup::Stack stack; + stack.Clear(); + { + auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { + return sqrt((Elab - m) * (Elab + m)); + }; + stack.AddParticle( + std::tuple<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, Point, + units::si::TimeType>{ + Code::NuMu, 1_GeV, + corsika::stack::MomentumVector( + rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::GetMass())}), + Point(rootCS, {1_m, 1_m, 10_m}), 0_ns}); + } + auto particle = stack.GetNextParticle(); + + SECTION("horizontal plane") { + + Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}), + Vector<dimensionless_d>(rootCS, {0., 0., 1.})); + ObservationPlane obs(obsPlane, "particles.dat"); + + obs.Init(); + [[maybe_unused]] const LengthType length = obs.MaxStepLength(particle, track); + [[maybe_unused]] const process::EProcessReturn ret = + obs.DoContinuous(particle, track); + + SECTION("steplength") { REQUIRE(length == 10_m); } + + /* + SECTION("horizontal plane") { + REQUIRE(true); // todo: we have to check content of output file... + + } + */ + } + + SECTION("inclined plane") { + Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}), + Vector<dimensionless_d>(rootCS, {1., 1., 0.5})); + ObservationPlane obs(obsPlane, "particles.dat"); + + obs.Init(); + [[maybe_unused]] const LengthType length = obs.MaxStepLength(particle, track); + [[maybe_unused]] const process::EProcessReturn ret = + obs.DoContinuous(particle, track); + + SECTION("steplength") { REQUIRE(length == 12_m); } + } +}