Merge branch '133-use-nuclear-stack-extension-in-nuclear-interaction' into 'master'

Resolve "use nuclear stack extension in nuclear interaction"

Closes #133

See merge request !72

.gitlab-ci.yml

@@ -18,8 +18,10 @@ build:
     - cmake ..
     - cmake --build .
     - ctest -j4 -V >& test.log
-    - gzip -9 -S .gz test.log
+  after_script:
+    - cd build
     - ls
+    - gzip -v -9 -S .gz test.log
     - pwd
   artifacts:
     expire_in: 1 week

Documentation/Examples/cascade_example.cc

@@ -71,12 +71,22 @@ public:
 
   template <typename Particle>
   bool isBelowEnergyCut(Particle& p) const {
-    // FOR NOW: center-of-mass energy hard coded
-    const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
-    if (p.GetEnergy() < fECut || Ecm < 10_GeV)
-      return true;
-    else
-      return false;
+    // nuclei
+    if (p.GetPID() == corsika::particles::Code::Nucleus) {
+      auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
+      auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
+      if (ElabNuc < fECut || EcmNN < 10_GeV)
+        return true;
+      else
+        return false;
+    } else {
+      // TODO: center-of-mass energy hard coded
+      const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+      if (p.GetEnergy() < fECut || Ecm < 10_GeV)
+        return true;
+      else
+        return false;
+    }
   }
 
   bool isEmParticle(Code pCode) const {
@@ -244,7 +254,6 @@ int main() {
 
   corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
   corsika::process::sibyll::Interaction sibyll(env);
-  //  corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
   corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
   corsika::process::sibyll::Decay decay;
   ProcessCut cut(20_GeV);
@@ -265,9 +274,14 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const Code beamCode = Code::Proton;
+  const Code beamCode = Code::Nucleus;
+  const int nuclA = 56;
+  const int nuclZ = int(nuclA / 2.15 + 0.7);
+  const HEPMassType mass = corsika::particles::Proton::GetMass() * nuclZ +
+                           (nuclA - nuclZ) * corsika::particles::Neutron::GetMass();
   const HEPEnergyType E0 =
-      100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
+      nuclA *
+      100_GeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
   double theta = 0.;
   double phi = 0.;
 
@@ -275,7 +289,7 @@ int main() {
     auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
       return sqrt(Elab * Elab - m * m);
     };
-    HEPMomentumType P0 = elab2plab(E0, corsika::particles::GetMass(beamCode));
+    HEPMomentumType P0 = elab2plab(E0, mass);
     auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
       return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
                              -ptot * cos(theta));
@@ -287,7 +301,7 @@ int main() {
     cout << "input angles: theta=" << theta << " phi=" << phi << endl;
     cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
     Point pos(rootCS, 0_m, 0_m, 0_m);
-    stack.AddParticle(beamCode, E0, plab, pos, 0_ns);
+    stack.AddParticle(beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ);
   }
 
   // define air shower object, run simulation

Framework/StackInterface/testStackInterface.cc

@@ -38,6 +38,7 @@ class TestStackData {
 public:
   // these functions are needed for the Stack interface
   void Init() {}
+  void Dump() const {}
   void Clear() { fData.clear(); }
   unsigned int GetSize() const { return fData.size(); }
   unsigned int GetCapacity() const { return fData.size(); }

Processes/Sibyll/Interaction.h

@@ -56,6 +56,11 @@ namespace corsika::process::sibyll {
 
     bool WasInitialized() { return fInitialized; }
 
+    bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+      using namespace corsika::units::si;
+      return (10_GeV < ecm) && (ecm < 1_PeV);
+    }
+
     std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
                int>
     GetCrossSection(const corsika::particles::Code BeamId,
@@ -124,6 +129,7 @@ namespace corsika::process::sibyll {
                 << " beam can interact:" << kInteraction << endl
                 << " beam pid:" << p.GetPID() << endl;
 
+      // TODO: move limits into variables
       if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {
 
         // get target from environment

Processes/Sibyll/NuclearInteraction.h

@@ -47,7 +47,7 @@ namespace corsika::process::sibyll {
       using std::endl;
 
       // initialize hadronic interaction module
-      // TODO: save to run multiple initializations?
+      // TODO: safe to run multiple initializations?
       if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
 
       // initialize nuclib
@@ -56,38 +56,40 @@ namespace corsika::process::sibyll {
     }
 
     // TODO: remove number of nucleons, avg target mass is available in environment
-    std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
-               int>
-    GetCrossSection(const corsika::particles::Code BeamId,
-                    const corsika::particles::Code TargetId,
-                    const corsika::units::si::HEPEnergyType CoMenergy) {
+    template <typename Particle>
+    std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
+    GetCrossSection(Particle const& p, const corsika::particles::Code TargetId) {
       using namespace corsika::units::si;
       double sigProd;
+      auto const pCode = p.GetPID();
+      if (pCode != corsika::particles::Code::Nucleus)
+        throw std::runtime_error(
+            "NuclearInteraction: GetCrossSection: particle not a nucleus!");
 
-      std::cout << "NuclearInteraction: GetCrossSection: called with: beamId= " << BeamId
-                << " TargetId= " << TargetId << " CoMenergy= " << CoMenergy / 1_GeV
-                << std::endl;
-
-      if (!corsika::particles::IsNucleus(BeamId)) {
-        // ask hadronic interaction model for hadron cross section
-        return fHadronicInteraction.GetCrossSection(BeamId, TargetId, CoMenergy);
-      }
+      auto const iBeam = p.GetNuclearA();
+      HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam;
+      std::cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= "
+                << iBeam << " TargetId= " << TargetId
+                << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV << std::endl;
 
       // use nuclib to calc. nuclear cross sections
       // TODO: for now assumes air with hard coded composition
       // extend to arbitrary mixtures, requires smarter initialization
-
-      const int iBeam = corsika::particles::GetNucleusA(BeamId);
+      // get nuclib projectile code: nucleon number
       if (iBeam > 56 || iBeam < 2) {
-        // std::cout<<"NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
-        // << std::endl;
+        std::cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
+                  << std::endl
+                  << "A=" << iBeam << std::endl;
         throw std::runtime_error(
             "NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
             "NUCLIB!");
       }
 
-      const double dEcm = CoMenergy / 1_GeV;
-      if (dEcm < 10.)
+      const double dElabNuc = LabEnergyPerNuc / 1_GeV;
+      // TODO: these limitations are still sibyll specific.
+      // available target nuclei depends on the hadronic interaction model and the
+      // initialization
+      if (dElabNuc < 10.)
         throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
 
       // TODO: these limitations are still sibyll specific.
@@ -99,12 +101,17 @@ namespace corsika::process::sibyll {
           throw std::runtime_error(
               "Sibyll target outside range. Only nuclei with A<18 are allowed.");
         std::cout << "NuclearInteraction: calling signuc.." << std::endl;
-        signuc_(iBeam, dEcm, sigProd);
+        std::cout << "WARNING: using hard coded cross section for Nucleus-Air with "
+                     "SIBYLL! (fix me!)"
+                  << std::endl;
+        // TODO: target id is not used because cross section is still hard coded and fixed
+        // to air.
+        signuc_(iBeam, dElabNuc, sigProd);
         std::cout << "cross section: " << sigProd << std::endl;
-        return std::make_tuple(sigProd * 1_mbarn, 0_mbarn, iTarget);
+        return std::make_tuple(sigProd * 1_mbarn, 0_mbarn);
       }
       return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
-                             std::numeric_limits<double>::infinity() * 1_mbarn, 0);
+                             std::numeric_limits<double>::infinity() * 1_mbarn);
     }
 
     template <typename Particle, typename Track>
@@ -121,11 +128,15 @@ namespace corsika::process::sibyll {
           RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
       const particles::Code corsikaBeamId = p.GetPID();
-
       if (!corsika::particles::IsNucleus(corsikaBeamId)) {
         // no nuclear interaction
         return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
       }
+      // check if target-style nucleus (enum)
+      if (corsikaBeamId != corsika::particles::Code::Nucleus)
+        throw std::runtime_error(
+            "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
+            "projectiles should use NuclearStackExtension!");
 
       // read from cross section code table
       const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
@@ -136,6 +147,9 @@ namespace corsika::process::sibyll {
 
       // total momentum and energy
       HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+      int const nuclA = p.GetNuclearA();
+      auto const ElabNuc = p.GetEnergy() / nuclA;
+
       MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
       pTotLab += p.GetMomentum();
       pTotLab += pTarget;
@@ -143,13 +157,19 @@ namespace corsika::process::sibyll {
       // calculate cm. energy
       const HEPEnergyType ECoM = sqrt(
           (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
-
+      auto const ECoMNN = sqrt(2. * ElabNuc * nucleon_mass);
       std::cout << "NuclearInteraction: LambdaInt: \n"
                 << " input energy: " << Elab / 1_GeV << endl
                 << " input energy CoM: " << ECoM / 1_GeV << endl
-                << " beam pid:" << corsikaBeamId << endl;
+                << " beam pid:" << corsikaBeamId << endl
+                << " beam A: " << nuclA << endl
+                << " input energy per nucleon: " << ElabNuc / 1_GeV << endl
+                << " input energy CoM per nucleon: " << ECoMNN / 1_GeV << endl;
+      //      throw std::runtime_error("stop here");
 
-      if (Elab >= 8.5_GeV && ECoM >= 10_GeV) {
+      // energy limits
+      // TODO: values depend on hadronic interaction model !! this is sibyll specific
+      if (ElabNuc >= 8.5_GeV && ECoMNN >= 10_GeV) {
 
         // get target from environment
         /*
@@ -164,7 +184,6 @@ namespace corsika::process::sibyll {
         // determine average interaction length
         // weighted sum
         int i = -1;
-        double avgTargetMassNumber = 0.;
         si::CrossSectionType weightedProdCrossSection = 0_mbarn;
         // get weights of components from environment/medium
         const auto w = mediumComposition.GetFractions();
@@ -173,24 +192,23 @@ namespace corsika::process::sibyll {
           i++;
           cout << "NuclearInteraction: get interaction length for target: " << targetId
                << endl;
-          // TODO: remove number of nucleons, avg target mass is available in environment
-          auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
-              GetCrossSection(corsikaBeamId, targetId, ECoM);
+          auto const [productionCrossSection, elaCrossSection] =
+              GetCrossSection(p, targetId);
           [[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
 
           std::cout << "NuclearInteraction: "
                     << "IntLength: nuclib return (mb): "
                     << productionCrossSection / 1_mbarn << std::endl;
           weightedProdCrossSection += w[i] * productionCrossSection;
-          avgTargetMassNumber += w[i] * numberOfNucleons;
         }
         cout << "NuclearInteraction: "
              << "IntLength: weighted CrossSection (mb): "
              << weightedProdCrossSection / 1_mbarn << endl;
 
         // calculate interaction length in medium
-        GrammageType const int_length =
-            avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
+        GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+                                        corsika::units::constants::u /
+                                        weightedProdCrossSection;
         std::cout << "NuclearInteraction: "
                   << "interaction length (g/cm2): "
                   << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
@@ -207,9 +225,6 @@ namespace corsika::process::sibyll {
       // this routine superimposes different nucleon-nucleon interactions
       // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
 
-      // this routine superimposes different nucleon-nucleon interactions
-      // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
-
       using namespace corsika::units;
       using namespace corsika::utl;
       using namespace corsika::units::si;
@@ -217,18 +232,29 @@ namespace corsika::process::sibyll {
       using std::cout;
       using std::endl;
 
-      const auto corsikaProjId = p.GetPID();
-      cout << "NuclearInteraction: DoInteraction: called with:" << corsikaProjId << endl
-           << "NuclearInteraction: projectile mass: "
-           << corsika::particles::GetMass(corsikaProjId) / 1_GeV << endl;
+      const auto ProjId = p.GetPID();
+      // TODO: calculate projectile mass in nuclearStackExtension
+      //      const auto ProjMass = p.GetMass();
+      cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
 
-      if (!IsNucleus(corsikaProjId)) {
+      if (!IsNucleus(ProjId)) {
         cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
         // this should not happen
         // throw std::runtime_error("Non nuclear projectile in NUCLIB!");
         return process::EProcessReturn::eOk;
       }
 
+      // check if target-style nucleus (enum)
+      if (ProjId != corsika::particles::Code::Nucleus)
+        throw std::runtime_error(
+            "NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
+            "should use NuclearStackExtension!");
+
+      auto const ProjMass =
+	p.GetNuclearZ() * corsika::particles::Proton::GetMass() +
+	(p.GetNuclearA() - p.GetNuclearZ()) * corsika::particles::Neutron::GetMass();
+      cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
+      
       fCount++;
 
       const CoordinateSystem& rootCS =
@@ -242,7 +268,7 @@ namespace corsika::process::sibyll {
       cout << "Interaction: time: " << tOrig << endl;
 
       // projectile nucleon number
-      const int kAProj = GetNucleusA(corsikaProjId);
+      const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
       if (kAProj > 56)
         throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
 
@@ -284,6 +310,13 @@ namespace corsika::process::sibyll {
       HEPEnergyType EcmNN = PtotNN4.GetNorm();
       cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
 
+      if (!fHadronicInteraction.ValidCoMEnergy(EcmNN)) {
+        cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
+                "interaction model!"
+             << endl;
+        throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
+      }
+
       // define boost to NUCLEON-NUCLEON frame
       COMBoost const boost(PprojNucLab, nucleon_mass);
       // boost projecticle
@@ -308,8 +341,9 @@ namespace corsika::process::sibyll {
       const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
       const auto& mediumComposition =
           currentNode->GetModelProperties().GetNuclearComposition();
-      cout << "get cross sections for target materials.." << endl;
+      cout << "get nucleon-nucleus cross sections for target materials.." << endl;
       // get cross sections for target materials
+      // using nucleon-target-nucleus cross section!!!
       /*
         Here we read the cross section from the interaction model again,
         should be passed from GetInteractionLength if possible
@@ -320,8 +354,9 @@ namespace corsika::process::sibyll {
       for (size_t i = 0; i < compVec.size(); ++i) {
         auto const targetId = compVec[i];
         cout << "target component: " << targetId << endl;
-        cout << "beam id: " << targetId << endl;
-        const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId, targetId, EcmNN);
+        cout << "beam id: " << beamId << endl;
+        const auto [sigProd, sigEla, nNuc] =
+            fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
         cross_section_of_components[i] = sigProd;
         [[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
         [[maybe_unused]] auto sigNucCopy = nNuc;   // ONLY TO AVOID COMPILER WARNINGS
@@ -396,97 +431,90 @@ namespace corsika::process::sibyll {
              << " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1]
              << " pz=" << fragments_.ppp[j][2] << endl;
 
-      auto Nucleus = [](int iA) {
-        //		       int znew = iA / 2.15 + 0.7;
-        corsika::particles::Code pCode;
-        switch (iA) {
-          case 2:
-            pCode = corsika::particles::Deuterium::GetCode();
-            break;
-
-          case 3:
-            pCode = corsika::particles::Tritium::GetCode();
-            break;
-
-          case 4:
-            pCode = corsika::particles::Helium::GetCode();
-            break;
-
-          case 12:
-            pCode = corsika::particles::Carbon::GetCode();
-            break;
-
-          case 14:
-            pCode = corsika::particles::Nitrogen::GetCode();
-            break;
-
-          case 16:
-            pCode = corsika::particles::Oxygen::GetCode();
-            break;
-
-          case 33:
-            pCode = corsika::particles::Sulphur::GetCode();
-            break;
-
-          case 40:
-            pCode = corsika::particles::Argon::GetCode();
-            break;
-
-          default:
-            pCode = corsika::particles::Proton::GetCode();
-        }
-        return pCode;
-      };
-
+      cout << "adding nuclear fragments to particle stack.." << endl;
       // put nuclear fragments on corsika stack
       for (int j = 0; j < nFragments; ++j) {
-        // here we need the nucleonNumber to corsika Id conversion
-        auto pCode = Nucleus(AFragments[j]);
+        corsika::particles::Code specCode;
+        const auto nuclA = AFragments[j];
+        // get Z from stability line
+        const auto nuclZ = int(nuclA / 2.15 + 0.7);
+
+        // TODO: do we need to catch single nucleons??
+        if (nuclA == 1)
+          // TODO: sample neutron or proton
+          specCode = corsika::particles::Code::Proton;
+        else
+          specCode = corsika::particles::Code::Nucleus;
+
+        // TODO: mass of nuclei?
+        const HEPMassType mass =
+            corsika::particles::Proton::GetMass() * nuclZ +
+            (nuclA - nuclZ) *
+                corsika::particles::Neutron::GetMass(); // this neglects binding energy
+
+        cout << "NuclearInteraction: adding fragment: " << specCode << endl;
+        cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << endl;
+        cout << "NuclearInteraction: mass: " << mass / 1_GeV << endl;
 
         // CORSIKA 7 way
         // spectators inherit momentum from original projectile
-        const double mass_ratio = corsika::particles::GetMass(pCode) /
-                                  corsika::particles::GetMass(corsikaProjId);
+        const double mass_ratio = mass / ProjMass;
+
+        cout << "NuclearInteraction: mass ratio " << mass_ratio << endl;
+
         auto const Plab = PprojLab * mass_ratio;
 
-        p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
-                       pOrig, tOrig);
+        cout << "NuclearInteraction: fragment momentum: "
+             << Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+
+        if (nuclA == 1)
+          // add nucleon
+          p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
+                         Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+        else
+          // add nucleus
+          p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
+                         Plab.GetSpaceLikeComponents(), pOrig, tOrig, nuclA, nuclZ);
       }
 
       // add elastic nucleons to corsika stack
       // TODO: the elastic interaction could be external like the inelastic interaction,
       // e.g. use existing ElasticModel
+      cout << "adding elastically scattered nucleons to particle stack.." << endl;
       for (int j = 0; j < nElasticNucleons; ++j) {
         // TODO: sample proton or neutron
-        auto pCode = corsika::particles::Proton::GetCode();
+        auto const elaNucCode = corsika::particles::Code::Proton;
 
         // CORSIKA 7 way
         // elastic nucleons inherit momentum from original projectile
         // neglecting momentum transfer in interaction
-        const double mass_ratio = corsika::particles::GetMass(pCode) /
-                                  corsika::particles::GetMass(corsikaProjId);
+        const double mass_ratio = corsika::particles::GetMass(elaNucCode) / ProjMass;
         auto const Plab = PprojLab * mass_ratio;
 
-        p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
-                       pOrig, tOrig);
+        p.AddSecondary(elaNucCode, Plab.GetTimeLikeComponent(),
+                       Plab.GetSpaceLikeComponents(), pOrig, tOrig);
       }
 
       // add inelastic interactions
+      cout << "calculate inelastic nucleon-nucleon interactions.." << endl;
       for (int j = 0; j < nInelNucleons; ++j) {
-
         // TODO: sample neutron or proton
         auto pCode = corsika::particles::Proton::GetCode();
         // temporarily add to stack, will be removed after interaction in DoInteraction
+        cout << "inelastic interaction no. " << j << endl;
         auto inelasticNucleon =
             p.AddSecondary(pCode, PprojNucLab.GetTimeLikeComponent(),
                            PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig);
         // create inelastic interaction
+        cout << "calling HadronicInteraction..." << endl;
         fHadronicInteraction.DoInteraction(inelasticNucleon, s);
       }
 
       // delete parent particle
       p.Delete();
 
+      cout << "NuclearInteraction: DoInteraction: done" << endl;
+
       return process::EProcessReturn::eOk;
     }
 

Processes/Sibyll/SibStack.h

@@ -27,6 +27,7 @@ namespace corsika::process::sibyll {
 
   public:
     void Init();
+    void Dump() const {}
 
     void Clear() { s_plist_.np = 0; }
     unsigned int GetSize() const { return s_plist_.np; }
@@ -65,8 +66,7 @@ namespace corsika::process::sibyll {
           RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
       QuantityVector<hepmomentum_d> components = {
           s_plist_.p[0][i] * 1_GeV, s_plist_.p[1][i] * 1_GeV, s_plist_.p[2][i] * 1_GeV};
-      MomentumVector v1(rootCS, components);
-      return v1;
+      return MomentumVector(rootCS, components);
     }
 
     void Copy(const unsigned int i1, const unsigned int i2) {

Processes/Sibyll/testSibyll.cc

@@ -136,13 +136,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
   SECTION("NuclearInteractionInterface") {
 
     setup::Stack stack;
-    const HEPEnergyType E0 = 100_GeV;
+    const HEPEnergyType E0 = 400_GeV;
     HEPMomentumType P0 =
         sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
     auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
     geometry::Point pos(cs, 0_m, 0_m, 0_m);
 
-    auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
+    auto particle = stack.AddParticle(particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2);
 
     Interaction hmodel(env);
     NuclearInteraction model(env, hmodel);

Processes/StackInspector/StackInspector.cc

@@ -51,7 +51,10 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
     auto pos = iterP.GetPosition().GetCoordinates(rootCS);
     cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
          << iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
-         << " pos=" << pos << endl;
+         << " pos=" << pos;
+    // if (iterP.GetPID()==Code::Nucleus)
+    // cout << " nuc_ref=" << iterP.GetNucleusRef();
+    cout << endl;
   }
   fCountStep++;
   cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()

Stack/NuclearStackExtension/NuclearStackExtension.h

@@ -24,11 +24,28 @@
 
 namespace corsika::stack {
 
+  /**
+   * @namespace nuclear_extension
+   *
+   * Add A and Z data to existing stack of particle properties.
+   *
+   * Only for Code::Nucleus particles A and Z are stored, not for all
+   * normal elementary particles.
+   *
+   * Thus in your code, make sure to always check <code>
+   * particle.GetPID()==Code::Nucleus </code> before attempting to
+   * read any nuclear information.
+   *
+   *
+   */
+
   namespace nuclear_extension {
 
     typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
 
     /**
+     * @class NuclearParticleInterface
+     *
      * Define ParticleInterface for NuclearStackExtension Stack derived from
      * ParticleInterface of Inner stack class
      */
@@ -58,12 +75,12 @@ namespace corsika::stack {
             err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
             throw std::runtime_error(err.str());
           }
-          SetNuclearRef(
+          SetNucleusRef(
               GetStackData().GetNucleusNextRef()); // store this nucleus data ref
           SetNuclearA(vA);
           SetNuclearZ(vZ);
         } else {
-          SetNuclearRef(-1); // this is not a nucleus
+          SetNucleusRef(-1); // this is not a nucleus
         }
         InnerParticleInterface<StackIteratorInterface>::SetParticleData(
             vDataPID, vDataE, vMomentum, vPosition, vTime);
@@ -99,28 +116,40 @@ namespace corsika::stack {
       int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); }
       /// @}
 
+      // int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
+
     protected:
-      void SetNuclearRef(const int vR) { GetStackData().SetNuclearRef(GetIndex(), vR); }
+      void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
     };
 
     /**
-     * Memory implementation of the most simple (stupid) particle stack object.
+     * @class NuclearStackExtension
+     *
+     * Memory implementation of adding nuclear inforamtion to the
+     * existing particle stack defined in class InnerStackImpl.
+     *
+     * Inside the NuclearStackExtension class there is a dedicated
+     * fNucleusRef index, where fNucleusRef[i] is referring to the
+     * correct A and Z for a specific particle index i. fNucleusRef[i]
+     * == -1 means that this is not a nucleus, and a subsequent call to
+     * GetNucleusA would produce an exception.
      */
     template <typename InnerStackImpl>
     class NuclearStackExtensionImpl : public InnerStackImpl {
 
     public:
       void Init() { InnerStackImpl::Init(); }
+      void Dump() { InnerStackImpl::Dump(); }
 
       void Clear() {
         InnerStackImpl::Clear();
-        fNuclearRef.clear();
+        fNucleusRef.clear();
         fNuclearA.clear();
         fNuclearZ.clear();
       }
 
-      unsigned int GetSize() const { return fNuclearRef.size(); }
-      unsigned int GetCapacity() const { return fNuclearRef.size(); }
+      unsigned int GetSize() const { return fNucleusRef.size(); }
+      unsigned int GetCapacity() const { return fNucleusRef.capacity(); }
 
       void SetNuclearA(const unsigned int i, const unsigned short vA) {
         fNuclearA[GetNucleusRef(i)] = vA;
@@ -128,7 +157,7 @@ namespace corsika::stack {
       void SetNuclearZ(const unsigned int i, const unsigned short vZ) {
         fNuclearZ[GetNucleusRef(i)] = vZ;
       }
-      void SetNuclearRef(const unsigned int i, const int v) { fNuclearRef[i] = v; }
+      void SetNucleusRef(const unsigned int i, const int v) { fNucleusRef[i] = v; }
 
       int GetNuclearA(const unsigned int i) const { return fNuclearA[GetNucleusRef(i)]; }
       int GetNuclearZ(const unsigned int i) const { return fNuclearZ[GetNucleusRef(i)]; }
@@ -141,7 +170,7 @@ namespace corsika::stack {
       }
 
       int GetNucleusRef(const unsigned int i) const {
-        if (fNuclearRef[i] >= 0) return fNuclearRef[i];
+        if (fNucleusRef[i] >= 0) return fNucleusRef[i];
         std::ostringstream err;
         err << "NuclearStackExtension: no nucleus at ref=" << i;
         throw std::runtime_error(err.str());
@@ -151,35 +180,41 @@ namespace corsika::stack {
        *   Function to copy particle at location i1 in stack to i2
        */
       void Copy(const unsigned int i1, const unsigned int i2) {
+        // index range check
         if (i1 >= GetSize() || i2 >= GetSize()) {
           std::ostringstream err;
           err << "NuclearStackExtension: trying to access data beyond size of stack!";
           throw std::runtime_error(err.str());
         }
+        // copy internal particle data p[i2] = p[i1]
         InnerStackImpl::Copy(i1, i2);
-        const int ref1 = fNuclearRef[i1];
-        const int ref2 = fNuclearRef[i2];
+        // check if any of p[i1] or p[i2] was a Code::Nucleus
+        const int ref1 = fNucleusRef[i1];
+        const int ref2 = fNucleusRef[i2];
         if (ref2 < 0) {
           if (ref1 >= 0) {
             // i1 is nucleus, i2 is not
-            fNuclearRef[i2] = GetNucleusNextRef();
-            fNuclearA[fNuclearRef[i2]] = fNuclearA[ref1];
-            fNuclearZ[fNuclearRef[i2]] = fNuclearZ[ref1];
+            fNucleusRef[i2] = GetNucleusNextRef();
+            fNuclearA[fNucleusRef[i2]] = fNuclearA[ref1];
+            fNuclearZ[fNucleusRef[i2]] = fNuclearZ[ref1];
           } else {
             // neither i1 nor i2 are nuclei
           }
         } else {
           if (ref1 >= 0) {
             // both are nuclei, i2 is overwritten with nucleus i1
+            // fNucleusRef stays the same, but A and Z data is overwritten
             fNuclearA[ref2] = fNuclearA[ref1];
             fNuclearZ[ref2] = fNuclearZ[ref1];
           } else {
             // i2 is overwritten with non-nucleus i1
-            fNuclearA.erase(fNuclearA.cbegin() + ref2);
-            fNuclearZ.erase(fNuclearZ.cbegin() + ref2);
-            const int n = fNuclearRef.size();
+            fNucleusRef[i2] = -1;                       // flag as non-nucleus
+            fNuclearA.erase(fNuclearA.cbegin() + ref2); // remove data for i2
+            fNuclearZ.erase(fNuclearZ.cbegin() + ref2); // remove data for i2
+            const int n = fNucleusRef.size(); // update fNucleusRef: indices above ref2
+                                              // must be decremented by 1
             for (int i = 0; i < n; ++i) {
-              if (fNuclearRef[i] >= ref2) { fNuclearRef[i] -= 1; }
+              if (fNucleusRef[i] > ref2) { fNucleusRef[i] -= 1; }
             }
           }
         }
@@ -189,48 +224,34 @@ namespace corsika::stack {
        *   Function to copy particle at location i2 in stack to i1
        */
       void Swap(const unsigned int i1, const unsigned int i2) {
+        // index range check
         if (i1 >= GetSize() || i2 >= GetSize()) {
           std::ostringstream err;
           err << "NuclearStackExtension: trying to access data beyond size of stack!";
           throw std::runtime_error(err.str());
         }
+        // swap original particle data
         InnerStackImpl::Swap(i1, i2);
-        const int ref1 = fNuclearRef[i1];
-        const int ref2 = fNuclearRef[i2];
-        if (ref2 < 0) {
-          if (ref1 >= 0) {
-            // i1 is nucleus, i2 is not
-            std::swap(fNuclearRef[i2], fNuclearRef[i1]);
-          } else {
-            // neither i1 nor i2 are nuclei
-          }
-        } else {
-          if (ref1 >= 0) {
-            // both are nuclei, i2 is overwritten with nucleus i1
-            std::swap(fNuclearRef[i2], fNuclearRef[i1]);
-          } else {
-            // i2 is overwritten with non-nucleus i1
-            std::swap(fNuclearRef[i2], fNuclearRef[i1]);
-          }
-        }
+        // swap corresponding nuclear reference data
+        std::swap(fNucleusRef[i2], fNucleusRef[i1]);
       }
 
       void IncrementSize() {
         InnerStackImpl::IncrementSize();
-        fNuclearRef.push_back(-1);
+        fNucleusRef.push_back(-1);
       }
 
       void DecrementSize() {
         InnerStackImpl::DecrementSize();
-        if (fNuclearRef.size() > 0) {
-          const int ref = fNuclearRef.back();
-          fNuclearRef.pop_back();
+        if (fNucleusRef.size() > 0) {
+          const int ref = fNucleusRef.back();
+          fNucleusRef.pop_back();
           if (ref >= 0) {
             fNuclearA.erase(fNuclearA.begin() + ref);
             fNuclearZ.erase(fNuclearZ.begin() + ref);
-            const int n = fNuclearRef.size();
+            const int n = fNucleusRef.size();
             for (int i = 0; i < n; ++i) {
-              if (fNuclearRef[i] >= ref) { fNuclearRef[i] -= 1; }
+              if (fNucleusRef[i] >= ref) { fNucleusRef[i] -= 1; }
             }
           }
         }
@@ -239,7 +260,7 @@ namespace corsika::stack {
     private:
       /// the actual memory to store particle data
 
-      std::vector<int> fNuclearRef;
+      std::vector<int> fNucleusRef;
       std::vector<unsigned short> fNuclearA;
       std::vector<unsigned short> fNuclearZ;
 

Stack/SuperStupidStack/SuperStupidStack.h

@@ -114,6 +114,7 @@ namespace corsika::stack {
 
     public:
       void Init() {}
+      void Dump() const {}
 
       void Clear() {
         fDataPID.clear();