qgsjetII coverage

README.rst

@@ -56,19 +56,16 @@ Get in contact
 
 Installation
 ------------
-
 CORSIKA 8 is tested regularly at least on gcc7.3.0 and clang-8.0.0. 
 
 Prerequisites
 ~~~~~~~~~~~~~
-
 You will also need:
-
-- Python 3 (supported versions are Python >= 3.6), with pip
-- conan (via pip)
-- cmake 
-- git
-- g++, gfortran, binutils, make
+  * Python 3 (supported versions are Python >= 3.6), with pip
+  * conan (via pip)
+  * cmake 
+  * git
+  * g++, gfortran, binutils, make
 
 On a bare Ubuntu 20.04, just add:
 ::

corsika/detail/modules/qgsjetII/Interaction.inl

@@ -66,15 +66,18 @@ namespace corsika::qgsjetII {
         iTarget = targetA;
         if (iTarget > maxMassNumber_ || iTarget <= 0) {
           std::ostringstream txt;
-          txt << "QgsjetII target outside range. iTarget=" << iTarget;
+          txt << "QgsjetII target outside range. Atarget=" << iTarget;
           throw std::runtime_error(txt.str().c_str());
         }
       }
       int iProjectile = 1;
       if (is_nucleus(beamId)) {
         iProjectile = Abeam;
-        if (iProjectile > maxMassNumber_ || iProjectile <= 0)
-          throw std::runtime_error("QgsjetII target outside range. ");
+        if (iProjectile > maxMassNumber_ || iProjectile <= 0) {
+          std::ostringstream txt;
+          txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
+          throw std::runtime_error(txt.str().c_str());
+	}
       }
 
       CORSIKA_LOG_DEBUG(
@@ -234,7 +237,7 @@ namespace corsika::qgsjetII {
 
       int targetMassNumber = 1;     // proton
       if (is_nucleus(targetCode)) { // nucleus
-        targetMassNumber = get_nucleus_A(targetCode);
+        targetMassNumber = get_nucleus_A(targetCode);       
         if (targetMassNumber > maxMassNumber_)
           throw std::runtime_error("QgsjetII target mass outside range.");
       } else {

tests/framework/testProcessSequence.cpp

@@ -696,7 +696,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
     checkInteract = 0;
     checkSec = 0;
     checkCont = 0;
-    particle.data_[0] = -100; // data negative
+    particle.data_[0] = -100; // data positive
     sequence4.selectInteraction(view, lambda_select);
     sequence4.doSecondaries(view);
     sequence4.selectDecay(view, time_select);

tests/modules/testQGSJetII.cpp

@@ -169,4 +169,56 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
     CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
           Approx(0).margin(1e-2));
   }
+  
+  SECTION("InteractionInterface Nuclei") {
+
+    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+        Code::Nucleus, 20, 10, 10100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+        *csPtr);
+    setup::StackView& view = *(secViewPtr.get());
+    auto particle = stackPtr->first();
+    auto projectile = secViewPtr->getProjectile();
+    auto const projectileMomentum = projectile.getMomentum();
+
+    corsika::qgsjetII::Interaction model;
+    model.doInteraction(view);
+    [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
+
+    CHECK(length / (1_g / square(1_cm)) == Approx(20.13).margin(0.1));
+  }
+
+  SECTION("Heavy nuclei") {
+    
+    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+        Code::Nucleus, 1000, 1000, 1100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+        *csPtr);
+    setup::StackView& view = *(secViewPtr.get());
+    auto particle = stackPtr->first();
+    auto projectile = secViewPtr->getProjectile();
+    auto const projectileMomentum = projectile.getMomentum();
+    
+    corsika::qgsjetII::Interaction model;
+        
+    CHECK_THROWS(model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 10., 1000.));
+    CHECK_THROWS(model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 1000., 10.));
+    CHECK_THROWS(model.doInteraction(view));
+    CHECK_THROWS(model.getInteractionLength(particle));
+  }
+
+  SECTION("Allowed Particles") {
+    
+    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+        Code::Electron, 0, 0, 1100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+        *csPtr);
+    setup::StackView& view = *(secViewPtr.get());
+    auto particle = stackPtr->first();
+    auto projectile = secViewPtr->getProjectile();
+    auto const projectileMomentum = projectile.getMomentum();
+    
+    corsika::qgsjetII::Interaction model;
+        
+    GrammageType const  length = model.getInteractionLength(particle);
+    CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());    
+  }
+
 }