HadronicElasticModel, first try

Processes/CMakeLists.txt

@@ -4,6 +4,7 @@ add_subdirectory (Sibyll)
 add_subdirectory (StackInspector)
 add_subdirectory (TrackWriter)
 add_subdirectory (TrackingLine)
+add_subdirectory (HadronicElasticModel)
 
 #add_custom_target(CORSIKAprocesses)
 add_library (CORSIKAprocesses INTERFACE)

Processes/HadronicElasticModel/CMakeLists.txt

+
+set (
+  MODEL_SOURCES
+  HadronicElasticModel.cc
+  )
+
+set (
+  MODEL_HEADERS
+  HadronicElasticModel.h
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/hadronic_elastic_model
+  )
+
+add_library (ProcessHadronicElasticModel STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessHadronicElasticModel ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+  ProcessHadronicElasticModel
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+#  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessHadronicElasticModel
+  CORSIKAunits
+  CORSIKAgeometry
+  CORSIKAsetup
+  CORSIKAutilities
+  )
+
+target_include_directories (
+  ProcessHadronicElasticModel 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS ProcessHadronicElasticModel
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+  )
+
+
+# --------------------
+# code unit testing
+#~ add_executable (ProcessHadronicElasticModel testProcessHadronicElasticModel.cc)
+
+#~ target_link_libraries (
+  #~ testProcessHadronicElasticModel
+  #~ ProcessHadronicElasticModel
+  #~ CORSIKAsetup
+  #~ CORSIKAgeometry
+  #~ CORSIKAunits
+  #~ CORSIKAthirdparty # for catch2
+  #~ )
+#~ CORSIKA_ADD_TEST(testNullModel)

Processes/HadronicElasticModel/HadronicElasticModel.cc

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
+
+namespace corsika::process::HadronicElasticModel {
+
+  void HadronicElasticInteraction::Init() {}
+
+  HadronicElasticInteraction::HadronicElasticInteraction(
+      corsika::environment::Environment const& env,
+      corsika::units::si::CrossSectionType x, corsika::units::si::CrossSectionType y)
+      : fX(x)
+      , fY(y)
+      , fEnvironment(env) {}
+
+} // namespace corsika::process::HadronicElasticModel

Processes/HadronicElasticModel/HadronicElasticModel.h

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_HadronicElasticInteraction_h
+#define _include_HadronicElasticInteraction_h
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/random/ExponentialDistribution.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalConstants.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+
+namespace corsika::process::HadronicElasticModel {
+  /**
+   * A simple model for elastic hadronic interactions based on the formulas
+   * in Gaisser, Engel, Resconi, Cosmic Rays and Particle Physics (Cambridge Univ. Press,
+   * 2016), section 4.2 and Donnachie, Landshoff, Phys. Lett. B 296, 227 (1992)
+   *
+   * Currently only \f$p\f$ projectiles are supported and cross-sections are assumed to be
+   * \f$pp\f$-like even for nuclei.
+   */
+  class HadronicElasticInteraction
+      : public corsika::process::InteractionProcess<HadronicElasticInteraction> {
+  private:
+    corsika::units::si::CrossSectionType const fX, fY;
+    static double constexpr fEpsilon = 0.0808;
+    static double constexpr fEta = 0.4525;
+
+    using SquaredHEPEnergyType = decltype(corsika::units::si::HEPEnergyType() *
+                                          corsika::units::si::HEPEnergyType());
+
+    corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream(
+            "HadronicElasticModel");
+    corsika::environment::Environment const& fEnvironment;
+
+    auto B(corsika::units::si::HEPEnergyType sqrtS) const {
+      using namespace corsika::units::si;
+      auto constexpr a =
+          2.5 / (1_GeV * 1_GeV); // read off by eye from Gaisser, Engel, Resconi, fig. 4.9
+      auto constexpr b = 10 / (1_GeV * 1_GeV);
+      return a * std::log10(sqrtS / 10_GeV) + b;
+    }
+
+    corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const {
+      using namespace corsika::units::si;
+      // according to Gaisser, Engel, Resconi, eq. (4.41)
+      // assuming every target behaves like a proton
+      CrossSectionType const sigmaTot = fX * pow(s * (1 / (1_GeV * 1_GeV)), fEpsilon) +
+                                        fY * pow(s * (1 / (1_GeV * 1_GeV)), -fEta);
+
+      // we ignore rho because rho^2 is just ~2 %
+      auto const sigmaElastic =
+          sigmaTot * sigmaTot /
+          (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(sqrt(s))));
+      return sigmaElastic;
+    }
+
+  public:
+    HadronicElasticInteraction(corsika::environment::Environment const&,
+                               units::si::CrossSectionType x = 2.17e-5 * units::si::barn,
+                               units::si::CrossSectionType y = 5.608e-5 *
+                                                               units::si::barn);
+    void Init();
+
+    template <typename Particle, typename Track>
+    corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
+      using namespace corsika::units::si;
+      if (p.GetPID() == corsika::particles::Code::Proton) {
+        auto const* currentNode =
+            fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+        auto const& mediumComposition =
+            currentNode->GetModelProperties().GetNuclearComposition();
+
+        auto const& components = mediumComposition.GetComponents();
+        auto const& fractions = mediumComposition.GetFractions();
+
+        auto const projectileMomentum = p.GetMomentum();
+        auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
+        auto const projectileEnergy = p.GetEnergy();
+
+        auto const avgCrossSection = [&]() {
+          CrossSectionType avgCrossSection = 0_barn;
+
+          for (size_t i = 0; i < fractions.size(); ++i) {
+            auto const targetMass = corsika::particles::GetMass(components[i]);
+            auto const s = detail::static_pow<2>(projectileEnergy + targetMass) -
+                           projectileMomentumSquaredNorm;
+            avgCrossSection += CrossSection(s) * fractions[i];
+          }
+
+          std::cout << "avgCrossSection: " << avgCrossSection / 1_mbarn << " mb"
+                    << std::endl;
+
+          return avgCrossSection;
+        }();
+
+        auto const avgTargetMassNumber = mediumComposition.GetAverageMassNumber();
+
+        GrammageType const interactionLength =
+            avgTargetMassNumber * corsika::units::constants::u / avgCrossSection;
+
+        return interactionLength;
+      } else {
+        return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+      }
+    }
+
+    template <typename Particle, typename Stack>
+    corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
+      if (p.GetPID() != corsika::particles::Code::Proton) {
+        return process::EProcessReturn::eOk;
+      }
+
+      using namespace units::si;
+
+      const auto* currentNode =
+          fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+      const auto& mediumComposition =
+          currentNode->GetModelProperties().GetNuclearComposition();
+      const auto& components = mediumComposition.GetComponents();
+
+      std::vector<units::si::CrossSectionType> cross_section_of_components(
+          mediumComposition.GetComponents().size());
+
+      auto const projectileMomentum = p.GetMomentum();
+      auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
+      auto const projectileEnergy = p.GetEnergy();
+
+      for (size_t i = 0; i < components.size(); ++i) {
+        auto const targetMass = corsika::particles::GetMass(components[i]);
+        auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) -
+                       projectileMomentumSquaredNorm;
+        cross_section_of_components[i] = CrossSection(s);
+      }
+
+      const auto targetCode = currentNode->GetModelProperties().SampleTarget(
+          cross_section_of_components, fRNG);
+
+      auto const targetMass = corsika::particles::GetMass(targetCode);
+
+      std::uniform_real_distribution phiDist(0., 2 * M_PI);
+
+      utl::COMBoost const boost(projectileEnergy, projectileMomentum, targetMass);
+      auto const [eProjectileCoM, pProjectileCoM] =
+          boost.toCoM(projectileEnergy, projectileMomentum);
+      auto const [eTargetCoM, pTargetCoM] = boost.toCoM(
+          targetMass, geometry::Vector<units::si::hepmomentum_d>(
+                          projectileMomentum.GetCoordinateSystem(), {0_eV, 0_eV, 0_eV}));
+
+      auto const pProjectileCoMSqNorm = pProjectileCoM.squaredNorm();
+      auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm);
+
+      auto const sqrtS = eProjectileCoM + eTargetCoM;
+      auto const s = units::si::detail::static_pow<2>(sqrtS);
+
+      auto const B = this->B(sqrtS);
+      std::cout << B << std::endl;
+
+      corsika::random::ExponentialDistribution tDist(1 / B);
+      auto const absT = [&]() {
+        decltype(tDist(fRNG)) absT;
+        auto const maxT = 4 * pProjectileCoMSqNorm;
+
+        do {
+          // |t| cannot become arbitrarily large, max. given by GEN eq. (4.16), so we just
+          // throw again until we have an acceptable value note that the formula holds in
+          // any frame despite of what is stated there
+          absT = tDist(fRNG);
+        } while (absT >= maxT);
+
+        return absT;
+      }();
+
+      auto constexpr GeVsq = 1_GeV * 1_GeV;
+      std::cout << "s = " << s / GeVsq << " GeV²; absT = " << absT / GeVsq
+                << " GeV² (max./GeV² = " << 4 * projectileMomentumSquaredNorm / GeVsq
+                << ')' << std::endl;
+
+      auto const theta =
+          2 *
+          asin(sqrt(absT / (4 * pProjectileCoMSqNorm))); // would work for lab frame, too
+      auto const phi = phiDist(fRNG);
+
+      geometry::QuantityVector<units::si::hepmomentum_d> const scatteredMomentum{
+          pProjectileCoMNorm * sin(theta) * cos(phi),
+          pProjectileCoMNorm * sin(theta) * cos(phi), pProjectileCoMNorm * cos(theta)};
+      auto const [eProjectileScattered, pProjectileScattered] =
+          boost.fromCoM(eProjectileCoM, scatteredMomentum);
+
+      p.SetMomentum(pProjectileScattered);
+      p.SetEnergy(sqrt(pProjectileScattered.squaredNorm() +
+                       units::si::detail::static_pow<2>(targetMass)));
+
+      return process::EProcessReturn::eOk;
+    }
+  };
+} // namespace corsika::process::HadronicElasticModel
+
+#endif