From feed1776fcc10eeba8d816224aad75c2a6c0ba22 Mon Sep 17 00:00:00 2001
From: Fan <fan_hu@pku.edu.cn>
Date: Sat, 27 Nov 2021 20:01:29 +0800
Subject: [PATCH] add cubic
---
corsika/detail/framework/geometry/Cubic.inl | 20 ++
corsika/detail/modules/ObservationCubic.inl | 148 ++++++++++++
.../modules/tracking/TrackingStraight.inl | 215 +++++++++++-------
.../writers/ObservationCubicWriterParquet.inl | 65 ++++++
corsika/framework/geometry/Cubic.hpp | 45 ++++
corsika/modules/ObservationCubic.hpp | 59 +++++
corsika/modules/tracking/TrackingStraight.hpp | 4 +
.../writers/ObservationCubicWriterParquet.hpp | 61 +++++
8 files changed, 537 insertions(+), 80 deletions(-)
create mode 100644 corsika/detail/framework/geometry/Cubic.inl
create mode 100644 corsika/detail/modules/ObservationCubic.inl
create mode 100644 corsika/detail/modules/writers/ObservationCubicWriterParquet.inl
create mode 100644 corsika/framework/geometry/Cubic.hpp
create mode 100644 corsika/modules/ObservationCubic.hpp
create mode 100644 corsika/modules/writers/ObservationCubicWriterParquet.hpp
diff --git a/corsika/detail/framework/geometry/Cubic.inl b/corsika/detail/framework/geometry/Cubic.inl
new file mode 100644
index 000000000..be944049e
--- /dev/null
+++ b/corsika/detail/framework/geometry/Cubic.inl
@@ -0,0 +1,20 @@
+
+
+namespace corsika {
+inline bool Cubic::contains(Point const &p) const {
+ if ((abs(p.getX(cs_)) < x_) && (abs(p.getY(cs_)) < y_) &&
+ (abs(p.getZ(cs_)) < z_))
+ return true;
+ else
+ return false;
+}
+
+inline std::string Cubic::asString() const {
+ std::ostringstream txt;
+ txt << "center=" << center_ << ", x-axis=" << DirectionVector{cs_, {1, 0, 0}}
+ << ", y-axis: " << DirectionVector{cs_, {0, 1, 0}}
+ << ", z-axis: " << DirectionVector{cs_, {0, 0, 1}};
+ return txt.str();
+}
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/ObservationCubic.inl b/corsika/detail/modules/ObservationCubic.inl
new file mode 100644
index 000000000..668150624
--- /dev/null
+++ b/corsika/detail/modules/ObservationCubic.inl
@@ -0,0 +1,148 @@
+/*
+ * (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.
+ */
+
+namespace corsika {
+
+template <typename TTracking, typename TOutput>
+ObservationCubic<TTracking, TOutput>::ObservationCubic(
+ Point const ¢er, LengthType const x, LengthType const y,
+ LengthType const z, bool deleteOnHit)
+ : Cubic(center, x, y, z), deleteOnHit_(deleteOnHit),
+ energy_(0_GeV), count_(0) {}
+
+template <typename TTracking, typename TOutput>
+template <typename TParticle, typename TTrajectory>
+inline ProcessReturn ObservationCubic<TTracking, TOutput>::doContinuous(
+ TParticle &particle, TTrajectory &, bool const stepLimit) {
+ /*
+ The current step did not yet reach the ObservationCubic, do nothing now and
+ wait:
+ */
+ if (!stepLimit) {
+ // @todo this is actually needed to fix small instabilities of the leap-frog
+ // tracking: Note, this is NOT a general solution and should be clearly
+ // revised with a more robust tracking. #ifdef DEBUG
+ if (deleteOnHit_) {
+ // since this is basically a bug, it cannot be tested LCOV_EXCL_START
+ LengthType const check = 1_m; // TODO, do I need to check?
+ if (check < 0_m) {
+ CORSIKA_LOG_WARN("PARTICLE AVOIDED ObservationCubic {}", check);
+ CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
+ } else
+ return ProcessReturn::Ok;
+ // LCOV_EXCL_STOP
+ } else
+ // #endif
+ return ProcessReturn::Ok;
+ }
+
+ HEPEnergyType const energy = particle.getEnergy();
+ Point const pointOfIntersection = particle.getPosition();
+ DirectionVector const dirction = particle.getDirection();
+
+ // add our particles to the output file stream
+ this->write(particle.getPID(), energy, pointOfIntersection.getX(cs_),
+ pointOfIntersection.getY(cs_), pointOfIntersection.getZ(cs_),
+ dirction.getX(cs_), dirction.getY(cs_), dirction.getZ(cs_),
+ particle.getTime());
+
+ CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
+
+ if (deleteOnHit_) {
+ count_++;
+ energy_ += energy;
+ return ProcessReturn::ParticleAbsorbed;
+ } else {
+ return ProcessReturn::Ok;
+ }
+}
+
+template <typename TTracking, typename TOutput>
+template <typename TParticle, typename TTrajectory>
+inline LengthType ObservationCubic<TTracking, TOutput>::getMaxStepLength(
+ TParticle const &particle, TTrajectory const &trajectory) {
+
+ CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, cubic={}",
+ particle.asString(), particle.getPosition(),
+ particle.getDirection(), cubic_.asString());
+
+ auto const intersection = TTracking::intersect(particle, *this);
+
+ TimeType const timeOfIntersection = intersection.getEntry();
+ CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+ if (timeOfIntersection > trajectory.getDuration()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+ double const fractionOfIntersection =
+ timeOfIntersection / trajectory.getDuration();
+ CORSIKA_LOG_TRACE("ObservationCubic: getMaxStepLength dist={} m, pos={}",
+ trajectory.getLength(fractionOfIntersection) / 1_m,
+ trajectory.getPosition(fractionOfIntersection));
+ return trajectory.getLength(fractionOfIntersection);
+}
+
+template <typename TTracking, typename TOutput>
+inline void ObservationCubic<TTracking, TOutput>::showResults() const {
+ CORSIKA_LOG_INFO(" ******************************\n"
+ " ObservationCubic: \n"
+ " energy at cubic (GeV) : {}\n"
+ " no. of particles at cubic : {}\n"
+ " ******************************",
+ energy_ / 1_GeV, count_);
+}
+
+template <typename TTracking, typename TOutput>
+inline YAML::Node ObservationCubic<TTracking, TOutput>::getConfig() const {
+ using namespace units::si;
+
+ // construct the top-level node
+ YAML::Node node;
+
+ // basic info
+ node["type"] = "ObservationCubic";
+ node["units"] = "m"; // add default units for values
+
+ // save each component in its native coordinate system
+ auto const root_cs = get_root_CoordinateSystem();
+ node["center"].push_back(center_.getX(root_cs) / 1_m);
+ node["center"].push_back(center_.getY(root_cs) / 1_m);
+ node["center"].push_back(center_.getZ(root_cs) / 1_m);
+
+ // the x-axis vector
+ DirectionVector const x_axis = DirectionVector{cs_, {1, 0, 0}};
+ node["x-axis"].push_back(x_axis.getX(root_cs));
+ node["x-axis"].push_back(x_axis.getY(root_cs));
+ node["x-axis"].push_back(x_axis.getZ(root_cs));
+
+ // the y-axis vector
+ DirectionVector const y_axis = DirectionVector{cs_, {0, 1, 0}};
+ node["y-axis"].push_back(y_axis.getX(root_cs));
+ node["y-axis"].push_back(y_axis.getY(root_cs));
+ node["y-axis"].push_back(y_axis.getZ(root_cs));
+
+ // the x-axis vector
+ DirectionVector const z_axis = DirectionVector{cs_, {0, 0, 1}};
+ node["z-axis"].push_back(z_axis.getX(root_cs));
+ node["z-axis"].push_back(z_axis.getY(root_cs));
+ node["z-axis"].push_back(z_axis.getZ(root_cs));
+
+ node["delete_on_hit"] = deleteOnHit_;
+
+ return node;
+}
+
+template <typename TTracking, typename TOutput>
+inline void ObservationCubic<TTracking, TOutput>::reset() {
+ energy_= 0_GeV;
+ count_= 0;
+}
+
+} // namespace corsika
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 3bdbb73df..a88aed654 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -8,103 +8,158 @@
#pragma once
+#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp>
-#include <corsika/framework/geometry/Intersections.hpp>
-#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/modules/tracking/Intersect.hpp>
+#include <cmath>
#include <type_traits>
#include <utility>
-#include <cmath>
namespace corsika::tracking_line {
- template <typename TParticle>
- inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
- if (particle.getMomentum().getNorm() == 0_GeV) {
- return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV);
- } // charge of the particle
- DirectionVector const dir = particle.getDirection();
- return std::make_tuple(particle.getPosition() + dir * steplength, dir.normalized());
+template <typename TParticle>
+inline auto Tracking::makeStep(TParticle const &particle,
+ LengthType steplength) {
+ if (particle.getMomentum().getNorm() == 0_GeV) {
+ return std::make_tuple(particle.getPosition(),
+ particle.getMomentum() / 1_GeV);
+ } // charge of the particle
+ DirectionVector const dir = particle.getDirection();
+ return std::make_tuple(particle.getPosition() + dir * steplength,
+ dir.normalized());
+}
+
+template <typename TParticle>
+inline auto Tracking::getTrack(TParticle const &particle) {
+ VelocityVector const initialVelocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+
+ auto const &initialPosition = particle.getPosition();
+ CORSIKA_LOG_DEBUG("TrackingStraight pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {}",
+ particle.getPID(), particle.getEnergy() / 1_GeV,
+ initialPosition.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV,
+ initialVelocity.getComponents());
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = nextIntersect(particle);
+
+ return std::make_tuple(
+ StraightTrajectory(Line(initialPosition, initialVelocity),
+ minTime), // trajectory
+ minNode); // next volume node
+}
+
+template <typename TParticle>
+inline Intersections Tracking::intersect(TParticle const &particle,
+ Sphere const &sphere) {
+ auto const &position = particle.getPosition();
+ auto const delta = position - sphere.getCenter();
+ auto const velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const vSqNorm = velocity.getSquaredNorm();
+ auto const R = sphere.getRadius();
+
+ auto const vDotDelta = velocity.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+
+ CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
+ return Intersections((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
}
-
- template <typename TParticle>
- inline auto Tracking::getTrack(TParticle const& particle) {
- VelocityVector const initialVelocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
-
- auto const& initialPosition = particle.getPosition();
- CORSIKA_LOG_DEBUG(
- "TrackingStraight pid: {}"
- " , E = {} GeV \n"
- "\tTracking pos: {} \n"
- "\tTracking p: {} GeV \n"
- "\tTracking v: {}",
- particle.getPID(), particle.getEnergy() / 1_GeV, initialPosition.getCoordinates(),
- particle.getMomentum().getComponents() / 1_GeV, initialVelocity.getComponents());
-
- // traverse the environment volume tree and find next
- // intersection
- auto [minTime, minNode] = nextIntersect(particle);
-
- return std::make_tuple(StraightTrajectory(Line(initialPosition, initialVelocity),
- minTime), // trajectory
- minNode); // next volume node
- }
-
- template <typename TParticle>
- inline Intersections Tracking::intersect(TParticle const& particle,
- Sphere const& sphere) {
- auto const& position = particle.getPosition();
- auto const delta = position - sphere.getCenter();
- auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const vSqNorm = velocity.getSquaredNorm();
- auto const R = sphere.getRadius();
-
- auto const vDotDelta = velocity.dot(delta);
- auto const discriminant =
- vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
-
- if (discriminant.magnitude() > 0) {
- auto const sqDisc = sqrt(discriminant);
- auto const invDenom = 1 / vSqNorm;
-
- CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
- return Intersections((-vDotDelta - sqDisc) * invDenom,
- (-vDotDelta + sqDisc) * invDenom);
- }
+ return Intersections();
+}
+
+template <typename TParticle>
+inline Intersections Tracking::intersect(TParticle const &particle,
+ Cubic const &cubic) {
+ Point const &position = particle.getPosition();
+ VelocityVector const velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ CoordinateSystemPtr const &cs = cubic.getCoordinateSystem();
+ LengthType x0 = position.getX(cs);
+ LengthType y0 = position.getY(cs);
+ LengthType z0 = position.getZ(cs);
+ SpeedType vx = velocity.getX(cs);
+ SpeedType vy = velocity.getY(cs);
+ SpeedType vz = velocity.getZ(cs);
+
+ auto get_intersect_min_max = [](LengthType x0, SpeedType v0, LengthType dx) {
+ auto t1 = (x0 - dx) / v0;
+ auto t2 = (x0 + dx) / v0;
+ if (t1 > t2)
+ return std::make_pair(t1, t2);
+ else
+ return std::make_pair(t2, t1);
+ };
+
+ auto [tx_max, tx_min] = get_intersect_min_max(x0, vx, cubic.getX());
+ auto [ty_max, ty_min] = get_intersect_min_max(y0, vy, cubic.getY());
+ auto [tz_max, tz_min] = get_intersect_min_max(z0, vz, cubic.getZ());
+
+ TimeType t_exit = std::min(std::min(tx_max, ty_max), tz_max);
+ TimeType t_enter = std::min(std::min(tx_min, ty_min), tz_min);
+
+ if ((t_exit > t_enter)) {
+ if (t_enter < 0_s && t_exit > 0_s)
+ CORSIKA_LOG_DEBUG("numericallyInside={}", cubic.contains(position));
+ else if (t_enter < 0_s && t_exit < 0_s)
+ CORSIKA_LOG_DEBUG("oppisite direction");
+ return Intersections(std::move(t_enter), std::move(t_exit));
+ } else
return Intersections();
+}
+
+template <typename TParticle, typename TBaseNodeType>
+inline Intersections Tracking::intersect(TParticle const &particle,
+ TBaseNodeType const &volumeNode) {
+ if (Sphere const *sphere =
+ dynamic_cast<Sphere const *>(&volumeNode.getVolume());
+ sphere) {
+ return Tracking::intersect<TParticle>(particle, *sphere);
+ } else if (Cubic const *cubic =
+ dynamic_cast<Cubic const *>(&volumeNode.getVolume());
+ cubic) {
+ return Tracking::intersect<TParticle>(particle, *cubic);
+ } else {
+ throw std::runtime_error("The Volume type provided is not supported in "
+ "Intersect(particle, node)");
}
+}
- template <typename TParticle, typename TBaseNodeType>
- inline Intersections Tracking::intersect(TParticle const& particle,
- TBaseNodeType const& volumeNode) {
- Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
- if (sphere) { return Tracking::intersect<TParticle>(particle, *sphere); }
- throw std::runtime_error(
- "The Volume type provided is not supported in Intersect(particle, node)");
- }
-
- template <typename TParticle>
- inline Intersections Tracking::intersect(TParticle const& particle,
- Plane const& plane) {
- auto const delta = plane.getCenter() - particle.getPosition();
- auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const n = plane.getNormal();
- auto const n_dot_v = n.dot(velocity);
+template <typename TParticle>
+inline Intersections Tracking::intersect(TParticle const &particle,
+ Plane const &plane) {
+ auto const delta = plane.getCenter() - particle.getPosition();
+ auto const velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const n = plane.getNormal();
+ auto const n_dot_v = n.dot(velocity);
- CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
- particle.getMomentum());
+ CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
+ particle.getMomentum());
- if (n_dot_v.magnitude() == 0)
- return Intersections();
- else
- return Intersections(n.dot(delta) / n_dot_v);
- }
+ if (n_dot_v.magnitude() == 0)
+ return Intersections();
+ else
+ return Intersections(n.dot(delta) / n_dot_v);
+}
} // namespace corsika::tracking_line
diff --git a/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl b/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl
new file mode 100644
index 000000000..bd85fd0df
--- /dev/null
+++ b/corsika/detail/modules/writers/ObservationCubicWriterParquet.inl
@@ -0,0 +1,65 @@
+/*
+ * (c) Copyright 2021 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
+
+namespace corsika {
+
+ inline ObservationCubicWriterParquet::ObservationCubicWriterParquet()
+ : output_() {}
+
+ inline void ObservationCubicWriterParquet::startOfLibrary(
+ boost::filesystem::path const& directory) {
+
+ // setup the streamer
+ output_.initStreamer((directory / "particles.parquet").string());
+
+ // enable compression with the default level
+ output_.enableCompression();
+
+ // build the schema
+ output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
+ parquet::ConvertedType::INT_32);
+ output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("nx", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("ny", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("nz", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+ output_.addField("t", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
+ parquet::ConvertedType::NONE);
+
+ // and build the streamer
+ output_.buildStreamer();
+ }
+
+ inline void ObservationCubicWriterParquet::endOfShower() { ++shower_; }
+
+ inline void ObservationCubicWriterParquet::endOfLibrary() { output_.closeStreamer(); }
+
+ inline void ObservationCubicWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
+ LengthType const& x, LengthType const& y, LengthType const& z,
+ double nx, double ny, double nz,
+ TimeType const& t) {
+ // write the next row - we must write `shower_` first.
+ *(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid))
+ << static_cast<float>(energy / 1_GeV)
+ << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m) << static_cast<float>(z / 1_m)
+ << static_cast<float>(nx) << static_cast<float>(ny) << static_cast<float>(nz)
+ << static_cast<float>(t / 1_ns) << parquet::EndRow;
+ }
+
+} // namespace corsika
diff --git a/corsika/framework/geometry/Cubic.hpp b/corsika/framework/geometry/Cubic.hpp
new file mode 100644
index 000000000..114965432
--- /dev/null
+++ b/corsika/framework/geometry/Cubic.hpp
@@ -0,0 +1,45 @@
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/IVolume.hpp>
+
+namespace corsika {
+ class Cubic : public IVolume {
+
+ public:
+ // a CoordinateSystemPtr to specify the orintation of coordinate
+ Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const x, LengthType const y, LengthType const z)
+ : center_(center)
+ , cs_(make_translation(cs, center.getCoordinates(center.getCoordinateSystem())))
+ , x_(x), y_(y), z_(z) {}
+
+ Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const side)
+ : center_(center)
+ , cs_(make_translation(cs, center.getCoordinates(center.getCoordinateSystem())))
+ , x_(side/2), y_(side/2), z_(side/2) {}
+
+
+ //! returns true if the Point p is within the sphere
+ bool contains(Point const& p) const override;
+
+ Point const& getCenter() const { return center_; };
+ CoordinateSystemPtr const getCoordinateSystem() const { return cs_; }
+ LengthType const getX() const { return x_; }
+ LengthType const getY() const { return y_; }
+ LengthType const getZ() const { return z_; }
+
+ std::string asString() const;
+
+ protected:
+ Point center_;
+ CoordinateSystemPtr cs_; // local coordinate system with center_ in coordinate (0, 0, 0) and user defined orientation
+ LengthType x_;
+ LengthType y_;
+ LengthType z_;
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/Cubic.inl>
diff --git a/corsika/modules/ObservationCubic.hpp b/corsika/modules/ObservationCubic.hpp
new file mode 100644
index 000000000..6095a0fd1
--- /dev/null
+++ b/corsika/modules/ObservationCubic.hpp
@@ -0,0 +1,59 @@
+#pragma once
+
+#include <corsika/framework/geometry/Cubic.hpp>
+
+#include <corsika/modules/writers/ObservationCubicWriterParquet.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+
+
+namespace corsika {
+
+ /**
+ @ingroup Modules
+ @{
+
+ The ObservationCubic 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.
+
+ **Note/Limitation:** as discussed in
+ https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397
+ you cannot put two ObservationCubics exactly on top of each
+ other. Even if one of them is "permeable". You have to put a
+ small gap in between the two plane in such a scenario, or develop
+ another more specialized output class.
+ */
+ template <typename TTracking, typename TOutputWriter = ObservationCubicWriterParquet>
+ class ObservationCubic
+ : public Cubic,
+ public ContinuousProcess<ObservationCubic<TTracking, TOutputWriter>>,
+ public TOutputWriter {
+
+ public:
+ ObservationCubic(Point const& center,
+ LengthType const x, LengthType const y, LengthType const z, bool = true);
+
+ template <typename TParticle, typename TTrajectory>
+ ProcessReturn doContinuous(TParticle& vParticle, TTrajectory& vTrajectory,
+ bool const stepLimit);
+
+ template <typename TParticle, typename TTrajectory>
+ LengthType getMaxStepLength(TParticle const&, TTrajectory const& vTrajectory);
+
+ void showResults() const;
+ void reset();
+ HEPEnergyType getEnergy() const { return energy_; }
+ YAML::Node getConfig() const;
+
+ private:
+ bool const deleteOnHit_;
+ HEPEnergyType energy_;
+ unsigned int count_;
+
+ };
+ //! @}
+} // namespace corsika
+
+
+
+#include <corsika/details/modules/ObservationCubic.inl>
\ No newline at end of file
diff --git a/corsika/modules/tracking/TrackingStraight.hpp b/corsika/modules/tracking/TrackingStraight.hpp
index 3ebc18b37..31b2fad4f 100644
--- a/corsika/modules/tracking/TrackingStraight.hpp
+++ b/corsika/modules/tracking/TrackingStraight.hpp
@@ -12,6 +12,7 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Cubic.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
@@ -50,6 +51,9 @@ namespace corsika::tracking_line {
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Sphere const& sphere);
+ template <typename TParticle>
+ static Intersections intersect(TParticle const& particle, Cubic const& cubic);
+
//! find intersection of Volume node with Track of particle
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(TParticle const& particle, TBaseNodeType const& node);
diff --git a/corsika/modules/writers/ObservationCubicWriterParquet.hpp b/corsika/modules/writers/ObservationCubicWriterParquet.hpp
new file mode 100644
index 000000000..86514abff
--- /dev/null
+++ b/corsika/modules/writers/ObservationCubicWriterParquet.hpp
@@ -0,0 +1,61 @@
+/*
+ * (c) Copyright 2021 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/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+ class ObservationCubicWriterParquet : public BaseOutput {
+
+ ParquetStreamer output_; ///< The primary output file.
+
+ public:
+ /**
+ * Construct an ObservationPlane.
+ *
+ * @param name The name of this output.
+ */
+ ObservationCubicWriterParquet();
+
+ /**
+ * Called at the start of each library.
+ */
+ void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+ /**
+ * Called at the end of each shower.
+ */
+ void endOfShower() final override;
+
+ /**
+ * Called at the end of each library.
+ *
+ * This must also increment the run number since we override
+ * the default behaviour of BaseOutput.
+ */
+ void endOfLibrary() final override;
+
+ protected:
+ /**
+ * Write a particle to the file.
+ */
+ void write(Code const& pid, HEPEnergyType const& energy,
+ LengthType const& x, LengthType const& y, LengthType const& z,
+ double nx, double ny, double nz,
+ TimeType const& t);
+
+ }; // class ObservationCubicWriterParquet
+
+} // namespace corsika
+
+#include <corsika/details/modules/writers/ObservationCubicWriterParquet.inl>
--
GitLab