From 709b74e04d2b1abbdc85bbb0f781a706231ae166 Mon Sep 17 00:00:00 2001
From: Alan Coleman <alanco@umich.edu>
Date: Wed, 13 Dec 2023 11:59:09 +0000
Subject: [PATCH] Resolve "Make EnergyLossWriter not require nbins"
---
applications/c8_air_shower.cpp | 5 ++---
.../modules/writers/EnergyLossWriter.inl | 8 ++++++++
corsika/modules/writers/EnergyLossWriter.hpp | 12 ++++++++++-
examples/cascade_examples/em_shower.cpp | 5 ++---
examples/cascade_examples/mars.cpp | 5 ++---
examples/cascade_examples/mc_conex.cpp | 5 ++---
examples/cascade_examples/radio_em_shower.cpp | 5 ++---
examples/cascade_examples/water.cpp | 5 ++---
tests/output/testWriterEnergyLoss.cpp | 20 +++++++++----------
9 files changed, 41 insertions(+), 29 deletions(-)
diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index fd7b02538..ba60fd009 100644
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -353,7 +353,6 @@ int main(int argc, char** argv) {
// profile will go beyond the core, depending on zenith angle
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
- uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
/* === END: CONSTRUCT GEOMETRY === */
double const emthinfrac = app["--emthin"]->as<double>();
@@ -372,7 +371,7 @@ int main(int argc, char** argv) {
OutputManager output(app["--filename"]->as<std::string>(), seed, args.str(), outputDir);
// register energy losses as output
- EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
+ EnergyLossWriter dEdX{showerAxis, dX};
output.add("energyloss", dEdX);
DynamicInteractionProcess<StackType> heModel;
@@ -434,7 +433,7 @@ int main(int argc, char** argv) {
auto emContinuous =
make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
- LongitudinalWriter profile{showerAxis, nAxisBins, dX};
+ LongitudinalWriter profile{showerAxis, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
diff --git a/corsika/detail/modules/writers/EnergyLossWriter.inl b/corsika/detail/modules/writers/EnergyLossWriter.inl
index 6f694ebf6..b999e862a 100644
--- a/corsika/detail/modules/writers/EnergyLossWriter.inl
+++ b/corsika/detail/modules/writers/EnergyLossWriter.inl
@@ -21,7 +21,15 @@ namespace corsika {
template <typename TOutput>
inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
GrammageType dX,
+ GrammageType dX_threshold)
+ : EnergyLossWriter<TOutput>{axis,
+ static_cast<unsigned int>(axis.getMaximumX() / dX) + 1,
+ dX, dX_threshold} {}
+
+ template <typename TOutput>
+ inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
unsigned int const nBins,
+ GrammageType dX,
GrammageType dX_threshold)
: TOutput(dEdX_output::ProfileIndexNames)
, showerAxis_(axis)
diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp
index c247b424e..6f42de944 100644
--- a/corsika/modules/writers/EnergyLossWriter.hpp
+++ b/corsika/modules/writers/EnergyLossWriter.hpp
@@ -111,9 +111,17 @@ namespace corsika {
/**
* Construct a new writer.
*/
+
+ // Number of bins defined explicitly
EnergyLossWriter(ShowerAxis const& axis,
GrammageType dX = 10_g / square(1_cm), // profile binning
- unsigned int const nBins = 200, // number of bins
+ GrammageType dX_threshold = 0.0001_g /
+ square(1_cm)); // ignore too short tracks
+
+ // Number of bins defined explicitly
+ EnergyLossWriter(ShowerAxis const& axis,
+ unsigned int const nBins, // number of bins
+ GrammageType dX = 10_g / square(1_cm), // profile binning
GrammageType dX_threshold = 0.0001_g /
square(1_cm)); // ignore too short tracks
@@ -141,6 +149,8 @@ namespace corsika {
void write(GrammageType const Xstart, GrammageType const Xend, Code const PID,
HEPEnergyType const dE);
+ auto GetNBins() const { return nBins_; }
+
/**
* Get total observed energy loss.
*
diff --git a/examples/cascade_examples/em_shower.cpp b/examples/cascade_examples/em_shower.cpp
index a4d6c3471..55a15844e 100644
--- a/examples/cascade_examples/em_shower.cpp
+++ b/examples/cascade_examples/em_shower.cpp
@@ -134,7 +134,6 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
- uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
@@ -145,7 +144,7 @@ int main(int argc, char** argv) {
(showerCore - injectionPos).getNorm() * 1.02);
// setup processes, decays and interactions
- EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
+ EnergyLossWriter energyloss{showerAxis, dX};
ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
energyloss);
@@ -167,7 +166,7 @@ int main(int argc, char** argv) {
TrackWriter tracks;
output.add("tracks", tracks);
- LongitudinalWriter profile{showerAxis, nAxisBins, dX};
+ LongitudinalWriter profile{showerAxis, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
diff --git a/examples/cascade_examples/mars.cpp b/examples/cascade_examples/mars.cpp
index 55739516b..26029b804 100644
--- a/examples/cascade_examples/mars.cpp
+++ b/examples/cascade_examples/mars.cpp
@@ -304,9 +304,8 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
- uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
- EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
+ EnergyLossWriter dEdX{showerAxis, dX};
output.add("energyloss", dEdX);
HEPEnergyType const emcut = 1_GeV;
@@ -338,7 +337,7 @@ int main(int argc, char** argv) {
auto emContinuous =
make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
- LongitudinalWriter longprof{showerAxis, nAxisBins, dX};
+ LongitudinalWriter longprof{showerAxis, dX};
output.add("profile", longprof);
LongitudinalProfile<SubWriter<decltype(longprof)>> profile{longprof};
diff --git a/examples/cascade_examples/mc_conex.cpp b/examples/cascade_examples/mc_conex.cpp
index cdee78eb4..3a6ce5970 100644
--- a/examples/cascade_examples/mc_conex.cpp
+++ b/examples/cascade_examples/mc_conex.cpp
@@ -207,7 +207,6 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
- uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
@@ -222,7 +221,7 @@ int main(int argc, char** argv) {
OutputManager output("hybrid_MC_outputs");
// register energy losses as output
- EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
+ EnergyLossWriter dEdX{showerAxis, dX};
output.add("energyloss", dEdX);
// create a track writer and register it with the output manager
@@ -232,7 +231,7 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX);
- LongitudinalWriter profile{showerAxis, nAxisBins, dX};
+ LongitudinalWriter profile{showerAxis, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
diff --git a/examples/cascade_examples/radio_em_shower.cpp b/examples/cascade_examples/radio_em_shower.cpp
index 65a0b1485..8bf95f83a 100644
--- a/examples/cascade_examples/radio_em_shower.cpp
+++ b/examples/cascade_examples/radio_em_shower.cpp
@@ -158,7 +158,6 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
- uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
// setup the radio antennas
TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
@@ -229,7 +228,7 @@ int main(int argc, char** argv) {
// }
// setup processes, decays and interactions
- EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
+ EnergyLossWriter energyloss{showerAxis, dX};
ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
energyloss);
@@ -251,7 +250,7 @@ int main(int argc, char** argv) {
TrackWriter tracks;
output.add("tracks", tracks);
- LongitudinalWriter profile{showerAxis, nAxisBins, dX};
+ LongitudinalWriter profile{showerAxis, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
diff --git a/examples/cascade_examples/water.cpp b/examples/cascade_examples/water.cpp
index 324e2e732..9943b634e 100644
--- a/examples/cascade_examples/water.cpp
+++ b/examples/cascade_examples/water.cpp
@@ -200,12 +200,11 @@ int main(int argc, char** argv) {
// * longitutional profile
ShowerAxis const showerAxis{injectionPos, 1.2 * injectorLength * downVec, env};
auto const dX = 1_g / square(1_cm); // Binning of the writers along the shower axis
- uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
- LongitudinalWriter longiWriter{showerAxis, nAxisBins, dX};
+ LongitudinalWriter longiWriter{showerAxis, dX};
LongitudinalProfile<SubWriter<decltype(longiWriter)>> longprof{longiWriter};
// * energy loss profile
- EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
+ EnergyLossWriter dEdX{showerAxis, dX};
// * physical process list
// particle production threshold
diff --git a/tests/output/testWriterEnergyLoss.cpp b/tests/output/testWriterEnergyLoss.cpp
index 100bf3b82..811a5284d 100644
--- a/tests/output/testWriterEnergyLoss.cpp
+++ b/tests/output/testWriterEnergyLoss.cpp
@@ -62,7 +62,7 @@ TEST_CASE("EnergyLossWriter") {
[[maybe_unused]] auto const& node_dummy = nodePtr;
auto const observationHeight = 0_km;
- auto const injectionHeight = 10_km;
+ auto const injectionHeight = 20_km;
auto const t = -observationHeight + injectionHeight;
Point const showerCore{cs, 0_m, 0_m, observationHeight};
Point const injectionPos = showerCore + DirectionVector{cs, {0, 0, 1}} * t;
@@ -71,14 +71,14 @@ TEST_CASE("EnergyLossWriter") {
false, // -> throw exceptions
1000}; // -> number of bins
+ std::string const outputDir = "./output_dir_eloss";
+
// preparation
- if (boost::filesystem::exists("./output_dir_eloss")) {
- boost::filesystem::remove_all("./output_dir_eloss");
- }
- boost::filesystem::create_directory("./output_dir_eloss");
+ if (boost::filesystem::exists(outputDir)) { boost::filesystem::remove_all(outputDir); }
+ boost::filesystem::create_directory(outputDir);
TestEnergyLoss test(showerAxis);
- test.startOfLibrary("./output_dir_eloss");
+ test.startOfLibrary(outputDir);
test.startOfShower(0);
CHECK(test.getEnergyLost() / 1_GeV == Approx(0));
@@ -116,7 +116,7 @@ TEST_CASE("EnergyLossWriter") {
test.endOfShower(0);
test.endOfLibrary();
- CHECK(boost::filesystem::exists("./output_dir_eloss/dEdX.parquet"));
+ CHECK(boost::filesystem::exists(outputDir + "/dEdX.parquet"));
auto const config = test.getConfig();
CHECK(config["type"].as<std::string>() == "EnergyLoss");
@@ -128,7 +128,7 @@ TEST_CASE("EnergyLossWriter") {
auto const summary = test.getSummary();
CHECK(summary["sum_dEdX"].as<double>() == 600);
- // makes not yet sense:
- // CHECK(summary["Xmax"].as<double>() == 200);
- // CHECK(summary["dEdXmax"].as<double>() == 200);
+
+ // clean up
+ if (boost::filesystem::exists(outputDir)) { boost::filesystem::remove_all(outputDir); }
}
--
GitLab