diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 5ceb0de05c8fd9d401289d30aa7a80858facabbf..042e64bc29f468afa32c45b654b9082adedc33dd 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -252,11 +252,8 @@ int main(int argc, char** argv) {
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
- hists.saveLab("inthist_lab.txt");
- hists.saveCMS("inthist_cms.txt");
-
- hists.saveLab("inthist_lab.txt");
- hists.saveCMS("inthist_cms.txt");
+ hists.saveLab("inthist_lab.npz");
+ hists.saveCMS("inthist_cms.npz");
longprof.save("longprof.txt");
diff --git a/Framework/Utilities/SaveBoostHistogram.hpp b/Framework/Utilities/SaveBoostHistogram.hpp
index 8418457f0ab8c9fe86412a77352cca278e919f35..813099c4c3a13f9ba3623dedf6e0590e8b12597e 100644
--- a/Framework/Utilities/SaveBoostHistogram.hpp
+++ b/Framework/Utilities/SaveBoostHistogram.hpp
@@ -1,3 +1,11 @@
+/*
+ * (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 <corsika/third_party/cnpy/cnpy.hpp>
#include <boost/histogram.hpp>
@@ -8,97 +16,96 @@
#include <utility>
#include <vector>
-namespace corsika {
- namespace utl {
-
- /**
- * 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
- * counts (i.e. a double or int) are supported, nothing fancy like profiles.
- *
- * 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)
- */
- template <class Axes, class Storage>
- void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
- std::string const& filename) {
- int const rank = h.rank();
-
- std::vector<size_t> axes_dims;
- axes_dims.reserve(rank);
+#pragma once
- auto overflow = std::make_unique<bool[]>(rank);
- auto underflow = std::make_unique<bool[]>(rank);
+namespace corsika::utl {
+ /**
+ * 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
+ * counts (i.e. a double or int) are supported, nothing fancy like profiles.
+ *
+ * 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)
+ */
+ template <class Axes, class Storage>
+ inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
+ std::string const& filename) {
+ int const rank = h.rank();
- std::vector<char> axis_types;
- axis_types.reserve(rank);
+ std::vector<size_t> axes_dims;
+ axes_dims.reserve(rank);
- for (int i = 0; i < rank; ++i) {
- auto const& ax = h.axis(i);
- unsigned const has_underflow = (ax.options() & 0x01) ? 1 : 0;
- unsigned const has_overflow = (ax.options() & 0x02) ? 1 : 0;
+ auto overflow = std::make_unique<bool[]>(rank);
+ auto underflow = std::make_unique<bool[]>(rank);
- underflow[i] = has_underflow;
- overflow[i] = has_overflow;
+ std::vector<char> axis_types;
+ axis_types.reserve(rank);
- axes_dims.emplace_back(ax.size() + has_underflow + has_overflow);
+ for (int i = 0; i < rank; ++i) {
+ auto const& ax = h.axis(i);
+ unsigned const has_underflow = (ax.options() & 0x01) ? 1 : 0;
+ unsigned const has_overflow = (ax.options() & 0x02) ? 1 : 0;
- if (ax.continuous()) {
- axis_types.push_back('c');
- std::vector<double> ax_edges;
- ax_edges.reserve(ax.size() + 1);
+ underflow[i] = has_underflow;
+ overflow[i] = has_overflow;
- for (int j = 0; j < ax.size(); ++j) { ax_edges.push_back(ax.bin(j).lower()); }
+ axes_dims.emplace_back(ax.size() + has_underflow + has_overflow);
- ax_edges.push_back(ax.bin(ax.size() - 1).upper());
+ if (ax.continuous()) {
+ axis_types.push_back('c');
+ std::vector<double> ax_edges;
+ ax_edges.reserve(ax.size() + 1);
- cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
- 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
- bins.reserve(ax.size());
+ for (int j = 0; j < ax.size(); ++j) { ax_edges.push_back(ax.bin(j).lower()); }
- for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
+ ax_edges.push_back(ax.bin(ax.size() - 1).upper());
- cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
- {bins.size()}, "a");
- }
+ cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
+ 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
+ bins.reserve(ax.size());
- cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(),
- {axis_types.size()}, "a");
+ for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
- cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(),
- {axis_types.size()}, "a");
- cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(),
- {axis_types.size()}, "a");
+ cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
+ {bins.size()}, "a");
}
- auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
- unsigned{1}, std::multiplies<>());
- auto temp = std::make_unique<float[]>(prod_axis_size);
+ cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(),
+ {axis_types.size()}, "a");
- assert(prod_axis_size == h.size());
+ cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(),
+ {axis_types.size()}, "a");
+ cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(),
+ {axis_types.size()}, "a");
+ }
+
+ auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
+ unsigned{1}, std::multiplies<>());
+ auto temp = std::make_unique<float[]>(prod_axis_size);
- // reduce multi-dim. to 1-dim, row-major (i.e., last axis index is contiguous in
- // memory) take special care of underflow bins which have -1 as index and thus need
- // to be shifted by +1
+ assert(prod_axis_size == h.size());
- for (auto&& x : indexed(h, boost::histogram::coverage::all)) {
- int p = 0;
+ // reduce multi-dim. to 1-dim, row-major (i.e., last axis index is contiguous in
+ // memory) take special care of underflow bins which have -1 as index and thus need
+ // to be shifted by +1
- for (int axis_index = 0; axis_index < rank; ++axis_index) {
- int const offset_underflow = (h.axis(axis_index).options() & 0x01) ? 1 : 0;
- auto k = x.index(axis_index) + offset_underflow;
- p = k + p * axes_dims.at(axis_index);
- }
+ for (auto&& x : indexed(h, boost::histogram::coverage::all)) {
+ int p = 0;
- temp[p] = *x;
+ for (int axis_index = 0; axis_index < rank; ++axis_index) {
+ int const offset_underflow = (h.axis(axis_index).options() & 0x01) ? 1 : 0;
+ auto k = x.index(axis_index) + offset_underflow;
+ p = k + p * axes_dims.at(axis_index);
}
- 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']
+ temp[p] = *x;
}
- } // namespace utl
-} // namespace corsika
+
+ 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::utl
diff --git a/Processes/InteractionCounter/CMakeLists.txt b/Processes/InteractionCounter/CMakeLists.txt
index 7602d4e7c8250d8960d1686290846a5e826d4283..e7f374dbf345b888659c7b035b03ec1f24cf7494 100644
--- a/Processes/InteractionCounter/CMakeLists.txt
+++ b/Processes/InteractionCounter/CMakeLists.txt
@@ -28,6 +28,7 @@ set_target_properties (
target_link_libraries (
ProcessInteractionCounter
CORSIKAunits
+ CORSIKAutilities
CORSIKAprocesssequence
CORSIKAthirdparty
C8::ext::boost
diff --git a/Processes/InteractionCounter/InteractionHistogram.cc b/Processes/InteractionCounter/InteractionHistogram.cc
index a54b63bd185cbc0d2892bc564f1c294889d5e9eb..8ae626ffe57d1b0dd1cc016667b409a26b697f9a 100644
--- a/Processes/InteractionCounter/InteractionHistogram.cc
+++ b/Processes/InteractionCounter/InteractionHistogram.cc
@@ -13,10 +13,15 @@
using namespace corsika::process::interaction_counter;
+InteractionHistogram::InteractionHistogram()
+ : inthist_cms_{detail::hist_factory(num_bins_cms, lower_edge_cms, upper_edge_cms)}
+ , inthist_lab_{detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab)} {}
+
void InteractionHistogram::fill(particles::Code projectile_id,
units::si::HEPEnergyType lab_energy,
units::si::HEPEnergyType mass_target, int A, int Z) {
using namespace units::si;
+ auto constexpr inv_eV = 1 / 1_eV;
if (projectile_id == particles::Code::Nucleus) {
auto const sqrtS =
sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) +
@@ -24,123 +29,37 @@ void InteractionHistogram::fill(particles::Code projectile_id,
int32_t const pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
- if (nuclear_inthist_cms_.count(pdg) == 0) {
- nuclear_inthist_lab_.emplace(
- pdg,
- boost::histogram::make_histogram(
- boost::histogram::axis::regular<double,
- boost::histogram::axis::transform::log>(
- num_bins_lab, lower_edge_lab, upper_edge_lab)));
-
- nuclear_inthist_cms_.emplace(
- pdg,
- boost::histogram::make_histogram(
- boost::histogram::axis::regular<double,
- boost::histogram::axis::transform::log>(
- num_bins_cms, lower_edge_cms, upper_edge_cms)));
- }
-
- nuclear_inthist_lab_[pdg](lab_energy / 1_GeV);
- nuclear_inthist_cms_[pdg](sqrtS / 1_GeV);
+ inthist_lab_(pdg, lab_energy * inv_eV);
+ inthist_cms_(pdg, sqrtS * inv_eV);
} else {
auto const projectile_mass = particles::GetMass(projectile_id);
auto const sqrtS = sqrt(projectile_mass * projectile_mass +
mass_target * mass_target + 2 * lab_energy * mass_target);
- interaction_histogram_cms_(
- static_cast<corsika::particles::CodeIntType>(projectile_id), sqrtS / 1_GeV);
- interaction_histogram_lab_(
- static_cast<corsika::particles::CodeIntType>(projectile_id), lab_energy / 1_GeV);
+ inthist_cms_(static_cast<int>(particles::GetPDG(projectile_id)), sqrtS * inv_eV);
+ inthist_lab_(static_cast<int>(particles::GetPDG(projectile_id)), lab_energy * inv_eV);
}
}
void InteractionHistogram::saveLab(std::string const& filename) const {
- save(interaction_histogram_lab_, nuclear_inthist_lab_, filename, "lab system");
+ corsika::utl::save_hist(inthist_lab_, filename);
}
void InteractionHistogram::saveCMS(std::string const& filename) const {
- save(interaction_histogram_lab_, nuclear_inthist_cms_, filename,
- "center-of-mass system");
+ corsika::utl::save_hist(inthist_cms_, filename);
}
InteractionHistogram& InteractionHistogram::operator+=(
InteractionHistogram const& other) {
- interaction_histogram_lab_ += other.interaction_histogram_lab_;
- interaction_histogram_cms_ += other.interaction_histogram_cms_;
-
- nuclear_inthist_lab_ += other.nuclear_inthist_lab_;
- nuclear_inthist_cms_ += other.nuclear_inthist_cms_;
+ inthist_lab_ += other.inthist_lab_;
+ inthist_cms_ += other.inthist_cms_;
return *this;
}
InteractionHistogram InteractionHistogram::operator+(InteractionHistogram other) const {
- other.nuclear_inthist_lab_ += nuclear_inthist_lab_;
- other.nuclear_inthist_cms_ += nuclear_inthist_cms_;
-
- other.interaction_histogram_lab_ += interaction_histogram_lab_;
- other.interaction_histogram_cms_ += interaction_histogram_cms_;
+ other.inthist_lab_ += inthist_lab_;
+ other.inthist_cms_ += inthist_cms_;
return other;
}
-
-void InteractionHistogram::saveHist(hist_type const& hist, std::string const& filename,
- std::string const& comment) {
- std::ofstream file;
- file.open(filename);
- saveHist(hist, file, comment);
-}
-
-void InteractionHistogram::saveHist(hist_type const& hist, std::ofstream& file,
- std::string const& comment) {
- auto const& energy_axis = hist.axis(1);
-
- file << "# interaction count histogram (" << comment << ")" << std::endl
- << "# " << energy_axis.size() << " bins between " << energy_axis.bin(0).lower()
- << " and " << energy_axis.bin(energy_axis.size() - 1).upper() << " GeV"
- << std::endl;
-
- for (particles::CodeIntType p = 0;
- p < static_cast<particles::CodeIntType>(particles::Code::LastParticle); ++p) {
-
- if (auto pdg = static_cast<particles::PDGCodeType>(
- particles::GetPDG(static_cast<particles::Code>(p)));
- pdg < 1'000'000'000l) {
- file << "# " << static_cast<particles::Code>(p) << std::endl;
- file << pdg;
- for (int i = 0; i < energy_axis.size(); ++i) { file << ' ' << hist.at(p, i); }
- file << std::endl;
- }
- }
-}
-
-void InteractionHistogram::saveHistMap(nuclear_hist_type const& histMap,
- std::ofstream& file) {
- file << "# nuclei" << std::endl;
- for (auto const& [pdg, hist] : histMap) {
- auto const num_ebins_nucl = hist.axis(0).size();
-
- file << pdg << ' ';
- for (int i = 0; i < num_ebins_nucl; ++i) { file << ' ' << hist.at(i); }
- file << std::endl;
- }
-}
-
-void InteractionHistogram::save(hist_type const& hist, nuclear_hist_type const& histMap,
- std::string const& filename, std::string const& comment) {
- std::ofstream file;
- file.open(filename);
-
- saveHist(hist, file, comment);
- saveHistMap(histMap, file);
-}
-
-InteractionHistogram::InteractionHistogram()
- : interaction_histogram_cms_{boost::histogram::make_histogram(
- boost::histogram::axis::integer<short>(0, numIds),
- boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
- num_bins_cms, lower_edge_cms, upper_edge_cms))}
- , interaction_histogram_lab_{boost::histogram::make_histogram(
- boost::histogram::axis::integer<short>(0, numIds),
- boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
- num_bins_lab, lower_edge_lab, upper_edge_lab))} {}
diff --git a/Processes/InteractionCounter/InteractionHistogram.hpp b/Processes/InteractionCounter/InteractionHistogram.hpp
index ec9124e491d2be28eab999561fde2c7c831e5224..96d243f1223505498d211d675fc8cba6ec620fb6 100644
--- a/Processes/InteractionCounter/InteractionHistogram.hpp
+++ b/Processes/InteractionCounter/InteractionHistogram.hpp
@@ -10,47 +10,38 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/SaveBoostHistogram.hpp>
#include <boost/histogram.hpp>
-#include <boost/histogram/ostream.hpp>
#include <fstream>
#include <functional>
#include <map>
#include <utility>
namespace corsika::process::interaction_counter {
- using hist_type = decltype(boost::histogram::make_histogram(
- std::declval<boost::histogram::axis::integer<particles::CodeIntType>>(),
- std::declval<boost::histogram::axis::regular<
- double, boost::histogram::axis::transform::log>>()));
-
- using nucl_hist_type = decltype(boost::histogram::make_histogram(
- std::declval<boost::histogram::axis::regular<
- double, boost::histogram::axis::transform::log>>()));
-
- using nuclear_hist_type = std::map<int32_t, nucl_hist_type>;
-
- template <typename TKey, typename TVal>
- auto operator+=(std::map<TKey, TVal>& a, std::map<TKey, TVal> b) {
- a.merge(b);
- for (auto const& [id, hist] : b) { a[id] += hist; }
- return a;
- }
+ namespace detail {
+ inline auto hist_factory(unsigned int bin_number, float e_low, float e_high) {
+ namespace bh = boost::histogram;
+ namespace bha = bh::axis;
+
+ auto h = bh::make_histogram(
+ bha::category<int, bh::use_default, bha::option::growth_t>{{2212, 2112},
+ "projectile PDG"},
+ bha::regular<float, bha::transform::log>{bin_number, e_low, e_high,
+ "energy/eV"});
+ return h;
+ }
+ } // namespace detail
class InteractionHistogram {
- static auto constexpr lower_edge_cms = 1., upper_edge_cms = 1e8; // GeV sqrt s
- static auto constexpr lower_edge_lab = 1., upper_edge_lab = 1e12; // GeV lab
- static auto constexpr num_bins_lab = 120, num_bins_cms = 80;
- static auto constexpr numIds =
- static_cast<particles::CodeIntType>(particles::Code::LastParticle);
+ static double constexpr lower_edge_cms = 1e3, upper_edge_cms = 1e17; // eV sqrt s
+ static double constexpr lower_edge_lab = 1e3, upper_edge_lab = 1e21; // eV lab
+ static unsigned int constexpr num_bins_lab = 18 * 10, num_bins_cms = 14 * 10;
- hist_type interaction_histogram_cms_;
- hist_type interaction_histogram_lab_;
+ using hist_type =
+ decltype(detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab));
- /*!
- * These maps map PDG nuclei codes to their corresponding interaction histograms
- */
- nuclear_hist_type nuclear_inthist_lab_, nuclear_inthist_cms_;
+ hist_type inthist_cms_, inthist_lab_;
public:
InteractionHistogram();
@@ -59,17 +50,9 @@ namespace corsika::process::interaction_counter {
void fill(particles::Code projectile_id, units::si::HEPEnergyType lab_energy,
units::si::HEPEnergyType mass_target, int A = 0, int Z = 0);
- auto CMSHists() const {
- return std::pair<decltype(interaction_histogram_cms_) const&,
- decltype(nuclear_inthist_cms_) const&>{interaction_histogram_cms_,
- nuclear_inthist_cms_};
- }
+ hist_type const& CMSHist() const { return inthist_cms_; }
- auto labHists() const {
- return std::pair<decltype(interaction_histogram_lab_) const&,
- decltype(nuclear_inthist_lab_) const&>(interaction_histogram_lab_,
- nuclear_inthist_lab_);
- }
+ hist_type const& labHist() const { return inthist_lab_; }
void saveLab(std::string const& filename) const;
@@ -78,30 +61,6 @@ namespace corsika::process::interaction_counter {
InteractionHistogram& operator+=(InteractionHistogram const& other);
InteractionHistogram operator+(InteractionHistogram other) const;
-
- private:
- /*!
- * Save a histogram into a text file. The \arg comment string is written
- * into the header of the text file.
- * This method is static so that you can sum the histograms of multiple
- * InteractionProcesses before saving them into the same file.
- */
-
- static void saveHist(hist_type const& hist, std::string const& filename,
- std::string const& comment = "");
-
- static void saveHist(hist_type const& hist, std::ofstream& file,
- std::string const& comment = "");
-
- static void saveHistMap(nuclear_hist_type const& histMap, std::ofstream& file);
-
- /*!
- * save both the "normal" particle histograms as well as the "nuclear" histograms
- * into the same file
- */
-
- static void save(hist_type const& hist, nuclear_hist_type const& histMap,
- std::string const& filename, std::string const& comment = "");
};
} // namespace corsika::process::interaction_counter
diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc
index 08fc95113c2955b5782247d55133336d15e01d5b..af43fe63eae05e941998fef7e180a47ca020e032 100644
--- a/Processes/InteractionCounter/testInteractionCounter.cc
+++ b/Processes/InteractionCounter/testInteractionCounter.cc
@@ -126,12 +126,13 @@ TEST_CASE("InteractionCounter") {
auto const ret = countedProcess.DoInteraction(*secViewPtr);
REQUIRE(ret == nullptr);
- auto const& h = countedProcess.GetHistogram().labHists().second.at(1'000'070'140);
- REQUIRE(h.at(50) == 1);
+ auto const& h = countedProcess.GetHistogram().labHist();
+ REQUIRE(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
- auto const& h2 = countedProcess.GetHistogram().CMSHists().second.at(1'000'070'140);
- REQUIRE(h2.at(32) == 1); // bin 1.584 .. 1.995 TeV √s
+ auto const& h2 = countedProcess.GetHistogram().CMSHist();
+ REQUIRE(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
+ // REQUIRE(h2.at(1'000'070'140, 92) == 1); // bin 1.584 .. 1.995 TeV √s
REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
}
@@ -145,12 +146,12 @@ TEST_CASE("InteractionCounter") {
auto const ret = countedProcess.DoInteraction(*secViewPtr);
REQUIRE(ret == nullptr);
- auto const& h = countedProcess.GetHistogram().labHists().first;
- REQUIRE(h.at(codeInt, 50) == 1);
+ auto const& h = countedProcess.GetHistogram().labHist();
+ REQUIRE(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
- auto const& h2 = countedProcess.GetHistogram().CMSHists().first;
- REQUIRE(h2.at(codeInt, 32) == 1); // bin 1.584 .. 1.995 TeV √s
+ auto const& h2 = countedProcess.GetHistogram().CMSHist();
+ REQUIRE(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
}
}