obs. plane

corsika/detail/modules/ObservationPlane.inl

@@ -16,49 +16,73 @@ namespace corsika::observation_plane {
                                      std::string const& filename, bool deleteOnHit)
       : plane_(obsPlane)
       , outputStream_(filename)
-      , deleteOnHit_(deleteOnHit) {
+      , deleteOnHit_(deleteOnHit)
+      , energy_ground_(0_GeV)
+      , count_ground_(0) {
     outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
   }
 
-  corsika::ProcessReturn ObservationPlane::doContinuous(
-      corsika::setup::Stack::ParticleType const& particle,
-      corsika::setup::Trajectory const& trajectory) {
+  ProcessReturn ObservationPlane::doContinuous(
+      corsika::setup::Stack::particle_type& particle,
+      corsika::setup::Trajectory& trajectory) {
     TimeType const timeOfIntersection =
-        (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-        trajectory.GetV0().dot(plane_.GetNormal());
+        (plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
+        trajectory.getVelocity().dot(plane_.getNormal());
 
-    if (timeOfIntersection < TimeType::zero()) { return corsika::ProcessReturn::Ok; }
+    if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
 
-    if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
-      return corsika::ProcessReturn::Ok;
+    if (plane_.isAbove(trajectory.getStartPoint()) ==
+        plane_.isAbove(trajectory.getPosition(1))) {
+      return ProcessReturn::Ok;
     }
 
-    outputStream_ << static_cast<int>(corsika::get_PDG(particle.GetPID())) << ' '
-                  << particle.GetEnergy() * (1 / 1_eV) << ' '
-                  << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
+    const auto energy = particle.getEnergy();
+    outputStream_ << static_cast<int>(corsika::get_PDG(particle.getPID())) << ' '
+                  << energy / 1_eV << ' '
+                  << (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
                   << std::endl;
 
     if (deleteOnHit_) {
-      return corsika::ProcessReturn::ParticleAbsorbed;
+      count_ground_++;
+      energy_ground_ += energy;
+      particle.erase();
+      return ProcessReturn::ParticleAbsorbed;
     } else {
-      return corsika::ProcessReturn::Ok;
+      return ProcessReturn::Ok;
     }
   }
 
-  corsika::LengthType ObservationPlane::MaxStepLength(
-      corsika::setup::Stack::ParticleType const&,
+  corsika::LengthType ObservationPlane::getMaxStepLength(
+      corsika::setup::Stack::particle_type const&,
       corsika::setup::Trajectory const& trajectory) {
 
     TimeType const timeOfIntersection =
-        (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-        trajectory.GetV0().dot(plane_.GetNormal());
+        (plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
+        trajectory.getVelocity().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;
+    auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
+    auto dist = (trajectory.getStartPoint() - pointOfIntersection).getNorm() * 1.0001;
+    CORSIKA_LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
+    return dist;
+  }
+
+  void ObservationPlane::showResults() const {
+    CORSIKA_LOG_INFO(
+        " ******************************\n"
+        " ObservationPlane: \n"
+        " energy in ground (GeV)     :  {}\n"
+        " no. of particles in ground :  {}\n"
+        " ******************************",
+        energy_ground_ / 1_GeV, count_ground_);
+  }
+
+  void ObservationPlane::reset() {
+    energy_ground_ = 0_GeV;
+    count_ground_ = 0;
   }
 
 } // namespace corsika::observation_plane

corsika/modules/ObservationPlane.hpp

@@ -15,6 +15,7 @@
 #include <corsika/setup/SetupTrajectory.hpp>
 
 #include <fstream>
+#include <string>
 
 namespace corsika::observation_plane {
 
@@ -29,17 +30,22 @@ namespace corsika::observation_plane {
     ObservationPlane(corsika::Plane const&, std::string const&, bool = true);
     void Init() {}
 
-    corsika::ProcessReturn doContinuous(
-        corsika::setup::Stack::ParticleType const& vParticle,
-        corsika::setup::Trajectory const& vTrajectory);
+    corsika::ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
+                                        corsika::setup::Trajectory& vTrajectory);
 
-    LengthType MaxStepLength(corsika::setup::Stack::ParticleType const&,
-                             corsika::setup::Trajectory const& vTrajectory);
+    LengthType getMaxStepLength(corsika::setup::Stack::particle_type const&,
+                                corsika::setup::Trajectory const& vTrajectory);
+
+    void showResults() const;
+    void reset();
+    HEPEnergyType getEnergyGround() const { return energy_ground_; }
 
   private:
-    corsika::Plane const plane_;
+    Plane const plane_;
     std::ofstream outputStream_;
     bool const deleteOnHit_;
+    HEPEnergyType energy_ground_;
+    unsigned int count_ground_;
   };
 } // namespace corsika::observation_plane
 

tests/modules/CMakeLists.txt

@@ -4,7 +4,7 @@ set (test_modules_sources
   testStackInspector.cpp
   #testTrackingLine.cpp
   #testExecTime.cpp
-  #testObservationPlane.cpp
+  testObservationPlane.cpp
   testQGSJetII.cpp
   #    testSwitchProcess.cpp  -> gone
   # testPythia8.cpp

tests/modules/testObservationPlane.cpp

@@ -22,7 +22,7 @@ using namespace corsika;
 
 TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
 
-  auto const& rootCS = RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+  auto const& rootCS = get_root_CoordinateSystem();
 
   /*
     Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
@@ -38,19 +38,17 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
 
   // setup particle stack, and add primary particle
   setup::Stack stack;
-  stack.Clear();
+  stack.clear();
   {
     auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
       return sqrt((Elab - m) * (Elab + m));
     };
-    stack.AddParticle(
-        std::tuple<Code, HEPEnergyType, corsika::MomentumVector, Point, TimeType>{
-            Code::NuMu, 1_GeV,
-            corsika::MomentumVector(rootCS,
-                                    {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::mass())}),
-            Point(rootCS, {1_m, 1_m, 10_m}), 0_ns});
+    stack.addParticle(std::make_tuple(
+        Code::NuMu, 1_GeV,
+        corsika::MomentumVector(rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::mass)}),
+        Point(rootCS, {1_m, 1_m, 10_m}), 0_ns));
   }
-  auto particle = stack.GetNextParticle();
+  auto particle = stack.getNextParticle();
 
   SECTION("horizontal plane") {
 
@@ -58,15 +56,15 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
                          Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
     ObservationPlane obs(obsPlane, "particles.dat", true);
 
-    const LengthType length = obs.MaxStepLength(particle, track);
+    const LengthType length = obs.getMaxStepLength(particle, track);
     const ProcessReturn ret = obs.doContinuous(particle, track);
 
-    REQUIRE(length / 10_m == Approx(1).margin(1e-4));
-    REQUIRE(ret == ProcessReturn::ParticleAbsorbed);
+    CHECK(length / 10_m == Approx(1).margin(1e-4));
+    CHECK(ret == ProcessReturn::ParticleAbsorbed);
 
     /*
     SECTION("horizontal plane") {
-      REQUIRE(true); // todo: we have to check content of output file...
+      CHECK(true); // todo: we have to check content of output file...
 
     }
     */
@@ -79,10 +77,10 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
                          Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
     ObservationPlane obs(obsPlane, "particles.dat", false);
 
-    const LengthType length = obs.MaxStepLength(particle, track);
+    const LengthType length = obs.getMaxStepLength(particle, track);
     const ProcessReturn ret = obs.doContinuous(particle, track);
 
-    REQUIRE(length / 10_m == Approx(1).margin(1e-4));
-    REQUIRE(ret == ProcessReturn::Ok);
+    CHECK(length / 10_m == Approx(1).margin(1e-4));
+    CHECK(ret == ProcessReturn::Ok);
   }
 }