adapted InteractionCounter

corsika/detail/framework/process/InteractionCounter.inl

@@ -18,22 +18,20 @@ namespace corsika {
 
   template <class TCountedProcess>
   template <typename TSecondaryView>
-  inline void InteractionCounter<TCountedProcess>::doInteraction(TSecondaryView& view) {
-    auto const projectile = view.getProjectile();
-    auto const massNumber = projectile.getNode()
-                                ->getModelProperties()
-                                .getNuclearComposition()
-                                .getAverageMassNumber();
+  inline void InteractionCounter<TCountedProcess>::doInteraction(
+      TSecondaryView& view, Code const projectileId, Code const targetId,
+      FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+    size_t const massNumber = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
     auto const massTarget = massNumber * constants::nucleonMass;
-    histogram_.fill(projectile.getPID(), projectile.getEnergy(), massTarget);
-    process_.doInteraction(view);
+    histogram_.fill(projectileId, projectileP4.getTimeLikeComponent(), massTarget);
+    process_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
   }
 
   template <class TCountedProcess>
-  template <typename TParticle>
-  inline GrammageType InteractionCounter<TCountedProcess>::getInteractionLength(
-      TParticle const& particle) const {
-    return process_.getInteractionLength(particle);
+  inline CrossSectionType InteractionCounter<TCountedProcess>::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+    return process_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
   }
 
   template <class TCountedProcess>

corsika/detail/modules/urqmd/UrQMD.inl

@@ -43,7 +43,7 @@ namespace corsika::urqmd {
       throw std::runtime_error("UrQMD projectile is not a compatible hadron.");
     }
     if (!is_nucleus(targetId)) {
-      throw std::runtime_error("UrQMD target is not a nucleus.");
+      throw std::runtime_error("UrQMD target is not a nucleus .");
     }
   }
 

corsika/framework/process/InteractionCounter.hpp

@@ -10,24 +10,25 @@
 
 #include <corsika/framework/process/InteractionHistogram.hpp>
 #include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
 
 namespace corsika {
 
-  /*!
-    @ingroup Processes
-    @{
-
+  /**
+   * @ingroup Processes
+   * @{
+   *
    * Wrapper around an InteractionProcess that fills histograms of the number
    * of calls to `doInteraction()` binned in projectile energy (both in
-   * lab and center-of-mass frame) and species
+   * lab and center-of-mass frame) and species.
    *
-   * Use by wrapping a normal InteractionProcess
+   * Use by wrapping a normal InteractionProcess:
    * @code{.cpp}
    * InteractionProcess collision1;
    * InteractionClounter<collision1> counted_collision1;
    * @endcode
-   *
    */
+
   template <class TCountedProcess>
   class InteractionCounter
       : public InteractionProcess<InteractionCounter<TCountedProcess>> {
@@ -35,17 +36,24 @@ namespace corsika {
   public:
     InteractionCounter(TCountedProcess& process);
 
-    //! wrapper around internall process doInteraction
+    /**
+     * Wrapper around internal process doInteraction.
+     */
     template <typename TSecondaryView>
-    void doInteraction(TSecondaryView& view);
+    void doInteraction(TSecondaryView& view, Code const, Code const, FourMomentum const&,
+                       FourMomentum const&);
 
-    ///! returns internal process getInteractionLength
-    template <typename TParticle>
-    GrammageType getInteractionLength(TParticle const& particle) const;
+    /**
+     * Wrapper around internal process getCrossSection.
+     */
+    CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+                                     FourMomentum const&) const;
 
-    /** returns the filles histograms
-        @return InteractionHistogram, which contains the histogram data
-    */
+    /**
+     * returns the filles histograms.
+     *
+     * @return InteractionHistogram, which contains the histogram data
+     */
     InteractionHistogram const& getHistogram() const;
 
   private:

examples/staticsequence_example.cpp

@@ -105,7 +105,6 @@ void modular() {
 int main() {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   std::cout << "staticsequence_example" << std::endl;
 

tests/framework/testInteractionCounter.cpp

@@ -7,17 +7,11 @@
  */
 
 #include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/NuclearComposition.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/geometry/RootCoordinateSystem.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
-#include <SetupTestStack.hpp>
-#include <SetupTestEnvironment.hpp>
-
 #include <catch2/catch.hpp>
 
 #include <numeric>
@@ -32,37 +26,46 @@ using namespace corsika;
 const std::string refDataDir = std::string(REFDATADIR); // from cmake
 
 struct DummyProcess {
-  template <typename TParticle>
-  GrammageType getInteractionLength(TParticle const&) {
-    return 100_g / 1_cm / 1_cm;
+
+  CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+                                   FourMomentum const&) {
+    return 100_mb;
   }
+
   template <typename TParticle>
-  void doInteraction(TParticle&) {}
+  void doInteraction(TParticle&, Code const, Code const, FourMomentum const&,
+                     FourMomentum const&) {}
 };
 
-TEST_CASE("InteractionCounter", "[process]") {
+struct DummyOutput {
+  /* can do nothing */
+};
+
+TEST_CASE("InteractionCounter", "process") {
 
   logging::set_level(logging::level::info);
 
   DummyProcess d;
   InteractionCounter countedProcess(d);
 
-  SECTION("getInteractionLength") {
-    CHECK(countedProcess.getInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
-  }
+  auto const rootCS = get_root_CoordinateSystem();
+  DummyOutput output;
 
-  auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
-  [[maybe_unused]] auto& env_dummy = env;
+  SECTION("cross section pass-through") {
+    CHECK(countedProcess.getCrossSection(
+              Code::Oxygen, Code::Proton, {10_GeV, {rootCS, {0_eV, 0_eV, 0_eV}}},
+              {10_GeV, {rootCS, {0_eV, 0_eV, 0_eV}}}) == 100_mb);
+  }
 
-  SECTION("DoInteraction nucleus") {
+  SECTION("doInteraction nucleus") {
     unsigned short constexpr A = 14, Z = 7;
-    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
-        get_nucleus_code(A, Z), 105_TeV, (setup::Environment::BaseNodeType* const)nodePtr,
-        *csPtr);
-    CHECK(stackPtr->getEntries() == 1);
-    CHECK(secViewPtr->getEntries() == 0);
+    Code const pid = get_nucleus_code(A, Z);
 
-    countedProcess.doInteraction(*secViewPtr);
+    countedProcess.doInteraction(
+        output, pid, Code::Oxygen,
+        {sqrt(static_pow<2>(105_TeV) + static_pow<2>(get_mass(pid))),
+         {rootCS, {105_TeV, 0_GeV, 0_GeV}}},
+        {Oxygen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
 
     auto const& h = countedProcess.getHistogram().labHist();
     CHECK(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
@@ -99,14 +102,14 @@ TEST_CASE("InteractionCounter", "[process]") {
     }
   }
 
-  SECTION("DoInteraction Lambda") {
-    auto constexpr code = Code::Lambda0;
-    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
-        code, 105_TeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
-    CHECK(stackPtr->getEntries() == 1);
-    CHECK(secViewPtr->getEntries() == 0);
+  SECTION("doInteraction Lambda") {
+    auto constexpr pid = Code::Lambda0;
 
-    countedProcess.doInteraction(*secViewPtr);
+    countedProcess.doInteraction(
+        output, pid, Code::Oxygen,
+        {sqrt(static_pow<2>(105_TeV) + static_pow<2>(get_mass(pid))),
+         {rootCS, {105_TeV, 0_GeV, 0_GeV}}},
+        {Oxygen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
 
     auto const& h = countedProcess.getHistogram().labHist();
     CHECK(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);