From d4fabf2c0c834793ace0886b85500d5abf4ef71b Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 13 Dec 2021 12:52:25 +0100
Subject: [PATCH] finalize MR
---
.../framework/core/ParticleProperties.inl | 4 +-
corsika/detail/modules/ObservationPlane.inl | 27 ++-
corsika/detail/modules/ParticleCut.inl | 40 ++++-
corsika/detail/modules/TrackWriter.inl | 5 +-
.../modules/energy_loss/BetheBlochPDG.inl | 33 ++--
.../modules/proposal/ContinuousProcess.inl | 4 +-
.../detail/modules/proposal/Interaction.inl | 2 +-
.../writers/BetheBlochPDGWriterParquet.inl | 160 ------------------
.../modules/writers/EnergyLossWriter.inl | 2 +-
.../writers/EnergyLossWriterParquet.inl | 2 +-
.../writers/ParticleCutWriterParquet.inl | 64 -------
.../modules/writers/TrackWriterParquet.inl | 6 +-
corsika/detail/setup/SetupStack.hpp | 31 +++-
corsika/framework/core/ParticleProperties.hpp | 4 +-
corsika/modules/ObservationPlane.hpp | 45 +++--
corsika/modules/ParticleCut.hpp | 75 +++++---
corsika/modules/TrackWriter.hpp | 17 +-
corsika/modules/conex/CONEXhybrid.hpp | 7 +-
corsika/modules/energy_loss/BetheBlochPDG.hpp | 6 +-
.../writers/BetheBlochPDGWriterParquet.hpp | 66 --------
corsika/modules/writers/EnergyLossWriter.hpp | 2 +-
.../writers/EnergyLossWriterParquet.hpp | 2 +-
.../writers/ParticleCutWriterParquet.hpp | 61 -------
.../modules/writers/TrackWriterParquet.hpp | 6 +-
corsika/setup/SetupStack.hpp | 4 +-
corsika/stack/WeightStackExtension.hpp | 4 +-
examples/CMakeLists.txt | 6 +-
examples/cascade_example.cpp | 17 +-
examples/cascade_proton_example.cpp | 14 +-
examples/corsika.cpp | 46 +++--
examples/em_shower.cpp | 47 +++--
examples/hybrid_MC.cpp | 24 ++-
examples/mars.cpp | 27 +--
examples/stopping_power.cpp | 9 +-
examples/vertical_EAS.cpp | 24 +--
tests/framework/testParticles.cpp | 9 +-
tests/modules/CMakeLists.txt | 2 +-
tests/modules/testCONEX.cpp | 2 +-
tests/modules/testParticleCut.cpp | 20 +--
tests/output/testWriterTracks.cpp | 3 +-
40 files changed, 338 insertions(+), 591 deletions(-)
delete mode 100644 corsika/detail/modules/writers/BetheBlochPDGWriterParquet.inl
delete mode 100644 corsika/detail/modules/writers/ParticleCutWriterParquet.inl
delete mode 100644 corsika/modules/writers/BetheBlochPDGWriterParquet.hpp
delete mode 100644 corsika/modules/writers/ParticleCutWriterParquet.hpp
diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 67cf78f64..fc9fbf636 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -14,7 +14,7 @@
namespace corsika {
- inline HEPEnergyType constexpr calculate_kinetic_energy_threshold(Code const code) {
+ inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const code) {
if (is_nucleus(code)) return particle::detail::threshold_nuclei;
return particle::detail::thresholds[static_cast<CodeIntType>(code)];
}
@@ -32,6 +32,8 @@ namespace corsika {
return particle::detail::masses[static_cast<CodeIntType>(code)];
}
+ inline bool constexpr is_charged(Code const c) { return get_charge_number(c) != 0; }
+
inline bool constexpr is_nucleus(Code const code) { return code >= Code::Nucleus; }
inline Code constexpr get_nucleus_code(size_t const A,
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 8fa3302e5..eba312b29 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -14,13 +14,13 @@ namespace corsika {
DirectionVector const& x_axis,
bool const deleteOnHit,
TArgs&&... args)
- : plane_(obsPlane)
- , deleteOnHit_(deleteOnHit)
- , energy_ground_(0_GeV)
- , count_ground_(0)
+ : TOutput(std::forward<TArgs>(args)...)
+ , plane_(obsPlane)
, xAxis_(x_axis.normalized())
, yAxis_(obsPlane.getNormal().cross(xAxis_))
- , TOutput(std::forward<TArgs>(args)...) {}
+ , deleteOnHit_(deleteOnHit)
+ , energy_ground_(0_GeV)
+ , count_ground_(0) {}
template <typename TTracking, typename TOutput>
template <typename TParticle, typename TTrajectory>
@@ -52,17 +52,10 @@ namespace corsika {
Point const pointOfIntersection = step.getPosition(1);
Vector const displacement = pointOfIntersection - plane_.getCenter();
- double const weight = 1.0;
- Code const pid = particle.getPID();
- if (pid == Code::Nucleus) {
- // add our particles to the output file stream
- this->write(particle.getNuclearA(), particle.getNuclearZ(), energy,
- displacement.dot(xAxis_), displacement.dot(yAxis_), 0_m, weight);
- } else {
- // add our particles to the output file stream
- this->write(particle.getPID(), energy, displacement.dot(xAxis_),
- displacement.dot(yAxis_), 0_m, weight);
- }
+ // add our particles to the output file stream
+ double const weight = 1.; // particle.getWeight()
+ this->write(particle.getPID(), energy, displacement.dot(xAxis_),
+ displacement.dot(yAxis_), 0_m, weight);
CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
@@ -73,7 +66,7 @@ namespace corsika {
} else {
return ProcessReturn::Ok;
}
- }
+ } // namespace corsika
template <typename TTracking, typename TOutput>
template <typename TParticle, typename TTrajectory>
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 647228b7e..058c2799c 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -18,14 +18,15 @@ namespace corsika {
HEPEnergyType const ePhoCut,
HEPEnergyType const eHadCut,
HEPEnergyType const eMuCut, bool const inv,
- TArgs&&... args)
+ bool const em, TArgs&&... args)
: TOutput(std::forward<TArgs>(args)...)
, doCutInv_(inv)
+ , doCutEm_(em)
, energy_cut_(0_GeV)
, energy_timecut_(0_GeV)
, energy_invcut_(0_GeV)
- , em_count_(0)
, inv_count_(0)
+ , em_count_(0)
, energy_count_() {
for (auto p : get_all_particles()) {
if (is_hadron(p)) // nuclei are also hadrons
@@ -46,14 +47,38 @@ namespace corsika {
}
template <typename TOutput>
- inline ParticleCut<TOutput>::ParticleCut(
- std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const inv)
- : doCutInv_(inv)
+ template <typename... TArgs>
+ inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eCut, bool const inv,
+ bool const em, TArgs&&... args)
+ : TOutput(std::forward<TArgs>(args)...)
+ , doCutInv_(inv)
+ , doCutEm_(em)
, energy_cut_(0_GeV)
, energy_timecut_(0_GeV)
, energy_invcut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , inv_count_(0)
, em_count_(0)
+ , energy_count_() {
+ for (auto p : get_all_particles()) { set_kinetic_energy_threshold(p, eCut); }
+ set_kinetic_energy_threshold(Code::Nucleus, eCut);
+ CORSIKA_LOG_DEBUG("setting kinetic energy threshold {} GeV", eCut / 1_GeV);
+ }
+
+ template <typename TOutput>
+ template <typename... TArgs>
+ inline ParticleCut<TOutput>::ParticleCut(
+ std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const inv,
+ bool const em, TArgs&&... args)
+ : TOutput(std::forward<TArgs>(args)...)
+ , doCutInv_(inv)
+ , doCutEm_(em)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_invcut_(0_GeV)
+ , energy_emcut_(0_GeV)
, inv_count_(0)
+ , em_count_(0)
, energy_count_(0) {
set_kinetic_energy_thresholds(eCuts);
CORSIKA_LOG_DEBUG("setting threshold particles individually");
@@ -68,9 +93,9 @@ namespace corsika {
if (is_nucleus(pid)) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / get_nucleus_A(pid);
- return (ElabNuc < calculate_kinetic_energy_threshold(pid));
+ return (ElabNuc < get_kinetic_energy_threshold(pid));
} else {
- return (energyLab < calculate_kinetic_energy_threshold(pid));
+ return (energyLab < get_kinetic_energy_threshold(pid));
}
}
@@ -120,7 +145,6 @@ namespace corsika {
return false; // this particle will not be removed/cut
}
-
template <typename TOutput>
template <typename TStackView>
inline void ParticleCut<TOutput>::doSecondaries(TStackView& vS) {
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index e44119069..39bb395ac 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -15,7 +15,6 @@
namespace corsika {
-
template <typename TOutput>
inline TrackWriter<TOutput>::TrackWriter() {}
@@ -28,8 +27,8 @@ namespace corsika {
auto const end = vT.getPosition(1).getCoordinates();
// write the track to the file
- output_.write(vP.getPID(), vP.getEnergy(), vP.getWeight(), start, vP.getTime() - vT.getDuration(), end,
- vP.getTime());
+ TOutput::write(vP.getPID(), vP.getEnergy(), vP.getWeight(), start,
+ vP.getTime() - vT.getDuration(), end, vP.getTime());
return ProcessReturn::Ok;
}
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 8423d1082..360487286 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -19,19 +19,15 @@
namespace corsika {
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
-
template <typename TOutput>
template <typename... TArgs>
inline BetheBlochPDG<TOutput>::BetheBlochPDG(TArgs&&... args)
: TOutput(std::forward<TArgs>(args)...) {}
- template <typename TParticle>
template <typename TOutput>
- inline HEPEnergyType BetheBlochPDG<TOutput>::getBetheBloch(
- TParticle const& p, GrammageType const dX) {
+ template <typename TParticle>
+ inline HEPEnergyType BetheBlochPDG<TOutput>::getBetheBloch(TParticle const& p,
+ GrammageType const dX) {
// all these are material constants and have to come through Environment
// right now: values for nitrogen_D
@@ -124,8 +120,8 @@ namespace corsika {
}
// radiation losses according to PDG 2018, ch. 33 ref. [5]
- template <typename TParticle>
template <typename TOutput>
+ template <typename TParticle>
inline HEPEnergyType BetheBlochPDG<TOutput>::getRadiationLosses(
TParticle const& vP, GrammageType const vDX) {
// simple-minded hard-coded value for b(E) inspired by data from
@@ -134,17 +130,18 @@ namespace corsika {
return -vP.getEnergy() * b * vDX;
}
- template <typename TParticle>
template <typename TOutput>
+ template <typename TParticle>
inline HEPEnergyType BetheBlochPDG<TOutput>::getTotalEnergyLoss(
TParticle const& vP, GrammageType const vDX) {
return getBetheBloch(vP, vDX) + getRadiationLosses(vP, vDX);
}
- template <typename TParticle, typename TTrajectory>
template <typename TOutput>
- inline ProcessReturn BetheBlochPDG<TOutput>::doContinuous(
- TParticle& particle, TTrajectory const& track, bool const) {
+ template <typename TParticle, typename TTrajectory>
+ inline ProcessReturn BetheBlochPDG<TOutput>::doContinuous(TParticle& particle,
+ TTrajectory const& track,
+ bool const) {
// if this step was limiting the CORSIKA stepping, the particle is lost
/* see Issue https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/389
@@ -158,8 +155,7 @@ namespace corsika {
if (particle.getChargeNumber() == 0) return ProcessReturn::Ok;
GrammageType const dX =
- particle.getNode()->getModelProperties().getIntegratedGrammage(track,
- track.getLength());
+ particle.getNode()->getModelProperties().getIntegratedGrammage(track);
CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", particle.getPID(),
particle.getChargeNumber(), dX / 1_g * square(1_cm));
HEPEnergyType const dE = getTotalEnergyLoss(particle, dX);
@@ -173,10 +169,10 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TParticle, typename TTrajectory>
template <typename TOutput>
- inline LengthType BetheBlochPDG<TOutput>::getMaxStepLength(TParticle const& vParticle,
- TTrajectory const& vTrack) const {
+ template <typename TParticle, typename TTrajectory>
+ inline LengthType BetheBlochPDG<TOutput>::getMaxStepLength(
+ TParticle const& vParticle, TTrajectory const& vTrack) const {
if (vParticle.getChargeNumber() == 0) {
return meter * std::numeric_limits<double>::infinity();
}
@@ -186,7 +182,7 @@ namespace corsika {
auto const energy = vParticle.getKineticEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- calculate_kinetic_energy_threshold(
+ get_kinetic_energy_threshold(
vParticle.getPID()) // energy thresholds globally defined
// for individual particles
* 0.99999 // need to go slightly below global e-cut to assure
@@ -199,7 +195,6 @@ namespace corsika {
vTrack, maxGrammage);
}
-
template <typename TOutput>
inline void BetheBlochPDG<TOutput>::showResults() const {
CORSIKA_LOG_INFO("energy lost dE (GeV) : {} ", energy_lost_ / 1_GeV);
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 89f571b69..9c51b815d 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -29,7 +29,7 @@ namespace corsika::proposal {
// take some minutes if you have to build the tables and cannot read the
// from disk
auto const emCut =
- calculate_kinetic_energy_threshold(code) +
+ get_kinetic_energy_threshold(code) +
get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
@@ -123,7 +123,7 @@ namespace corsika::proposal {
auto const energy = vP.getEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- calculate_kinetic_energy_threshold(
+ get_kinetic_energy_threshold(
code) // energy thresholds globally defined for individual particles
* 0.9999 // need to go slightly below global e-cut to assure removal
// in ParticleCut. This does not matter since at cut-time
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index af6d603f6..d8916cbe9 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -33,7 +33,7 @@ namespace corsika::proposal {
// take some minutes if you have to build the tables and cannot read the
// from disk
auto const emCut =
- calculate_kinetic_energy_threshold(code) +
+ get_kinetic_energy_threshold(code) +
get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
diff --git a/corsika/detail/modules/writers/BetheBlochPDGWriterParquet.inl b/corsika/detail/modules/writers/BetheBlochPDGWriterParquet.inl
deleted file mode 100644
index bb3943005..000000000
--- a/corsika/detail/modules/writers/BetheBlochPDGWriterParquet.inl
+++ /dev/null
@@ -1,160 +0,0 @@
-/*
- * (c) Copyright 2021 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.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <string>
-
-namespace corsika {
-
- inline BetheBlochPDGWriterParquet::BetheBlochPDGWriterParquet()
- : output_()
- , showerId_(0.) {}
-
- inline void BetheBlochPDGWriterParquet::startOfLibrary(
- boost::filesystem::path const& directory) {
- // setup the streamer
- output_.initStreamer((directory / "energyloss.parquet").string());
-
- // enable compression with the default level
- // output_.enableCompression();
-
- // build the schema
- output_.addField("start_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("start_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("start_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("end_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("end_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("end_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
- parquet::ConvertedType::INT_32);
- output_.addField("dE", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
-
- // and build the streamer
- output_.buildStreamer();
- }
-
- inline void BetheBlochPDGWriterParquet::startOfShower(unsigned int const showerId) {
- showerId_ = showerId;
- }
-
- inline void BetheBlochPDGWriterParquet::endOfShower(unsigned int const showerId) {}
-
- inline void BetheBlochPDGWriterParquet::endOfLibrary() { output_.closeStreamer(); }
-
- template <typename TTrack>
- inline void BetheBlochPDGWriterParquet::write(TTrack const& track, Code const pid,
- HEPEnergyType dE) {
-
- auto const start = track.getPosition(0).getCoordinates();
- auto const end = track.getPosition(1).getCoordinates();
-
- // clang-format off
- *(output_.getWriter()) << showerId_
- << static_cast<float>(start[0] / 1_m)
- << static_cast<float>(start[1] / 1_m)
- << static_cast<float>(start[2] / 1_m)
- << static_cast<float>(end[0] / 1_m)
- << static_cast<float>(end[1] / 1_m)
- << static_cast<float>(end[2] / 1_m)
- << static_cast<int>(get_PDG(pid))
- << static_cast<float>(dE / 1_GeV)
- << parquet::EndRow;
- // clang-format on
- }
-
- // inline void BetheBlochPDGWriterParquet::write(GrammageType const Xstart,
- // GrammageType const Xend, Code const
- // pid, double const weight) {
- // double const bstart = Xstart / dX_;
- // double const bend = Xend / dX_;
-
- // if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
- // abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
- // CORSIKA_LOGGER_ERROR(getLogger(),
- // "CONEX and Corsika8 dX grammage binning are not the same! "
- // "Xstart={} Xend={} dX={}",
- // Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
- // dX_ / 1_g * square(1_cm));
- // throw std::runtime_error(
- // "CONEX and Corsika8 dX grammage binning are not the same!");
- // }
-
- // int const bin = int((bend + bstart) / 2);
-
- // if (pid == Code::Photon) {
- // profile_.at(bin)[static_cast<int>(ProfileIndex::Photon)] += weight;
- // } else if (pid == Code::Positron) {
- // profile_.at(bin)[static_cast<int>(ProfileIndex::Positron)] += weight;
- // } else if (pid == Code::Electron) {
- // profile_.at(bin)[static_cast<int>(ProfileIndex::Electron)] += weight;
- // } else if (pid == Code::MuPlus) {
- // profile_.at(bin)[static_cast<int>(ProfileIndex::MuPlus)] += weight;
- // } else if (pid == Code::MuMinus) {
- // profile_.at(bin)[static_cast<int>(ProfileIndex::MuMinus)] += weight;
- // } else if (is_hadron(pid)) {
- // profile_.at(bin)[static_cast<int>(ProfileIndex::Hadron)] += weight;
- // }
- // if (is_charged(pid)) {
- // profile_[bin][static_cast<int>(ProfileIndex::Charged)] += weight;
- // }
- // }
-
- // inline YAML::Node BetheBlochPDGWriterParquet::getSummary() const {
- // // determined Xmax and dEdXmax from quadratic interpolation
-
- // YAML::Node summary;
-
- // for (int index = 0; index < static_cast<int>(ProfileIndex::Entries); ++index) {
- // // first find highest 3-boxcar sum
- // double maximum = 0;
- // unsigned int iMaximum = 0;
- // for (unsigned int i = 0; i < profile_.size() - 3; ++i) {
- // double value = profile_[i + 0].at(index) + profile_[i + 1].at(index) +
- // profile_[i + 2].at(index);
- // if (value > maximum) {
- // maximum = value;
- // iMaximum = i;
- // }
- // }
-
- // double const dX = dX_ / 1_g * square(1_cm);
-
- // // quadratic interpolation of maximum in 3 highest points
- // auto [Xmax, Nmax] = FindXmax::interpolateProfile(
- // dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum),
- // profile_[iMaximum + 0].at(index), profile_[iMaximum + 1].at(index),
- // profile_[iMaximum + 2].at(index));
-
- // std::string const name = ProfileIndexNames[index];
- // summary[name]["Xmax"] = Xmax;
- // summary[name]["Nmax"] = Nmax;
- // }
- // return summary;
- // }
-
- // inline YAML::Node BetheBlochPDGWriterParquet::getConfig() const {
- // // determined Xmax and dEdXmax from quadratic interpolation
-
- // YAML::Node node;
- // node["type"] = "BetheBlochPDG";
-
- // return node;
- // }
-
-} // namespace corsika
diff --git a/corsika/detail/modules/writers/EnergyLossWriter.inl b/corsika/detail/modules/writers/EnergyLossWriter.inl
index 1fbf3b5ba..29f07e5fc 100644
--- a/corsika/detail/modules/writers/EnergyLossWriter.inl
+++ b/corsika/detail/modules/writers/EnergyLossWriter.inl
@@ -178,7 +178,7 @@ namespace corsika {
node["nbins"] = nBins_;
node["grammage_threshold"] = dX_threshold_ / (1_g / square(1_cm));
- // TODO: add shower axis to config
+ //! \todo add shower axis to config
return node;
}
diff --git a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
index b7e82d669..1680e421a 100644
--- a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
+++ b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl
@@ -52,7 +52,7 @@ namespace corsika {
inline void EnergyLossWriterParquet::startOfShower(unsigned int const) {}
- inline void EnergyLossWriterParquet::endOfShower(unsigned int const showerId) {}
+ inline void EnergyLossWriterParquet::endOfShower(unsigned int const) {}
inline void EnergyLossWriterParquet::endOfLibrary() { output_.closeStreamer(); }
diff --git a/corsika/detail/modules/writers/ParticleCutWriterParquet.inl b/corsika/detail/modules/writers/ParticleCutWriterParquet.inl
deleted file mode 100644
index 1405e9416..000000000
--- a/corsika/detail/modules/writers/ParticleCutWriterParquet.inl
+++ /dev/null
@@ -1,64 +0,0 @@
-/*
- * (c) Copyright 2021 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.
- */
-
-#pragma once
-
-namespace corsika {
-
- inline ParticleCutWriterParquet::ParticleCutWriterParquet()
- : output_() {}
-
- inline void ParticleCutWriterParquet::startOfLibrary(
- boost::filesystem::path const& directory) {
-
- // setup the streamer
- output_.initStreamer((directory / "energyloss.parquet").string());
-
- // build the schema
- output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
- parquet::ConvertedType::INT_32);
- output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
- output_.addField("z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
- parquet::ConvertedType::NONE);
-
- // and build the streamer
- output_.buildStreamer();
- }
-
- inline void ParticleCutWriterParquet::startOfShower(unsigned int const showerId) {
- showerId_ = showerId;
- }
-
- inline void ParticleCutWriterParquet::endOfShower(unsigned int const) {}
-
- inline void ParticleCutWriterParquet::endOfLibrary() { output_.closeStreamer(); }
-
- inline void ParticleCutWriterParquet::write(Point const& point, Code const pid,
- HEPEnergyType energy) {
-
- auto location{point.getCoordinates()};
-
- // write the next row - we must write `shower_` first.
- // clang-format off
- *(output_.getWriter())
- << showerId_
- << static_cast<int>(get_PDG(pid))
- << static_cast<float>(energy / 1_GeV)
- << static_cast<float>(location[0] / 1_m)
- << static_cast<float>(location[1] / 1_m)
- << static_cast<float>(location[2] / 1_m)
- << parquet::EndRow;
- // clang-format on
- }
-
-} // namespace corsika
diff --git a/corsika/detail/modules/writers/TrackWriterParquet.inl b/corsika/detail/modules/writers/TrackWriterParquet.inl
index f233ca318..88dfa8e03 100644
--- a/corsika/detail/modules/writers/TrackWriterParquet.inl
+++ b/corsika/detail/modules/writers/TrackWriterParquet.inl
@@ -57,12 +57,12 @@ namespace corsika {
inline void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); }
- inline void TrackWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
+ inline void TrackWriterParquet::write(Code const pid, HEPEnergyType const energy,
double const weight,
QuantityVector<length_d> const& start,
- TimeType const& t_start,
+ TimeType const t_start,
QuantityVector<length_d> const& end,
- TimeType const& t_end) {
+ TimeType const t_end) {
// write the next row - we must write `shower_` first.
// clang-format off
diff --git a/corsika/detail/setup/SetupStack.hpp b/corsika/detail/setup/SetupStack.hpp
index 7a1327b2e..2b40ccbdb 100644
--- a/corsika/detail/setup/SetupStack.hpp
+++ b/corsika/detail/setup/SetupStack.hpp
@@ -22,7 +22,8 @@ namespace corsika {
namespace setup::detail {
// ------------------------------------------
- // add geometry node tracking data to stack:
+ // add geometry node data to stack. This is fundamentally needed
+ // for robust tracking through multiple volumes.
// the GeometryNode stack needs to know the type of geometry-nodes from the
// environment:
@@ -42,16 +43,36 @@ namespace corsika {
DefaultSecondaryProducer>;
// ------------------------------------------
- // Add [optional] history data to stack, too:
+ // add weight data to stack. This is fundamentally needed
+ // for thinning.
+
+ // the "pure" weight stack (interface)
+ template <typename TStackIter>
+ using SetupWeightDataInterface =
+ typename weights::MakeWeightDataInterface<TStackIter>::type;
+
+ // combine geometry-node-vector data stack with weight information for tracking
+ template <typename TStackIter>
+ using StackWithWeightInterface =
+ CombinedParticleInterface<StackWithGeometry::pi_type, SetupWeightDataInterface,
+ TStackIter>;
+
+ // the combined stack data: particle + geometry + weight
+ using StackWithWeight =
+ CombinedStack<typename StackWithGeometry::stack_data_type, weights::WeightData,
+ StackWithWeightInterface, DefaultSecondaryProducer>;
+
+ // ------------------------------------------
+ // Add [OPTIONAL] history data to stack, too.
+ // This keeps the entire lineage of particles in memory.
- // combine dummy stack with geometry information for tracking
template <typename TStackIter>
using StackWithHistoryInterface =
- CombinedParticleInterface<StackWithGeometry::pi_type,
+ CombinedParticleInterface<StackWithWeight::pi_type,
history::HistoryEventDataInterface, TStackIter>;
using StackWithHistory =
- CombinedStack<typename StackWithGeometry::stack_data_type,
+ CombinedStack<typename StackWithWeight::stack_data_type,
history::HistoryEventData, StackWithHistoryInterface,
history::HistorySecondaryProducer>;
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index 267c67838..dd03aa106 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -92,7 +92,7 @@ namespace corsika {
int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e
ElectricChargeType constexpr get_charge(Code const); //!< electric charge
HEPMassType constexpr get_mass(Code const); //!< mass
- HEPEnergyType constexpr calculate_kinetic_energy_threshold(
+ HEPEnergyType constexpr get_kinetic_energy_threshold(
Code const); //!< get kinetic energy threshold below which the particle is
//!< discarded, by default set to zero
void constexpr set_kinetic_energy_threshold(
@@ -117,7 +117,7 @@ namespace corsika {
bool constexpr is_em(Code const); //!< true if particle is electron, positron or photon
bool constexpr is_muon(Code const); //!< true if particle is mu+ or mu-
bool constexpr is_neutrino(Code const); //!< true if particle is (anti-) neutrino
- bool constexpr is_charged(Code const); //!< true if particle is charged
+ bool constexpr is_charged(Code const); //!< true if particle is charged
/**
* @brief Creates the Code for a nucleus of type 10LZZZAAAI.
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index 43b9d74be..ae8f72af9 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -17,28 +17,39 @@
namespace corsika {
/**
- @ingroup Modules
- @{
-
- The ObservationPlane writes PDG codes, energies, and distances of particles to the
- central point of the plane into its output file. The particles are considered
- "absorbed" afterwards.
-
- **Note/Limitation:** as discussed in
- https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397
- you cannot put two ObservationPlanes exactly on top of each
- other. Even if one of them is "permeable". You have to put a
- small gap in between the two plane in such a scenario, or develop
- another more specialized output class.
+ * @ingroup Modules
+ * @{
+ *
+ * The ObservationPlane writes PDG codes, energies, and distances of particles to the
+ * central point of the plane into its output file. The particles are considered
+ * "absorbed" afterwards.
+ *
+ * **Note/Limitation:** as discussed in
+ * https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397
+ * you cannot put two ObservationPlanes exactly on top of each
+ * other. Even if one of them is "permeable". You have to put a
+ * small gap in between the two plane in such a scenario, or develop
+ * another more specialized output class.
*/
template <typename TTracking, typename TOutput = WriterOff>
class ObservationPlane : public ContinuousProcess<ObservationPlane<TTracking, TOutput>>,
public TOutput {
+ using TOutput::write;
+
public:
+ /**
+ * Construct a new Observation Plane object.
+ *
+ * @tparam TArgs
+ * @param plane The plane.
+ * @param x_dir The x-direction/axis.
+ * @param absorbing Flag to make the plane absorbing.
+ * @param args
+ */
template <typename... TArgs>
- ObservationPlane(Plane const&, DirectionVector const&, bool const = true,
- TArgs&&... args);
+ ObservationPlane(Plane const& plane, DirectionVector const& x_dir,
+ bool const absorbing = true, TArgs&&... args);
~ObservationPlane() {}
@@ -56,11 +67,11 @@ namespace corsika {
private:
Plane const plane_;
+ DirectionVector const xAxis_;
+ DirectionVector const yAxis_;
bool const deleteOnHit_;
HEPEnergyType energy_ground_;
unsigned int count_ground_;
- DirectionVector const xAxis_;
- DirectionVector const yAxis_;
};
//! @}
} // namespace corsika
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 7ce7139ed..5727ab3cb 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -19,13 +19,15 @@
namespace corsika {
/**
- simple ParticleCut process. Goes through the secondaries of an interaction and
- removes particles according to their kinetic energy. Particles with a time delay of
- more than 10ms are removed as well. Invisible particles (neutrinos) can be removed if
- selected. The threshold value is set to 0 by default but in principle can be configured
- for each particle. Special constructors for cuts by the following groups are
- implemented: (electrons,positrons), photons, hadrons and muons.
- **/
+ * ParticleCut process to kill particles.
+ *
+ * Goes through the secondaries of an interaction and
+ * removes particles according to their kinetic energy. Particles with a time delay of
+ * more than 10ms are removed as well. Invisible particles (neutrinos) can be removed if
+ * selected. The threshold value is set to 0 by default but in principle can be
+ * configured for each particle. Special constructors for cuts by the following groups
+ * are implemented: (electrons,positrons), photons, hadrons and muons.
+ */
template <typename TOutput = WriterOff>
class ParticleCut : public SecondariesProcess<ParticleCut<TOutput>>,
public ContinuousProcess<ParticleCut<TOutput>>,
@@ -34,29 +36,58 @@ namespace corsika {
public:
/**
* particle cut with kinetic energy thresholds for electrons, photons,
- * hadrons (including nuclei with energy per nucleon) and muons
- * invisible particles (neutrinos) can be cut or not
+ * hadrons (including nuclei with energy per nucleon) and muons
+ * invisible particles (neutrinos) can be cut or not.
*/
template <typename... TArgs>
ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv,
- TArgs&&... args);
+ bool const em, TArgs&&... args);
+
+ /**
+ * particle cut with kinetic energy thresholds for all particles.
+ */
+ template <typename... TArgs>
+ ParticleCut(HEPEnergyType const eCut, bool const inv, bool const em, TArgs&&... args);
/**
* Threshold for specific particles redefined. EM and invisible particles can be set
* to be discarded altogether.
*/
+ template <typename... TArgs>
ParticleCut(std::unordered_map<Code const, HEPEnergyType const> const& eCuts,
- bool const inv);
+ bool const inv, bool const em, TArgs&&... args);
+ /**
+ * Cut particles which are secondaries from discrete processes.
+ *
+ * @tparam TStackView
+ */
template <typename TStackView>
void doSecondaries(TStackView&);
+ /**
+ * Cut particles during contunuous processes (energy losses etc).
+ *
+ * @tparam TParticle
+ * @tparam TTrajectory
+ * @param vParticle
+ * @param vTrajectory
+ * @param limitFlag
+ * @return ProcessReturn
+ */
template <typename TParticle, typename TTrajectory>
ProcessReturn doContinuous(
TParticle& vParticle, TTrajectory const& vTrajectory,
const bool limitFlag = false); // this is not used for ParticleCut
+ /**
+ * Limit on continuous step length imposed by ParticleCut: none.
+ *
+ * @tparam TParticle
+ * @tparam TTrajectory
+ * @return LengthType
+ */
template <typename TParticle, typename TTrajectory>
LengthType getMaxStepLength(TParticle const&, TTrajectory const&) {
return meter * std::numeric_limits<double>::infinity();
@@ -67,27 +98,30 @@ namespace corsika {
void reset();
HEPEnergyType getElectronKineticECut() const {
- return calculate_kinetic_energy_threshold(Code::Electron);
+ return get_kinetic_energy_threshold(Code::Electron);
}
HEPEnergyType getPhotonKineticECut() const {
- return calculate_kinetic_energy_threshold(Code::Photon);
+ return get_kinetic_energy_threshold(Code::Photon);
}
HEPEnergyType getMuonKineticECut() const {
- return calculate_kinetic_energy_threshold(Code::MuPlus);
+ return get_kinetic_energy_threshold(Code::MuPlus);
}
HEPEnergyType getHadronKineticECut() const {
- return calculate_kinetic_energy_threshold(Code::Proton);
+ return get_kinetic_energy_threshold(Code::Proton);
}
//! returns total energy of particles that were removed by cut for invisible particles
HEPEnergyType getInvEnergy() const { return energy_invcut_; }
+ //! returns total energy of particles that were removed by cut for invisible particles
+ HEPEnergyType getEmEnergy() const { return energy_emcut_; }
//! returns total energy of particles that were removed by cut in time
HEPEnergyType getTimeCutEnergy() const { return energy_timecut_; }
//! returns total energy of particles that were removed by cut in kinetic energy
HEPEnergyType getCutEnergy() const { return energy_cut_; }
//! returns number of invisible particles
- unsigned int getNumberInvParticles() const { return inv_count_; }
+ size_t getNumberInvParticles() const { return inv_count_; }
+ size_t getNumberEmParticles() const { return em_count_; }
- // get configuration of this node
+ //! get configuration of this node, for output
YAML::Node getConfig() const override;
private:
@@ -102,11 +136,14 @@ namespace corsika {
private:
bool doCutInv_;
+ bool doCutEm_;
HEPEnergyType energy_cut_ = 0 * electronvolt;
HEPEnergyType energy_timecut_ = 0 * electronvolt;
HEPEnergyType energy_invcut_ = 0 * electronvolt;
- unsigned int inv_count_ = 0;
- unsigned int energy_count_ = 0;
+ HEPEnergyType energy_emcut_ = 0 * electronvolt;
+ size_t inv_count_ = 0;
+ size_t em_count_ = 0;
+ size_t energy_count_ = 0;
HEPEnergyType energy_event_; // per event sum
};
diff --git a/corsika/modules/TrackWriter.hpp b/corsika/modules/TrackWriter.hpp
index 4d1a956bb..646688411 100644
--- a/corsika/modules/TrackWriter.hpp
+++ b/corsika/modules/TrackWriter.hpp
@@ -10,12 +10,22 @@
#include <corsika/framework/process/ContinuousProcess.hpp>
#include <corsika/modules/writers/TrackWriterParquet.hpp>
+#include <corsika/modules/writers/WriterOff.hpp>
namespace corsika {
- template <typename TOutput = WriterOff>
- class TrackWriter : public ContinuousProcess<TrackWriter<TOutput>>,
- public TOutput {
+ /**
+ * To write 3D track data to disk.
+ *
+ * Since the only sole purpose of this module is to generate track
+ * output on disk, the default output mode is not "WriterOff" but
+ * directly TrackWriterParquet. It can of course be changed.
+ *
+ * @tparam TOutput with default TrackWriterParquet
+ */
+
+ template <typename TOutput = TrackWriterParquet>
+ class TrackWriter : public ContinuousProcess<TrackWriter<TOutput>>, public TOutput {
public:
TrackWriter();
@@ -29,7 +39,6 @@ namespace corsika {
YAML::Node getConfig() const;
private:
-
};
} // namespace corsika
diff --git a/corsika/modules/conex/CONEXhybrid.hpp b/corsika/modules/conex/CONEXhybrid.hpp
index a4942727f..990189b37 100644
--- a/corsika/modules/conex/CONEXhybrid.hpp
+++ b/corsika/modules/conex/CONEXhybrid.hpp
@@ -30,8 +30,9 @@ namespace corsika {
} // namespace conex
template <typename TOutput = WriterOff, typename TProfileOutput = WriterOff>
- class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid<TOutput, TProfileOutput>>,
- public SecondariesProcess<CONEXhybrid<TOutput, TProfileOutput>> {
+ class CONEXhybrid
+ : public CascadeEquationsProcess<CONEXhybrid<TOutput, TProfileOutput>>,
+ public SecondariesProcess<CONEXhybrid<TOutput, TProfileOutput>> {
public:
/**
@@ -78,7 +79,7 @@ namespace corsika {
void initCascadeEquations();
/**
- * Cascade equations are solved basoned on the data in the tables
+ * Cascade equations are solved basoned on the data in the tables.
*/
template <typename TStack>
void doCascadeEquations(TStack& stack);
diff --git a/corsika/modules/energy_loss/BetheBlochPDG.hpp b/corsika/modules/energy_loss/BetheBlochPDG.hpp
index b0c847ce5..2add0e2d9 100644
--- a/corsika/modules/energy_loss/BetheBlochPDG.hpp
+++ b/corsika/modules/energy_loss/BetheBlochPDG.hpp
@@ -13,7 +13,6 @@
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
-
#include <corsika/modules/writers/WriterOff.hpp>
#include <map>
@@ -58,9 +57,8 @@ namespace corsika {
bool const limitFlag);
template <typename TParticle, typename TTrajectory>
- LengthType getMaxStepLength(TParticle const&,
- TTrajectory const&);
-
+ LengthType getMaxStepLength(TParticle const&, TTrajectory const&) const;
+
template <typename TParticle>
static HEPEnergyType getBetheBloch(TParticle const&, const GrammageType);
diff --git a/corsika/modules/writers/BetheBlochPDGWriterParquet.hpp b/corsika/modules/writers/BetheBlochPDGWriterParquet.hpp
deleted file mode 100644
index 1f7003c99..000000000
--- a/corsika/modules/writers/BetheBlochPDGWriterParquet.hpp
+++ /dev/null
@@ -1,66 +0,0 @@
-/*
- * (c) Copyright 2021 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.
- */
-
-#pragma once
-
-#include <corsika/output/BaseOutput.hpp>
-#include <corsika/output/ParquetStreamer.hpp>
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/geometry/QuantityVector.hpp>
-
-namespace corsika {
-
- class BetheBlochPDGWriterParquet : public BaseOutput {
-
- public:
- /**
- * Construct a new writer.
- *
- * @param name The name of this output.
- */
- BetheBlochPDGWriterParquet();
-
- /**
- * Called at the start of each library.
- */
- void startOfLibrary(boost::filesystem::path const& directory) final override;
-
- /**
- * Called at the start of each shower.
- */
- void startOfShower(unsigned int const showerId) final override;
-
- /**
- * Called at the end of each shower.
- */
- void endOfShower(unsigned int const showerId) final override;
-
- /**
- * Called at the end of each library.
- *
- * This must also increment the run number since we override
- * the default behaviour of BaseOutput.
- */
- void endOfLibrary() final override;
-
- /**
- * Write a track to the file.
- */
- template <typename TTrack>
- void write(TTrack const& track, Code const pid, HEPEnergyType dE);
-
- private:
- ParquetStreamer output_; ///< The primary output file.
- unsigned int showerId_; ///< event Id counter
-
- }; // class BetheBlochPDGWriterParquet
-
-} // namespace corsika
-
-#include <corsika/detail/modules/writers/BetheBlochPDGWriterParquet.inl>
diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp
index 8806acd8c..50175b49d 100644
--- a/corsika/modules/writers/EnergyLossWriter.hpp
+++ b/corsika/modules/writers/EnergyLossWriter.hpp
@@ -81,7 +81,7 @@ namespace corsika {
private:
ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage
GrammageType dX_; ///< binning of profile.
- unsigned int nBins_; ///< number of profile bins.
+ size_t nBins_; ///< number of profile bins.
GrammageType dX_threshold_; ///< too short tracks are discarded.
std::vector<Profile> profile_; // longitudinal profile
diff --git a/corsika/modules/writers/EnergyLossWriterParquet.hpp b/corsika/modules/writers/EnergyLossWriterParquet.hpp
index 5f2cbac28..3edc4cdd6 100644
--- a/corsika/modules/writers/EnergyLossWriterParquet.hpp
+++ b/corsika/modules/writers/EnergyLossWriterParquet.hpp
@@ -51,7 +51,7 @@ namespace corsika {
void endOfLibrary() override;
/**
- * Write energy lost to the file
+ * Write energy lost to the file.
*/
void write(unsigned int const showerId, GrammageType const grammage,
HEPEnergyType const total);
diff --git a/corsika/modules/writers/ParticleCutWriterParquet.hpp b/corsika/modules/writers/ParticleCutWriterParquet.hpp
deleted file mode 100644
index 466bf5571..000000000
--- a/corsika/modules/writers/ParticleCutWriterParquet.hpp
+++ /dev/null
@@ -1,61 +0,0 @@
-/*
- * (c) Copyright 2021 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.
- */
-
-#pragma once
-
-#include <corsika/output/BaseOutput.hpp>
-#include <corsika/output/ParquetStreamer.hpp>
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-namespace corsika {
-
- class ParticleCutWriterParquet : public BaseOutput {
-
- public:
- /**
- * Construct a new writer.
- */
- ParticleCutWriterParquet();
-
- /**
- * Called at the start of each library.
- */
- void startOfLibrary(boost::filesystem::path const& directory) final override;
-
- /**
- * Called at the start of each shower.
- */
- void startOfShower(unsigned int const showerId) final override;
-
- /**
- * Called at the end of each shower.
- */
- void endOfShower(unsigned int const showerId) final override;
-
- /**
- * Called at the end of each library.
- *
- * This must also increment the run number since we override
- * the default behaviour of BaseOutput.
- */
- void endOfLibrary() final override;
-
- /**
- * Write a particle cut.
- */
- void write(Point const& point, Code const PID, HEPEnergyType const energy);
-
- private:
- ParquetStreamer output_; ///< The parquet streamer for this process.
- unsigned int showerId_{0};
- }; // namespace corsika
-
-} // namespace corsika
-
-#include <corsika/detail/modules/writers/ParticleCutWriterParquet.inl>
diff --git a/corsika/modules/writers/TrackWriterParquet.hpp b/corsika/modules/writers/TrackWriterParquet.hpp
index 91fea4e54..2d7ac7b28 100644
--- a/corsika/modules/writers/TrackWriterParquet.hpp
+++ b/corsika/modules/writers/TrackWriterParquet.hpp
@@ -52,9 +52,9 @@ namespace corsika {
/**
* Write a track to the file.
*/
- void write(Code const& pid, units::si::HEPEnergyType const& energy,
- QuantityVector<length_d> const& start, TimeType const& t_start,
- QuantityVector<length_d> const& end, TimeType const& t_end);
+ void write(Code const pid, HEPEnergyType const energy, double const weight,
+ QuantityVector<length_d> const& start, TimeType const t_start,
+ QuantityVector<length_d> const& end, TimeType const t_end);
private:
ParquetStreamer output_; ///< The primary output file.
diff --git a/corsika/setup/SetupStack.hpp b/corsika/setup/SetupStack.hpp
index 305089ec5..9028c2ab1 100644
--- a/corsika/setup/SetupStack.hpp
+++ b/corsika/setup/SetupStack.hpp
@@ -29,9 +29,9 @@ namespace corsika::setup {
#else // WITH_HISTORY
/*
- * the version without history
+ * the version without history (and geometry data and weights)
*/
- using Stack = detail::StackWithGeometry;
+ using Stack = detail::StackWithWeight;
#endif
diff --git a/corsika/stack/WeightStackExtension.hpp b/corsika/stack/WeightStackExtension.hpp
index b2e6081f2..bbe636db1 100644
--- a/corsika/stack/WeightStackExtension.hpp
+++ b/corsika/stack/WeightStackExtension.hpp
@@ -23,7 +23,7 @@ namespace corsika::weights {
* Corresponding defintion of a stack-readout object, the iteractor
* dereference operator will deliver access to these function
* defintion of a stack-readout object, the iteractor dereference
- * operator will deliver access to these function
+ * operator will deliver access to these function.
*/
/**
@@ -53,7 +53,7 @@ namespace corsika::weights {
/**
* @class WeightData
*
- * definition of stack-data object to store geometry information
+ * definition of stack-data object to store geometry information.
*/
class WeightData {
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 70816c0d3..aa48469e3 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -53,9 +53,9 @@ add_executable (em_shower em_shower.cpp)
target_link_libraries (em_shower CORSIKA8::CORSIKA8)
CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS "100.")
-add_executable (hybrid_MC hybrid_MC.cpp)
-target_link_libraries (hybrid_MC CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
+#add_executable (hybrid_MC hybrid_MC.cpp)
+#target_link_libraries (hybrid_MC CORSIKA8::CORSIKA8)
+#CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
add_executable (corsika corsika.cpp)
target_link_libraries (corsika CORSIKA8::CORSIKA8)
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index 2f0d25447..05b989f15 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -17,6 +17,9 @@
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/output/OutputManager.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -110,6 +113,8 @@ int main() {
OutputManager output("cascade_outputs");
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ output.add("energyloss", dEdX);
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
@@ -140,14 +145,14 @@ int main() {
corsika::sibyll::Interaction sibyll;
corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
corsika::sibyll::Decay decay;
+
// cascade with only HE model ==> HE cut
- ParticleCut cut(80_GeV, true, true);
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(80_GeV, true, true, dEdX);
+ BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
TrackWriter trackWriter;
output.add("tracks", trackWriter); // register TrackWriter
- BetheBlochPDG eLoss{showerAxis};
-
// assemble all processes into an ordered process list
auto sequence = make_sequence(stackInspect, make_sequence(sibyllNuc, sibyll), decay,
eLoss, cut, trackWriter);
@@ -159,15 +164,13 @@ int main() {
EAS.run();
output.endOfShower();
- eLoss.printProfile(); // print longitudinal profile
-
cut.showResults();
const HEPEnergyType Efinal =
cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
- cout << "total dEdX energy (GeV): " << eLoss.getTotal() / 1_GeV << endl
- << "relative difference (%): " << eLoss.getTotal() / E0 * 100 << endl;
+ cout << "total dEdX energy (GeV): " << eLoss.getEnergyLost() / 1_GeV << endl
+ << "relative difference (%): " << eLoss.getEnergyLost() / E0 * 100 << endl;
cut.reset();
output.endOfLibrary();
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 9934211b7..856484c88 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -17,6 +17,9 @@
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/output/OutputManager.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -124,7 +127,13 @@ int main() {
RNGManager<>::getInstance().registerRandomStream("pythia");
corsika::pythia8::Interaction pythia;
corsika::pythia8::Decay decay;
- ParticleCut cut(60_GeV, true, true);
+
+ ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ output.add("energyloss", dEdX);
+
+ BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, true, true, dEdX);
cut.printThresholds();
// RNGManager::getInstance().registerRandomStream("HadronicElasticModel");
@@ -134,9 +143,6 @@ int main() {
TrackWriter trackWriter;
output.add("tracks", trackWriter); // register TrackWriter
- ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
- BetheBlochPDG eLoss{showerAxis};
-
// assemble all processes into an ordered process list
auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index e8e688553..e64ff7f7a 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -27,7 +27,9 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
-#include <corsika/output/NoOutput.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -278,17 +280,12 @@ int main(int argc, char** argv) {
// create the output manager that we then register outputs with
OutputManager output(app["--filename"]->as<std::string>());
- EnergyLossWriterParquet dEdX_output{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output
- output.add("dEdX", dEdX_output);
- // register profile output
- LongitudinalProfileWriterParquet profile{showerAxis};
- output.add("profile", profile);
- // register ground particle output
- ParticleWriterParquet particles;
- output.add("particles", particles);
- // register TrackWriter
- TrackWriterParquet tracks;
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis, 10_g / square(1_cm), 200};
+ output.add("energyloss", dEdX);
+
+ // create a track writer and register it with the output manager
+ TrackWriter<TrackWriterParquet> tracks;
output.add("tracks", tracks);
corsika::sibyll::Interaction sibyll;
@@ -322,20 +319,20 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
- HEPEnergyType const emcut = 1_GeV;
- HEPEnergyType const hadcut = 1_GeV;
- ParticleCut cut(emcut, emcut, hadcut, hadcut, true);
+ HEPEnergyType const emcut = 50_GeV;
+ HEPEnergyType const hadcut = 50_GeV;
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, true,
+ dEdX);
- ParticleCut cut{dEdX_output, 50_GeV, 50_GeV, 50_GeV, 50_GeV, false};
corsika::proposal::Interaction emCascade(env);
// NOT available for PROPOSAL due to interface trouble:
// InteractionCounter emCascadeCounted(emCascade);
- // corsika::proposal::ContinuousProcess emContinuous(env);
- BetheBlochPDG emContinuous(showerAxis);
+ // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env);
+ BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
- // cut.printThresholds();
-
- LongitudinalProfile longprof{profile};
+ LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> longprof{
+ showerAxis, 10_g / square(1_cm), 200};
+ output.add("profile", longprof);
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted(urqmd);
@@ -355,8 +352,10 @@ int main(int argc, char** argv) {
// observation plane
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane observationLevel{obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
- particles};
+ ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
+ // register ground particle output
+ output.add("particles", observationLevel);
// assemble the final process sequence
auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, cut,
@@ -417,8 +416,7 @@ int main(int argc, char** argv) {
cut.showResults();
// emContinuous.showResults();
observationLevel.showResults();
- HEPEnergyType const Efinal =
- dEdX_output.getTotal() + particleOutput.getEnergyGround();
+ HEPEnergyType const Efinal = dEdX.getTotal() + observationLevel.getTotalEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 60bb71c72..57bb1fb94 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -22,6 +22,9 @@
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -107,10 +110,7 @@ int main(int argc, char** argv) {
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
@@ -137,47 +137,44 @@ int main(int argc, char** argv) {
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
- std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
- << std::endl;
+ CORSIKA_LOG_INFO("shower axis length: {} ",
+ (showerCore - injectionPos).getNorm() * 1.02);
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
OutputManager output("em_shower_outputs");
- EnergyLossWriterParquet dEdX_output{showerAxis, 10_g / square(1_cm), 200};
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output
- output.add("dEdX", dEdX_output);
- // register profile output
- LongitudinalProfileWriterParquet profile{showerAxis};
- output.add("profile", profile);
- // register ground particle output
- ParticleWriterParquet particles;
- output.add("particles", particles);
- // register TrackWriter
- TrackWriterParquet tracks;
- output.add("tracks", tracks);
+ output.add("dEdX", dEdX);
// setup processes, decays and interactions
- ParticleCut cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true);
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true, true,
+ dEdX);
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
+ // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade);
- TrackWriter trackWriter{tracks};
+ TrackWriter<TrackWriterParquet> tracks;
+ output.add("tracks", tracks);
// long. profile; columns for photon, e+, e- still need to be added
- LongitudinalProfile longprof(showerAxis);
+ LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> longprof{
+ showerAxis, 10_g / square(1_cm), 200};
+ output.add("profile", longprof);
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking> observationLevel{obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
- "particles.dat", particles};
+ ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
+ output.add("particles", observationLevel);
- auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, observationLevel,
- trackWriter);
+ auto sequence =
+ make_sequence(emCascade, emContinuous, longprof, cut, observationLevel, tracks);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, output, stack);
@@ -202,7 +199,5 @@ int main(int argc, char** argv) {
cut.reset();
emContinuous.reset();
- longprof.save("longprof_emShower.txt");
-
output.endOfLibrary();
}
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 2de965757..904c0d082 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -27,6 +27,9 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -157,10 +160,7 @@ int main(int argc, char** argv) {
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
@@ -228,13 +228,23 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
- ParticleCut cut(3_GeV, false, true);
- BetheBlochPDG eLoss(showerAxis);
+ // register energy losses as output
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis, 10_g / square(1_cm), 200};
+ output.add("energyloss", dEdX);
+
+ // create a track writer and register it with the output manager
+ TrackWriter<TrackWriterParquet> tracks;
+ output.add("tracks", tracks);
+
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, true, dEdX);
+ BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX);
CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
get_PDG(Code::Proton));
- OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+ LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> longprof{
+ showerAxis, 10_g / square(1_cm), 200};
+ output.add("profile", longprof);
LongitudinalProfile longprof(showerAxis);
diff --git a/examples/mars.cpp b/examples/mars.cpp
index 12a4039bb..28b9efabd 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -29,7 +29,9 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
-#include <corsika/output/NoOutput.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -305,12 +307,21 @@ int main(int argc, char** argv) {
// we make the axis much longer than the inj-core distance since the
// profile will go beyond the core, depending on zenith angle
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
/* === END: CONSTRUCT GEOMETRY === */
// create the output manager that we then register outputs with
OutputManager output(app["--filename"]->as<std::string>());
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
+
+ EnergyLossWriter<EnergyLossWriterParquet> dEdX{showerAxis};
+ output.add("energyloss", dEdX);
+
+ HEPEnergyType const emcut = 1_GeV;
+ HEPEnergyType const hadcut = 1_GeV;
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, true,
+ dEdX);
+
/* === START: SETUP PROCESS LIST === */
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
@@ -343,16 +354,14 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
- HEPEnergyType const emcut = 1_GeV;
- HEPEnergyType const hadcut = 1_GeV;
- ParticleCut cut(emcut, emcut, hadcut, hadcut, true);
corsika::proposal::Interaction emCascade(env);
// NOT possible right now, due to interface difference for PROPOSAL:
// InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess emContinuous(env);
- BetheBlochPDG emContinuous(showerAxis);
+ BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
- LongitudinalProfile longprof{showerAxis, 1_g / square(1_cm)};
+ LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> profile{showerAxis};
+ output.add("profile", profile);
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
@@ -382,7 +391,7 @@ int main(int argc, char** argv) {
// assemble the final process sequence
auto sequence =
make_sequence(stackInspect, hadronSequence, decaySequence, emCascade, emContinuous,
- cut, trackWriter, observationLevel, longprof);
+ cut, trackWriter, observationLevel, profile);
/* === END: SETUP PROCESS LIST === */
// create the cascade object using the default stack and tracking implementation
@@ -415,7 +424,6 @@ int main(int argc, char** argv) {
string const outdir(app["--filename"]->as<std::string>());
string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
- string const longprof_file = outdir + "/longprof_" + to_string(i_shower) + ".txt";
// setup particle stack, and add primary particle
stack.clear();
@@ -445,7 +453,6 @@ int main(int argc, char** argv) {
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
- longprof.save(longprof_file);
// trigger the output manager to save this shower to disk
output.endOfShower();
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index 8ca094f6e..8d0641960 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -48,9 +48,7 @@ int main() {
rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
- ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env, false,
- 100};
- BetheBlochPDG eLoss{showerAxis};
+ BetheBlochPDG eLoss;
setup::Stack stack;
@@ -64,10 +62,7 @@ int main() {
double theta = 0.;
double phi = 0.;
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index a8ee6af6f..13b335ab3 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -29,7 +29,9 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
-#include <corsika/output/NoOutput.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
@@ -41,17 +43,12 @@
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriterParquet.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
-#include <corsika/modules/writers/BetheBlochPDGWriterParquet.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
-#include <corsika/modules/writers/ParticleCutWriterParquet.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/UrQMD.hpp>
@@ -209,11 +206,10 @@ int main(int argc, char** argv) {
// construct the continuous energy loss model
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
- // output.add("bethebloch", emContinuous);
// construct a particle cut
- ParticleCut<SubWriter<decltype(dEdX)>> cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true, dEdX};
- // output.add("cuts", cut);
+ ParticleCut<SubWriter<decltype(dEdX)>> cut{60_GeV, 60_GeV, 60_GeV, 60_GeV,
+ true, false, dEdX};
// setup longitudinal profile
LongitudinalProfile<corsika::LongitudinalProfileWriterParquet> profile{showerAxis};
@@ -284,8 +280,8 @@ int main(int argc, char** argv) {
plab.normalized(), injectionPos, 0_ns));
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<<setup::Tracking,ParticleWriterParquet> observationLevel{
- obsPlane, DirectionVector(rootCS, {1., 0., 0.}), particleOutput};
+ ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel);
// auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence,
@@ -301,14 +297,8 @@ int main(int argc, char** argv) {
EAS.run();
output.endOfShower();
- // cut.showResults();
- // emContinuous.showResults();
observationLevel.showResults();
- // HEPEnergyType const Efinal = dEdX_output.getTotal() +
- // particleOutput.getTotalEnergy();
observationLevel.reset();
- // cut.reset();
- // emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index 3ec511b62..eaa229222 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -41,6 +41,9 @@ TEST_CASE("ParticleProperties", "[Particles]") {
CHECK(Positron::charge / constants::e == Approx(+1));
CHECK(get_charge(Positron::anti_code) / constants::e == Approx(-1));
CHECK(Photon::charge / constants::e == 0.);
+ CHECK_FALSE(is_charged(Code::Photon));
+ CHECK(is_charged(Code::Iron));
+ CHECK(is_charged(Code::PiPlus));
}
SECTION("Names") {
@@ -87,11 +90,11 @@ TEST_CASE("ParticleProperties", "[Particles]") {
SECTION("Energy threshold") {
//! by default energy thresholds are set to zero
- CHECK(calculate_kinetic_energy_threshold(Electron::code) == 0_GeV);
+ CHECK(get_kinetic_energy_threshold(Electron::code) == 0_GeV);
set_kinetic_energy_threshold(Electron::code, 10_GeV);
- CHECK_FALSE(calculate_kinetic_energy_threshold(Code::Electron) == 1_GeV);
- CHECK(calculate_kinetic_energy_threshold(Code::Electron) == 10_GeV);
+ CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == 1_GeV);
+ CHECK(get_kinetic_energy_threshold(Code::Electron) == 10_GeV);
}
SECTION("Particle groups: electromagnetic") {
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 00f28c41c..a5bd79ed9 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -7,7 +7,7 @@ set (test_modules_sources
testQGSJetII.cpp
testPythia8.cpp
testUrQMD.cpp
- testCONEX.cpp
+# testCONEX.cpp
testParticleCut.cpp
testSibyll.cpp
testEpos.cpp
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index d04f3ac80..703bf6dc1 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -48,7 +48,7 @@ using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>;
struct DummyStack {};
-TEST_CASE("CONEXSourceCut") {
+TEST_CASE("CONEX") {
logging::set_level(logging::level::info);
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index a1803af75..33ada16d4 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -52,7 +52,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut on particle type: inv") {
// particle cut with 20GeV threshold for all, also cut invisible
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true, false);
CHECK(cut.getHadronKineticECut() == 20_GeV);
// add primary particle to stack
@@ -81,7 +81,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut on particle type: em") {
- ParticleCut cut(20_GeV, true, false);
+ ParticleCut cut(20_GeV, false, true);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
@@ -105,7 +105,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
}
SECTION("cut low energy") {
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true, false);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
@@ -134,7 +134,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
}
SECTION("cut low energy: electrons, photons, hadrons and muons") {
- ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, true);
+ ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, true, false);
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
@@ -168,15 +168,15 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut low energy: reset thresholds of arbitrary set of particles") {
ParticleCut cut({{Code::Electron, 5_MeV}, {Code::Positron, 50_MeV}}, false, true);
- CHECK(calculate_kinetic_energy_threshold(Code::Electron) !=
- calculate_kinetic_energy_threshold(Code::Positron));
- CHECK_FALSE(calculate_kinetic_energy_threshold(Code::Electron) == Electron::mass);
+ CHECK(get_kinetic_energy_threshold(Code::Electron) !=
+ get_kinetic_energy_threshold(Code::Positron));
+ CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == Electron::mass);
// test default values still correct
- CHECK(calculate_kinetic_energy_threshold(Code::Proton) == 5_GeV);
+ CHECK(get_kinetic_energy_threshold(Code::Proton) == 5_GeV);
}
SECTION("cut on time") {
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, false);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, false, false);
const TimeType too_late = 1_s;
// add primary particle to stack
@@ -208,7 +208,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
SECTION("cut on doContinous, just invisibles") {
- ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+ ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true, false);
// add particles, all with energies above the threshold
// only cut is by species
diff --git a/tests/output/testWriterTracks.cpp b/tests/output/testWriterTracks.cpp
index ea1cbcda8..f1d109fbe 100644
--- a/tests/output/testWriterTracks.cpp
+++ b/tests/output/testWriterTracks.cpp
@@ -21,7 +21,8 @@ struct TestWriterTrack : public TrackWriterParquet {
YAML::Node getConfig() const { return YAML::Node(); }
void checkWrite() {
- TrackWriterParquet::write(Code::Unknown, 1_eV, 1.0, {2_m, 3_m, 4_m}, {5_m, 6_m, 7_m});
+ TrackWriterParquet::write(Code::Unknown, 1_eV, 1.0, {2_m, 3_m, 4_m}, 1_ns,
+ {5_m, 6_m, 7_m}, 2_ns);
}
};
--
GitLab