diff --git a/applications/c8_air_shower.cpp b/applications/c8_air_shower.cpp
index fd7b025388f64009910e2f95983b3d1a8d2752fe..ba60fd009c180d620cc461056226ebe59cf9eede 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 6f694ebf605c2ea032de8fc88bd55f2e22087b1c..b999e862af7c6da85d5ec94da1670b73f63726b2 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 c247b424ec150b78d81f23d8d78e82cc94229df5..6f42de944526bcc32b6f33f73c83401c8e66df09 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 a4d6c347186a80444e7893ffa6dfc7fa3424d69d..55a15844e70bcb0d4d0a93573b05c09ac473ba29 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 55739516b55c1a55c7102960b705ff2f6c485aff..26029b80422cc9a28086ee7c265389245e83d004 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 cdee78eb4eba362168b1defecde73fd91c72cbd1..3a6ce5970d887aa3a1872bbe1ce8ab209ea3736f 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 65a0b1485fe72d0cf181f31c59a8f50927c5c3cb..8bf95f83a24861b851941f9759f53e2c5b5047d2 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 324e2e7328fde642704b5e37a00ca90f435d01c4..9943b634e48466faf04a227953671aec49df6c07 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 100bf3b82b9f0b3aa81dab7c64a649dcaec2cd0a..811a5284d22f6987dde25caf0f94ff000798fc5a 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); }
}