more unit tests

corsika/detail/modules/proposal/ContinuousProcess.inl

@@ -48,9 +48,11 @@ namespace corsika::proposal {
     calc[std::make_pair(comp.getHash(), code)] = std::move(calculator);
   }
 
-  template <typename TEnvironment>
-  inline ContinuousProcess::ContinuousProcess(TEnvironment const& _env)
-      : ProposalProcessBase(_env) {}
+  template <typename TEnvironment, typename... TOutputArgs>
+  inline ContinuousProcess::ContinuousProcess(TEnvironment const& _env,
+                                              TOutputArgs&&... args)
+      : ProposalProcessBase(_env)
+      , TOutput(args) {}
 
   template <typename TParticle>
   inline void ContinuousProcess::scatter(TParticle& vP, HEPEnergyType const& loss,
@@ -103,11 +105,14 @@ namespace corsika::proposal {
                             vP.getEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
                         1_MeV;
     auto dE = vP.getEnergy() - final_energy;
-    energy_lost_ += dE;
 
     // if the particle has a charge take multiple scattering into account
     if (vP.getChargeNumber() != 0) scatter(vP, dE, dX);
     vP.setEnergy(final_energy); // on the stack, this is just kinetic energy, E-m
+
+    // also send to output
+    TOutput::write(vT, vP.getPID(), dE);
+
     return ProcessReturn::Ok;
   }
 
@@ -143,14 +148,4 @@ namespace corsika::proposal {
     return dist;
   }
 
-  inline void ContinuousProcess::showResults() const {
-    CORSIKA_LOG_DEBUG(
-        " ******************************\n"
-        " PROCESS::ContinuousProcess: \n"
-        " energy lost dE (GeV)      :  {}",
-        energy_lost_ / 1_GeV);
-  }
-
-  inline void ContinuousProcess::reset() { energy_lost_ = 0_GeV; }
-
 } // namespace corsika::proposal

corsika/modules/energy_loss/BetheBlochPDG.hpp

@@ -40,8 +40,8 @@ namespace corsika {
     using MeVgcm2 = decltype(1e6 * electronvolt / gram * square(1e-2 * meter));
 
   public:
-    template <typename... TArgs>
-    BetheBlochPDG(TArgs&&... args);
+    template <typename... TOutputArgs>
+    BetheBlochPDG(TOutputArgs&&... args);
 
     /**
      * Interface function of ContinuousProcess.

corsika/modules/proposal/ContinuousProcess.hpp

@@ -29,9 +29,11 @@ namespace corsika::proposal {
   //! use of interpolation tables which are runtime intensive calculation, but can be
   //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
   //!
+  template <typename TOutput = WriterOff>
   class ContinuousProcess
-      : public corsika::ContinuousProcess<proposal::ContinuousProcess>,
-        ProposalProcessBase {
+      : public corsika::ContinuousProcess<proposal::ContinuousProcess<TOutput>>,
+        public ProposalProcessBase,
+        public TOutput {
 
     struct Calculator {
       std::unique_ptr<PROPOSAL::Displacement> disp;
@@ -41,8 +43,6 @@ namespace corsika::proposal {
     std::unordered_map<calc_key_t, Calculator, hash>
         calc; //!< Stores the displacement and scattering calculators.
 
-    HEPEnergyType energy_lost_ = 0 * electronvolt;
-
     //!
     //! Build the displacement and scattering calculators and add it to calc.
     //!
@@ -53,8 +53,8 @@ namespace corsika::proposal {
     //! Produces the continuous loss calculator for leptons based on nuclear
     //! compositions and stochastic description limited by the particle cut.
     //!
-    template <typename TEnvironment>
-    ContinuousProcess(TEnvironment const&);
+    template <typename TEnvironment, typename... TOutputArgs>
+    ContinuousProcess(TEnvironment const&, TOutputArgs&&...);
 
     //!
     //! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
@@ -80,10 +80,6 @@ namespace corsika::proposal {
     //!
     template <typename TParticle, typename TTrack>
     LengthType getMaxStepLength(TParticle const&, TTrack const&);
-
-    void showResults() const;
-    void reset();
-    HEPEnergyType getEnergyLost() const { return energy_lost_; }
   };
 } // namespace corsika::proposal
 

tests/output/CMakeLists.txt

@@ -7,6 +7,7 @@ set (test_output_sources
   testWriterTracks.cpp
   testWriterEnergyLoss.cpp
   testWriterLongitudinal.cpp
+  testWriterOff.cpp
   )
 
 CORSIKA_ADD_TEST (testOutput SOURCES ${test_output_sources})

tests/output/testWriterEnergyLoss.cpp

@@ -49,8 +49,6 @@ class TestEnergyLoss : public corsika::EnergyLossWriter<> {
 public:
   TestEnergyLoss(corsika::ShowerAxis const& axis)
       : EnergyLossWriter(axis) {}
-
-  YAML::Node getConfig() const { return YAML::Node(); }
 };
 
 TEST_CASE("EnergyLossWriter") {
@@ -90,21 +88,37 @@ TEST_CASE("EnergyLossWriter") {
 
   // generate straight simple track
   CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
-  Point r0(rootCS, {0_m, 0_m, 5_km});
+  Point r0(rootCS, {0_m, 0_m, 7_km});
   SpeedType const V0 = constants::c;
   VelocityVector v0(rootCS, {V0, 0_m / second, 0_m / second});
   Line const line(r0, v0);
-  auto const time = 10_ns;
+  auto const time = 1000_ns;
   StraightTrajectory track(line, time);
   // test write
   test.write(track, Code::Proton, 100_GeV);
 
   // incompatible binning
-  CHECK_THROWS(
-      test.write(100_g / square(1_cm), 120_g / square(1_cm), Code::Photon, 100_GeV));
+  CHECK_THROWS(test.write(100_g / square(1_cm), // extra line break by purpose
+                          120_g / square(1_cm), Code::Photon, 100_GeV));
+  test.write(100000_g / square(1_cm), 100010_g / square(1_cm), Code::Photon,
+             100_GeV); // this doesn't throw, but it skips
 
   test.endOfShower(0);
   test.endOfLibrary();
 
   CHECK(boost::filesystem::exists("./output_dir_eloss/dEdX.parquet"));
+
+  auto const config = test.getConfig();
+  CHECK(config["type"].as<std::string>() == "EnergyLoss");
+  CHECK(config["units"]["energy"].as<std::string>() == "GeV");
+  CHECK(config["units"]["grammage"].as<std::string>() == "g/cm^2");
+  CHECK(config["bin-size"].as<double>() == 10.);
+  CHECK(config["nbins"].as<int>() == 200);
+  CHECK(config["grammage_threshold"].as<double>() == Approx(0.0001));
+
+  auto const summary = test.getSummary();
+  CHECK(summary["sum_dEdX"].as<double>() == 300);
+  // makes not yet sense:
+  // CHECK(summary["Xmax"].as<double>() == 200);
+  // CHECK(summary["dEdXmax"].as<double>() == 200);
 }

tests/output/testWriterLongitudinal.cpp

@@ -23,6 +23,8 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/core/Logging.hpp>
 
+#include <string>
+
 using namespace corsika;
 
 const auto density = 1_kg / (1_m * 1_m * 1_m);
@@ -50,8 +52,6 @@ class TestLongitudinal : public corsika::LongitudinalWriter<> {
 public:
   TestLongitudinal(corsika::ShowerAxis const& axis)
       : LongitudinalWriter(axis) {}
-
-  YAML::Node getConfig() const { return YAML::Node(); }
 };
 
 TEST_CASE("LongitudinalWriter") {
@@ -85,17 +85,42 @@ TEST_CASE("LongitudinalWriter") {
 
   // generate straight simple track
   CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
-  Point r0(rootCS, {0_km, 0_m, 5_m});
+  Point r0(rootCS, {0_km, 0_m, 7_m});
   SpeedType const V0 = constants::c;
   VelocityVector v0(rootCS, {V0, 0_m / second, 0_m / second});
   Line const line(r0, v0);
-  auto const time = 10_ns;
+  auto const time = 1000_ns;
   StraightTrajectory track(line, time);
   // test write
   test.write(track, Code::Proton, 1.0);
+  test.write(track, Code::Photon, 1.0);
+  test.write(track, Code::Electron, 1.0);
+  test.write(track, Code::Positron, 1.0);
+  test.write(track, Code::MuPlus, 1.0);
+  test.write(track, Code::MuMinus, 1.0);
+
+  test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::PiPlus, 1.0);
+  test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::Electron, 1.0);
+  test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::Positron, 1.0);
+  test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::Photon, 1.0);
+  test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::MuPlus, 1.0);
+  test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::MuMinus, 1.0);
+
+  // wrong binning
+  CHECK_THROWS(test.write(10_g / square(1_cm), 10.1_g / square(1_cm), Code::PiPlus, 1.0));
+  test.write(100000_g / square(1_cm), 100010_g / square(1_cm), Code::PiPlus,
+             1.0); // this doesn't throw, it just skips
 
   test.endOfShower(0);
   test.endOfLibrary();
 
   CHECK(boost::filesystem::exists("./output_dir_long/profile.parquet"));
+
+  auto const config = test.getConfig();
+  CHECK(config["type"].as<std::string>() == "LongitudinalProfile");
+  CHECK(config["units"]["grammage"].as<std::string>() == "g/cm^2");
+  CHECK(config["bin-size"].as<double>() == 10.);
+  CHECK(config["nbins"].as<int>() == 200);
+
+  auto const summary = test.getSummary(); // nothing to check yet
 }

tests/output/testWriterObservationPlane.cpp

@@ -22,7 +22,11 @@ struct TestWriterPlane : public ParticleWriterParquet {
   YAML::Node getConfig() const { return YAML::Node(); }
 
   void checkWrite() {
-    ParticleWriterParquet::write(Code::Unknown, 1_eV, 2_m, 3_m, 0_m, 1.0);
+    ParticleWriterParquet::write(Code::Unknown, 1_GeV, 2_m, 3_m, 0_m, 1.0);
+    ParticleWriterParquet::write(Code::Proton, 1_GeV, 2_m, 3_m, 0_m, 1.0);
+    ParticleWriterParquet::write(Code::MuPlus, 1_GeV, 2_m, 3_m, 0_m, 1.0);
+    ParticleWriterParquet::write(Code::MuMinus, 1_GeV, 2_m, 3_m, 0_m, 1.0);
+    ParticleWriterParquet::write(Code::Photon, 1_GeV, 2_m, 3_m, 0_m, 1.0);
   }
 };
 
@@ -41,10 +45,21 @@ TEST_CASE("ObservationPlaneWriterParquet") {
     TestWriterPlane test;
     test.startOfLibrary("./output_dir");
     test.startOfShower(0);
+
+    // write a few particles
     test.checkWrite();
+
     test.endOfShower(0);
     test.endOfLibrary();
 
     CHECK(boost::filesystem::exists("./output_dir/particles.parquet"));
+
+    auto const summary = test.getSummary();
+
+    CHECK(summary["Eground"].as<double>() == Approx(5));
+    CHECK(summary["hadrons"].as<int>() == Approx(1));
+    CHECK(summary["muons"].as<int>() == Approx(2));
+    CHECK(summary["em"].as<int>() == Approx(1));
+    CHECK(summary["others"].as<int>() == Approx(1));
   }
 }

tests/output/testWriterOff.cpp

+/*
+ * (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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <corsika/modules/writers/WriterOff.hpp>
+
+#include <corsika/framework/core/Logging.hpp>
+
+using namespace corsika;
+
+/*
+class TestEnergyLoss : public corsika::EnergyLossWriter<> {
+public:
+  TestEnergyLoss(corsika::ShowerAxis const& axis)
+      : EnergyLossWriter(axis) {}
+
+  YAML::Node getConfig() const { return YAML::Node(); }
+};
+*/
+
+TEST_CASE("WriterOff") {
+
+  logging::set_level(logging::level::info);
+
+  WriterOff test("irrelevant", 3);
+  WriterOff test2();
+
+  test.startOfLibrary("./output_dir_eloss");
+  test.startOfShower(0);
+
+  test.endOfShower(0);
+
+  test.endOfLibrary();
+
+  auto const config = test.getConfig();
+
+  auto const summary = test.getSummary();
+}