Merge branch 'interaction_counter' into 'master'

New output format for interaction histogram and simplifications

See merge request AirShowerPhysics/corsika!268

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");
 

Framework/Utilities/SaveBoostHistogram.hpp

+/*
+ * (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

Processes/InteractionCounter/CMakeLists.txt

@@ -28,6 +28,7 @@ set_target_properties (
 target_link_libraries (
   ProcessInteractionCounter
   CORSIKAunits
+  CORSIKAutilities
   CORSIKAprocesssequence
   CORSIKAthirdparty
   C8::ext::boost

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))} {}

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

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);
   }
 }