diff --git a/corsika/detail/modules/writers/EnergyLossWriter.inl b/corsika/detail/modules/writers/EnergyLossWriter.inl index e4101afeb2145f70e9f3069bd43b1a7c380b4a6e..f6a0cbd6eebd80ebd13ab96bd37f62d3ce0fc6a4 100644 --- a/corsika/detail/modules/writers/EnergyLossWriter.inl +++ b/corsika/detail/modules/writers/EnergyLossWriter.inl @@ -152,10 +152,16 @@ namespace corsika { "CONEX and Corsika8 dX grammage binning are not the same!"); } - int const bin = int((bend + bstart) / 2); + size_t const bin = size_t((bend + bstart) / 2); CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "add binned energy loss {} {} bin={} dE={} GeV ", bstart, bend, bin, dE / 1_GeV); + if (bin >= profile_.size()) { + CORSIKA_LOGGER_WARN(TOutput::getLogger(), + "Grammage bin {} outside of profile {}. skipping.", bin, + profile_.size()); + return; + } profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE; } diff --git a/corsika/detail/modules/writers/LongitudinalWriter.inl b/corsika/detail/modules/writers/LongitudinalWriter.inl index f971e8b33e1adfd83b111e5e43b6074321b07a61..5c8f8cd4420fe7402fb3b5241654900c7d304aae 100644 --- a/corsika/detail/modules/writers/LongitudinalWriter.inl +++ b/corsika/detail/modules/writers/LongitudinalWriter.inl @@ -113,7 +113,13 @@ namespace corsika { "CONEX and Corsika8 dX grammage binning are not the same!"); } - int const bin = int((bend + bstart) / 2); + size_t const bin = size_t((bend + bstart) / 2); + if (bin >= profile_.size()) { + CORSIKA_LOGGER_WARN(TOutput::getLogger(), + "Grammage bin {} outside of profile {}. skipping.", bin, + profile_.size()); + return; + } if (pid == Code::Photon) { profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] += weight; diff --git a/corsika/modules/conex/CONEXhybrid.hpp b/corsika/modules/conex/CONEXhybrid.hpp index 11ae68dc1a0bf3acd67bfa11b51863ab8feec862..4e5119e5cddf31f060a0547c2454f625f82a4d47 100644 --- a/corsika/modules/conex/CONEXhybrid.hpp +++ b/corsika/modules/conex/CONEXhybrid.hpp @@ -17,10 +17,7 @@ #include <corsika/framework/geometry/Vector.hpp> #include <corsika/media/ShowerAxis.hpp> -#include <corsika/modules/writers/WriterOff.hpp> #include <corsika/modules/writers/SubWriter.hpp> -#include <corsika/modules/writers/EnergyLossWriterParquet.hpp> -#include <corsika/modules/writers/LongitudinalProfileWriterParquet.hpp> #include <corsika/modules/writers/EnergyLossWriter.hpp> #include <corsika/modules/writers/LongitudinalWriter.hpp> @@ -52,7 +49,7 @@ namespace corsika { * @tparam TOutputE -- Output writer for dEdX data. * @tparam TOutputN -- Output writer for particle number profile data. */ - template <typename TOutputE = WriterOff, typename TOutputN = WriterOff2> + template <typename TOutputE, typename TOutputN> class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid<TOutputE, TOutputN>>, public SecondariesProcess<CONEXhybrid<TOutputE, TOutputN>>, public SubWriter<TOutputE>, diff --git a/corsika/modules/writers/LongitudinalWriter.hpp b/corsika/modules/writers/LongitudinalWriter.hpp index e1073e64eced5fc5147a16c185d998f1f715db56..ac00dc3f0b25a98d1e138f7ee75bd8adf3033869 100644 --- a/corsika/modules/writers/LongitudinalWriter.hpp +++ b/corsika/modules/writers/LongitudinalWriter.hpp @@ -128,6 +128,11 @@ namespace corsika { */ YAML::Node getConfig() const; + number_profile::ProfileData const& getProfile( + number_profile::ProfileIndex index) const { + return profile_.at(static_cast<int>(index)); + } + private: ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage GrammageType dX_; ///< binning of profile. diff --git a/corsika/modules/writers/WriterOff.hpp b/corsika/modules/writers/WriterOff.hpp index 0aa01006ad9fb72a4b386a6c4290cee4df7e457a..3d752cc9a80de7a85b0705c3d6dac9dac9d8722d 100644 --- a/corsika/modules/writers/WriterOff.hpp +++ b/corsika/modules/writers/WriterOff.hpp @@ -21,38 +21,29 @@ namespace corsika { class WriterOff : public BaseOutput { public: - WriterOff() {} - virtual ~WriterOff() {} - - void startOfLibrary(boost::filesystem::path const&) final override {} - - void endOfShower(unsigned int const) final override {} - - void endOfLibrary() final override {} - + /** + * The purpose of this catch-all constructor is to discard all parameters. + * @tparam TArgs + */ template <typename... TArgs> - void write(TArgs&&...) {} + WriterOff(TArgs&&...) {} - virtual YAML::Node getConfig() const override { return YAML::Node(); } - - }; // class WriterOff - - class WriterOff2 : public BaseOutput { - - public: - WriterOff2() {} - virtual ~WriterOff2() {} + virtual ~WriterOff() {} - void startOfLibrary(boost::filesystem::path const&) final override {} + void startOfLibrary(boost::filesystem::path const&) override {} - void endOfShower(unsigned int const) final override {} + void endOfShower(unsigned int const) override {} - void endOfLibrary() final override {} + void endOfLibrary() override {} + /** + * The purpose of this catch-all method is to discard all data. + * @tparam TArgs + */ template <typename... TArgs> void write(TArgs&&...) {} - virtual YAML::Node getConfig() const override { return YAML::Node(); } + virtual YAML::Node getConfig() const { return YAML::Node(); } }; // class WriterOff diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp index 03e4640881387367e24126d8ddc6e7426587ad80..8d8ce52efe95717594f58c3e965d67f19bc8eae9 100644 --- a/tests/modules/testCONEX.cpp +++ b/tests/modules/testCONEX.cpp @@ -102,10 +102,16 @@ TEST_CASE("CONEX") { corsika::sibyll::Interaction sibyll; [[maybe_unused]] corsika::sibyll::NuclearInteractionModel sibyllNuc(sibyll, env); - WriterOff w1; - WriterOff2 w2; - CONEXhybrid<WriterOff, WriterOff2> conex(center, showerAxis, t, injectionHeight, E0, - get_PDG(Code::Proton), w1, w2); + EnergyLossWriter<WriterOff> w1(showerAxis); + LongitudinalWriter<WriterOff> w2(showerAxis); + CONEXhybrid<decltype(w1), decltype(w2)> conex(center, showerAxis, t, injectionHeight, + E0, get_PDG(Code::Proton), w1, w2); + // initialize writers + w1.startOfLibrary("test"); + w1.startOfShower(0); + w2.startOfLibrary("test"); + w2.startOfShower(0); + // init conex conex.initCascadeEquations(); HEPEnergyType const Eem{1_PeV}; @@ -123,45 +129,14 @@ TEST_CASE("CONEX") { emPosition.getCoordinates(conex.getObserverCS()), emPosition.getCoordinates(rootCS)); - conex.addParticle(Code::Proton, Eem, 0_eV, emPosition, momentum.normalized(), 0_s); + conex.addParticle(Code::Proton, Eem, Proton::mass, emPosition, momentum.normalized(), + 0_s); // supperimpose a photon auto const momentumPhoton = showerAxis.getDirection() * 1_TeV; conex.addParticle(Code::Photon, 1_TeV, 0_eV, emPosition, momentumPhoton.normalized(), 0_s); DummyStack stack; conex.doCascadeEquations(stack); -} - -#include <algorithm> -#include <iterator> -#include <string> -#include <fstream> - -TEST_CASE("ConexOutput", "[output validation]") { - - logging::set_level(logging::level::info); - - auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output"); - - SECTION(std::string("check saved data, ") + file + ".txt") { - - // compare to reference data - std::ifstream file1(file + ".txt"); - std::ifstream file1ref(refDataDir + "/" + file + "_REF.txt"); - - std::istreambuf_iterator<char> begin1(file1); - std::istreambuf_iterator<char> begin1ref(file1ref); - - std::istreambuf_iterator<char> end; - while (begin1 != end && begin1ref != end) { - CHECK(*begin1 == *begin1ref); - ++begin1; - ++begin1ref; - } - CHECK(begin1 == end); - CHECK(begin1ref == end); - file1.close(); - file1ref.close(); - } + CHECK(w1.getEnergyLost() / 1_TeV == Approx(1.0).epsilon(0.1)); }