corsika/detail/modules/radio/LimitedRadioProcess.inl

+/*
+ * (c) Copyright 2024 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
+
+namespace corsika {
+
+  template <class TUnderlyingProcess>
+  inline LimitedRadioProcess<TUnderlyingProcess>::LimitedRadioProcess(
+      TUnderlyingProcess&& process,
+      std::function<LimitedRadioProcess::functor_signature> runThis)
+      : process_{std::forward<TUnderlyingProcess>(process)}
+      , runThis_{std::move(runThis)} {}
+
+  template <class TUnderlyingProcess>
+  template <typename Particle>
+  inline ProcessReturn LimitedRadioProcess<TUnderlyingProcess>::doContinuous(
+      const Step<Particle>& step, const bool arg) {
+
+    if ((step.getParticlePre().getPID() != Code::Electron) &&
+        (step.getParticlePre().getPID() != Code::Positron)) {
+      CORSIKA_LOG_DEBUG("Not e+/-, will not run RadioProcess");
+      return ProcessReturn::Ok;
+    }
+
+    // Check for running the wrapped RadioProcess function
+    if (runThis_(step.getPositionPre())) {
+      CORSIKA_LOG_TRACE("Will run RadioProcess for particle at {}",
+                        step.getPositionPre());
+      return process_.doContinuous(step, arg);
+    }
+    CORSIKA_LOG_DEBUG("Particle at {} is not selected by LimitedRadioProcess",
+                      step.getPositionPre());
+    return ProcessReturn::Ok;
+  }
+
+  template <class TUnderlyingProcess>
+  template <typename Particle, typename Track>
+  inline LengthType LimitedRadioProcess<TUnderlyingProcess>::getMaxStepLength(
+      Particle const& vParticle, Track const& vTrack) const {
+    return process_.getMaxStepLength(vParticle, vTrack);
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/radio/RadioProcess.inl

+/*
+ * (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/radio/RadioProcess.hpp>
+
+namespace corsika {
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline TRadioImpl&
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::implementation() {
+    return static_cast<TRadioImpl&>(*this);
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline TRadioImpl const&
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::implementation() const {
+    return static_cast<TRadioImpl const&>(*this);
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::RadioProcess(
+      TObserverCollection& observers, TPropagator& propagator)
+      : observers_(observers)
+      , propagator_(propagator) {}
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn RadioProcess<TObserverCollection, TRadioImpl,
+                                    TPropagator>::doContinuous(const Step<Particle>& step,
+                                                               const bool) {
+
+    // return immediately if radio process does not have any observers
+    if (observers_.size() == 0) return ProcessReturn::Ok;
+
+    // we want the following particles:
+    // Code::Electron & Code::Positron
+
+    // we wrap Simulate() in doContinuous as the plan is to add particle level
+    // filtering or thinning for calculation of the radio emission. This is
+    // important for controlling the runtime of radio (by ignoring particles
+    // that aren't going to contribute i.e. heavy hadrons)
+    // if (valid(step)) {
+    auto const particleID_{step.getParticlePre().getPID()};
+    if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
+      CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
+      return this->implementation().simulate(step);
+    } else {
+      CORSIKA_LOG_DEBUG("Particle {} is irrelevant for radio", particleID_);
+      return ProcessReturn::Ok;
+    }
+    //}
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  template <typename Particle, typename Track>
+  inline LengthType
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::getMaxStepLength(
+      [[maybe_unused]] const Particle& vParticle,
+      [[maybe_unused]] const Track& vTrack) const {
+    return meter * std::numeric_limits<double>::infinity();
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline void RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::startOfLibrary(
+      const boost::filesystem::path& directory) {
+
+    // setup the streamer
+    output_.initStreamer((directory / ("observers.parquet")).string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // LCOV_EXCL_START
+    // build the schema
+    output_.addField("Time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ex", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ey", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ez", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+    // LCOV_EXCL_STOP
+    // and build the streamer
+    output_.buildStreamer();
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline void RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::endOfShower(
+      const unsigned int) {
+
+    // loop over every observer and instruct them to
+    // flush data to disk, and then reset the observer
+    // before the next event
+    for (auto& observer : observers_.getObservers()) {
+
+      auto const sampleRate = observer.getSampleRate() * 1_s;
+      auto const radioImplementation =
+          static_cast<std::string>(this->implementation().algorithm);
+
+      // get the axis labels for this observer and write the first row.
+      axistype axis = observer.implementation().getAxis();
+
+      // get the copy of the waveform data for this event
+      std::vector<double> const& dataX = observer.implementation().getWaveformX();
+      std::vector<double> const& dataY = observer.implementation().getWaveformY();
+      std::vector<double> const& dataZ = observer.implementation().getWaveformZ();
+
+      // check for the axis name
+      std::string label = "Unknown";
+      if (observer.getDomainLabel() == "Time") {
+        label = "Time";
+      }
+      // LCOV_EXCL_START
+      else if (observer.getDomainLabel() == "Frequency") {
+        label = "Frequency";
+      }
+      // LCOV_EXCL_STOP
+      if (radioImplementation == "ZHS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          auto time = (axis.at(i + 1) + axis.at(i)) / 2.;
+          auto Ex = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
+          auto Ey = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
+          auto Ez = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
+
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(time) << static_cast<double>(Ex)
+              << static_cast<double>(Ey) << static_cast<double>(Ez) << parquet::EndRow;
+        }
+      } else if (radioImplementation == "CoREAS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(axis[i])
+              << static_cast<double>(dataX[i]) << static_cast<double>(dataY[i])
+              << static_cast<double>(dataZ[i]) << parquet::EndRow;
+        }
+      }
+
+      observer.reset();
+    }
+    output_.closeStreamer();
+
+    // increment our event counter
+    showerId_++;
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline YAML::Node
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::getConfig() const {
+
+    // top-level YAML node
+    YAML::Node config;
+
+    // fill in some basics
+    config["type"] = "RadioProcess";
+    config["algorithm"] = this->implementation().algorithm;
+    config["units"]["time"] = "ns";
+    config["units"]["frequency"] = "GHz";
+    config["units"]["electric field"] = "V/m";
+    config["units"]["distance"] = "m";
+
+    for (auto& observer : observers_.getObservers()) {
+      // get the name/location of this observer
+      auto name = observer.getName();
+      auto location = observer.getLocation().getCoordinates();
+
+      // get the observers config
+      config["observers"][name] = observer.getConfig();
+
+      // write the location of this observer
+      config["observers"][name]["location"].push_back(location.getX() / 1_m);
+      config["observers"][name]["location"].push_back(location.getY() / 1_m);
+      config["observers"][name]["location"].push_back(location.getZ() / 1_m);
+    }
+
+    return config;
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/ZHS.inl

+/*
+ * (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/radio/ZHS.hpp>
+
+namespace corsika {
+
+  template <typename TRadioDetector, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn ZHS<TRadioDetector, TPropagator>::simulate(
+      Step<Particle> const& step) {
+    auto const startTime{step.getTimePre()};
+    auto const endTime{step.getTimePost()};
+
+    // LCOV_EXCL_START
+    if (startTime == endTime) {
+      return ProcessReturn::Ok;
+      // LCOV_EXCL_STOP
+    } else {
+
+      auto const startPoint{step.getPositionPre()};
+      auto const endPoint{step.getPositionPost()};
+      LengthType const trackLength{(startPoint - endPoint).getNorm()};
+
+      auto const betaModule{(endPoint - startPoint).getNorm() /
+                            (constants::c * (endTime - startTime))};
+      auto const beta{(endPoint - startPoint).normalized() * betaModule};
+
+      auto const charge{get_charge(step.getParticlePre().getPID())};
+
+      // // get "mid" position of the track geometrically
+
+      auto const halfVector{(startPoint - endPoint) / 2};
+      auto const midPoint{endPoint + halfVector};
+
+      // get thinning weight
+      auto const thinningWeight{step.getParticlePre().getWeight()};
+
+      auto const constants{charge * emConstant_ * thinningWeight};
+
+      // we loop over each observer in the collection
+      for (auto& observer : observers_.getObservers()) {
+        auto midPaths{this->propagator_.propagate(step.getParticlePre(), midPoint,
+                                                  observer.getLocation())};
+        // Loop over midPaths, first check Fraunhoffer limit
+        for (size_t i{0}; i < midPaths.size(); i++) {
+          double const uTimesK{beta.dot(midPaths[i].emit_) / betaModule};
+          double const sinTheta2{1. - uTimesK * uTimesK};
+          LengthType const lambda{constants::c / observer.getSampleRate()};
+          double const fraunhLimit{sinTheta2 * trackLength * trackLength /
+                                   midPaths[i].R_distance_ / lambda * 2 * M_PI};
+          // Checks if we are in fraunhoffer domain
+          if (fraunhLimit > 1.0) {
+            /// code for dividing track and calculating field.
+            double const nSubTracks{sqrt(fraunhLimit) + 1};
+            auto const step_{(endPoint - startPoint) / nSubTracks};
+            TimeType const timeStep{(endTime - startTime) / nSubTracks};
+            // energy should be divided up when it is possible to get the energy at end of
+            // track!!!!
+            auto point1{startPoint};
+            TimeType time1{startTime};
+            for (int j{0}; j < nSubTracks; j++) {
+              auto const point2{point1 + step_};
+              TimeType const time2{time1 + timeStep};
+              auto const newHalfVector{(point1 - point2) / 2.};
+              auto const newMidPoint{point2 + newHalfVector};
+              auto const newMidPaths{this->propagator_.propagate(
+                  step.getParticlePre(), newMidPoint, observer.getLocation())};
+              // A function for calculating the field should be made since it is repeated
+              // later
+              for (size_t k{0}; k < newMidPaths.size(); k++) {
+
+                double const n_source{newMidPaths[k].refractive_index_source_};
+
+                double const betaTimesK{beta.dot(newMidPaths[k].emit_)};
+                TimeType const midTime{(time1 + time2) / 2.};
+                TimeType detectionTime1{time1 + newMidPaths[k].propagation_time_ -
+                                        n_source * betaTimesK * (time1 - midTime)};
+                TimeType detectionTime2{time2 + newMidPaths[k].propagation_time_ -
+                                        n_source * betaTimesK * (time2 - midTime)};
+
+                // make detectionTime1_ be the smallest time =>
+                // changes step function order so sign is changed to account for it
+                double sign = 1.;
+                if (detectionTime1 > detectionTime2) {
+                  detectionTime1 = time2 + newMidPaths[k].propagation_time_ -
+                                   n_source * betaTimesK * (time2 - midTime);
+                  detectionTime2 = time1 + newMidPaths[k].propagation_time_ -
+                                   n_source * betaTimesK * (time1 - midTime);
+                  sign = -1.;
+                } // end if statement for time structure
+
+                double const startBin{
+                    std::floor((detectionTime1 - observer.getStartTime()) *
+                                   observer.getSampleRate() +
+                               0.5)};
+                double const endBin{
+                    std::floor((detectionTime2 - observer.getStartTime()) *
+                                   observer.getSampleRate() +
+                               0.5)};
+
+                auto const betaPerp{
+                    newMidPaths[k].emit_.cross(beta.cross(newMidPaths[k].emit_))};
+                double const denominator{1. - n_source * betaTimesK};
+
+                if (startBin == endBin) {
+                  // track contained in bin
+                  // if not in Cerenkov angle then
+                  if (std::fabs(denominator) > 1.e-15) {
+                    double const f{
+                        std::fabs((detectionTime2 * observer.getSampleRate() -
+                                   detectionTime1 * observer.getSampleRate()))};
+                    VectorPotential const Vp = betaPerp * sign * constants * f /
+                                               denominator / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime2, betaPerp, Vp);
+                  } else { // If emission in Cerenkov angle => approximation
+                    double const f{time2 * observer.getSampleRate() -
+                                   time1 * observer.getSampleRate()};
+                    VectorPotential const Vp =
+                        betaPerp * sign * constants * f / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime2, betaPerp, Vp);
+                  } // end if Cerenkov angle approx
+                } else {
+                  /*Track is contained in more than one bin*/
+                  int const numberOfBins{static_cast<int>(endBin - startBin)};
+                  // first contribution/ plus 1 bin minus 0.5 from new observer ruonding
+                  double f{std::fabs(startBin + 0.5 -
+                                     (detectionTime1 - observer.getStartTime()) *
+                                         observer.getSampleRate())};
+                  VectorPotential Vp = betaPerp * sign * constants * f / denominator /
+                                       newMidPaths[k].R_distance_;
+                  observer.receive(detectionTime1, betaPerp, Vp);
+                  // intermidiate contributions
+                  for (int it{1}; it < numberOfBins; ++it) {
+                    Vp = betaPerp * constants / denominator / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime1 + static_cast<double>(it) /
+                                                          observer.getSampleRate(),
+                                     betaPerp, Vp);
+                  } // end loop over bins in which potential vector is not zero
+                  // final contribution// f +0.5 from new observer rounding
+                  f = std::fabs((detectionTime2 - observer.getStartTime()) *
+                                    observer.getSampleRate() +
+                                0.5 - endBin);
+                  Vp = betaPerp * sign * constants * f / denominator /
+                       newMidPaths[k].R_distance_;
+                  observer.receive(detectionTime2, betaPerp, Vp);
+                } // end if statement for track in multiple bins
+
+              } // end of loop over newMidPaths
+              // update points for next sub track
+              point1 = point1 + step_;
+              time1 = time1 + timeStep;
+            }
+
+          } else // Calculate vector potential of whole track
+          {
+            double const n_source{midPaths[i].refractive_index_source_};
+
+            double const betaTimesK{beta.dot(midPaths[i].emit_)};
+            TimeType const midTime{(startTime + endTime) / 2};
+            TimeType detectionTime1{startTime + midPaths[i].propagation_time_ -
+                                    n_source * betaTimesK * (startTime - midTime)};
+            TimeType detectionTime2{endTime + midPaths[i].propagation_time_ -
+                                    n_source * betaTimesK * (endTime - midTime)};
+
+            // make detectionTime1_ be the smallest time =>
+            // changes step function order so sign is changed to account for it
+            double sign = 1.;
+            if (detectionTime1 > detectionTime2) {
+              detectionTime1 = endTime + midPaths[i].propagation_time_ -
+                               n_source * betaTimesK * (endTime - midTime);
+              detectionTime2 = startTime + midPaths[i].propagation_time_ -
+                               n_source * betaTimesK * (startTime - midTime);
+              sign = -1.;
+            } // end if statement for time structure
+
+            double const startBin{std::floor((detectionTime1 - observer.getStartTime()) *
+                                                 observer.getSampleRate() +
+                                             0.5)};
+            double const endBin{std::floor((detectionTime2 - observer.getStartTime()) *
+                                               observer.getSampleRate() +
+                                           0.5)};
+
+            auto const betaPerp{midPaths[i].emit_.cross(beta.cross(midPaths[i].emit_))};
+            double const denominator{1. -
+                                     midPaths[i].refractive_index_source_ * betaTimesK};
+
+            if (startBin == endBin) {
+              // track contained in bin
+              // if not in Cerenkov angle then
+              if (std::fabs(denominator) > 1.e-15) {
+                double const f{std::fabs((detectionTime2 * observer.getSampleRate() -
+                                          detectionTime1 * observer.getSampleRate()))};
+
+                VectorPotential const Vp = betaPerp * sign * constants * f / denominator /
+                                           midPaths[i].R_distance_;
+                observer.receive(detectionTime2, betaPerp, Vp);
+              } else { // If emission in Cerenkov angle => approximation
+                double const f{endTime * observer.getSampleRate() -
+                               startTime * observer.getSampleRate()};
+                VectorPotential const Vp =
+                    betaPerp * sign * constants * f / midPaths[i].R_distance_;
+                observer.receive(detectionTime2, betaPerp, Vp);
+              } // end if Cerenkov angle approx
+            } else {
+              /*Track is contained in more than one bin*/
+              int const numberOfBins{static_cast<int>(endBin - startBin)};
+              // TODO: should we check for Cerenkov angle?
+              // first contribution
+              double f{std::fabs(startBin + 0.5 -
+                                 (detectionTime1 - observer.getStartTime()) *
+                                     observer.getSampleRate())};
+              VectorPotential Vp =
+                  betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+              observer.receive(detectionTime1, betaPerp, Vp);
+              // intermediate contributions
+              for (int it{1}; it < numberOfBins; ++it) {
+                Vp = betaPerp * sign * constants / denominator / midPaths[i].R_distance_;
+                observer.receive(
+                    detectionTime1 + static_cast<double>(it) / observer.getSampleRate(),
+                    betaPerp, Vp);
+              } // end loop over bins in which potential vector is not zero
+              // final contribution
+              f = std::fabs((detectionTime2 - observer.getStartTime()) *
+                                observer.getSampleRate() +
+                            0.5 - endBin);
+              Vp =
+                  betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+              observer.receive(detectionTime2, betaPerp, Vp);
+            } // end if statement for track in multiple bins
+
+          } // finish if statement of track in fraunhoffer or not
+
+        } // end loop over mid paths
+
+      } // END: loop over observers
+      return ProcessReturn::Ok;
+    }
+  } // end simulate
+
+} // namespace corsika

corsika/detail/modules/radio/detectors/ObserverCollection.inl

+/*
+ * (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/radio/detectors/ObserverCollection.hpp>
+
+namespace corsika {
+
+  template <typename TObserverImpl>
+  inline void ObserverCollection<TObserverImpl>::addObserver(
+      TObserverImpl const& observer) {
+    observers_.push_back(observer);
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl& ObserverCollection<TObserverImpl>::at(std::size_t const i) {
+    return observers_.at(i);
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl const& ObserverCollection<TObserverImpl>::at(
+      std::size_t const i) const {
+    return observers_.at(i);
+  }
+
+  template <typename TObserverImpl>
+  inline int ObserverCollection<TObserverImpl>::size() const {
+    return observers_.size();
+  }
+
+  template <typename TObserverImpl>
+  inline std::vector<TObserverImpl>& ObserverCollection<TObserverImpl>::getObservers() {
+    return observers_;
+  }
+
+  template <typename TObserverImpl>
+  inline void ObserverCollection<TObserverImpl>::reset() {
+    std::for_each(observers_.begin(), observers_.end(),
+                  std::mem_fn(&TObserverImpl::reset));
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/observers/Observer.inl

+/*
+ * (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/radio/observers/Observer.hpp>
+
+namespace corsika {
+
+  template <typename TObserverImpl>
+  inline Observer<TObserverImpl>::Observer(std::string const& name, Point const& location,
+                                           CoordinateSystemPtr const& coordinateSystem)
+      : name_(name)
+      , location_(location)
+      , coordinateSystem_(coordinateSystem) {}
+
+  template <typename TObserverImpl>
+  inline Point const& Observer<TObserverImpl>::getLocation() const {
+    return location_;
+  }
+
+  template <typename TObserverImpl>
+  inline std::string const& Observer<TObserverImpl>::getName() const {
+    return name_;
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl& Observer<TObserverImpl>::implementation() {
+    return static_cast<TObserverImpl&>(*this);
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/observers/TimeDomainObserver.inl

+/*
+ * (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/radio/observers/TimeDomainObserver.hpp>
+
+namespace corsika {
+
+  inline TimeDomainObserver::TimeDomainObserver(std::string const& name,
+                                                Point const& location,
+                                                CoordinateSystemPtr coordinateSystem,
+                                                TimeType const& start_time,
+                                                TimeType const& duration,
+                                                InverseTimeType const& sample_rate,
+                                                TimeType const ground_hit_time)
+      : Observer(name, location, coordinateSystem)
+      , start_time_(start_time)
+      , duration_(std::abs(duration / 1_s) * 1_s)
+      , sample_rate_(std::abs(sample_rate / 1_Hz) * 1_Hz)
+      , ground_hit_time_(ground_hit_time)
+      , num_bins_(static_cast<std::size_t>(duration_ * sample_rate_ + 1.5l))
+      , waveformEX_(num_bins_, 0)
+      , waveformEY_(num_bins_, 0)
+      , waveformEZ_(num_bins_, 0)
+      , time_axis_(createTimeAxis()) {
+    if (0_s == duration_) {
+      CORSIKA_LOG_WARN(
+          "Observer: \"{}\" has a duration of zero. Nothing will be injected into it",
+          name);
+    } else if (duration_ != duration) {
+      CORSIKA_LOG_WARN(
+          "Observer: \"{}\" was given a negative duration. Set to absolute value.", name);
+    }
+
+    if (sample_rate_ != sample_rate) {
+      CORSIKA_LOG_WARN(
+          "Observer: \"{}\" was given a negative sampling rate. Set to absolute value.",
+          name);
+    }
+  }
+
+  inline void TimeDomainObserver::receive(
+      const TimeType time, [[maybe_unused]] const Vector<dimensionless_d>& receive_vector,
+      const ElectricFieldVector& efield) {
+
+    if (time < start_time_ || time > (start_time_ + duration_)) {
+      return;
+    } else {
+      // figure out the correct timebin to store the E-field value.
+      // NOTE: static cast is implicitly flooring
+      auto timebin{static_cast<std::size_t>(
+          std::floor((time - start_time_) * sample_rate_ + 0.5l))};
+
+      // store the x,y,z electric field components.
+      auto const& Electric_field_components{efield.getComponents(coordinateSystem_)};
+      waveformEX_.at(timebin) += (Electric_field_components.getX() * (1_m / 1_V));
+      waveformEY_.at(timebin) += (Electric_field_components.getY() * (1_m / 1_V));
+      waveformEZ_.at(timebin) += (Electric_field_components.getZ() * (1_m / 1_V));
+      // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+      // a 3D object
+    }
+  }
+
+  inline void TimeDomainObserver::receive(
+      const TimeType time, [[maybe_unused]] const Vector<dimensionless_d>& receive_vector,
+      const VectorPotential& vectorP) {
+
+    if (time < start_time_ || time > (start_time_ + duration_)) {
+      return;
+    } else {
+      // figure out the correct timebin to store the E-field value.
+      // NOTE: static cast is implicitly flooring
+      auto timebin{static_cast<std::size_t>(
+          std::floor((time - start_time_) * sample_rate_ + 0.5l))};
+
+      // store the x,y,z electric field components.
+      auto const& Vector_potential_components{vectorP.getComponents(coordinateSystem_)};
+      waveformEX_.at(timebin) +=
+          (Vector_potential_components.getX() * (1_m / (1_V * 1_s)));
+      waveformEY_.at(timebin) +=
+          (Vector_potential_components.getY() * (1_m / (1_V * 1_s)));
+      waveformEZ_.at(timebin) +=
+          (Vector_potential_components.getZ() * (1_m / (1_V * 1_s)));
+      // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+      // a 3D object
+    }
+  }
+
+  inline auto const& TimeDomainObserver::getWaveformX() const { return waveformEX_; }
+
+  inline auto const& TimeDomainObserver::getWaveformY() const { return waveformEY_; }
+
+  inline auto const& TimeDomainObserver::getWaveformZ() const { return waveformEZ_; }
+
+  inline std::string const TimeDomainObserver::getDomainLabel() { return "Time"; }
+
+  inline std::vector<long double> TimeDomainObserver::createTimeAxis() const {
+
+    // create a 1-D xtensor to store time values so we can print them later.
+    std::vector<long double> times(num_bins_, 0);
+
+    // calculate the sample_period
+    auto sample_period{1 / sample_rate_};
+
+    // fill in every time-value
+    for (uint64_t i = 0; i < num_bins_; i++) {
+      // create the current time in nanoseconds
+      times.at(i) = static_cast<long double>(
+          ((start_time_ - ground_hit_time_) + i * sample_period) / 1_ns);
+    }
+
+    return times;
+  }
+
+  inline auto const TimeDomainObserver::getAxis() const { return time_axis_; }
+
+  inline InverseTimeType const& TimeDomainObserver::getSampleRate() const {
+    return sample_rate_;
+  }
+
+  inline TimeType const& TimeDomainObserver::getStartTime() const { return start_time_; }
+
+  inline void TimeDomainObserver::reset() {
+    std::fill(waveformEX_.begin(), waveformEX_.end(), 0);
+    std::fill(waveformEY_.begin(), waveformEY_.end(), 0);
+    std::fill(waveformEZ_.begin(), waveformEZ_.end(), 0);
+  }
+
+  inline YAML::Node TimeDomainObserver::getConfig() const {
+
+    // top-level config
+    YAML::Node config;
+
+    config["type"] = "TimeDomainObserver";
+    config["start time"] = start_time_ / 1_ns;
+    config["duration"] = duration_ / 1_ns;
+    config["number of bins"] = duration_ * sample_rate_;
+    config["sampling frequency"] = sample_rate_ / 1_GHz;
+
+    return config;
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/DummyTestPropagator.inl

+/*
+ * (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/radio/propagators/DummyTestPropagator.hpp>
+
+namespace corsika {
+
+  template <typename TEnvironment>
+  inline DummyTestPropagator<TEnvironment>::DummyTestPropagator(TEnvironment const& env)
+      : RadioPropagator<DummyTestPropagator, TEnvironment>(env) {
+    rindex.reserve(2);
+  }
+
+  template <typename TEnvironment>
+  template <typename Particle>
+  inline typename DummyTestPropagator<TEnvironment>::SignalPathCollection
+  DummyTestPropagator<TEnvironment>::propagate(Particle const& particle,
+                                               Point const& source,
+                                               Point const& destination) {
+
+    /**
+     * This is the simplest case of straight propagator
+     * where no integration takes place.
+     * This can be used for fast tests and checks of the radio module.
+     *
+     */
+
+    // these are used for the direction of emission and reception of signal at the
+    // observer
+    auto const emit_{(destination - source).normalized()};
+    auto const receive_{-emit_};
+
+    // the geometrical distance from the point of emission to an observer
+    auto const distance_{(destination - source).getNorm()};
+
+    // get the universe for this environment
+    auto const* const universe{Base::env_.getUniverse().get()};
+
+    // clear the refractive index vector and points deque for this signal propagation.
+    rindex.clear();
+    points.clear();
+
+    // get and store the refractive index of the first point 'source'.
+    auto const* const nodeSource{particle.getNode()};
+    auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
+    rindex.push_back(ri_source);
+    points.push_back(source);
+
+    // add the refractive index of last point 'destination' and store it.
+    auto const* const node{universe->getContainingNode(destination)};
+    auto const ri_destination{node->getModelProperties().getRefractiveIndex(destination)};
+    rindex.push_back(ri_destination);
+    points.push_back(destination);
+
+    // compute the average refractive index.
+    auto const averageRefractiveIndex_ = (ri_source + ri_destination) * 0.5;
+
+    // compute the total time delay.
+    TimeType const time = averageRefractiveIndex_ * (distance_ / constants::c);
+
+    return std::vector<SignalPath>(
+        1, SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
+                      receive_, distance_, points));
+
+  } // END: propagate()
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/NumericalIntegratingPropagator.inl

+/*
+ * (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/radio/propagators/NumericalIntegratingPropagator.hpp>
+
+namespace corsika {
+
+  template <typename TEnvironment>
+  // TODO: maybe the constructor doesn't take any arguments for the environment (?)
+  inline NumericalIntegratingPropagator<TEnvironment>::NumericalIntegratingPropagator(
+      TEnvironment const& env, LengthType const stepsize)
+      : RadioPropagator<NumericalIntegratingPropagator, TEnvironment>(env)
+      , stepsize_(stepsize) {}
+
+  template <typename TEnvironment>
+  template <typename Particle>
+  inline typename NumericalIntegratingPropagator<TEnvironment>::SignalPathCollection
+  NumericalIntegratingPropagator<TEnvironment>::propagate(
+      Particle const& particle, Point const& source, Point const& destination) const {
+
+    /*
+     * get the normalized (unit) vector from `source` to `destination'.
+     * this is also the `emit` and `receive` vectors in the SignalPath class.
+     * in this case emit and receive unit vectors should be the same
+     * so they are both called direction
+     */
+
+    // these are used for the direction of emission and reception of signal at the
+    // observer
+    auto const emit{(destination - source).normalized()};
+    auto const receive{-emit};
+
+    // the distance from the point of emission to an observer
+    auto const distance{(destination - source).getNorm()};
+
+    try {
+      if (stepsize_ <= 0.5 * distance) {
+
+        // "step" is the direction vector with length `stepsize`
+        auto const step{emit * stepsize_};
+
+        // calculate the number of points (roughly) for the numerical integration
+        auto const n_points{(destination - source).getNorm() / stepsize_};
+
+        // get the universe for this environment
+        auto const* const universe{Base::env_.getUniverse().get()};
+
+        // the points that we build along the way for the numerical integration
+        std::deque<Point> points;
+
+        // store value of the refractive index at points
+        std::vector<double> rindex;
+        rindex.reserve(n_points);
+
+        // get and store the refractive index of the first point 'source'
+        auto const* const nodeSource{particle.getNode()};
+        auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
+        rindex.push_back(ri_source);
+        points.push_back(source);
+
+        // loop from `source` to `destination` to store values before Simpson's rule.
+        // this loop skips the last point 'destination' and "misses" the extra point
+        for (auto point = source + step;
+             (point - destination).getNorm() > 0.6 * stepsize_; point = point + step) {
+
+          // get the environment node at this specific 'point'
+          auto const* const node{universe->getContainingNode(point)};
+
+          // get the associated refractivity at 'point'
+          auto const refractive_index{
+              node->getModelProperties().getRefractiveIndex(point)};
+          rindex.push_back(refractive_index);
+
+          // add this 'point' to our deque collection
+          points.push_back(point);
+        }
+
+        // Get the extra points that the for loop misses until the destination
+        auto const extrapoint_{points.back() + step};
+
+        // add the refractive index of last point 'destination' and store it
+        auto const* const node{universe->getContainingNode(destination)};
+        auto const ri_destination{
+            node->getModelProperties().getRefractiveIndex(destination)};
+        rindex.push_back(ri_destination);
+        points.push_back(destination);
+
+        auto N = rindex.size();
+        std::size_t index = 0;
+        double sum = rindex.at(index);
+        auto refra = rindex.at(index);
+        TimeType time{0_s};
+
+        if ((N - 1) % 2 == 0) {
+          // Apply the standard Simpson's rule
+          auto const h = ((destination - source).getNorm()) / (N - 1);
+
+          for (std::size_t index = 1; index < (N - 1); index += 2) {
+            sum += 4 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          for (std::size_t index = 2; index < (N - 1); index += 2) {
+            sum += 2 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          index = N - 1;
+          sum = sum + rindex.at(index);
+          refra += rindex.at(index);
+
+          // compute the total time delay.
+          time = sum * (h / (3 * constants::c));
+        } else {
+          // Apply Simpson's rule for one "extra" point and then subtract the difference
+          points.pop_back();
+          rindex.pop_back();
+          auto const* const node{universe->getContainingNode(extrapoint_)};
+          auto const ri_extrapoint{
+              node->getModelProperties().getRefractiveIndex(extrapoint_)};
+          rindex.push_back(ri_extrapoint);
+          points.push_back(extrapoint_);
+          auto const extrapoint2_{extrapoint_ + step};
+          auto const* const node2{universe->getContainingNode(extrapoint2_)};
+          auto const ri_extrapoint2{
+              node2->getModelProperties().getRefractiveIndex(extrapoint2_)};
+          rindex.push_back(ri_extrapoint2);
+          points.push_back(extrapoint2_);
+          N = rindex.size();
+          auto const h = ((extrapoint2_ - source).getNorm()) / (N - 1);
+          for (std::size_t index = 1; index < (N - 1); index += 2) {
+            sum += 4 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          for (std::size_t index = 2; index < (N - 1); index += 2) {
+            sum += 2 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          index = N - 1;
+          sum = sum + rindex.at(index);
+          refra += rindex.at(index);
+
+          // compute the total time delay including the correction
+          time =
+              sum * (h / (3 * constants::c)) -
+              (ri_extrapoint2 * ((extrapoint2_ - destination).getNorm()) / constants::c);
+        }
+
+        // compute the average refractive index.
+        auto const averageRefractiveIndex = refra / N;
+
+        return std::vector<SignalPath>(
+            1, SignalPath(time, averageRefractiveIndex, ri_source, ri_destination, emit,
+                          receive, distance, points));
+      } else {
+        throw stepsize_;
+      }
+    } catch (const LengthType& s) {
+      CORSIKA_LOG_ERROR("Please choose a smaller stepsize for the numerical integration");
+    }
+
+    std::deque<Point> const defaultPoints;
+    return std::vector<SignalPath>(
+        1, SignalPath(0_s, 0, 0, 0, emit, receive, distance,
+                      defaultPoints)); // Dummy return that is never called (combat
+                                       // compile warnings)
+  }                                    // END: propagate()
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/RadioPropagator.inl

+/*
+ * (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/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+  template <typename TImpl, typename TEnvironment>
+  inline RadioPropagator<TImpl, TEnvironment>::RadioPropagator(TEnvironment const& env)
+      : env_(env) {}
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/SignalPath.inl

+/*
+ * (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
+
+namespace corsika {
+
+  inline SignalPath::SignalPath(
+      TimeType const propagation_time, double const average_refractive_index,
+      double const refractive_index_source, double const refractive_index_destination,
+      Vector<dimensionless_d> const& emit, Vector<dimensionless_d> const& receive,
+      LengthType const R_distance, std::deque<Point> const& points)
+      : Path(points)
+      , propagation_time_(propagation_time)
+      , average_refractive_index_(average_refractive_index)
+      , refractive_index_source_(refractive_index_source)
+      , refractive_index_destination_(refractive_index_destination)
+      , emit_(emit)
+      , receive_(receive)
+      , R_distance_(R_distance) {}
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl

+/*
+ * (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 <corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp>
+
+namespace corsika {
+
+  template <typename TEnvironment>
+  inline TabulatedFlatAtmospherePropagator<
+      TEnvironment>::TabulatedFlatAtmospherePropagator(TEnvironment const& env,
+                                                       Point const& upperLimit,
+                                                       Point const& lowerLimit,
+                                                       LengthType const step)
+      : RadioPropagator<TabulatedFlatAtmospherePropagator, TEnvironment>(env)
+      , upperLimit_(upperLimit)
+      , lowerLimit_(lowerLimit)
+      , step_(step)
+      , inverseStep_(1 / step)
+      , maxHeight_(upperLimit_.getCoordinates().getZ())
+      , minHeight_(lowerLimit_.getCoordinates().getZ() - 1_km) {
+    auto const minX_ = lowerLimit_.getCoordinates().getX();
+    auto const minY_ = lowerLimit_.getCoordinates().getY();
+    std::size_t const nBins_ = (maxHeight_ - minHeight_) * inverseStep_ + 1;
+    refractivityTable_.reserve(nBins_);
+    heightTable_.reserve(nBins_);
+    integratedRefractivityTable_.reserve(nBins_);
+
+    // get the root coordinate system of this environment
+    CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+    // get the universe for this environment
+    auto const* const universe{Base::env_.getUniverse().get()};
+
+    for (std::size_t i = 0; i < nBins_; i++) {
+      Point point_{rootCS, minX_, minY_, minHeight_ + i * step_};
+      auto const* const node{universe->getContainingNode(point_)};
+      auto const ri_ = node->getModelProperties().getRefractiveIndex(point_);
+      refractivityTable_.push_back(ri_ - 1);
+      auto const height_ = minHeight_ + i * step_;
+      heightTable_.push_back(height_);
+    }
+
+    double const stepOverMeter_{inverseStep_ * 1_m};
+    auto const intRerfZero_ = refractivityTable_.at(0) * stepOverMeter_;
+    integratedRefractivityTable_.push_back(intRerfZero_);
+    for (std::size_t i = 1; i < nBins_; i++) {
+      auto const intRefrI_ =
+          integratedRefractivityTable_[i - 1] + refractivityTable_[i] * stepOverMeter_;
+      integratedRefractivityTable_.push_back(intRefrI_);
+    }
+
+    // this is used to interpolate refractivity and integrated refractivity for particles
+    // below the lower limit
+    slopeRefrLower_ = (refractivityTable_.at(10) - refractivityTable_.front()) / 10.;
+    slopeIntRefrLower_ =
+        (integratedRefractivityTable_.at(10) - integratedRefractivityTable_.front()) /
+        10.;
+
+    // this is used to interpolate refractivity and integrated refractivity for particles
+    // above the upper limit
+    lastElement_ = integratedRefractivityTable_.size() - 1;
+    slopeRefrUpper_ =
+        (refractivityTable_.at(lastElement_) - refractivityTable_.at(lastElement_ - 10)) /
+        10.;
+    slopeIntRefrUpper_ = (integratedRefractivityTable_.at(lastElement_) -
+                          integratedRefractivityTable_.at(lastElement_ - 10)) /
+                         10.;
+  }
+
+  template <typename TEnvironment>
+  template <typename Particle>
+  inline typename TabulatedFlatAtmospherePropagator<TEnvironment>::SignalPathCollection
+  TabulatedFlatAtmospherePropagator<TEnvironment>::propagate(
+      [[maybe_unused]] Particle const& particle, Point const& source,
+      Point const& destination) {
+
+    /**
+     * This is a simple case of straight propagator where
+     * tabulated values of refractive index are called assuming
+     * a flat atmosphere.
+     *
+     */
+
+    // these are used for the direction of emission and reception of signal at the
+    // observer
+    auto const emit_{(destination - source).normalized()};
+    auto const receive_{-emit_};
+
+    // the geometrical distance from the point of emission to an observer
+    auto const distance_{(destination - source).getNorm()};
+
+    // clear the refractive index vector and points deque for this signal propagation.
+    rindex.clear();
+    points.clear();
+    auto const sourceZ_{source.getCoordinates().getZ()};
+    auto const checkMaxZ_{(maxHeight_ - sourceZ_) / 1_m};
+    auto ri_source{1.};
+    auto intRerfr_source{1.};
+    auto ri_destination{1.};
+    auto intRerfr_destination{1.};
+    auto height_{1.};
+
+    double const sourceHeight_{(sourceZ_ - heightTable_.front()) * inverseStep_};
+    double const destinationHeight_{
+        (destination.getCoordinates().getZ() - heightTable_.front()) * inverseStep_};
+
+    if ((sourceHeight_ + 0.5) >=
+        lastElement_) { // source particle is above maximum height
+      ri_source =
+          refractivityTable_.back() +
+          slopeRefrUpper_ * std::abs((sourceZ_ - heightTable_.back()) * inverseStep_) + 1;
+      intRerfr_source =
+          integratedRefractivityTable_.back() + slopeIntRefrUpper_ * std::abs(checkMaxZ_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      ri_destination = refractivityTable_.at(indexDestination_) + 1;
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = sourceHeight_ - destinationHeight_;
+
+    } else if ((sourceHeight_ + 0.5 < lastElement_) &&
+               sourceHeight_ > 0) { // source particle in the table
+      // get and store the refractive index of the first point 'source'.
+      std::size_t const indexSource_{static_cast<std::size_t>(sourceHeight_ + 0.5)};
+      ri_source = refractivityTable_.at(indexSource_) + 1;
+      intRerfr_source = integratedRefractivityTable_.at(indexSource_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      ri_destination = refractivityTable_.at(indexDestination_) + 1;
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ =
+          (heightTable_.at(indexSource_) - heightTable_.at(indexDestination_)) / 1_m;
+      if (height_ == 0) { height_ = 1.; }
+    } else if (sourceHeight_ ==
+               0) { // source particle is exactly at the lower edge of the table
+      // get and store the refractive index of the first point 'source'.
+      std::size_t const indexSource_{static_cast<std::size_t>(sourceHeight_ + 0.5)};
+      auto const ri_source{refractivityTable_.at(indexSource_) + 1};
+      intRerfr_source = integratedRefractivityTable_.at(indexSource_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = destinationHeight_ - sourceHeight_;
+    } else if (sourceHeight_ < 0) { // source particle is in the ground.
+      // get and store the refractive index of the first point 'source'.
+      ri_source =
+          refractivityTable_.front() + slopeRefrLower_ * std::abs(sourceHeight_) + 1;
+      intRerfr_source = integratedRefractivityTable_.front() +
+                        slopeIntRefrLower_ * std::abs(sourceHeight_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = destinationHeight_ - std::abs(sourceHeight_);
+    }
+
+    // compute the average refractive index.
+    auto const averageRefractiveIndex_ = (ri_source + ri_destination) * 0.5;
+
+    // compute the total time delay.
+    TimeType const time = (1 + (intRerfr_source - intRerfr_destination) / height_) *
+                          (distance_ / constants::c);
+
+    return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
+                       receive_, distance_, points)};
+
+  } // END: propagate()
+
+} // namespace corsika

corsika/detail/modules/sibyll/Decay.inl

 /*
  * (c) Copyright 2018 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.
+ * 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
@@ -11,6 +10,7 @@
 #include <corsika/modules/sibyll/Decay.hpp>
 #include <corsika/modules/sibyll/ParticleConversion.hpp>
 #include <corsika/modules/sibyll/SibStack.hpp>
+#include <corsika/modules/Random.hpp>
 
 #include <iostream>
 #include <vector>
@@ -19,6 +19,7 @@ namespace corsika::sibyll {
 
   inline Decay::Decay(const bool sibyll_printout_on)
       : sibyll_listing_(sibyll_printout_on) {
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
     // switch off decays to avoid internal decay chains
     setAllStable();
     // handle all decays by default
@@ -31,7 +32,9 @@ namespace corsika::sibyll {
     setAllStable();
   }
 
-  inline Decay::~Decay() { CORSIKA_LOG_TRACE("Sibyll::Decay n={}", count_); }
+  inline Decay::~Decay() {
+    CORSIKA_LOGGER_TRACE(logger_, "Total number of Sibyll decays n={}", count_);
+  }
 
   inline bool Decay::canHandleDecay(const Code vParticleCode) {
     // if known to sibyll and not proton or neutrino it can decay
@@ -51,7 +54,7 @@ namespace corsika::sibyll {
 
   inline void Decay::setHandleDecay(const Code vParticleCode) {
     handleAllDecays_ = false;
-    CORSIKA_LOG_DEBUG("Sibyll::Decay: set to handle decay of {}", vParticleCode);
+    CORSIKA_LOGGER_DEBUG(logger_, "set to handle decay of {}", vParticleCode);
     if (Decay::canHandleDecay(vParticleCode))
       handledDecays_.insert(vParticleCode);
     else
@@ -60,7 +63,7 @@ namespace corsika::sibyll {
 
   inline void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
     handleAllDecays_ = false;
-    for (auto p : vParticleList) Decay::setHandleDecay(p);
+    for (auto const p : vParticleList) Decay::setHandleDecay(p);
   }
 
   inline bool Decay::isDecayHandled(corsika::Code const vParticleCode) {
@@ -73,11 +76,11 @@ namespace corsika::sibyll {
   }
 
   inline void Decay::setStable(std::vector<Code> const& vParticleList) {
-    for (auto p : vParticleList) Decay::setStable(p);
+    for (auto const p : vParticleList) Decay::setStable(p);
   }
 
   inline void Decay::setUnstable(std::vector<Code> const& vParticleList) {
-    for (auto p : vParticleList) Decay::setUnstable(p);
+    for (auto const p : vParticleList) Decay::setUnstable(p);
   }
 
   inline bool Decay::isStable(Code const vCode) {
@@ -93,14 +96,14 @@ namespace corsika::sibyll {
   }
 
   inline void Decay::setUnstable(Code const vCode) {
-    CORSIKA_LOG_DEBUG("Sibyll::Decay: setting {} unstable. ", vCode);
+    CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"unstable\". ", vCode);
 
     const int s_id = abs(sibyll::convertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
   }
 
   inline void Decay::setStable(Code const vCode) {
-    CORSIKA_LOG_DEBUG("Sibyll::Decay: setting {} stable. ", vCode);
+    CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"stable\". ", vCode);
 
     const int s_id = abs(sibyll::convertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
@@ -117,18 +120,18 @@ namespace corsika::sibyll {
   inline void Decay::printDecayConfig([[maybe_unused]] const Code vCode) {
     [[maybe_unused]] const int sibCode = corsika::sibyll::convertToSibyllRaw(vCode);
     [[maybe_unused]] const int absSibCode = abs(sibCode);
-    CORSIKA_LOG_DEBUG("Decay: Sibyll decay configuration: {} is {}", vCode,
-                      (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable");
+    CORSIKA_LOGGER_DEBUG(logger_, "decay configuration: {} is \"{}\"", vCode,
+                         (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable");
   }
   inline void Decay::printDecayConfig() {
-    CORSIKA_LOG_DEBUG("Sibyll::Decay: decay configuration:");
+    CORSIKA_LOGGER_DEBUG(logger_, "decay configuration:");
     if (handleAllDecays_) {
-      CORSIKA_LOG_DEBUG(
-          "     all particles known to Sibyll are handled by Sibyll::Decay!");
+      CORSIKA_LOGGER_DEBUG(
+          logger_, "     all particles known to Sibyll are handled by Sibyll::Decay!");
 
     } else {
-      for ([[maybe_unused]] auto& pCode : handledDecays_) {
-        CORSIKA_LOG_DEBUG("      Decay of {}  is handled by Sibyll!", pCode);
+      for ([[maybe_unused]] auto const pCode : handledDecays_) {
+        CORSIKA_LOGGER_DEBUG(logger_, "      Decay of {} is handled by Sibyll!", pCode);
       }
     }
   }
@@ -146,18 +149,18 @@ namespace corsika::sibyll {
       [[maybe_unused]] const auto mkin =
           +(E * E - projectile.getMomentum()
                         .getSquaredNorm()); // delta_mass(projectile.getMomentum(), E, m);
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: code: {} ", projectile.getPID());
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: MinStep: t0: {} ", t0);
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: MinStep: energy: {} GeV ", E / 1_GeV);
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: momentum: {} GeV ",
-                        projectile.getMomentum().getComponents() / 1_GeV);
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: momentum: shell mass-kin. inv. mass {} {}",
-                        mkin / 1_GeV / 1_GeV, m / 1_GeV * m / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "code: {} ", projectile.getPID());
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: t0: {} ", t0);
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: energy: {} GeV ", E / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "momentum: {} GeV ",
+                           projectile.getMomentum().getComponents() / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "momentum: shell mass-kin. inv. mass {} {}",
+                           mkin / 1_GeV / 1_GeV, m / 1_GeV * m / 1_GeV);
       [[maybe_unused]] auto sib_id =
           corsika::sibyll::convertToSibyllRaw(projectile.getPID());
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: sib mass: {}", get_sibyll_mass2(sib_id));
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: MinStep: gamma:  {}", gamma);
-      CORSIKA_LOG_DEBUG("Sibyll::Decay: MinStep: tau {} s: ", lifetime / 1_s);
+      CORSIKA_LOGGER_DEBUG(logger_, "sib mass: {}", get_sibyll_mass2(sib_id));
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: gamma:  {}", gamma);
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: tau {} s: ", lifetime / 1_s);
       return lifetime;
     }
     return std::numeric_limits<double>::infinity() * 1_s;
@@ -179,7 +182,7 @@ namespace corsika::sibyll {
     printDecayConfig(pCode);
 
     // call sibyll decay
-    CORSIKA_LOG_DEBUG("Decay: calling Sibyll decay routine..");
+    CORSIKA_LOGGER_DEBUG(logger_, "calling Sibyll decay routine..");
 
     // particle to pass to sibyll decay
     int inputSibPID = sibyll::convertToSibyllRaw(pCode);
@@ -198,11 +201,10 @@ namespace corsika::sibyll {
     int outputSibPID[10];
 
     // run decay routine
-    decpar_(inputSibPID, inputMomentum, nFinalParticles, outputSibPID,
-            &outputMomentum[0]);
+    decpar_sib_(inputSibPID, inputMomentum, nFinalParticles, outputSibPID,
+                &outputMomentum[0]);
 
-    CORSIKA_LOG_TRACE("Sibyll::Decay: number of final state particles: {}",
-                      nFinalParticles);
+    CORSIKA_LOGGER_TRACE(logger_, "number of final state particles: {}", nFinalParticles);
 
     // reset to stable
     setStable(pCode);
@@ -218,7 +220,8 @@ namespace corsika::sibyll {
       auto const pid = sibyll::convertFromSibyll(
           static_cast<corsika::sibyll::SibyllCode>(outputSibPID[i]));
 
-      CORSIKA_LOG_TRACE("Sibyll::Decay: i={} id={} p={} GeV", i, pid, components / 1_GeV);
+      CORSIKA_LOGGER_TRACE(logger_, "final particle i={} id={} p={} GeV", i, pid,
+                           components / 1_GeV);
 
       auto const p3 = MomentumVector(rootCS, components);
       HEPEnergyType const mass = get_mass(pid);

corsika/detail/modules/sibyll/HadronInteractionModel.inl

+/*
+ * (c) Copyright 2018 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/geometry/Point.hpp>
+
+#include <corsika/modules/sibyll/ParticleConversion.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/sibyll/SibStack.hpp>
+#include <corsika/modules/Random.hpp>
+
+#include <sibyll2.3d.hpp>
+
+#include <tuple>
+
+namespace corsika::sibyll {
+
+  inline void HadronInteractionModel::setVerbose(bool const flag) {
+    sibyll_listing_ = flag;
+  }
+
+  inline HadronInteractionModel::HadronInteractionModel()
+      : sibyll_listing_(false)
+      , internal_decays_(false)
+      , stable_particles_{} {
+    // initialize Sibyll
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+    static bool initialized = false;
+    if (!initialized) {
+      sibyll_ini_();
+      initialized = true;
+    }
+  }
+
+  inline HadronInteractionModel::HadronInteractionModel(std::set<Code> vlist)
+      : sibyll_listing_(false)
+      , internal_decays_(true)
+      , stable_particles_(vlist) {
+    // initialize Sibyll
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+    static bool initialized = false;
+    if (!initialized) {
+      sibyll_ini_();
+      initialized = true;
+    }
+  }
+
+  inline void HadronInteractionModel::setAllParticlesUnstable() {
+    // activate all decays in SIBYLL
+    CORSIKA_LOGGER_DEBUG(logger_, "setting all particles \"unstable\".");
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+  }
+
+  inline void HadronInteractionModel::setParticleListStable(std::set<Code> vList) {
+    // de-activate specific decays in SIBYLL
+    for (auto const p : vList) {
+      CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"stable\". ", p);
+      auto const sib_code = sibyll::convertToSibyll(p);
+      if (sib_code != corsika::sibyll::SibyllCode::Unknown) {
+        int const s_id = abs(static_cast<int>(sib_code));
+        s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+      } else {
+        CORSIKA_LOGGER_WARN(logger_,
+                            "particle {} not known to SIBYLL! Cannot set stable!", p);
+      }
+    }
+  }
+
+  inline HadronInteractionModel::~HadronInteractionModel() {
+    CORSIKA_LOGGER_DEBUG(logger_, "Sibyll::Model n={}, Nnuc={}", count_, nucCount_);
+  }
+
+  inline bool constexpr HadronInteractionModel::isValid(
+      Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn) const {
+    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
+
+    if (is_nucleus(targetId)) {
+      size_t const targA = get_nucleus_A(targetId);
+      if (targA != 1 && (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) {
+        return false;
+      }
+    } else if (targetId != Code::Proton && targetId != Code::Neutron &&
+               targetId != Code::Hydrogen) {
+      return false;
+    }
+    if (is_nucleus(projectileId) || !corsika::sibyll::canInteract(projectileId)) {
+      return false;
+    }
+    return true;
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  HadronInteractionModel::getCrossSectionInelEla(Code const projectileId,
+                                                 Code const targetId,
+                                                 FourMomentum const& projectileP4,
+                                                 FourMomentum const& targetP4) const {
+
+    int targetSibCode = 1; // nucleon or particle count
+    if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
+
+    if (!isValid(projectileId, targetId, sqrtSnn)) {
+      return {CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+
+    double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
+    int const iBeam = corsika::sibyll::getSibyllXSCode(
+        projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like,
+                       // 3:kaon-like)
+    if (!iBeam) {
+      CORSIKA_LOGGER_TRACE(
+          logger_, "The cross section for projectile {} and target {} is zero, E={}",
+          projectileId, targetId, sqrtSnn);
+      return {CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+
+    double const dEcm = sqrtSnn / 1_GeV;
+    // single nucleon target (p,n, hydrogen) or 4<=A<=18
+    double sigProd = 0;
+    double sigEla = 0;
+    if (targetId == Code::Proton || targetId == Code::Hydrogen ||
+        targetId == Code::Neutron) {
+      // single nucleon target
+      sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+    } else {
+      // nuclear target
+      int const iTarget = get_nucleus_A(targetId);
+      sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+    }
+    return {sigProd * 1_mb, sigEla * 1_mb};
+  }
+
+  /**
+   * In this function SIBYLL is called to produce one event. The
+   * event is copied (and boosted) into the shower lab frame.
+   */
+
+  template <typename TSecondaryView>
+  inline void HadronInteractionModel::doInteraction(TSecondaryView& secondaries,
+                                                    Code const projectileId,
+                                                    Code const targetId,
+                                                    FourMomentum const& projectileP4,
+                                                    FourMomentum const& targetP4) {
+
+    int targetSibCode = 1; // nucleon or particle count
+    if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
+    CORSIKA_LOGGER_DEBUG(logger_, "sibyll code: {} (nucleon/particle count)",
+                         targetSibCode);
+
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
+    COMBoost const boost(projectileP4, targetP4 / targetSibCode);
+
+    if (!isValid(projectileId, targetId, sqrtSnn)) {
+      throw std::runtime_error("Invalid target/projectile/energy combination");
+    }
+
+    CORSIKA_LOGGER_DEBUG(logger_, "pId={} tId={} sqrtSnn={}GeV", projectileId, targetId,
+                         sqrtSnn);
+
+    // beam id for sibyll
+    int const projectileSibyllCode = corsika::sibyll::convertToSibyllRaw(projectileId);
+
+    count_++;
+    // Sibyll does not know about units..
+    double const sqs = sqrtSnn / 1_GeV;
+    if (internal_decays_) {
+      setAllParticlesUnstable();
+      setParticleListStable(stable_particles_);
+    }
+    // running sibyll, filling stack
+    sibyll_(projectileSibyllCode, targetSibCode, sqs);
+    if (internal_decays_) decsib_();
+
+    if (sibyll_listing_) {
+      // print final state
+      int print_unit = 6;
+      sib_list_(print_unit);
+      nucCount_ += get_nwounded() - 1;
+    }
+
+    // ------ output and particle readout -----
+    auto const& csPrime = boost.getRotatedCS();
+
+    // add particles from sibyll to stack
+
+    // link to sibyll stack
+    SibStack ss;
+
+    auto const& originalCS = boost.getOriginalCS();
+    MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
+    for (auto const& psib : ss) {
+      // abort on particles that have decayed in Sibyll. Should not happen!
+      if (psib.hasDecayed()) { // LCOV_EXCL_START
+        if (internal_decays_) {
+          continue;
+        } else {
+          throw std::runtime_error(
+              "internal decayse switched off! found particle that decayed in SIBYLL!");
+        } // LCOV_EXCL_STOP
+      }
+      // transform 4-momentum to lab. frame
+      // note that the momentum needs to be rotated back
+      auto const tmp = psib.getMomentum().getComponents();
+      auto const pCoM = MomentumVector(csPrime, tmp);
+      HEPEnergyType const eCoM = psib.getEnergy();
+      auto const P4lab = boost.fromCoM(FourVector{eCoM, pCoM});
+      auto const p3lab = P4lab.getSpaceLikeComponents();
+
+      Code const pid = corsika::sibyll::convertFromSibyll(psib.getPID());
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
+      // add to corsika stack
+      auto pnew =
+          secondaries.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
+
+      Plab_final += pnew.getMomentum();
+      Elab_final += pnew.getEnergy();
+      Ecm_final += psib.getEnergy();
+    }
+    { // just output
+      [[maybe_unused]] HEPEnergyType const Elab_initial =
+          static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
+      CORSIKA_LOGGER_DEBUG(logger_,
+                           "conservation (all GeV): "
+                           "sqrtSnn={}, sqrtSnn_final={}, "
+                           "Elab_initial={}, Elab_final={}, "
+                           "diff(%)={}, "
+                           "E in nucleons={}, "
+                           "Plab_final={} ",
+                           sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV,
+                           Elab_initial, Elab_final / 1_GeV,
+                           (Elab_final - Elab_initial) / Elab_initial * 100,
+                           constants::nucleonMass * get_nwounded() / 1_GeV,
+                           (Plab_final / 1_GeV).getComponents());
+    }
+  }
+} // namespace corsika::sibyll

corsika/detail/modules/sibyll/InteractionModel.inl

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2022 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.
+ * 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/geometry/Point.hpp>
-
-#include <corsika/modules/sibyll/ParticleConversion.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/modules/sibyll/SibStack.hpp>
-
-#include <sibyll2.3d.hpp>
-
-#include <tuple>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/modules/sibyll/HadronInteractionModel.hpp>
+#include <corsika/modules/sibyll/NuclearInteractionModel.hpp>
 
 namespace corsika::sibyll {
 
-  inline void InteractionModel::setVerbose(bool const flag) { sibyll_listing_ = flag; }
+  inline InteractionModel::InteractionModel(std::set<Code> const& nuccomp,
+                                            std::set<Code> const& list)
+      : hadronSibyll_{list}
+      , nuclearSibyll_{hadronSibyll_, nuccomp} {}
 
-  inline InteractionModel::InteractionModel()
-      : sibyll_listing_(false) {
-    // initialize Sibyll
-    static bool initialized = false;
-    if (!initialized) {
-      sibyll_ini_();
-      initialized = true;
-    }
+  inline HadronInteractionModel& InteractionModel::getHadronInteractionModel() {
+    return hadronSibyll_;
   }
 
-  inline InteractionModel::~InteractionModel() {
-    CORSIKA_LOG_DEBUG("Sibyll::Model n={}, Nnuc={}", count_, nucCount_);
+  inline HadronInteractionModel const& InteractionModel::getHadronInteractionModel()
+      const {
+    return hadronSibyll_;
   }
 
-  inline bool constexpr InteractionModel::isValid(Code const projectileId,
-                                                  Code const targetId,
-                                                  HEPEnergyType const sqrtSnn) const {
-    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
-
-    if (is_nucleus(targetId)) {
-      size_t const targA = get_nucleus_A(targetId);
-      if (targA != 1 && (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) {
-        return false;
-      }
-    } else if (targetId != Code::Proton && targetId != Code::Neutron &&
-               targetId != Code::Hydrogen) {
-      return false;
-    }
-    if (is_nucleus(projectileId) || !corsika::sibyll::canInteract(projectileId)) {
-      return false;
-    }
-    return true;
+  inline typename InteractionModel::nuclear_model_type&
+  InteractionModel::getNuclearInteractionModel() {
+    return nuclearSibyll_;
   }
 
-  inline std::tuple<CrossSectionType, CrossSectionType>
-  InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
-                                           FourMomentum const& projectileP4,
-                                           FourMomentum const& targetP4) const {
-
-    int targetSibCode = 1; // nucleon or particle count
-    if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
-    // sqrtS per target nucleon
-    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
-
-    if (!isValid(projectileId, targetId, sqrtSnn)) {
-      return {CrossSectionType::zero(), CrossSectionType::zero()};
-    }
-
-    double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
-    int const iBeam = corsika::sibyll::getSibyllXSCode(
-        projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like,
-                       // 3:kaon-like)
-
-    double const dEcm = sqrtSnn / 1_GeV;
-    // single nucleon target (p,n, hydrogen) or 4<=A<=18
-    double sigProd = 0;
-    double sigEla = 0;
-    if (targetId == Code::Proton || targetId == Code::Hydrogen ||
-        targetId == Code::Neutron) {
-      // single nucleon target
-      sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
-    } else {
-      // nuclear target
-      int const iTarget = get_nucleus_A(targetId);
-      sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
-    }
-    return {sigProd * 1_mb, sigEla * 1_mb};
-  } // namespace corsika::sibyll
-
-  /**
-   * In this function SIBYLL is called to produce one event. The
-   * event is copied (and boosted) into the shower lab frame.
-   */
-
-  template <typename TSecondaryView>
-  inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
-                                              Code const projectileId,
-                                              Code const targetId,
-                                              FourMomentum const& projectileP4,
-                                              FourMomentum const& targetP4) {
-
-    int targetSibCode = 1; // nucleon or particle count
-    if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
-    CORSIKA_LOG_DEBUG("sibyll code: {} (nucleon/particle count)", targetSibCode);
-
-    // sqrtS per target nucleon
-    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
-    COMBoost const boost(projectileP4, targetP4 / targetSibCode);
-
-    if (!isValid(projectileId, targetId, sqrtSnn)) {
-      throw std::runtime_error("Invalid target/projectile/energy combination");
-    }
-
-    CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnn={}GeV", projectileId, targetId, sqrtSnn);
-
-    // beam id for sibyll
-    int const projectileSibyllCode = corsika::sibyll::convertToSibyllRaw(projectileId);
-
-    count_++;
-    // Sibyll does not know about units..
-    double const sqs = sqrtSnn / 1_GeV;
-    // running sibyll, filling stack
-    sibyll_(projectileSibyllCode, targetSibCode, sqs);
-
-    if (sibyll_listing_) {
-      // print final state
-      int print_unit = 6;
-      sib_list_(print_unit);
-      nucCount_ += get_nwounded() - 1;
-    }
-
-    // ------ output and particle readout -----
-    auto const& csPrime = boost.getRotatedCS();
-
-    // add particles from sibyll to stack
-
-    // link to sibyll stack
-    SibStack ss;
-
-    auto const& originalCS = boost.getOriginalCS();
-    MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-    HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
-    for (auto& psib : ss) {
-      // abort on particles that have decayed in Sibyll. Should not happen!
-      if (psib.hasDecayed()) { // LCOV_EXCL_START
-        throw std::runtime_error("found particle that decayed in SIBYLL!");
-      } // LCOV_EXCL_STOP
-
-      // transform 4-momentum to lab. frame
-      // note that the momentum needs to be rotated back
-      auto const tmp = psib.getMomentum().getComponents();
-      auto const pCoM = MomentumVector(csPrime, tmp);
-      HEPEnergyType const eCoM = psib.getEnergy();
-      auto const P4lab = boost.fromCoM(FourVector{eCoM, pCoM});
-      auto const p3lab = P4lab.getSpaceLikeComponents();
+  inline typename InteractionModel::nuclear_model_type const&
+  InteractionModel::getNuclearInteractionModel() const {
+    return nuclearSibyll_;
+  }
 
-      Code const pid = corsika::sibyll::convertFromSibyll(psib.getPID());
-      HEPEnergyType const mass = get_mass(pid);
-      HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code projCode, Code targetCode, FourMomentum const& proj4mom,
+      FourMomentum const& target4mom) const {
+    if (is_nucleus(projCode))
+      return getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                          target4mom);
+    else
+      return getHadronInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                         target4mom);
+  }
 
-      // add to corsika stack
-      auto pnew =
-          secondaries.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code projCode, Code targetCode,
+                                           FourMomentum const& proj4mom,
+                                           FourMomentum const& target4mom) const {
+    if (is_nucleus(projCode))
+      return {getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                           target4mom),
+              CrossSectionType::zero()};
+    else
+      return getHadronInteractionModel().getCrossSectionInelEla(projCode, targetCode,
+                                                                proj4mom, target4mom);
+  }
 
-      Plab_final += pnew.getMomentum();
-      Elab_final += pnew.getEnergy();
-      Ecm_final += psib.getEnergy();
-    }
-    { // just output
-      HEPEnergyType const Elab_initial =
-          static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
-      CORSIKA_LOG_DEBUG(
-          "conservation (all GeV): "
-          "sqrtSnn={}, sqrtSnn_final={}, "
-          "Elab_initial={}, Elab_final={}, "
-          "diff(%)={}, "
-          "E in nucleons={}, "
-          "Plab_final={} ",
-          sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Elab_initial,
-          Elab_final / 1_GeV, (Elab_final - Elab_initial) / Elab_initial * 100,
-          constants::nucleonMass * get_nwounded() / 1_GeV,
-          (Plab_final / 1_GeV).getComponents());
-    }
+  template <typename TSecondaries>
+  inline void InteractionModel::doInteraction(TSecondaries& view, Code projCode,
+                                              Code targetCode,
+                                              FourMomentum const& proj4mom,
+                                              FourMomentum const& target4mom) {
+    if (is_nucleus(projCode))
+      return getNuclearInteractionModel().doInteraction(view, projCode, targetCode,
+                                                        proj4mom, target4mom);
+    else
+      return getHadronInteractionModel().doInteraction(view, projCode, targetCode,
+                                                       proj4mom, target4mom);
   }
 } // namespace corsika::sibyll

corsika/detail/modules/sibyll/NuclearInteractionModel.inl

 /*
  * (c) Copyright 2018 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.
+ * 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/media/Environment.hpp>
 #include <corsika/media/NuclearComposition.hpp>
-
+#include <corsika/modules/Random.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
 #include <corsika/framework/core/Logging.hpp>
 
@@ -19,28 +19,27 @@
 
 namespace corsika::sibyll {
 
-  template <typename TEnvironment, typename TNucleonModel>
-  inline NuclearInteractionModel<TEnvironment, TNucleonModel>::NuclearInteractionModel(
-      TNucleonModel& hadint, TEnvironment const& env)
-      : environment_(env)
-      , hadronicInteraction_(hadint) {
+  template <typename TNucleonModel>
+  inline NuclearInteractionModel<TNucleonModel>::NuclearInteractionModel(
+      TNucleonModel& hadint, std::set<Code> const& nuccomp)
+      : hadronicInteraction_(hadint) {
 
     // initialize nuclib
     // TODO: make sure this does not overlap with sibyll
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
     nuc_nuc_ini_();
 
     // initialize cross sections
-    initializeNuclearCrossSections();
+    initializeNuclearCrossSections(nuccomp);
   }
 
-  template <typename TEnvironment, typename TNucleonModel>
-  inline NuclearInteractionModel<TEnvironment,
-                                 TNucleonModel>::~NuclearInteractionModel() {
-    CORSIKA_LOG_DEBUG("Nuclib::NuclearInteractionModel n={} Nnuc={}", count_, nucCount_);
+  template <typename TNucleonModel>
+  inline NuclearInteractionModel<TNucleonModel>::~NuclearInteractionModel() {
+    CORSIKA_LOGGER_DEBUG(logger_, "nuclear interactions handled by Sibyll n={}", count_);
   }
 
-  template <typename TEnvironment, typename TNucleonModel>
-  inline bool constexpr NuclearInteractionModel<TEnvironment, TNucleonModel>::isValid(
+  template <typename TNucleonModel>
+  inline bool constexpr NuclearInteractionModel<TNucleonModel>::isValid(
       Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn) const {
 
     // also depends on underlying model, for Proton/Neutron projectile
@@ -53,12 +52,12 @@ namespace corsika::sibyll {
     return true;
   } // namespace corsika::sibyll
 
-  template <typename TEnvironment, typename TNucleonModel>
-  inline void
-  NuclearInteractionModel<TEnvironment, TNucleonModel>::printCrossSectionTable(
+  template <typename TNucleonModel>
+  inline void NuclearInteractionModel<TNucleonModel>::printCrossSectionTable(
       Code const pCode) const {
     if (!hadronicInteraction_.isValid(Code::Proton, pCode, 100_GeV)) { // LCOV_EXCL_START
-      CORSIKA_LOG_ERROR("Invalid target type {} for hadron interaction model.", pCode);
+      CORSIKA_LOGGER_WARN(logger_, "Invalid target type {} for hadron interaction model.",
+                          pCode);
       return;
     } // LCOV_EXCL_STOP
 
@@ -82,39 +81,25 @@ namespace corsika::sibyll {
       }
       table << "\n";
     }
-    CORSIKA_LOG_DEBUG(table.str());
+    CORSIKA_LOGGER_DEBUG(logger_, table.str());
   }
 
-  template <typename TEnvironment, typename TNucleonModel>
-  inline void
-  NuclearInteractionModel<TEnvironment, TNucleonModel>::initializeNuclearCrossSections() {
-
-    auto& universe = *(environment_.getUniverse());
-    // generate complete list of all nuclei types in universe
-
-    auto const allElementsInUniverse = std::invoke([&]() {
-      std::set<Code> allElementsInUniverse;
-      auto collectElements = [&](auto& vtn) {
-        if (vtn.hasModelProperties()) {
-          auto const& comp =
-              vtn.getModelProperties().getNuclearComposition().getComponents();
-          for (auto const c : comp) allElementsInUniverse.insert(c);
-        }
-      };
-      universe.walk(collectElements);
-      return allElementsInUniverse;
-    });
+  template <typename TNucleonModel>
+  inline void NuclearInteractionModel<TNucleonModel>::initializeNuclearCrossSections(
+      std::set<Code> const& allElementsInUniverse) {
 
-    CORSIKA_LOG_DEBUG("initializing nuclear cross sections...");
+    CORSIKA_LOGGER_DEBUG(logger_, "initializing nuclear cross sections...");
 
     // loop over target components, at most 4!!
     int k = -1;
     for (Code const ptarg : allElementsInUniverse) {
       ++k;
-      CORSIKA_LOG_DEBUG("init target component: {} A={}", ptarg, get_nucleus_A(ptarg));
+      CORSIKA_LOGGER_DEBUG(logger_, "init target component: {} A={}", ptarg,
+                           get_nucleus_A(ptarg));
       int const ib = get_nucleus_A(ptarg);
       if (!hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV)) {
-        CORSIKA_LOG_ERROR("Invalid target type {} for hadron interaction model.", ptarg);
+        CORSIKA_LOGGER_WARN(
+            logger_, "Invalid target type {} for hadron interaction model.", ptarg);
         continue;
       }
       targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
@@ -132,14 +117,15 @@ namespace corsika::sibyll {
         FourMomentum targetP4(EcmHalve, {cs, -pcm, 0_eV, 0_eV});
         // get p-p cross sections
         if (!hadronicInteraction_.isValid(Code::Proton, Code::Proton, Ecm)) {
-          throw std::runtime_error("invalid projectile,target,ecm combination");
+          throw std::runtime_error("invalid (projectile,target,ecm) combination");
         }
         auto const [siginel, sigela] = hadronicInteraction_.getCrossSectionInelEla(
             Code::Proton, Code::Proton, projectileP4, targetP4);
         double const dsig = siginel / 1_mb;
         double const dsigela = sigela / 1_mb;
         // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
-        CORSIKA_LOG_TRACE("Ecm={} siginel={} sigela={}", Ecm / 1_GeV, dsig, dsigela);
+        CORSIKA_LOGGER_TRACE(logger_, "Ecm={} siginel={} sigela={}", Ecm / 1_GeV, dsig,
+                             dsigela);
         for (size_t j = 1; j < gMaxNucleusAProjectile_; ++j) {
           const int jj = j + 1;
           double sig_out, dsig_out, sigqe_out, dsigqe_out;
@@ -148,53 +134,76 @@ namespace corsika::sibyll {
           // write to table
           cnucsignuc_.sigma[j][k][i] = sig_out;
           cnucsignuc_.sigqe[j][k][i] = sigqe_out;
-          CORSIKA_LOG_TRACE("nuc A={} sig={} qe={}", j, sig_out, sigqe_out);
+          CORSIKA_LOGGER_TRACE(logger_, "nuc A={} sig={} qe={}", j, sig_out, sigqe_out);
         }
       }
     }
-    CORSIKA_LOG_DEBUG("cross sections for {} components initialized!",
-                      targetComponentsIndex_.size());
+    CORSIKA_LOGGER_DEBUG(logger_, "cross sections for {} components initialized!",
+                         targetComponentsIndex_.size());
     for (auto& ptarg : allElementsInUniverse) { printCrossSectionTable(ptarg); }
   }
 
-  template <typename TEnvironment, typename TNucleonModel>
-  inline CrossSectionType
-  NuclearInteractionModel<TEnvironment, TNucleonModel>::readCrossSectionTable(
+  template <typename TNucleonModel>
+  inline CrossSectionType NuclearInteractionModel<TNucleonModel>::readCrossSectionTable(
       int const ia, Code const pTarget, HEPEnergyType const elabnuc) const {
-
+    CORSIKA_LOGGER_DEBUG(logger_, "ia={}, target={}, ElabNuc={} GeV", ia, pTarget,
+                         elabnuc / 1_GeV);
     int const ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
     auto const ECoMNuc = sqrt(2. * constants::nucleonMass * elabnuc);
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtSnn= {} GeV", ECoMNuc / 1_GeV);
     if (ECoMNuc < getMinEnergyPerNucleonCoM() || ECoMNuc > getMaxEnergyPerNucleonCoM()) {
-      throw std::runtime_error("energy outside tabulated range!");
+      CORSIKA_LOGGER_WARN(
+          logger_,
+          "nucleon-nucleon energy outside range! sqrtSnn={}GeV (limits: {} .. {} GeV)",
+          ECoMNuc / 1_GeV, getMinEnergyPerNucleonCoM() / 1_GeV,
+          getMaxEnergyPerNucleonCoM() / 1_GeV);
+      // throw std::runtime_error("energy outside tabulated range!");
     }
     double const e0 = elabnuc / 1_GeV;
     double sig;
-    CORSIKA_LOG_DEBUG("ReadCrossSectionTable: {} {} {}", ia, ib, e0);
+    CORSIKA_LOGGER_DEBUG(logger_, "ReadCrossSectionTable: {} {} {}", ia, ib, e0);
     signuc2_(ia, ib, e0, sig);
-    CORSIKA_LOG_DEBUG("ReadCrossSectionTable: sig={}", sig);
+    CORSIKA_LOGGER_DEBUG(logger_, "ReadCrossSectionTable: sig={}", sig);
     return sig * 1_mb;
   }
 
-  template <typename TEnvironment, typename TNucleonModel>
-  CrossSectionType inline NuclearInteractionModel<
-      TEnvironment, TNucleonModel>::getCrossSection(Code const projectileId,
-                                                    Code const targetId,
-                                                    FourMomentum const& projectileP4,
-                                                    FourMomentum const& targetP4) const {
+  template <typename TNucleonModel>
+  CrossSectionType inline NuclearInteractionModel<TNucleonModel>::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+    CORSIKA_LOGGER_DEBUG(logger_, "projectile: E={}, p3={} \n target: E={}, p3={}",
+                         projectileP4.getTimeLikeComponent() / 1_GeV,
+                         projectileP4.getSpaceLikeComponents() / 1_GeV,
+                         targetP4.getTimeLikeComponent() / 1_GeV,
+                         targetP4.getSpaceLikeComponents() / 1_GeV);
+    // check if projectile and target are nuclei!
+    if (!is_nucleus(projectileId) || !is_nucleus(targetId)) {
+      return CrossSectionType::zero();
+    }
+
+    // calculate sqrt(Snn) (only works if projectile and target are nuclei)
+    HEPEnergyType const sqrtSnn =
+        (projectileP4 / get_nucleus_A(projectileId) + targetP4 / get_nucleus_A(targetId))
+            .getNorm();
 
-    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4).getNorm();
+    CORSIKA_LOG_DEBUG("proj={}, targ={}, sqrtSNN={}GeV", projectileId, targetId,
+                      sqrtSnn / 1_GeV);
     if (!isValid(projectileId, targetId, sqrtSnn)) { return CrossSectionType::zero(); }
-    HEPEnergyType const LabEnergyPerNuc =
-        static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
+
+    // lab-frame energy per projectile nucleon as required by signuc2()
+    HEPEnergyType const LabEnergyPerNuc = calculate_lab_energy(
+        static_pow<2>(sqrtSnn), get_mass(projectileId) / get_nucleus_A(projectileId),
+        get_mass(targetId) / get_nucleus_A(targetId));
     auto const sigProd =
         readCrossSectionTable(get_nucleus_A(projectileId), targetId, LabEnergyPerNuc);
-    CORSIKA_LOG_DEBUG("cross section (mb): {}", sigProd / 1_mb);
+    CORSIKA_LOGGER_DEBUG(logger_, "cross section (mb): sqrtSnn={} sig={}",
+                         sqrtSnn / 1_GeV, sigProd / 1_mb);
     return sigProd;
   }
 
-  template <typename TEnvironment, typename TNucleonModel>
+  template <typename TNucleonModel>
   template <typename TSecondaryView>
-  inline void NuclearInteractionModel<TEnvironment, TNucleonModel>::doInteraction(
+  inline void NuclearInteractionModel<TNucleonModel>::doInteraction(
       TSecondaryView& view, Code const projectileId, Code const targetId,
       FourMomentum const& projectileP4, FourMomentum const& targetP4) {
 
@@ -215,8 +224,8 @@ namespace corsika::sibyll {
     // Elab corresponding to sqrtSnucleon -> fixed target projectile
     COMBoost const boost(nucleonP4, targetP4);
 
-    CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnucleon={}GeV Aproj={}", projectileId, targetId,
-                      sqrtSnucleon / 1_GeV, projectileA);
+    CORSIKA_LOGGER_DEBUG(logger_, "pId={} tId={} sqrtSnucleon={}GeV Aproj={}",
+                         projectileId, targetId, sqrtSnucleon / 1_GeV, projectileA);
     count_++;
 
     // lab. momentum per projectile nucleon
@@ -238,12 +247,12 @@ namespace corsika::sibyll {
                targetId == Code::Hydrogen) {
       kATarget = 1;
     }
-    CORSIKA_LOG_DEBUG("nuclib target code: {}", kATarget);
+    CORSIKA_LOGGER_DEBUG(logger_, "nuclib target code: {}", kATarget);
 
     // end of target sampling
 
     // superposition
-    CORSIKA_LOG_DEBUG("sampling nuc. multiple interaction structure.. ");
+    CORSIKA_LOGGER_DEBUG(logger_, "sampling nuc. multiple interaction structure.. ");
     // get nucleon-nucleon cross section
     // (needed to determine number of nucleon-nucleon scatterings)
     auto const protonId = Code::Proton;
@@ -257,20 +266,20 @@ namespace corsika::sibyll {
     // nuclear multiple scattering according to glauber (r.i.p.)
     int_nuc_(kATarget, projectileA, sigProd, sigEla);
 
-    CORSIKA_LOG_DEBUG(
-        "number of nucleons in target           : {}\n"
-        "number of wounded nucleons in target   : {}\n"
-        "number of nucleons in projectile       : {}\n"
-        "number of wounded nucleons in project. : {}\n"
-        "number of inel. nuc.-nuc. interactions : {}\n"
-        "number of elastic nucleons in target   : {}\n"
-        "number of elastic nucleons in project. : {}\n"
-        "impact parameter: {}",
-        kATarget, cnucms_.na, projectileA, cnucms_.nb, cnucms_.ni, cnucms_.nael,
-        cnucms_.nbel, cnucms_.b);
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "number of nucleons in target           : {}\n"
+                         "number of wounded nucleons in target   : {}\n"
+                         "number of nucleons in projectile       : {}\n"
+                         "number of wounded nucleons in project. : {}\n"
+                         "number of inel. nuc.-nuc. interactions : {}\n"
+                         "number of elastic nucleons in target   : {}\n"
+                         "number of elastic nucleons in project. : {}\n"
+                         "impact parameter: {}",
+                         kATarget, cnucms_.na, projectileA, cnucms_.nb, cnucms_.ni,
+                         cnucms_.nael, cnucms_.nbel, cnucms_.b);
 
     // calculate fragmentation
-    CORSIKA_LOG_DEBUG("calculating nuclear fragments..");
+    CORSIKA_LOGGER_DEBUG(logger_, "calculating nuclear fragments..");
     // number of interactions
     // include elastic
     int const nElasticNucleons = cnucms_.nbel;
@@ -292,12 +301,13 @@ namespace corsika::sibyll {
     }
     // (LCOV_EXCL_STOP)
 
-    CORSIKA_LOG_DEBUG("number of fragments: {}", nFragments);
-    CORSIKA_LOG_DEBUG("adding nuclear fragments to particle stack..");
+    CORSIKA_LOGGER_DEBUG(logger_, "number of fragments: {}", nFragments);
+    CORSIKA_LOGGER_DEBUG(logger_, "adding nuclear fragments to particle stack..");
     // put nuclear fragments on corsika stack
     for (int j = 0; j < nFragments; ++j) {
-      CORSIKA_LOG_DEBUG("fragment {}: A={} px={} py={} pz={}", j, AFragments[j],
-                        fragments_.ppp[j][0], fragments_.ppp[j][1], fragments_.ppp[j][2]);
+      CORSIKA_LOGGER_DEBUG(logger_, "fragment {}: A={} px={} py={} pz={}", j,
+                           AFragments[j], fragments_.ppp[j][0], fragments_.ppp[j][1],
+                           fragments_.ppp[j][2]);
       auto const nuclA = AFragments[j];
       // get Z from stability line
       auto const nuclZ = int(nuclA / 2.15 + 0.7);
@@ -309,9 +319,9 @@ namespace corsika::sibyll {
                                         : get_nucleus_code(nuclA, nuclZ));
       HEPMassType const mass = get_mass(specCode);
 
-      CORSIKA_LOG_DEBUG("adding fragment: {}", get_name(specCode));
-      CORSIKA_LOG_DEBUG("A,Z: {}, {}", nuclA, nuclZ);
-      CORSIKA_LOG_DEBUG("mass: {} GeV", mass / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "adding fragment: {}", get_name(specCode));
+      CORSIKA_LOGGER_DEBUG(logger_, "A,Z: {}, {}", nuclA, nuclZ);
+      CORSIKA_LOGGER_DEBUG(logger_, "mass: {} GeV", mass / 1_GeV);
 
       // CORSIKA 7 way
       // spectators inherit momentum from original projectile
@@ -319,14 +329,16 @@ namespace corsika::sibyll {
 
       HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
 
-      CORSIKA_LOG_DEBUG("fragment momentum {}", p3lab.getComponents() / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "fragment momentum {}",
+                           p3lab.getComponents() / 1_GeV);
       view.addSecondary(std::make_tuple(specCode, Ekin, p3lab.normalized()));
     }
 
     // add elastic nucleons to corsika stack
     // TODO: the elastic interaction could be external like the inelastic interaction,
     // e.g. use existing ElasticModel
-    CORSIKA_LOG_DEBUG("adding elastically scattered nucleons to particle stack..");
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "adding elastically scattered nucleons to particle stack..");
     for (int j = 0; j < nElasticNucleons; ++j) {
       // TODO: sample proton or neutron
       Code const elaNucCode = Code::Proton;
@@ -343,7 +355,7 @@ namespace corsika::sibyll {
     }
 
     // add inelastic interactions
-    CORSIKA_LOG_DEBUG("calculate inelastic nucleon-nucleon interactions..");
+    CORSIKA_LOGGER_DEBUG(logger_, "calculate inelastic nucleon-nucleon interactions..");
     for (int j = 0; j < nInelNucleons; ++j) {
       // TODO: sample neutron or proton
       auto const pCode = Code::Proton;
@@ -351,7 +363,7 @@ namespace corsika::sibyll {
       HEPEnergyType const Ekin = sqrt(p3NucleonLab.getSquaredNorm() + mass * mass) - mass;
 
       // temporarily add to stack, will be removed after interaction in DoInteraction
-      CORSIKA_LOG_DEBUG("inelastic interaction no. {}", j);
+      CORSIKA_LOGGER_DEBUG(logger_, "inelastic interaction no. {}", j);
       typename TSecondaryView::inner_stack_value_type nucleonStack;
       Point const pDummy(boost.getOriginalCS(), {0_m, 0_m, 0_m});
       TimeType const tDummy = 0_ns;
@@ -360,7 +372,7 @@ namespace corsika::sibyll {
       inelasticNucleon.setNode(view.getProjectile().getNode());
 
       // create inelastic interaction for each nucleon
-      CORSIKA_LOG_TRACE("calling HadronicInteraction...");
+      CORSIKA_LOGGER_TRACE(logger_, "calling HadronicInteraction...");
       // create new StackView for each of the nucleons
       TSecondaryView nucleon_secondaries(inelasticNucleon);
       // all inner hadronic event generator

corsika/detail/modules/sibyll/ParticleConversion.inl

 /*
  * (c) Copyright 2018 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.
+ * 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

corsika/detail/modules/sophia/InteractionModel.inl

+/*
+ * (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/framework/geometry/Point.hpp>
+#include <corsika/modules/sophia/ParticleConversion.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/sophia/SophiaStack.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline void InteractionModel::setVerbose(bool const flag) { sophia_listing_ = flag; }
+
+  inline InteractionModel::InteractionModel()
+      : sophia_listing_(false) {
+    corsika::connect_random_stream(RNG_, ::sophia::set_rng_function);
+    // set all particles stable in SOPHIA
+    for (int i = 0; i < 49; ++i) so_csydec_.idb[i] = -abs(so_csydec_.idb[i]);
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOG_DEBUG("Sophia::Model n={}", count_);
+  }
+
+  inline bool constexpr InteractionModel::isValid(Code const projectileId,
+                                                  Code const targetId,
+                                                  HEPEnergyType const sqrtSnn) const {
+    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
+
+    if (!(targetId == Code::Proton || targetId == Code::Neutron ||
+          targetId == Code::Hydrogen))
+      return false;
+
+    if (projectileId != Code::Photon) return false;
+
+    return true;
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "projectile: Id={}, E={} GeV, p3={} GeV", projectileId,
+                         projectileP4.getTimeLikeComponent() / 1_GeV,
+                         projectileP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+    CORSIKA_LOGGER_DEBUG(logger_, "target: Id={}, E={} GeV, p3={} GeV", targetId,
+                         targetP4.getTimeLikeComponent() / 1_GeV,
+                         targetP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={}GeV", sqrtS / 1_GeV);
+
+    // accepts only photon-nucleon interactions
+    if (!isValid(projectileId, targetId, sqrtS)) {
+      CORSIKA_LOGGER_ERROR(logger_,
+                           "Invalid target/projectile/energy combination: {},{},{} GeV",
+                           projectileId, targetId, sqrtS / 1_GeV);
+      throw std::runtime_error("SOPHIA: Invalid target/projectile/energy combination");
+    }
+
+    COMBoost const boost(projectileP4, targetP4);
+
+    int nucleonSophiaCode = convertToSophiaRaw(targetId); // either proton or neutron
+    // initialize resonance spectrum
+    initial_(nucleonSophiaCode);
+    // Sophia does sqrt(1 - mass_sophia / mass_c8), so we need to make sure that E0 >=
+    // m_sophia
+    double Enucleon = std::max(corsika::sophia::getSophiaMass(targetId) / 1_GeV,
+                               targetP4.getTimeLikeComponent() / 1_GeV);
+    double Ephoton = projectileP4.getTimeLikeComponent() / 1_GeV;
+    double theta = 0.0; // set nucleon at rest in collision
+    int Imode = -1;     // overwritten inside SOPHIA
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
+                         nucleonSophiaCode, Enucleon, Ephoton, theta);
+    count_++;
+    // call sophia
+    eventgen_(nucleonSophiaCode, Enucleon, Ephoton, theta, Imode);
+
+    if (sophia_listing_) {
+      int arg = 3;
+      print_event_(arg);
+    }
+
+    auto const& originalCS = boost.getOriginalCS();
+    // SOPHIA has photon along -z  and nucleon along +z (GZK calc..)
+    COMBoost const boostInternal(targetP4, projectileP4);
+    auto const& csPrime = boost.getRotatedCS();
+    CoordinateSystemPtr csPrimePrime =
+        make_rotation(csPrime, QuantityVector<length_d>{1_m, 0_m, 0_m}, M_PI);
+
+    SophiaStack ss;
+
+    MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType E_final = 0_GeV;
+    for (auto& psop : ss) {
+      // abort on particles that have decayed in SOPHIA. Should not happen!
+      if (psop.hasDecayed()) { // LCOV_EXCL_START
+        throw std::runtime_error("found particle that decayed in SOPHIA!");
+      } // LCOV_EXCL_STOP
+
+      auto momentumSophia = psop.getMomentum(csPrimePrime);
+      momentumSophia.rebase(csPrime);
+      auto const energySophia = psop.getEnergy();
+      auto const P4com = boostInternal.toCoM(FourVector{energySophia, momentumSophia});
+      auto const P4lab = boost.fromCoM(P4com);
+      SophiaCode const pidSophia = psop.getPID();
+      Code const pid = convertFromSophia(pidSophia);
+      auto momentum = P4lab.getSpaceLikeComponents();
+      momentum.rebase(originalCS);
+      HEPEnergyType const Ekin =
+          calculate_kinetic_energy(momentum.getNorm(), get_mass(pid));
+
+      CORSIKA_LOGGER_TRACE(logger_, "SOPHIA: pid={}, p={} GeV", pidSophia,
+                           momentumSophia.getComponents() / 1_GeV);
+
+      CORSIKA_LOGGER_TRACE(logger_, "CORSIKA: pid={}, p={} GeV", pid,
+                           momentum.getComponents() / 1_GeV);
+
+      auto pnew =
+          secondaries.addSecondary(std::make_tuple(pid, Ekin, momentum.normalized()));
+
+      P_final += pnew.getMomentum();
+      E_final += pnew.getEnergy();
+    }
+    CORSIKA_LOGGER_TRACE(logger_, "Efinal={} GeV,Pfinal={} GeV", E_final / 1_GeV,
+                         P_final.getComponents() / 1_GeV);
+  }
+
+} // namespace corsika::sophia

corsika/detail/modules/sophia/ParticleConversion.inl

+/*
+ * (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/framework/core/ParticleProperties.hpp>
+
+#include <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline HEPMassType getSophiaMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
+    auto sCode = convertToSophiaRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getSophiaMass: unknown particle!");
+    else
+      return sqrt(get_sophia_mass2(sCode)) * 1_GeV;
+  }
+} // namespace corsika::sophia
\ No newline at end of file

corsika/detail/modules/tauola/Decay.inl

+/*
+ * (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 <Tauola/Tauola.h>
+
+#include <corsika/framework/utility/COMBoost.hpp>
+
+namespace Tpp = Tauolapp;
+
+namespace corsika::tauola {
+
+  // default constructor
+  Decay::Decay(const Helicity helicity, const bool radiative, const bool spincorr)
+      : helicity_(helicity) {
+
+    // if TAUOLA has not already been initialized,
+    if (!Tpp::Tauola::getIsTauolaIni()) {
+
+      // set the default units explicitly.
+      // we use GeV from momentum and cm for length
+      Tpp::Tauola::setUnits(Tpp::Tauola::GEV, Tpp::Tauola::CM);
+
+      // we set the lifetime as zero so that the tau decays
+      // at the vertex point given by CORSIKA since
+      // CORSIKA handles proapgating the tau until the decay point.
+      Tpp::Tauola::setTauLifetime(0.0);
+
+      // don't decay any intermediaries - add these to the stack
+      // in particular, eta's, K0, Pi-0's.
+      // See Pg. 46 in the TAUOLA manual.
+      Tpp::Tauola::setEtaK0sPi(1, 1, 0);
+
+      // use the new currents (Pg. 40 in the user manual).
+      Tpp::Tauola::setNewCurrents(1);
+
+      // whether to use QCD radiative corrections in calculating
+      // leptonic tau decay. If this is enabled, TAUOLA generates
+      // an order of magnitude more photons than Pythia (for example)
+      Tpp::Tauola::setRadiation(radiative);
+
+      // we have to set the spin correlation state after
+      // we have initialized TAUOLA
+      Tpp::Tauola::spin_correlation.setAll(spincorr);
+
+      // initialize the TAUOLA library
+      Tpp::Tauola::initialize();
+    }
+  }
+
+  Decay::~Decay() {
+    CORSIKA_LOGGER_INFO(logger_, "Number of particles decayed using TAUOLA: {}", count_);
+  }
+
+  bool Decay::isDecayHandled(const Code code) {
+    // if known to tauola, it can decay
+    if (code == Code::WPlus || code == Code::WMinus || code == Code::Z0 ||
+        code == Code::TauPlus || code == Code::TauMinus || code == Code::Eta ||
+        code == Code::K0Short || code == Code::RhoPlus || code == Code::RhoMinus ||
+        code == Code::KStarPlus || code == Code::KStarMinus)
+      return true;
+    else
+      return false;
+  }
+
+  // perform the decay
+  template <typename TSecondaryView>
+  void Decay::doDecay(TSecondaryView& view) {
+    CORSIKA_LOGGER_DEBUG(logger_, "Decaying primary particle using TAUOLA...");
+    count_++;
+
+    // get the coordinate system of the decay
+    auto projectile = view.getProjectile();
+
+    const auto& cs{projectile.getMomentum().getCoordinateSystem()};
+
+    // get the TAUOLA particle from this projectile
+    auto particle{TauolaInterfaceParticle::from_projectile(projectile)};
+
+    CORSIKA_LOGGER_TRACE(
+        logger_, "[decay primary] code: {}, momentum: ({}, {}, {}) GeV, energy: {} GeV",
+        convert_from_PDG(static_cast<PDGCode>(particle->getPdgID())), particle->getPx(),
+        particle->getPy(), particle->getPz(), particle->getE());
+
+    // decay this particle using TAUOLA
+    decay(particle.get());
+
+    std::vector<TauolaInterfaceParticle*> daughters;
+
+    // we iterate over the daughter particles and check for
+    // short-lived secondaries that we have to decay again using TAUOLA.
+    // We only ever have to go one layer down in the hierarchy.
+    for (const auto& daughter : particle->getDaughters()) {
+
+      // check if this is decay can be handled by TAUOLA
+      if (abs(daughter->getPdgID()) == 20213 ||
+          isDecayHandled(convert_from_PDG(static_cast<PDGCode>(daughter->getPdgID())))) {
+        // decay this secondary
+        decay(daughter);
+        // and now add all the daughters of this particle to our list
+        for (const auto& child : daughter->getDaughters()) {
+
+          // check if this is further decayable (K0 and Eta's typically)
+          if (abs(child->getPdgID()) == 20213 ||
+              isDecayHandled(convert_from_PDG(static_cast<PDGCode>(child->getPdgID())))) {
+            // decay the cild
+            decay(child);
+
+            // and add it's daughters to the top-level daughters vector
+            for (const auto& grandchild : child->getDaughters()) {
+              daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(grandchild));
+            }
+          } else {
+            // if we are here, the original daughter is non-decayable so add it to our
+            // struct.
+            daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(child));
+          }
+        }
+
+      } else {
+        // if we get here, `daughter` is long-lived and doesn't require a decay.
+        daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(daughter));
+      }
+    }
+
+    // we now have a complete list of the daughter particles of this decay
+    // we now iterate over them and add them to CORSIKA stack.
+    for (const auto& daughter : daughters) {
+
+      // get particle's PDG ID
+      const auto pdg{daughter->getPdgID()};
+
+      // get the CORSIKA PID from our PDG ID
+      const auto pid{convert_from_PDG(static_cast<PDGCode>(pdg))};
+
+      // construct the momentum in the ROOT coordinate system
+      const MomentumVector P(cs, {daughter->getPx() * 1_GeV, daughter->getPy() * 1_GeV,
+                                  daughter->getPz() * 1_GeV});
+
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(P.getSquaredNorm() + mass * mass) - mass;
+
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "[decay result] code: {}, momentum: {} GeV, energy: {} GeV",
+                           pid, P.getComponents() / 1_GeV, daughter->getE());
+
+      // and add the particle to the CORSIKA stack
+      projectile.addSecondary(std::make_tuple(pid, Ekin, P.normalized()));
+    }
+  }
+
+  template <typename InterfaceParticle>
+  void Decay::decay(InterfaceParticle* particle) {
+
+    // return if we get a nullptr here
+    // this should never happen (but to be safe!)
+    // LCOV_EXCL_START
+    if (particle == nullptr) {
+      CORSIKA_LOGGER_WARN(logger_, "Got nullptr in decay routine. Skipping decay.");
+      return;
+    }
+    // LCOV_EXCL_STOP
+
+    // get the particle's PDG
+    const auto pdg{static_cast<PDGCode>(particle->getPdgID())};
+
+    // if this is a tau- or tau+, use the helicity in the decay
+    if (pdg == PDGCode::TauMinus || pdg == PDGCode::TauPlus) {
+
+      // get helicity from enum class
+      double helic;
+      switch (helicity_) {
+        case Helicity::LeftHanded:
+          helic = -1.;
+          break;
+        case Helicity::RightHanded:
+          helic = 1.;
+          break;
+        case Helicity::Unpolarized:
+          helic = 0.;
+          break;
+          // NOT testing undefined
+        default: // LCOV_EXCL_START
+          CORSIKA_LOGGER_ERROR(logger_, "Unknown helicity setting: {}", helicity_);
+          throw std::runtime_error("Unknown helicity setting!");
+          // LCOV_EXCL_STOP
+      }
+
+      double Polx{0};
+      double Poly{0};
+      double Polz{-helic};
+
+      // and call the decay routine with the polarization
+      Tpp::Tauola::decayOne(particle, false, Polx, Poly, Polz);
+
+    } else { // this is not a tau, so just do a standard decay
+      Tpp::Tauola::decayOne(particle);
+    }
+  }
+
+  template <typename TParticle>
+  TimeType Decay::getLifetime(TParticle const& particle) {
+    // Can be moved into Cascade.inl, see Issue #271
+    auto const pid = particle.getPID();
+    if ((pid == Code::TauMinus) || (pid == Code::TauPlus)) {
+      HEPEnergyType E = particle.getEnergy();
+      HEPMassType m = particle.getMass();
+      double const gamma = E / m;
+      TimeType const t0 = get_lifetime(pid);
+      auto const lifetime = gamma * t0;
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "[decay time] code: {}, t0: {} ns, momentum: {} GeV, energy: "
+                           "{} GeV, gamma: {}, tau: {} ns",
+                           particle.getPID(), t0 / 1_ns,
+                           particle.getMomentum().getComponents() / 1_GeV, E / 1_GeV,
+                           gamma, lifetime / 1_ns);
+      return lifetime;
+    } else
+      return std::numeric_limits<double>::infinity() * 1_s;
+  }
+
+  void Decay::PrintSummary() const { Tpp::Tauola::summary(); }
+
+} // namespace corsika::tauola

corsika/detail/modules/thinning/EMThinning.inl

+/*
+ * (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/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/SecondariesProcess.hpp>
+
+namespace corsika {
+
+  EMThinning::EMThinning(HEPEnergyType threshold, double maxWeight,
+                         bool const eraseParticles)
+      : threshold_{threshold}
+      , maxWeight_{maxWeight}
+      , eraseParticles_{eraseParticles} {}
+
+  template <typename TStackView>
+  void EMThinning::doSecondaries(TStackView& view) {
+    if (view.getSize() != 2) return;
+
+    auto projectile = view.getProjectile();
+    if (!is_em(projectile.getPID())) return;
+
+    double const parentWeight = projectile.getWeight();
+    if (parentWeight >= maxWeight_) { return; }
+
+    auto particle1 = view.begin();
+    auto particle2 = ++(view.begin());
+
+    // skip non-EM interactions
+    if (!is_em(particle1.getPID()) || !is_em(particle2.getPID())) { return; }
+
+    /* skip high-energy interactions
+     * We consider only the primary energy. A more sophisticated algorithm might
+     * sort out above-threshold secondaries and apply thinning only on the subset of
+     * those below the threshold.
+     */
+    if (projectile.getEnergy() > threshold_) { return; }
+
+    corsika::HEPEnergyType const Esum = particle1.getEnergy() + particle2.getEnergy();
+
+    double const p1 = particle1.getEnergy() / Esum;
+    double const p2 = particle2.getEnergy() / Esum;
+
+    // weight factors
+    auto const f1 = 1 / p1;
+    auto const f2 = 1 / p2;
+
+    // max. allowed weight factor
+    double const fMax = maxWeight_ / parentWeight;
+
+    if (f1 <= fMax && f2 <= fMax) {
+      /* cannot possibly reach wmax in this vertex; apply
+        Hillas thinning; select one of the two */
+      if (uniform_(rng_) <= p1) { // keep 1st with probability p1
+        particle2.setWeight(0);
+        particle1.setWeight(parentWeight * f1);
+      } else { // keep 2nd
+        particle1.setWeight(0);
+        particle2.setWeight(parentWeight * f2);
+      }
+    } else { // weight limitation kicks in, do statistical thinning
+      double const f1prime = std::min(f1, fMax);
+      double const f2prime = std::min(f2, fMax);
+
+      if (uniform_(rng_) * f1prime <= 1) { // keep 1st
+        particle1.setWeight(parentWeight * f1prime);
+      } else {
+        particle1.setWeight(0);
+      }
+      if (uniform_(rng_) * f2prime <= 1) { // keep 2nd
+        particle2.setWeight(parentWeight * f2prime);
+      } else {
+        particle2.setWeight(0);
+      }
+    }
+
+    // erase discared particles in case of multithinning
+    if (eraseParticles_) {
+      for (auto& p : view) {
+        if (auto const w = p.getWeight(); w == 0) { p.erase(); }
+      }
+    } else {
+      return;
+    }
+  }
+
+} // namespace corsika