diff --git a/Processes/InteractionCounter/CMakeLists.txt b/Processes/InteractionCounter/CMakeLists.txt
index 30b8815312afe511b787c84a5dbee1efff274cf6..2323207a37e42a6723df1433553f80974606b732 100644
--- a/Processes/InteractionCounter/CMakeLists.txt
+++ b/Processes/InteractionCounter/CMakeLists.txt
@@ -1,24 +1,36 @@
set (
MODEL_HEADERS
InteractionCounter.h
+ InteractionHistogram.h
)
+set (
+ MODEL_SOURCES
+ InteractionHistogram.cc
+)
+
set (
MODEL_NAMESPACE
corsika/process/interaction_counter
)
-FIND_PACKAGE( Boost 1.70 REQUIRED )
-
-add_library (ProcessInteractionCounter INTERFACE)
+add_library (ProcessInteractionCounter STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessInteractionCounter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+set_target_properties (
+ ProcessInteractionCounter
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+# PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessInteractionCounter
- INTERFACE
CORSIKAunits
CORSIKAprocesssequence
+ CORSIKAthirdparty
)
target_include_directories (
@@ -26,7 +38,6 @@ target_include_directories (
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
- ${Boost_INCLUDE_DIR}
)
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
diff --git a/Processes/InteractionCounter/InteractionCounter.h b/Processes/InteractionCounter/InteractionCounter.h
index a6fb6b43b8f8bb245697b935d0d54810ad2a1506..a03c528881ce6e3a67f530c1fe170312cd90ce7f 100644
--- a/Processes/InteractionCounter/InteractionCounter.h
+++ b/Processes/InteractionCounter/InteractionCounter.h
@@ -13,207 +13,10 @@
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/interaction_counter/InteractionHistogram.h>
#include <corsika/setup/SetupStack.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <boost/histogram.hpp>
-#include <boost/histogram/ostream.hpp>
-#include <fstream>
-#include <functional>
-#include <map>
-#include <utility>
-
-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 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>;
-
- 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);
-
- hist_type interaction_histogram_cms_;
- hist_type interaction_histogram_lab_;
-
- /*!
- * These maps map PDG nuclei codes to their corresponding interaction histograms
- */
- nuclear_hist_type nuclear_inthist_lab_, nuclear_inthist_cms_;
-
- public:
- 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))} {}
-
- void fill(particles::Code projectile_id, units::si::HEPEnergyType lab_energy,
- units::si::HEPEnergyType mass_target, int A = 0, int Z = 0) {
- using namespace units::si;
-
- if (projectile_id == particles::Code::Nucleus) {
- auto const sqrtS =
- sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) +
- mass_target * mass_target + 2 * lab_energy * mass_target);
-
- 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);
- } 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);
- }
- }
-
- auto CMSHists() const {
- return std::pair<decltype(interaction_histogram_cms_) const&,
- decltype(nuclear_inthist_cms_) const&>{interaction_histogram_cms_,
- nuclear_inthist_cms_};
- }
-
- auto labHists() const {
- return std::pair<decltype(interaction_histogram_lab_) const&,
- decltype(nuclear_inthist_lab_) const&>(interaction_histogram_lab_,
- nuclear_inthist_lab_);
- }
-
- void saveLab(std::string const& filename) const {
- save(interaction_histogram_lab_, nuclear_inthist_lab_, filename, "lab system");
- }
-
- void saveCMS(std::string const& filename) const {
- save(interaction_histogram_lab_, nuclear_inthist_cms_, filename,
- "center-of-mass system");
- }
-
- auto& 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_;
-
- return *this;
- }
-
- auto operator+(InteractionHistogram other) const {
- other.nuclear_inthist_lab_ += nuclear_inthist_lab_;
- other.nuclear_inthist_cms_ += nuclear_inthist_cms_;
-
- return other;
- }
-
- 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 = "") {
- std::ofstream file;
- file.open(filename);
- saveHist(hist, file, comment);
- }
-
- static void 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;
- }
- }
- }
-
- static void 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;
- }
- }
-
- /*!
- * 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 = "") {
- std::ofstream file;
- file.open(filename);
-
- saveHist(hist, file, comment);
- saveHistMap(histMap, file);
- }
- };
-
/*!
* Wrapper around an InteractionProcess that fills histograms of the number
* of calls to DoInteraction() binned in projectile energy (both in
diff --git a/Processes/InteractionCounter/InteractionHistogram.cc b/Processes/InteractionCounter/InteractionHistogram.cc
new file mode 100644
index 0000000000000000000000000000000000000000..556b58f2788f1063a594bb41fb06a552ed5c563d
--- /dev/null
+++ b/Processes/InteractionCounter/InteractionHistogram.cc
@@ -0,0 +1,147 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/process/interaction_counter/InteractionHistogram.h>
+
+#include <fstream>
+#include <string>
+
+using namespace corsika::process::interaction_counter;
+
+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;
+
+ if (projectile_id == particles::Code::Nucleus) {
+ auto const sqrtS =
+ sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) +
+ mass_target * mass_target + 2 * lab_energy * mass_target);
+
+ 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);
+ } 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);
+ }
+}
+
+void InteractionHistogram::saveLab(std::string const& filename) const {
+ save(interaction_histogram_lab_, nuclear_inthist_lab_, filename, "lab system");
+}
+
+void InteractionHistogram::saveCMS(std::string const& filename) const {
+ save(interaction_histogram_lab_, nuclear_inthist_cms_, filename,
+ "center-of-mass system");
+}
+
+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_;
+
+ return *this;
+}
+
+InteractionHistogram InteractionHistogram::operator+(InteractionHistogram other) const {
+ other.nuclear_inthist_lab_ += nuclear_inthist_lab_;
+ other.nuclear_inthist_cms_ += nuclear_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.h b/Processes/InteractionCounter/InteractionHistogram.h
new file mode 100644
index 0000000000000000000000000000000000000000..98fedae84ee5b101e3904b2d1d18e5dc26e29d32
--- /dev/null
+++ b/Processes/InteractionCounter/InteractionHistogram.h
@@ -0,0 +1,107 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#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;
+ }
+
+ 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);
+
+ hist_type interaction_histogram_cms_;
+ hist_type interaction_histogram_lab_;
+
+ /*!
+ * These maps map PDG nuclei codes to their corresponding interaction histograms
+ */
+ nuclear_hist_type nuclear_inthist_lab_, nuclear_inthist_cms_;
+
+ public:
+ InteractionHistogram();
+
+ //! fill both CMS and lab histograms at the same time
+ 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_};
+ }
+
+ auto labHists() const {
+ return std::pair<decltype(interaction_histogram_lab_) const&,
+ decltype(nuclear_inthist_lab_) const&>(interaction_histogram_lab_,
+ nuclear_inthist_lab_);
+ }
+
+ void saveLab(std::string const& filename) const;
+
+ void saveCMS(std::string const& filename) const;
+
+ 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