diff --git a/corsika/detail/framework/process/InteractionHistogram.inl b/corsika/detail/framework/process/InteractionHistogram.inl
index ecb582368df13fc4efdb93347878148773b578f6..e6cfae869ea9ba6a9c8adc36d4e817a872822609 100644
--- a/corsika/detail/framework/process/InteractionHistogram.inl
+++ b/corsika/detail/framework/process/InteractionHistogram.inl
@@ -42,14 +42,6 @@ namespace corsika {
}
}
- void InteractionHistogram::saveLab(std::string const& filename, SaveMode mode) const {
- corsika::save_hist(inthist_lab_, filename, mode);
- }
-
- void InteractionHistogram::saveCMS(std::string const& filename, SaveMode mode) const {
- corsika::save_hist(inthist_cms_, filename, mode);
- }
-
InteractionHistogram& InteractionHistogram::operator+=(
InteractionHistogram const& other) {
inthist_lab_ += other.inthist_lab_;
@@ -64,4 +56,4 @@ namespace corsika {
return other;
}
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/SaveBoostHistogram.inl b/corsika/detail/framework/utility/SaveBoostHistogram.inl
index 94e5e1bbc8a383f99ef6a70776e622e6f0f2ad28..817fd6714a0a00ba492c31e9a45fa7c89ac2b5ff 100644
--- a/corsika/detail/framework/utility/SaveBoostHistogram.inl
+++ b/corsika/detail/framework/utility/SaveBoostHistogram.inl
@@ -11,6 +11,7 @@
#include <cnpy.hpp>
#include <boost/histogram.hpp>
+#include <boost/filesystem.hpp> // can be changed to std::filesystem if compiler supports it
#include <functional>
#include <memory>
@@ -23,11 +24,19 @@ namespace corsika {
template <class Axes, class Storage>
inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
- std::string const& filename, SaveMode mode) {
- unsigned const rank = h.rank();
+ std::string const& filename, bool overwrite) {
+ if (boost::filesystem::status(filename).type() !=
+ boost::filesystem::file_type::file_not_found) {
+ if (overwrite) {
+ boost::filesystem::remove(filename);
+ } else {
+ using namespace std::literals;
+ throw std::runtime_error(
+ ("save_hist(): "s + filename + " already exists"s).c_str());
+ }
+ }
- // append vs. overwrite
- const std::string mode_str = (mode == SaveMode::append ? "a" : "w");
+ unsigned const rank = h.rank();
std::vector<size_t> axes_dims;
axes_dims.reserve(rank);
@@ -58,7 +67,7 @@ namespace corsika {
ax_edges.push_back(ax.bin(ax.size() - 1).upper());
cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
- ax_edges.data(), {ax_edges.size()}, mode_str);
+ ax_edges.data(), {ax_edges.size()}, "a");
} else {
axis_types.push_back('d');
std::vector<int64_t> bins; // we assume that discrete axes have integer bins
@@ -67,14 +76,13 @@ namespace corsika {
for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
- {bins.size()}, mode_str);
+ {bins.size()}, "a");
}
}
- cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank},
- mode_str);
- cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, mode_str);
- cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, mode_str);
+ cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank}, "a");
+ cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, "a");
+ cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, "a");
auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
unsigned{1}, std::multiplies<>());
@@ -98,9 +106,9 @@ namespace corsika {
temp[p] = *x;
}
- cnpy::npz_save(filename, "data", temp.get(), axes_dims, mode_str);
+ cnpy::npz_save(filename, "data", temp.get(), axes_dims, "a");
// In Python this array can directly be assigned to a histogram view if that
// histogram has its axes correspondingly: hist.view(flow=True)[:] = file['data']
} // namespace corsika
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/framework/process/InteractionHistogram.hpp b/corsika/framework/process/InteractionHistogram.hpp
index 667ba3b9820d29bfdbe32b7d7a5643a4699f1a95..ab0fee8fd299a36d3f150db1e8d0f1213a615315 100644
--- a/corsika/framework/process/InteractionHistogram.hpp
+++ b/corsika/framework/process/InteractionHistogram.hpp
@@ -43,9 +43,6 @@ namespace corsika {
hist_type const& CMSHist() const { return inthist_cms_; }
hist_type const& labHist() const { return inthist_lab_; }
- void saveLab(std::string const& filename, SaveMode mode = SaveMode::append) const;
- void saveCMS(std::string const& filename, SaveMode mode = SaveMode::append) const;
-
InteractionHistogram& operator+=(InteractionHistogram const& other);
InteractionHistogram operator+(InteractionHistogram other) const;
};
diff --git a/corsika/framework/utility/SaveBoostHistogram.hpp b/corsika/framework/utility/SaveBoostHistogram.hpp
index 13a079410e7c4735e88eb0ec56f5346f96ad2407..90a5621bd54715964a117b041380d20470cf4af3 100644
--- a/corsika/framework/utility/SaveBoostHistogram.hpp
+++ b/corsika/framework/utility/SaveBoostHistogram.hpp
@@ -12,8 +12,6 @@
namespace corsika {
- enum class SaveMode { overwrite, append };
-
/**
* This functions saves a boost::histogram into a numpy file. Only rather basic axis
* types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
@@ -21,10 +19,13 @@ namespace corsika {
*
* Note that this function makes a temporary, dense copy of the histogram, which could
* be an issue for huge sizes (e.g. for high dimensions)
+ *
+ * @param overwrite silently overwrite existing files if true, otherwise throw
+ * runtime_error
*/
template <class Axes, class Storage>
inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
- std::string const& filename, SaveMode mode = SaveMode::append);
+ std::string const& filename, bool overwrite = true);
} // namespace corsika
#include <corsika/detail/framework/utility/SaveBoostHistogram.inl>
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 855fc5a388ed59e7e1b037cf077c6ca580feac5e..edb377953427cb89120828193319b41ce152d8b4 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -7,6 +7,7 @@
*/
#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
@@ -180,7 +181,7 @@ int main(int argc, char** argv) {
em_continuous.reset();
auto const hists = proposalCounted.getHistogram();
- hists.saveLab("inthist_lab_emShower.npz");
- hists.saveCMS("inthist_cms_emShower.npz");
+ save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
+ save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
longprof.save("longprof_emShower.txt");
}
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index b4eab698c94feb5a182ea7c068272eb050b88289..71d8cf82875f39c85a7f6b11dd39f7ce2faac5bf 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -13,6 +13,7 @@
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
@@ -267,11 +268,8 @@ int main(int argc, char** argv) {
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
- hists.saveLab("inthist_lab.txt");
- hists.saveCMS("inthist_cms.txt");
-
- hists.saveLab("inthist_lab.txt");
- hists.saveCMS("inthist_cms.txt");
+ save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true);
+ save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true);
longprof.save("longprof.txt");
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 7d0a97d062940867454349e88f9520769f4ad0fd..28488a1bdb53a8cc9a6b8f234152acf289919dfa 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -15,6 +15,7 @@
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
@@ -281,7 +282,7 @@ int main(int argc, char** argv) {
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram() + proposalCounted.getHistogram();
- hists.saveLab("inthist_lab_verticalEAS.npz");
- hists.saveCMS("inthist_cms_verticalEAS.npz");
+ save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
+ save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
longprof.save("longprof_verticalEAS.txt");
}
diff --git a/tests/framework/testInteractionCounter.cpp b/tests/framework/testInteractionCounter.cpp
index b9cb8b013bbc94dd883367c6bafd4c41b91ae7f1..bf73fc578316eb3c7e948c22671f6a1aad220ddc 100644
--- a/tests/framework/testInteractionCounter.cpp
+++ b/tests/framework/testInteractionCounter.cpp
@@ -21,6 +21,11 @@
#include <catch2/catch.hpp>
#include <numeric>
+#include <algorithm>
+#include <iterator>
+#include <string>
+#include <fstream>
+#include <cstdio>
using namespace corsika;
@@ -69,10 +74,31 @@ TEST_CASE("InteractionCounter", "[process]") {
CHECK(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
- countedProcess.getHistogram().saveLab("testInteractionCounter_file1.npz",
- SaveMode::overwrite);
- countedProcess.getHistogram().saveCMS("testInteractionCounter_file2.npz",
- SaveMode::overwrite);
+ save_hist(countedProcess.getHistogram().labHist(), "testInteractionCounter_file1.npz",
+ true);
+ save_hist(countedProcess.getHistogram().CMSHist(), "testInteractionCounter_file2.npz",
+ true);
+
+ SECTION("output validation") {
+ auto const file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
+ "testInteractionCounter_file2");
+
+ std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
+
+ // compare to binary reference data
+ // note that this currenly compares the whole files byte by byte. If the new
+ // or the reference file are compressed this would be a false-negative outcome
+ // of this test
+ std::ifstream file1(file + ".npz");
+ std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
+
+ std::istreambuf_iterator<char> begin1(file1);
+ std::istreambuf_iterator<char> begin1ref(file1ref);
+
+ std::istreambuf_iterator<char> end;
+
+ CHECK(std::equal(begin1, end, begin1ref));
+ }
}
SECTION("DoInteraction Lambda") {
@@ -93,41 +119,3 @@ TEST_CASE("InteractionCounter", "[process]") {
CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
}
}
-
-#include <algorithm>
-#include <iterator>
-#include <string>
-#include <fstream>
-
-TEST_CASE("InteractionCounterOutput", "[output validation]") {
-
- logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
- "testInteractionCounter_file2");
-
- SECTION(std::string("check saved data, ") + file + ".npz") {
-
- std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
-
- // compare to binary reference data
- std::ifstream file1(file + ".npz");
- std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
-
- std::istreambuf_iterator<char> begin1(file1);
- std::istreambuf_iterator<char> begin1ref(file1ref);
-
- std::istreambuf_iterator<char> end;
-
- while (begin1 != end && begin1ref != end) {
- CHECK(*begin1 == *begin1ref);
- ++begin1;
- ++begin1ref;
- }
- CHECK(begin1 == end);
- CHECK(begin1ref == end);
- file1.close();
- file1ref.close();
- }
-}
diff --git a/tests/framework/testInteractionCounter_file1_REF.npz b/tests/framework/testInteractionCounter_file1_REF.npz
index 13e2d5d49e16b6a2e7eb6b95dc1a26bccab18ae8..ca75c99e4df0292da9db9f05f5a7b3319c482653 100644
Binary files a/tests/framework/testInteractionCounter_file1_REF.npz and b/tests/framework/testInteractionCounter_file1_REF.npz differ
diff --git a/tests/framework/testInteractionCounter_file2_REF.npz b/tests/framework/testInteractionCounter_file2_REF.npz
index 8509b27e4d690a386c2710a24ceda6a1143bf757..2aaa26ed6240a2de96ae24cc566d695d7f4467b1 100644
Binary files a/tests/framework/testInteractionCounter_file2_REF.npz and b/tests/framework/testInteractionCounter_file2_REF.npz differ
diff --git a/tests/framework/testSaveBoostHistogram.cpp b/tests/framework/testSaveBoostHistogram.cpp
index 91bac8b82120b0ce2595a0f8d91a5f503706403c..d5c16438a2ab6de1595393205d4eff4320a43d08 100644
--- a/tests/framework/testSaveBoostHistogram.cpp
+++ b/tests/framework/testSaveBoostHistogram.cpp
@@ -40,5 +40,6 @@ TEST_CASE("SaveHistogram") {
h(a, b, c, d);
}
- corsika::save_hist(h, "hist.npz");
+ REQUIRE_NOTHROW(corsika::save_hist(h, "hist.npz", true));
+ REQUIRE_THROWS(corsika::save_hist(h, "hist.npz", false));
}
diff --git a/tests/modules/testInteractionCounter.cpp b/tests/modules/testInteractionCounter.cpp
deleted file mode 100644
index 350e1dc1de38a96a69d91941a37f0fb371eaccc3..0000000000000000000000000000000000000000
--- a/tests/modules/testInteractionCounter.cpp
+++ /dev/null
@@ -1,136 +0,0 @@
-/*
- * (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.
- */
-
-#include <corsika/modules/InteractionCounter.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-
-#include <catch2/catch.hpp>
-
-#include <numeric>
-
-using namespace corsika;
-
-const std::string refDataDir = std::string(REFDATADIR); // from cmake
-
-struct DummyProcess {
- template <typename TParticle>
- GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
- return 100_g / 1_cm / 1_cm;
- }
-
- template <typename TParticle>
- auto DoInteraction([[maybe_unused]] TParticle& projectile) {
- return nullptr;
- }
-};
-
-TEST_CASE("InteractionCounter", "[process]") {
-
- logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- DummyProcess d;
- InteractionCounter countedProcess(d);
-
- SECTION("GetInteractionLength") {
- CHECK(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
- }
-
- auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
- [[maybe_unused]] auto& env_dummy = env;
-
- SECTION("DoInteraction nucleus") {
- unsigned short constexpr A = 14, Z = 7;
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
- Z, 105_TeV, nodePtr, *csPtr);
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
-
- auto const ret = countedProcess.DoInteraction(*secViewPtr);
- CHECK(ret == nullptr);
-
- auto const& h = countedProcess.GetHistogram().labHist();
- CHECK(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
- CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
- auto const& h2 = countedProcess.GetHistogram().CMSHist();
- CHECK(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
- CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-
- countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
- utl::SaveMode::overwrite);
- countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
- utl::SaveMode::overwrite);
- }
-
- SECTION("DoInteraction Lambda") {
- auto constexpr code = particles::Code::Lambda0;
- auto [stackPtr, secViewPtr] =
- setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
-
- auto const ret = countedProcess.DoInteraction(*secViewPtr);
- CHECK(ret == nullptr);
-
- auto const& h = countedProcess.GetHistogram().labHist();
- CHECK(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
- CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
- auto const& h2 = countedProcess.GetHistogram().CMSHist();
- CHECK(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
- CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
- }
-}
-
-#include <algorithm>
-#include <iterator>
-#include <string>
-#include <fstream>
-
-TEST_CASE("InteractionCounterOutput", "[output validation]") {
-
- logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
- "testInteractionCounter_file2");
-
- SECTION(std::string("check saved data, ") + file + ".npz") {
-
- std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
-
- // compare to binary reference data
- std::ifstream file1(file + ".npz");
- std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
-
- std::istreambuf_iterator<char> begin1(file1);
- std::istreambuf_iterator<char> begin1ref(file1ref);
-
- std::istreambuf_iterator<char> end;
-
- while (begin1 != end && begin1ref != end) {
- CHECK(*begin1 == *begin1ref);
- ++begin1;
- ++begin1ref;
- }
- CHECK(begin1 == end);
- CHECK(begin1ref == end);
- file1.close();
- file1ref.close();
- }
-}
diff --git a/tools/read_hist.py b/tools/read_hist.py
new file mode 100755
index 0000000000000000000000000000000000000000..184ba5e2e72e3d01a8e1e109731793f8cfc9c665
--- /dev/null
+++ b/tools/read_hist.py
@@ -0,0 +1,27 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import boost_histogram as bh
+
+def read_hist(filename):
+ """
+ read numpy file produced with CORSIKA 8's save_hist() function into
+ boost-histogram object.
+ """
+
+ d = np.load(filename)
+ axistypes = d['axistypes'].view('c')
+ overflow = d['overflow']
+ underflow = d['underflow']
+
+ axes = []
+ for i, (at, has_overflow, has_underflow) in enumerate(zip(axistypes, overflow, underflow)):
+ if at == b'c':
+ axes.append(bh.axis.Variable(d[f'binedges_{i}'], overflow=has_overflow, underflow=has_underflow))
+ elif at == b'd':
+ axes.append(bh.axis.IntCategory(d[f'bins_{i}'], growth=(not has_overflow)))
+
+ h = bh.Histogram(*axes)
+ h.view(flow=True)[:] = d['data']
+
+ return h
+