diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 3cce0678741508b8de9453dd482e7626cd1cd7f9..936e54ed20ba87e53a25a1bc020fbadc382c773e 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -87,7 +87,7 @@ if (Pythia8_FOUND)
install (TARGETS cascade_proton_example DESTINATION share/examples)
endif()
-CORSIKA_ADD_TEST(vertical_EAS)
+add_executable(vertical_EAS vertical_EAS.cc)
target_link_libraries (vertical_EAS
SuperStupidStack
CORSIKAunits
@@ -100,10 +100,10 @@ target_link_libraries (vertical_EAS
CORSIKAcascade
ProcessEnergyLoss
ProcessObservationPlane
- ProcessTrackWriter
ProcessTrackingLine
ProcessParticleCut
ProcessStackInspector
+ ProcessInteractionCounter
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 8dd6f827df46e28c98c2354c2e122734250d87e1..2cdb101b631f592111452c88738a43715916d7a9 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -9,45 +9,35 @@
*/
#include <corsika/cascade/Cascade.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/StackProcess.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/process/observation_plane/ObservationPlane.h>
-#include <corsika/process/particle_cut/ParticleCut.h>
-#include <corsika/process/switch_process/SwitchProcess.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/NuclearComposition.h>
-
#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
-
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/StackProcess.h>
+#include <corsika/process/energy_loss/EnergyLoss.h>
+#include <corsika/process/interaction_counter/InteractionCounter.h>
+#include <corsika/process/observation_plane/ObservationPlane.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
-
-#include <corsika/process/pythia/Decay.h>
-
+#include <corsika/process/switch_process/SwitchProcess.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/process/urqmd/UrQMD.h>
-
-#include <corsika/process/particle_cut/ParticleCut.h>
-#include <corsika/process/track_writer/TrackWriter.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
#include <corsika/random/RNGManager.h>
-
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iomanip>
#include <iostream>
#include <limits>
+#include <string>
#include <typeinfo>
using namespace corsika;
@@ -64,6 +54,7 @@ using namespace corsika::units::si;
void registerRandomStreams() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("qgran");
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
@@ -71,7 +62,11 @@ void registerRandomStreams() {
random::RNGManager::GetInstance().SeedAll();
}
-int main() {
+int main(int argc, char** argv) {
+ if (argc != 4) {
+ std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV>" << std::endl;
+ return 1;
+ }
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams();
@@ -97,15 +92,17 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const Code beamCode = Code::Proton;
- auto const mass = particles::GetMass(beamCode);
- const HEPEnergyType E0 = 0.1_PeV;
+ const Code beamCode = Code::Nucleus;
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ unsigned short Z = std::stoi(std::string(argv[2]));
+ auto const mass = particles::GetNucleusMass(A, Z);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
double phi = 0.;
Point const injectionPos(
rootCS, 0_m, 0_m,
- 112.8_km * 0.999 +
+ 112.7_km +
builder.earthRadius); // this is the CORSIKA 7 start of atmosphere/universe
// {
@@ -124,23 +121,33 @@ int main() {
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
- stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- beamCode, E0, plab, injectionPos, 0_ns});
- // }
+ if (A != 1) {
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ beamCode, E0, plab, injectionPos, 0_ns, A, Z});
+
+ } else {
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, E0, plab, injectionPos, 0_ns});
+ }
Line const line(injectionPos, plab.normalized() * 1_m * 1_Hz);
auto const velocity = line.GetV0().norm();
- auto const observationHeight = 1.425_km + builder.earthRadius;
+ auto const observationHeight = 1.4_km + builder.earthRadius;
- setup::Trajectory const showerAxis(line, (112.8_km - observationHeight) / velocity);
+ setup::Trajectory const showerAxis(line, (112.7_km - observationHeight) / velocity);
// setup processes, decays and interactions
process::sibyll::Interaction sibyll;
+ process::interaction_counter::InteractionCounter sibyllCounted(sibyll);
+
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ process::interaction_counter::InteractionCounter sibyllNucCounted(sibyllNuc);
process::pythia::Decay decayPythia;
@@ -166,9 +173,8 @@ int main() {
});
decaySibyll.PrintDecayConfig();
- process::particle_cut::ParticleCut cut(5_GeV);
+ process::particle_cut::ParticleCut cut(100_GeV);
- process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss(showerAxis);
Plane const obsPlane(Point(rootCS, 0_m, 0_m, observationHeight),
@@ -180,16 +186,20 @@ int main() {
process::UrQMD::UrQMD urqmd;
- auto sibyllSequence = sibyll << sibyllNuc;
+ auto sibyllSequence = sibyllNucCounted << sibyllCounted;
process::switch_process::SwitchProcess switchProcess(urqmd, sibyllSequence, 55_GeV);
auto decaySequence = decayPythia << decaySibyll;
- auto sequence = switchProcess << decaySequence << eLoss << cut << observationLevel
- << trackWriter;
+ auto sequence = switchProcess << decaySequence << eLoss << cut << observationLevel;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.SetNodes();
+ // EAS.forceInteraction();
+
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
@@ -202,6 +212,10 @@ int main() {
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+ auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram();
+ hists.saveLab("inthist_lab.txt");
+ hists.saveCMS("inthist_cms.txt");
+
std::ofstream finish("finished");
finish << "run completed without error" << std::endl;
}
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 9ca109dee4bbef7442de84944047730ffa0ba98c..6e93793eceed7a20cdaf7363dd585cc93b1afbac 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -10,7 +10,6 @@ if (PYTHIA8_FOUND)
endif (PYTHIA8_FOUND)
add_subdirectory (HadronicElasticModel)
add_subdirectory (UrQMD)
-add_subdirectory (SwitchProcess)
# continuous physics
add_subdirectory (EnergyLoss)
@@ -22,6 +21,10 @@ add_subdirectory (StackInspector)
# cuts, thinning, etc.
add_subdirectory (ParticleCut)
+# meta-processes
+add_subdirectory (InteractionCounter)
+add_subdirectory (SwitchProcess)
+
##########################################
# add_custom_target(CORSIKAprocesses)
add_library (CORSIKAprocesses INTERFACE)
diff --git a/Processes/InteractionCounter/CMakeLists.txt b/Processes/InteractionCounter/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2323207a37e42a6723df1433553f80974606b732
--- /dev/null
+++ b/Processes/InteractionCounter/CMakeLists.txt
@@ -0,0 +1,58 @@
+set (
+ MODEL_HEADERS
+ InteractionCounter.h
+ InteractionHistogram.h
+ )
+
+set (
+ MODEL_SOURCES
+ InteractionHistogram.cc
+)
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/interaction_counter
+ )
+
+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
+ CORSIKAunits
+ CORSIKAprocesssequence
+ CORSIKAthirdparty
+ )
+
+target_include_directories (
+ ProcessInteractionCounter
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+# --------------------
+# code unit testing
+CORSIKA_ADD_TEST(testInteractionCounter SOURCES
+ testInteractionCounter.cc
+ ${MODEL_HEADERS}
+)
+
+target_link_libraries (
+ testInteractionCounter
+ ProcessInteractionCounter
+ CORSIKAunits
+ CORSIKAenvironment
+ CORSIKAtesting
+ )
diff --git a/Processes/InteractionCounter/InteractionCounter.h b/Processes/InteractionCounter/InteractionCounter.h
new file mode 100644
index 0000000000000000000000000000000000000000..a03c528881ce6e3a67f530c1fe170312cd90ce7f
--- /dev/null
+++ b/Processes/InteractionCounter/InteractionCounter.h
@@ -0,0 +1,67 @@
+/*
+ * (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.
+ */
+
+#ifndef _corsika_InteractionCounter_h
+#define _corsika_InteractionCounter_h
+
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/interaction_counter/InteractionHistogram.h>
+#include <corsika/setup/SetupStack.h>
+
+namespace corsika::process::interaction_counter {
+ /*!
+ * Wrapper around an InteractionProcess that fills histograms of the number
+ * of calls to DoInteraction() binned in projectile energy (both in
+ * lab and center-of-mass frame) and species
+ */
+ template <class TCountedProcess>
+ class InteractionCounter
+ : public InteractionProcess<InteractionCounter<TCountedProcess>> {
+
+ TCountedProcess& process_;
+ InteractionHistogram histogram_;
+
+ public:
+ InteractionCounter(TCountedProcess& process)
+ : process_(process) {}
+
+ template <typename TProjectile>
+ auto DoInteraction(TProjectile& projectile) {
+ auto const massNumber = projectile.GetNode()
+ ->GetModelProperties()
+ .GetNuclearComposition()
+ .GetAverageMassNumber();
+ auto const massTarget = massNumber * units::constants::nucleonMass;
+
+ if (auto const projectile_id = projectile.GetPID();
+ projectile_id == particles::Code::Nucleus) {
+ auto const A = projectile.GetNuclearA();
+ auto const Z = projectile.GetNuclearZ();
+ histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget, A, Z);
+ } else {
+ histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget);
+ }
+ return process_.DoInteraction(projectile);
+ }
+
+ void Init() { process_.Init(); }
+
+ template <typename TParticle>
+ auto GetInteractionLength(TParticle const& particle) const {
+ return process_.GetInteractionLength(particle);
+ }
+
+ auto const& GetHistogram() const { return histogram_; }
+ };
+
+} // namespace corsika::process::interaction_counter
+
+#endif
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
diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc
new file mode 100644
index 0000000000000000000000000000000000000000..797d0f87d1205fb0db28ae66b1ca21c50e709d3b
--- /dev/null
+++ b/Processes/InteractionCounter/testInteractionCounter.cc
@@ -0,0 +1,156 @@
+
+/*
+ * (c) Copyright 2018 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/InteractionCounter.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupStack.h>
+
+#include <catch2/catch.hpp>
+
+#include <numeric>
+
+using namespace corsika;
+using namespace corsika::process::interaction_counter;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+auto setupEnvironment(particles::Code target_code) {
+ // setup environment, geometry
+ auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
+ auto& universe = *(env->GetUniverse());
+ const geometry::CoordinateSystem& cs = env->GetCoordinateSystem();
+
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{cs, 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(std::vector<particles::Code>{target_code},
+ std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
+
+ return std::make_tuple(std::move(env), &cs, nodePtr);
+}
+
+template <typename TNodeType>
+auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
+ geometry::CoordinateSystem const& cs) {
+ auto stack = std::make_unique<setup::Stack>();
+ auto constexpr mN = corsika::units::constants::nucleonMass;
+
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+
+ HEPEnergyType const E0 =
+ sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
+ auto particle =
+ stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType, unsigned short, unsigned short>{
+ particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ});
+
+ particle.SetNode(vNodePtr);
+ return std::make_tuple(
+ std::move(stack),
+ std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
+}
+
+template <typename TNodeType>
+auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
+ TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
+ auto stack = std::make_unique<setup::Stack>();
+
+ geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+ corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+
+ HEPEnergyType const E0 =
+ sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
+ pLab.squaredNorm());
+ auto particle = stack->AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ vProjectileType, E0, pLab, origin, 0_ns});
+
+ particle.SetNode(vNodePtr);
+ return std::make_tuple(
+ std::move(stack),
+ std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
+}
+
+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") {
+ DummyProcess d;
+ InteractionCounter countedProcess(d);
+
+ SECTION("GetInteractionLength") {
+ REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
+ }
+
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ SECTION("DoInteraction nucleus") {
+ unsigned short constexpr A = 14, Z = 7;
+ auto [stackPtr, secViewPtr] = setupStack(A, Z, 105_TeV, nodePtr, *csPtr);
+
+ auto projectile = secViewPtr->GetProjectile();
+ auto const ret = countedProcess.DoInteraction(projectile);
+ REQUIRE(ret == nullptr);
+
+ auto const& h = countedProcess.GetHistogram().labHists().second.at(1'000'070'140);
+ REQUIRE(h.at(50) == 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
+ REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
+ }
+
+ SECTION("DoInteraction Lambda") {
+ auto constexpr code = particles::Code::Lambda0;
+ auto constexpr codeInt = static_cast<particles::CodeIntType>(code);
+ auto [stackPtr, secViewPtr] = setupStack(code, 105_TeV, nodePtr, *csPtr);
+
+ auto projectile = secViewPtr->GetProjectile();
+ auto const ret = countedProcess.DoInteraction(projectile);
+ REQUIRE(ret == nullptr);
+
+ auto const& h = countedProcess.GetHistogram().labHists().first;
+ REQUIRE(h.at(codeInt, 50) == 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
+ REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
+ }
+}
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index a79fd80c2a5bc751bf247b4a01cc725986f9501f..d557f28662be5ea754074d2ebf3f25828146ca6b 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -177,7 +177,7 @@ namespace corsika::process::sibyll {
template <>
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle& vP,
+ NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle const& vP,
const particles::Code TargetId) {
using namespace units::si;
if (vP.GetPID() != particles::Code::Nucleus)
@@ -216,7 +216,7 @@ namespace corsika::process::sibyll {
template <>
template <>
units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
- Particle& vP) {
+ Particle const& vP) {
using namespace units;
using namespace units::si;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 199901cb8cd3b1061b6b2c2a30ea1f37e77d0e03..1be24afccb741d5f4c4d96dc59c1be9105c87507 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -50,10 +50,10 @@ namespace corsika::process::sibyll {
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
- GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
+ GetCrossSection(Particle const& p, const corsika::particles::Code TargetId);
template <typename Particle>
- corsika::units::si::GrammageType GetInteractionLength(Particle&);
+ corsika::units::si::GrammageType GetInteractionLength(Particle const&);
template <typename Projectile>
corsika::process::EProcessReturn DoInteraction(Projectile&);