corsika/modules/Random.hpp

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <iterator>
+#include <functional>
+#include <random>
+#include <string_view>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+namespace corsika {
+  using rng_function_type = std::function<void(double*, std::size_t)>;
+
+  namespace detail {
+    inline void rng_func(corsika::default_prng_type& rng, double* dest, std::size_t N) {
+      std::uniform_real_distribution<double> udist(0.0, 1.0);
+      std::generate(dest, std::next(dest, N), std::bind(udist, std::ref(rng)));
+    }
+  } // namespace detail
+
+  inline void connect_random_stream(corsika::default_prng_type& rng,
+                                    void (*injection_func)(rng_function_type)) {
+    using namespace std::placeholders;
+    injection_func(std::bind(detail::rng_func, rng, _1, _2));
+  }
+
+  inline void connect_random_stream(std::string_view stream_name,
+                                    void (*injection_func)(rng_function_type)) {
+    auto& rng = RNGManager<>::getInstance().getRandomStream(std::string{stream_name});
+    connect_random_stream(rng, injection_func);
+  }
+} // namespace corsika

corsika/modules/Sibyll.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/sibyll/ParticleConversion.hpp>
+#include <corsika/modules/sibyll/HadronInteractionModel.hpp>
+#include <corsika/modules/sibyll/Decay.hpp>
+#include <corsika/modules/sibyll/NuclearInteractionModel.hpp>
+
+#include <corsika/modules/sibyll/InteractionModel.hpp>
+
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file Sibyll.hpp
+ *
+ * Includes all the parts of the Sibyll model. Defines the InteractionProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::sibyll {
+  /**
+   * @brief sibyll::Interaction is the process for ProcessSequence.
+   *
+   * The sibyll::InteractionModel is wrapped as an InteractionProcess here in order
+   * to provide all the functions for ProcessSequence.
+   */
+  struct Interaction : public InteractionModel, public InteractionProcess<Interaction> {
+    Interaction(std::set<Code> const& nuccomp, std::set<Code> const& stablehad)
+        : InteractionModel{nuccomp, stablehad} {}
+  };
+
+  /**
+   * @brief sibyll::NuclearInteraction is the process for ProcessSequence.
+   *
+   * The sibyll::NuclearInteractionModel is wrapped as an InteractionProcess here in order
+   * to provide all the functions for ProcessSequence.
+   */
+  template <class TNucleonModel>
+  class NuclearInteraction
+      : public NuclearInteractionModel<TNucleonModel>,
+        public InteractionProcess<NuclearInteraction<TNucleonModel>> {
+  public:
+    NuclearInteraction(TNucleonModel& model, std::set<Code> const& nuccomp)
+        : NuclearInteractionModel<TNucleonModel>{model, nuccomp} {}
+  };
+
+} // namespace corsika::sibyll

corsika/modules/Sophia.hpp

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/sophia/ParticleConversion.hpp>
+
+#include <corsika/modules/sophia/InteractionModel.hpp>
+
+#include <corsika/framework/process/InteractionProcess.hpp>

corsika/modules/StackInspector.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/StackProcess.hpp>
+
+#include <chrono>
+
+namespace corsika {
+
+  /**
+   * StackProcess that will act each @f$n_{step}@f$ steps to perform diagnostics on the
+   * full stack.
+   *
+   * The StackInspector can dump the entrie stack content for debugging, or also just
+   * determine the total energy remaining on the stack. From the decrease of energy on the
+   * stack an ETA for the completion of the simulation is determined.
+   *
+   * @tparam TStack Is the type of the particle stack.
+   */
+
+  template <typename TStack>
+  class StackInspector : public StackProcess<StackInspector<TStack>> {
+
+    typedef typename TStack::particle_type Particle;
+
+    using StackProcess<StackInspector<TStack>>::getStep;
+
+  public:
+    StackInspector(int const nStep, bool const reportStack, HEPEnergyType const vE0);
+    ~StackInspector();
+
+    void doStack(TStack const&);
+
+    /**
+     * To set a new E0, for example when a new shower event is started.
+     */
+    void setE0(HEPEnergyType const E0) { E0_ = E0; }
+
+  private:
+    bool ReportStack_;
+    int PrintoutCounter_ = 0;
+    const int MaxNumberOfPrintouts_ = 10;
+    HEPEnergyType E0_;
+    const HEPEnergyType dE_threshold_ = 1_eV;
+    std::chrono::system_clock::time_point StartTime_;
+    HEPEnergyType energyPostInit_;
+    std::chrono::system_clock::time_point timePostInit_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/modules/StackInspector.inl>

corsika/modules/TAUOLA.hpp

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/modules/tauola/Decay.hpp>

corsika/modules/TimeCut.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * 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/process/ContinuousProcess.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+namespace corsika {
+
+  /**
+   * Simple TimeCut process, removes particles older than specified cut time.
+   */
+  class TimeCut : public ContinuousProcess<TimeCut> {
+
+  public:
+    TimeCut(TimeType const time);
+
+    template <typename Particle>
+    ProcessReturn doContinuous(
+        Step<Particle> const& step,
+        const bool limitFlag = false); // this is not used for TimeCut
+    template <typename Particle, typename Track>
+    LengthType getMaxStepLength(Particle const&, Track const&) {
+      return meter * std::numeric_limits<double>::infinity();
+    }
+
+  private:
+    TimeType time_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/modules/TimeCut.inl>

corsika/modules/TrackWriter.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/modules/writers/TrackWriterParquet.hpp>
+#include <corsika/modules/writers/WriterOff.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+namespace corsika {
+
+  /**
+   * @ingroup Modules
+   * @{
+   *
+   * To write 3D track data to disk.
+   *
+   * Since the only sole purpose of this module is to generate track
+   * output on disk, the default output mode is not "WriterOff" but
+   * directly TrackWriterParquet. It can of course be changed.
+   *
+   * @tparam TOutput with default TrackWriterParquet
+   */
+
+  template <typename TOutput = TrackWriterParquet>
+  class TrackWriter : public ContinuousProcess<TrackWriter<TOutput>>, public TOutput {
+
+  public:
+    TrackWriter();
+
+    template <typename TParticle>
+    ProcessReturn doContinuous(Step<TParticle> const&, bool const limitFlag);
+
+    template <typename TParticle, typename TTrack>
+    LengthType getMaxStepLength(TParticle const&, TTrack const&);
+
+    YAML::Node getConfig() const;
+
+  private:
+  };
+
+  //! @}
+
+} // namespace corsika
+
+#include <corsika/detail/modules/TrackWriter.inl>

corsika/modules/Tracking.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/tracking/TrackingStraight.hpp>
+#include <corsika/modules/tracking/TrackingLeapFrogStraight.hpp> // simple leap-frog implementation with two straight lines
+#include <corsika/modules/tracking/TrackingLeapFrogCurved.hpp> // // more complete, curved, leap-frog
+
+/**
+ * \todo add TrackingCurved, with simple circular trajectories
+   \todo add Boris algorithm
+ */

corsika/modules/UrQMD.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/urqmd/UrQMD.hpp>

corsika/modules/conex/CONEX_f.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+//----------------------------------------------
+//  C++ interface for the CONEX fortran code
+//----------------------------------------------
+// wrapper
+
+#include <conexConfig.h>
+#include <conexHEModels.h>
+#include <array>
+
+namespace conex {
+  // the CONEX fortran interface
+
+  extern "C" {
+  extern struct { std::array<double, 16> dptl; } cxoptl_;
+
+  //! common block for atmosphere composition
+  extern struct {
+    std::array<double, 3> airz, aira, airw; //!< nuclear Z, A, composition fraction
+    double airavz, airava;                  //!< average Z, A
+    std::array<double, 3> airi;             //!< ionization potential, not used in cxroot
+  } cxair_;
+
+  void cegs4_(int&, int&);
+
+  void initconex_(int&, int*, int&, int&,
+#ifdef CONEX_EXTENSIONS
+                  int&,
+#endif
+                  const char*, int);
+  void conexrun_(int& ipart, double& energy, double& theta, double& phi,
+                 double& injectionHeight, double& dimpact, int ioseed[3]);
+  void conexcascade_();
+  void hadroncascade_(int&, int&, int&, int&);
+  void solvemomentequations_(int&);
+  void show_(int& iqi, double& ei, double& xmi, double& ymi, double& zmi, double& dmi,
+             double& xi, double& yi, double& zi, double& tmi, double& ui, double& vi,
+             double& wi, int& iri, double& wti, int& latchi);
+
+  int get_number_of_depth_bins_();
+
+  void get_shower_data_(const int&, const int&, const int&, float&, float&, float&,
+                        float&, float&);
+  void get_shower_edep_(const int&, const int&, float&, float&);
+  void get_shower_muon_(const int&, const int&, float&, float&);
+  void get_shower_gamma_(const int&, const int&, float&);
+  void get_shower_electron_(const int&, const int&, float&);
+  void get_shower_positron_(const int&, const int&, float&);
+  void get_shower_hadron_(const int&, const int&, float&);
+  }
+
+} // namespace conex

corsika/modules/conex/CONEXhybrid.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/SecondariesProcess.hpp>
+#include <corsika/framework/process/CascadeEquationsProcess.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
+
+#include <corsika/modules/conex/CONEX_f.hpp>
+
+/**
+ * @file CONEXhybrid.hpp
+ */
+
+namespace corsika {
+
+  namespace conex {
+    LengthType constexpr earthRadius{6371315 * meter};
+  } // namespace conex
+
+  /**
+   * Access to the CONEX model.
+   *
+   * The fortran version of CONEX is interfaced. This is a SecondariesProcess since it
+   * skims the Stack for input particles (specific energies, or types) in order to fill
+   * internal CONEX source tables. CONEX is the called after the Cascade is finished to do
+   * a CascadeEquations step. This typically produces output data, and in addition can
+   * even generate new secondaries on the main Stack.
+   *
+   * Note that the output is processed by `SubWriter<TOutputE/N>`. The SubWriters have to
+   * be initialized with valid objects of type TOutputE/N at construction time.If no
+   * output is wished, `WriterOff` can be used for this purpose.
+   *
+   * @tparam TOutputE -- Output writer for dEdX data.
+   * @tparam TOutputN -- Output writer for particle number profile data.
+   */
+  template <typename TOutputE, typename TOutputN>
+  class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid<TOutputE, TOutputN>>,
+                      public SecondariesProcess<CONEXhybrid<TOutputE, TOutputN>>,
+                      public SubWriter<TOutputE>,
+                      public SubWriter<TOutputN> {
+
+  public:
+    /**
+     * Constructor with output sink specified.
+     *
+     * The output will be generated by the output sink.
+     *
+     * @param center center of earth.
+     * @param showerAxis shower axis to convert geometry to grammage.
+     * @param groundDist distance to ground.
+     * @param injectionHeight height of injection.
+     * @param primaryEnergy energy of primary particle
+     * @param pdg type of primary particle.
+     * @param outputE object to initialized SubWriter<TOutputE>
+     * @param outputN object to initialized SubWriter<TOutputN>
+     */
+    CONEXhybrid(Point const& center, ShowerAxis const& showerAxis, LengthType groundDist,
+                LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode pdg,
+                TOutputE& outputE, TOutputN& outputN);
+
+    template <typename TStackView>
+    void doSecondaries(TStackView&);
+
+    /**
+     * init currently needs to be called to initializa a new
+     * event. All tables are cleared, etc.
+     */
+    void initCascadeEquations();
+
+    /**
+     * Cascade equations are solved basoned on the data in the tables.
+     */
+    template <typename TStack>
+    void doCascadeEquations(TStack& stack);
+
+    /**
+     * Internal function to fill particle data inside CONEX
+     * tables. Only e.m. particles are selected right now.
+     */
+    bool addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass,
+                     Point const& position, Vector<dimensionless_d> const& direction,
+                     TimeType t, double weight = 1);
+
+    CoordinateSystemPtr const& getObserverCS() const { return conexObservationCS_; }
+
+    YAML::Node getConfig() const final override;
+
+  private:
+    // data members
+    //! CONEX e.m. particle codes
+    static std::array<std::pair<Code, int>, 3> constexpr egs_em_codes_{
+        {{Code::Photon, 0}, {Code::Electron, -1}, {Code::Positron, -1}}};
+
+    Point const center_; //!< center of CONEX Earth
+    ShowerAxis const& showerAxis_;
+    LengthType groundDist_;       //!< length from injection point to shower core
+    LengthType injectionHeight_;  //!< starting height of primary particle
+    HEPEnergyType primaryEnergy_; //!< primary particle energy
+    PDGCode primaryPDG_;          //!< primary particle PDG
+    Point const showerCore_;      //!< shower core
+    CoordinateSystemPtr const conexObservationCS_; //!< CONEX observation frame
+    DirectionVector const x_sf_,
+        y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction
+  };
+} // namespace corsika
+
+#include <corsika/detail/modules/conex/CONEXhybrid.inl>

corsika/modules/conex/CONEXrandom.hpp

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <functional>
+
+namespace conex {
+  void set_rng_function(std::function<void(double*, std::size_t)>);
+}

corsika/modules/energy_loss/BetheBlochPDG.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/core/Step.hpp>
+
+#include <corsika/modules/writers/WriterOff.hpp>
+
+#include <map>
+
+namespace corsika {
+
+  /**
+   *   PDG2018, passage of particles through matter
+   *
+   * Note, that \f$I_{\mathrm{eff}}\f$ of composite media a determined from \f$ \ln I =
+   * \sum_i a_i \ln(I_i) \f$ where \f$ a_i \f$ is the fraction of the electron population
+   * (\f$\sim Z_i\f$) of the \f$i\f$-th element. This can also be used for shell
+   * corrections or density effects.
+   *
+   * The \f$I_{\mathrm{eff}}\f$ of compounds is not better than a few percent, if not
+   * measured explicitly.
+   *
+   * For shell correction, see Sec 6 of https://www.nap.edu/read/20066/chapter/8#115
+   *
+   */
+
+  template <typename TOutput = WriterOff>
+  class BetheBlochPDG : public ContinuousProcess<BetheBlochPDG<TOutput>>, public TOutput {
+
+    using MeVgcm2 = decltype(1e6 * electronvolt / gram * square(1e-2 * meter));
+
+  public:
+    template <typename... TOutputArgs>
+    BetheBlochPDG(TOutputArgs&&... args);
+
+    /**
+     * Interface function of ContinuousProcess.
+     *
+     * @param particle The particle to process in its current state
+     * @param track The trajectory in space of this particle, on which doContinuous
+     *should act
+     * @param limitFlag flag to identify, if BetheBlochPDG::getMaxStepLength is the
+     *        globally limiting factor (or not)
+     clang-format-on **/
+    template <typename TParticle>
+    ProcessReturn doContinuous(Step<TParticle>& step, bool const limitFlag);
+
+    template <typename TParticle, typename TTrajectory>
+    LengthType getMaxStepLength(TParticle const&, TTrajectory const&) const;
+
+    template <typename TParticle>
+    static HEPEnergyType getBetheBloch(TParticle const&, const GrammageType);
+
+    template <typename TParticle>
+    static HEPEnergyType getRadiationLosses(TParticle const&, const GrammageType);
+
+    template <typename TParticle>
+    static HEPEnergyType getTotalEnergyLoss(TParticle const&, const GrammageType);
+
+    YAML::Node getConfig() const override;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/modules/energy_loss/BetheBlochPDG.inl>