diff --git a/corsika/detail/modules/writers/LongitudinalWriter.inl b/corsika/detail/modules/writers/LongitudinalWriter.inl new file mode 100644 index 0000000000000000000000000000000000000000..f971e8b33e1adfd83b111e5e43b6074321b07a61 --- /dev/null +++ b/corsika/detail/modules/writers/LongitudinalWriter.inl @@ -0,0 +1,157 @@ +/* + * (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 <corsika/framework/utility/FindXmax.hpp> + +#include <corsika/media/ShowerAxis.hpp> + +#include <exception> + +namespace corsika { + + template <typename TOutput> + inline LongitudinalWriter<TOutput>::LongitudinalWriter(ShowerAxis const& axis, + GrammageType dX, + size_t const nBins) + : TOutput(number_profile::ProfileIndexNames) + , showerAxis_(axis) + , dX_(dX) + , nBins_(nBins) {} + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::startOfLibrary( + boost::filesystem::path const& directory) { + TOutput::startOfLibrary(directory); + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::startOfShower(unsigned int const showerId) { + TOutput::startOfShower(showerId); + // reset profile + profile_.clear(); + profile_.resize(nBins_); + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::endOfLibrary() { + TOutput::endOfLibrary(); + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::endOfShower(unsigned int const showerId) { + + int iRow{0}; + for (number_profile::ProfileData const& row : profile_) { + // here: write to underlying writer (e.g. parquet) + TOutput::write(showerId, iRow * dX_, row); + iRow++; + } + TOutput::endOfShower(showerId); + } + + template <typename TOutput> + template <typename TTrack> + inline void LongitudinalWriter<TOutput>::write(TTrack const& track, Code const pid, + double const weight) { + GrammageType const grammageStart = showerAxis_.getProjectedX(track.getPosition(0)); + GrammageType const grammageEnd = showerAxis_.getProjectedX(track.getPosition(1)); + + // Note: particle may go also "upward", thus, grammageEnd<grammageStart + int const binStart = std::ceil(grammageStart / dX_); + int const binEnd = std::floor(grammageEnd / dX_); + + for (int bin = binStart; bin <= binEnd; ++bin) { + if (pid == Code::Photon) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] += + weight; + } else if (pid == Code::Positron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] += + weight; + } else if (pid == Code::Electron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] += + weight; + } else if (pid == Code::MuPlus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] += + weight; + } else if (pid == Code::MuMinus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] += + weight; + } else if (is_hadron(pid)) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] += + weight; + } + if (is_charged(pid)) { + profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight; + } + } + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::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(TOutput::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>(number_profile::ProfileIndex::Photon)] += weight; + } else if (pid == Code::Positron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] += + weight; + } else if (pid == Code::Electron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] += + weight; + } else if (pid == Code::MuPlus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] += weight; + } else if (pid == Code::MuMinus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] += weight; + } else if (is_hadron(pid)) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] += weight; + } + if (is_charged(pid)) { + profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight; + } + } + + template <typename TOutput> + inline YAML::Node LongitudinalWriter<TOutput>::getConfig() const { + + YAML::Node node; + + node["type"] = "LongitudinalProfile"; + node["units"]["grammage"] = "g/cm^2"; + node["bin-size"] = dX_ / (1_g / square(1_cm)); + node["nbins"] = nBins_; + + return node; + } + + template <typename TOutput> + inline YAML::Node LongitudinalWriter<TOutput>::getSummary() const { + YAML::Node summary; + return summary; + } + +} // namespace corsika diff --git a/corsika/modules/writers/LongitudinalWriter.hpp b/corsika/modules/writers/LongitudinalWriter.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e1073e64eced5fc5147a16c185d998f1f715db56 --- /dev/null +++ b/corsika/modules/writers/LongitudinalWriter.hpp @@ -0,0 +1,140 @@ +/* + * (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/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/media/ShowerAxis.hpp> +#include <corsika/modules/writers/WriterOff.hpp> +#include <corsika/modules/writers/LongitudinalProfileWriterParquet.hpp> + +#include <vector> +#include <array> + +/** + * @file LongitudinalWriter.hpp + */ + +namespace corsika { + + /** + * Local helper namespace to store number and names of particle number profile columns. + */ + namespace number_profile { + + /** + * Definition of longitudinal profile columns. + */ + enum class ProfileIndex { + Charged, + Hadron, + Photon, + Electron, + Positron, + MuPlus, + MuMinus, + Entries + }; + + /** + * Number of columns (static). + */ + size_t constexpr NColumns = static_cast<size_t>(ProfileIndex::Entries); + + /** + * Names of columns in output. + */ + static std::array<char const*, NColumns> constexpr ProfileIndexNames{ + {"charged", "hadron", "photon", "electron", "positron", "muplus", "muminus"}}; + + /** + * Data type to store column data. + */ + typedef std::array<double, NColumns> ProfileData; + } // namespace number_profile + + // clang-format-off + /** + * The LongitudinalWriter can be used to pool the particle counts of several + * longitudinal profile processes into one output file/stream. + * + * Typically several processes/modules can lead to particle counts along the shower axis + * in the shower. The LongitudianalWriter can be used in combination with the SubWriter + * class to collect all of them into a single output stream: + * + * \code {.cpp} + * # showerAxis must be a ShowerAxis object + * # the X binning can be specified. + * LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200}; + * # add to OutputManager: + * output.add("profile", profile); + * # add SubWriters, e.g. LongitudinalProfile, CONEX: + * LongitudinalProfile<SubWriter<decltype(profile)>> long{profile}; + * CONEXhybrid<SubWriter<decltype(profile)>> conex{..., profile}; + * ... + * \endcode + * + * The default output option is parquet format. + * + * @tparam TOutput + */ + // clang-format-on + + template <typename TOutput = LongitudinalProfileWriterParquet<number_profile::NColumns>> + class LongitudinalWriter : public TOutput { + + public: + /** + * Construct a new writer. + */ + LongitudinalWriter(ShowerAxis const& axis, + GrammageType dX = 10_g / square(1_cm), // profile binning + size_t const nBins = 200); // number of bins + + void startOfLibrary(boost::filesystem::path const& directory) final override; + + void startOfShower(unsigned int const showerId) final override; + + void endOfShower(unsigned int const showerId) final override; + + void endOfLibrary() final override; + + /** + * Add continuous profile. + */ + template <typename TTrack> + void write(TTrack const& track, Code const pid, double const weight); + + /** + * Add binned profile. + */ + void write(GrammageType const Xstart, GrammageType const Xend, Code const pid, + double const weight); + + /** + * Return a summary. + */ + YAML::Node getSummary() const; + + /** + * Return the configuration of this output. + */ + YAML::Node getConfig() const; + + private: + ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage + GrammageType dX_; ///< binning of profile. + size_t nBins_; ///< number of profile bins. + std::vector<number_profile::ProfileData> profile_; // longitudinal profile + }; + +} // namespace corsika + +#include <corsika/detail/modules/writers/LongitudinalWriter.inl>