missing new class

corsika/detail/modules/writers/LongitudinalWriter.inl

+/*
+ * (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

corsika/modules/writers/LongitudinalWriter.hpp

+/*
+ * (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>