diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 3cce0678741508b8de9453dd482e7626cd1cd7f9..13848d99808f1e9ec83f7727674e76ae82b4c2f6 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -87,7 +87,8 @@ if (Pythia8_FOUND)
install (TARGETS cascade_proton_example DESTINATION share/examples)
endif()
-CORSIKA_ADD_TEST(vertical_EAS)
+#CORSIKA_ADD_TEST(vertical_EAS)
+add_executable(vertical_EAS vertical_EAS.cc)
target_link_libraries (vertical_EAS
SuperStupidStack
CORSIKAunits
@@ -104,6 +105,7 @@ target_link_libraries (vertical_EAS
ProcessTrackingLine
ProcessParticleCut
ProcessStackInspector
+ ProcessInteractionCounter
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 8dd6f827df46e28c98c2354c2e122734250d87e1..ae48ca770681eeba48e7f6626e0e4a7dd7afbfa5 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -12,6 +12,7 @@
#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/switch_process/SwitchProcess.h>
@@ -48,6 +49,7 @@
#include <iomanip>
#include <iostream>
#include <limits>
+#include <string>
#include <typeinfo>
using namespace corsika;
@@ -64,6 +66,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 +74,11 @@ void registerRandomStreams() {
random::RNGManager::GetInstance().SeedAll();
}
-int main() {
+int main(int argc, char** argv) {
+ if (argc != 4) {
+ std::cerr << "must provide A, Z, energy" << std::endl;
+ return 1;
+ }
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams();
@@ -97,15 +104,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 +133,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,30 +185,31 @@ 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::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss(showerAxis);
Plane const obsPlane(Point(rootCS, 0_m, 0_m, observationHeight),
Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- process::observation_plane::ObservationPlane observationLevel(obsPlane,
- "particles.dat");
+ process::observation_plane::ObservationPlane observationLevel(obsPlane, "/dev/null");
// assemble all processes into an ordered process list
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;
+ // << trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
+// EAS.SetNodes();
+// EAS.forceInteraction();
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
@@ -201,6 +221,13 @@ int main() {
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+
+ auto const cms_hists = *std::get<0>(sibyllCounted.CMSHists()) + *std::get<0>(sibyllNucCounted.CMSHists());
+ auto const lab_hists = *std::get<0>(sibyllCounted.labHists()) + *std::get<0>(sibyllNucCounted.labHists());
+
+ process::interaction_counter::saveHist(cms_hists, *std::get<1>(sibyllNucCounted.CMSHists()), "intcount_hist_cms.txt", "center-of-mass system");
+ process::interaction_counter::saveHist(lab_hists, *std::get<1>(sibyllNucCounted.labHists()), "intcount_hist_lab.txt", "lab system");
+
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..4ce877826cc206270220754cad4bcb206bd3fe97
--- /dev/null
+++ b/Processes/InteractionCounter/CMakeLists.txt
@@ -0,0 +1,43 @@
+set (
+ MODEL_HEADERS
+ InteractionCounter.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/interaction_counter
+ )
+
+FIND_PACKAGE( Boost 1.70 REQUIRED )
+
+add_library (ProcessInteractionCounter INTERFACE)
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessInteractionCounter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessInteractionCounter
+ INTERFACE
+ CORSIKAunits
+ CORSIKAprocesssequence
+ )
+
+target_include_directories (
+ ProcessInteractionCounter
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ ${Boost_INCLUDE_DIR}
+ )
+
+install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
+
+# --------------------
+# code unit testing
+#~ CORSIKA_ADD_TEST(testLuminosityCounter)
+#~ target_link_libraries (
+ #~ testLuminosityCounter
+ #~ ProcessLuminosityCounter
+ #~ CORSIKAstackinterface
+ #~ CORSIKAtesting
+#~ )
+
diff --git a/Processes/InteractionCounter/InteractionCounter.h b/Processes/InteractionCounter/InteractionCounter.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c5de97a8944ea0ec28292471d30a8509048824f
--- /dev/null
+++ b/Processes/InteractionCounter/InteractionCounter.h
@@ -0,0 +1,172 @@
+/*
+ * (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/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>
+
+namespace corsika::process::interaction_counter {
+
+ template <class TCountedProcess>
+ class InteractionCounter
+ : public InteractionProcess<InteractionCounter<TCountedProcess>> {
+ 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);
+
+ using hist_type_count_cms = decltype(boost::histogram::make_histogram(
+ boost::histogram::axis::integer<particles::CodeIntType>(0, numIds),
+ boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
+ num_bins_cms, lower_edge_cms, upper_edge_cms)));
+
+ using hist_type_count_lab = decltype(boost::histogram::make_histogram(
+ boost::histogram::axis::integer<particles::CodeIntType>(0, numIds),
+ boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
+ num_bins_lab, lower_edge_lab, upper_edge_lab)));
+
+ static_assert(std::is_same_v<hist_type_count_cms, hist_type_count_lab>);
+
+ using nucl_hist_type = decltype(boost::histogram::make_histogram(
+ boost::histogram::axis::regular<double, boost::histogram::axis::transform::log>(
+ num_bins_lab, lower_edge_lab, upper_edge_lab)));
+
+ TCountedProcess& process_;
+ hist_type_count_cms interaction_histogram_cms_;
+ hist_type_count_lab interaction_histogram_lab_;
+
+ std::map<int32_t, nucl_hist_type> nuclIntHistLab, nuclIntHistCMS;
+
+ public:
+ InteractionCounter(TCountedProcess& process)
+ : process_(process)
+ , 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))} {}
+
+ template <typename TProjectile>
+ auto DoInteraction(TProjectile& projectile) {
+ using namespace units::si;
+ auto constexpr massAir = 15.994_GeV * (1.f - 0.7847) + 14.003_GeV * .7847;
+ auto const sqrtS = sqrt(projectile.GetMass() * projectile.GetMass() +
+ massAir * massAir + 2 * projectile.GetEnergy() * massAir);
+
+ if (projectile.GetPID() == corsika::particles::Code::Nucleus) {
+ std::cerr << "NUCLEUS " << projectile.GetEnergy() << " "
+ << projectile.GetNuclearA() << " " << projectile.GetNuclearZ()
+ << std::endl;
+ int const A = projectile.GetNuclearA();
+ int const Z = projectile.GetNuclearZ();
+ int32_t pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
+
+ if (nuclIntHistCMS.count(pdg) == 0) {
+
+ nuclIntHistLab.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)));
+
+ nuclIntHistCMS.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)));
+ }
+
+ nuclIntHistLab[pdg](projectile.GetEnergy() / 1_GeV);
+ nuclIntHistCMS[pdg](sqrtS / 1_GeV);
+ } else {
+ interaction_histogram_cms_(
+ static_cast<corsika::particles::CodeIntType>(projectile.GetPID()),
+ sqrtS / 1_GeV);
+
+ interaction_histogram_lab_(
+ static_cast<corsika::particles::CodeIntType>(projectile.GetPID()),
+ projectile.GetEnergy() / 1_GeV);
+ }
+
+ return process_.DoInteraction(projectile);
+ }
+
+ void Init() { process_.Init(); }
+
+ template <typename TParticle>
+ auto GetInteractionLength(TParticle const& particle) const {
+ return process_.GetInteractionLength(particle);
+ }
+
+ auto CMSHists() const {
+ return std::make_pair(&interaction_histogram_cms_, &nuclIntHistCMS);
+ }
+
+ auto labHists() const {
+ return std::make_pair(&interaction_histogram_lab_, &nuclIntHistLab);
+ }
+ };
+
+ template <typename T1, typename T2>
+ static void saveHist(T1 const& hist, T2 const& histMap, std::string const& filename,
+ std::string const& comment = "") {
+ auto const& energy_axis = hist.axis(1);
+ std::ofstream myfile;
+ myfile.open(filename);
+ myfile << "# 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) {
+ myfile << "# " << static_cast<particles::Code>(p) << std::endl;
+ myfile << pdg;
+ for (int i = 0; i < energy_axis.size(); ++i) { myfile << ' ' << hist.at(p, i); }
+ myfile << std::endl;
+ }
+ }
+
+ myfile << "# nuclei" << std::endl;
+ for (auto const& [pdg, hist] : histMap) {
+ auto const num_ebins_nucl = hist.axis(0).size();
+ assert(energy_axis.size() == num_ebins_nucl);
+
+ myfile << pdg << " ";
+ for (int i = 0; i < num_ebins_nucl; ++i) { myfile << ' ' << hist.at(i); }
+ myfile << std::endl;
+ }
+ }
+} // namespace corsika::process::interaction_counter
+#endif
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&);