fixed sibyll tests

corsika/detail/modules/sibyll/InteractionModel.inl

@@ -51,14 +51,14 @@ namespace corsika::sibyll {
       targA = get_nucleus_A(targetId);
     }
 
-    if ((is_nucleus(targetId) &&
-         (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) &&
-        targetId != Code::Proton && targetId != Code::Hydrogen &&
-        targetId != Code::Neutron) {
+    if (is_nucleus(targetId)) {
+      if (targA != 1 && (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) {
+        throw std::runtime_error("Target outside of allowed range for SIBYLL");
+      }
+    } else if (targetId != Code::Proton && targetId != Code::Neutron) {
       throw std::runtime_error("Target cannot be handled by SIBYLL");
     }
-
-    if (is_nucleus(projectileId) && !corsika::sibyll::canInteract(projectileId)) {
+    if (is_nucleus(projectileId) || !corsika::sibyll::canInteract(projectileId)) {
       throw std::runtime_error("Projectile cannot be handled by SIBYLL");
     }
   }
@@ -73,7 +73,8 @@ namespace corsika::sibyll {
 
     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)
+        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

corsika/detail/modules/sibyll/NuclearInteractionModel.inl

@@ -49,7 +49,7 @@ namespace corsika::sibyll {
     hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn, 1, targetA); // throws
 
     // projectile limits:
-    if (is_nucleus(projectileId)) {
+    if (!is_nucleus(projectileId)) {
       throw std::runtime_error("can only handle nuclear projectile");
     }
     if (projectileA >= getMaxNucleusAProjectile() || projectileA < 2) {

corsika/modules/sibyll/InteractionModel.hpp

+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <tuple>
+
+namespace corsika::sibyll {
+
+  /**
+   * @brief sibyll::InteractionModel provides the SIBYLL proton-nucleus interaction model.
+   *
+   * This is a TModel argument for InteractionProcess<TModel>.
+   */
+
+  class InteractionModel {
+
+  public:
+    InteractionModel();
+    ~InteractionModel();
+
+    /**
+     * @brief Set the Verbose flag.
+     *
+     * If flag is true, SIBYLL will printout additional secondary particle information
+     * lists, etc.
+     *
+     * @param flag to switch.
+     */
+    void setVerbose(bool const flag);
+
+    /**
+     * @brief evaluated validity of collision system.
+     *
+     * sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
+     * neutrons (p,n == nucleon).
+     */
+    void constexpr isValid(Code const projectileId, Code const targetId,
+                           HEPEnergyType const sqrtSnn, unsigned int const projectileA,
+                           unsigned int const targetA) const;
+
+    /**
+     * @brief returns inelastic AND elastic cross sections.
+     *
+     * These cross sections must correspond to the process described in doInteraction
+     * AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
+     * nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
+     * Sibyll must be useful inside the NuclearInteraction model, which requires that.
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+     * @param Aprojectil is the mass number of the projectils, if it is a nucleus
+     * @param Atarget is the mass number of the target, if it is a nucleus
+     *
+     * @return a tuple of: inelastic cross section, elastic cross section
+     */
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code const projectile, Code const target, HEPEnergyType const sqrtSnn,
+        unsigned int const Aprojectile = 1, unsigned int const Atarget = 1) const;
+
+    /**
+     * @brief returns inelastic (production) cross section.
+     *
+     * This cross section must correspond to the process described in doInteraction.
+     * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+     * @param Aprojectil is the mass number of the projectils, if it is a nucleus
+     * @param Atarget is the mass number of the target, if it is a nucleus
+     *
+     * @return inelastic cross section
+     * elastic cross section
+     */
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     HEPEnergyType const sqrtSnn,
+                                     unsigned int const Aprojectile = 1,
+                                     unsigned int const Atarget = 1) const {
+      return std::get<0>(
+          getCrossSectionInelEla(projectile, target, sqrtSnn, Aprojectile, Atarget));
+    }
+
+    /**
+     * In this function SIBYLL is called to produce one event. The
+     * event is copied (and boosted) into the shower lab frame.
+     */
+
+    template <typename TSecondaries>
+    void doInteraction(TSecondaries&, COMBoost const& boost, Code const projectile,
+                       Code const target, HEPEnergyType const sqrtSnn,
+                       unsigned int const Aprojectile = 1,
+                       unsigned int const Atarget = 1);
+
+  private:
+    HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
+    HEPEnergyType constexpr getMaxEnergyCoM() const { return maxEnergyCoM_; }
+
+    // hard model limits
+    static HEPEnergyType constexpr minEnergyCoM_ = 10. * 1e9 * electronvolt;
+    static HEPEnergyType constexpr maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
+    static unsigned int constexpr maxTargetMassNumber_ = 18;
+    static unsigned int constexpr minNuclearTargetA_ = 4;
+
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
+
+    // data members
+    int count_ = 0;
+    int nucCount_ = 0;
+    bool sibyll_listing_;
+  };
+
+} // namespace corsika::sibyll
+
+#include <corsika/detail/modules/sibyll/InteractionModel.inl>

corsika/modules/sibyll/NuclearInteractionModel.hpp

+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
+namespace corsika::sibyll {
+
+  /**
+   * The sibyll::NuclearInteractionModel provides the SIBYLL semi superposition model.
+   *
+   * can transform a proton-nucleus interaction model into a nucleus-nucleus interaction
+   * model.
+   *
+   * @tparam TNucleonModel
+   */
+  template <class TEnvironment, class TNucleonModel>
+  class NuclearInteractionModel {
+
+  public:
+    NuclearInteractionModel(TNucleonModel&, TEnvironment const&);
+    ~NuclearInteractionModel();
+
+    void constexpr isValid(Code const projectileId, Code const targetId,
+                           HEPEnergyType const sqrtSnn, unsigned int const projectileA,
+                           unsigned int const targetA) const;
+
+    void initializeNuclearCrossSections();
+    void printCrossSectionTable(Code) const;
+    CrossSectionType readCrossSectionTable(int const, Code const,
+                                           HEPEnergyType const) const;
+    HEPEnergyType getMinEnergyPerNucleonCoM() const { return gMinEnergyPerNucleonCoM_; }
+    HEPEnergyType getMaxEnergyPerNucleonCoM() const { return gMaxEnergyPerNucleonCoM_; }
+    unsigned int constexpr getMaxNucleusAProjectile() const {
+      return gMaxNucleusAProjectile_;
+    }
+    unsigned int constexpr getMaxNFragments() const { return gMaxNFragments_; }
+    unsigned int constexpr getNEnergyBins() const { return gNEnBins_; }
+
+    CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const,
+                                     unsigned int const projectileA = 1,
+                                     unsigned int const targetA = 1) const;
+
+    template <typename TSecondaryView>
+    void doInteraction(TSecondaryView&, COMBoost const& boost, Code const, Code const,
+                       HEPEnergyType const, unsigned int const projectileA = 1,
+                       unsigned int const targetA = 1);
+
+  private:
+    int count_ = 0;
+    int nucCount_ = 0;
+
+    TEnvironment const& environment_;
+    TNucleonModel& hadronicInteraction_;
+    std::map<Code, int> targetComponentsIndex_;
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
+    static unsigned int constexpr gNSample_ =
+        500; // number of samples in MC estimation of cross section
+    static unsigned int constexpr gMaxNucleusAProjectile_ = 56;
+    static unsigned int constexpr gNEnBins_ = 6;
+    static unsigned int constexpr gMaxNFragments_ = 60;
+    // energy limits defined by table used for cross section in signuc.f
+    // 10**1 GeV to 10**6 GeV
+    static HEPEnergyType constexpr gMinEnergyPerNucleonCoM_ = 10. * 1e9 * electronvolt;
+    static HEPEnergyType constexpr gMaxEnergyPerNucleonCoM_ = 1.e6 * 1e9 * electronvolt;
+  };
+
+} // namespace corsika::sibyll
+
+#include <corsika/detail/modules/sibyll/NuclearInteractionModel.inl>

tests/modules/testSibyll.cpp

@@ -146,8 +146,8 @@ TEST_CASE("SibyllInterface", "modules") {
     CHECK_THROWS(model.isValid(Code::Proton, Code::Iron, 100_GeV, 1, 56));
     CHECK_NOTHROW(model.isValid(Code::Proton, Code::Oxygen, 100_GeV, 1, 16));
     // beam particles
-    CHECK_NOTHROW(model.isValid(Code::Electron, Code::Oxygen, 100_GeV, 1, 1));
-    CHECK_NOTHROW(model.isValid(Code::Nucleus, Code::Oxygen, 100_GeV, 1, 20));
+    CHECK_THROWS(model.isValid(Code::Electron, Code::Oxygen, 100_GeV, 1, 1));
+    CHECK_THROWS(model.isValid(Code::Nucleus, Code::Oxygen, 100_GeV, 1, 20));
     // energy too low
     CHECK_THROWS(model.isValid(Code::Proton, Code::Proton, 9_GeV, 1, 1));
     CHECK_NOTHROW(model.isValid(Code::Proton, Code::Proton, 11_GeV, 1, 1));