Merge branch '188-sibyll-decay-of-sigmabar0-yields-sigmabar0' into 'master'

Resolve "SIBYLL: decay of Sigmabar0 yields Sigmabar0"

Closes #188

See merge request AirShowerPhysics/corsika!129

Documentation/Examples/boundary_example.cc

@@ -121,9 +121,7 @@ int main() {
 
   random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
   process::sibyll::Interaction sibyll;
-  process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
-                                particles::Code::KPlus, particles::Code::KMinus,
-                                particles::Code::K0Long, particles::Code::K0Short}};
+  process::sibyll::Decay decay;
 
   process::particle_cut::ParticleCut cut(20_GeV);
 
@@ -180,4 +178,8 @@ int main() {
       cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
   cout << "total energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+
+  // basic check for unit-tests
+  assert(cut.GetNumberEmParticles() == 29785);
+  assert(cut.GetNumberInvParticles() == 26697);
 }

Documentation/Examples/cascade_example.cc

@@ -130,15 +130,11 @@ int main() {
   tracking_line::TrackingLine tracking;
   stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
 
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
   random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
   random::RNGManager::GetInstance().RegisterRandomStream("pythia");
   process::sibyll::Interaction sibyll;
   process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
+  process::sibyll::Decay decay;
   process::particle_cut::ParticleCut cut(20_GeV);
 
   process::track_writer::TrackWriter trackWriter("tracks.dat");
@@ -162,4 +158,8 @@ int main() {
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
   cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
        << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+
+  // basic check for unit-tests
+  assert(cut.GetNumberEmParticles() == 127);
+  assert(cut.GetNumberInvParticles() == 116);
 }

Documentation/Examples/vertical_EAS.cc

@@ -125,14 +125,10 @@ int main() {
   tracking_line::TrackingLine tracking;
   stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
 
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
   random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
   process::sibyll::Interaction sibyll;
   process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
+  process::sibyll::Decay decay;
   process::particle_cut::ParticleCut cut(20_GeV);
 
   process::track_writer::TrackWriter trackWriter("tracks.dat");
@@ -155,4 +151,8 @@ int main() {
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
   cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
        << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+
+  // basic check for unit-tests
+  assert(cut.GetNumberEmParticles() == 526);
+  assert(cut.GetNumberInvParticles() == 645);
 }

Framework/Cascade/Cascade.h

@@ -263,6 +263,10 @@ namespace corsika::cascade {
             InverseTimeType inv_decay_count = 0 / second;
             fProcessSequence.SelectDecay(vParticle, projectile, sample_process,
                                          inv_decay_count);
+            // make sure particle actually did decay if it should have done so
+            if (secondaries.GetSize() == 1 &&
+                projectile.GetPID() == secondaries.GetNextParticle().GetPID())
+              throw std::runtime_error("Cascade::Step: Particle decays into itself!");
           }
 
           fProcessSequence.DoSecondaries(secondaries);
@@ -273,6 +277,7 @@ namespace corsika::cascade {
                               // be "protected" and not accessible to physics
 
         } else { // step-length limitation within volume
+
           std::cout << "step-length limitation" << std::endl;
           fProcessSequence.DoSecondaries(secondaries);
         }

Processes/ParticleCut/ParticleCut.h

@@ -24,9 +24,9 @@ namespace corsika::process {
 
       units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
       units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
-      int fEmCount = 0;
+      unsigned int fEmCount = 0;
       units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
-      int fInvCount = 0;
+      unsigned int fInvCount = 0;
 
     public:
       ParticleCut(const units::si::HEPEnergyType vCut)
@@ -46,6 +46,8 @@ namespace corsika::process {
       units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
       units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
       units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+      unsigned int GetNumberEmParticles() const { return fEmCount; }
+      unsigned int GetNumberInvParticles() const { return fInvCount; }
     };
   } // namespace particle_cut
 } // namespace corsika::process

Processes/Sibyll/Decay.cc

@@ -23,82 +23,69 @@ using std::vector;
 
 using namespace corsika;
 using namespace corsika::setup;
+
+using SetupView = corsika::setup::StackView;
 using SetupProjectile = corsika::setup::StackView::ParticleType;
 using SetupParticle = corsika::setup::Stack::ParticleType;
 
 namespace corsika::process::sibyll {
 
-  Decay::Decay(vector<particles::Code> pParticles)
-      : fTrackedParticles(pParticles) {}
+  Decay::Decay() {}
   Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
   void Decay::Init() {
-    SetHadronsUnstable();
-    SetParticleListStable(fTrackedParticles);
+    // switch off decays to avoid internal decay chains
+    SetAllStable();
+  }
+
+  void Decay::SetStable(const vector<particles::Code> vParticleList) {
+    for (auto p : vParticleList) Decay::SetStable(p);
+  }
+
+  void Decay::SetUnstable(const vector<particles::Code> vParticleList) {
+    for (auto p : vParticleList) Decay::SetUnstable(p);
+  }
+
+  bool Decay::IsStable(const particles::Code vCode) {
+    return abs(process::sibyll::ConvertToSibyllRaw(vCode)) <= 0 ? true : false;
+  }
+
+  bool Decay::IsUnstable(const particles::Code vCode) {
+    return abs(process::sibyll::ConvertToSibyllRaw(vCode)) > 0 ? true : false;
   }
 
-  void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
-    /*
-       Sibyll is hadronic generator
-       only hadrons decay
-     */
-    // set particles unstable
-    SetHadronsUnstable();
-    cout << "Interaction: setting tracked hadrons stable.." << endl;
-    for (auto p : particleList) Decay::SetStable(p);
+  void Decay::SetDecay(const particles::Code vCode, const bool vMakeUnstable) {
+    vMakeUnstable ? SetUnstable(vCode) : SetStable(vCode);
   }
 
-  void Decay::SetUnstable(const particles::Code pCode) {
-    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+  void Decay::SetUnstable(const particles::Code vCode) {
+    cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
+    const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
   }
 
-  void Decay::SetStable(const particles::Code pCode) {
-    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+  void Decay::SetStable(const particles::Code vCode) {
+    cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
+    const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
   }
 
   void Decay::SetAllStable() {
-    // name? also makes EM particles stable
-
-    cout << "Decay: setting all particles stable.." << endl;
-
-    // loop over all particles in sibyll
-    // should be changed to loop over human readable list
-    // i.e. particles::ListOfParticles()
-    for (auto& p : corsika2sibyll) {
-      // cout << (int)p << endl;
-      const int sibCode = static_cast<int>(p);
-      // skip unknown and antiparticles
-      if (sibCode < 1) continue;
-      s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
-    }
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
   }
 
-  void Decay::SetHadronsUnstable() {
-
-    // name? also makes EM particles stable
+  void Decay::SetAllUnstable() {
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+  }
 
-    // loop over all particles in sibyll
-    // should be changed to loop over human readable list
-    // i.e. particles::ListOfParticles()
-    cout << "Sibyll: setting hadrons unstable.." << endl;
-    // make ALL particles unstable, then set EM stable
-    for (int sibCode : corsika2sibyll) {
-      if (sibCode < 1) continue;
-      s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
-    }
-    // set Leptons and Proton and Neutron stable
-    // use stack to loop over particles
-    constexpr particles::Code particleList[] = {
-        particles::Code::Proton,   particles::Code::Neutron, particles::Code::Electron,
-        particles::Code::Positron, particles::Code::NuE,     particles::Code::NuEBar,
-        particles::Code::MuMinus,  particles::Code::MuPlus,  particles::Code::NuMu,
-        particles::Code::NuMuBar};
-
-    for (auto p : particleList) {
-      const int sibid = process::sibyll::ConvertToSibyllRaw(p);
-      s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
-    }
+  void Decay::PrintDecayConfig(const particles::Code vCode) {
+    cout << "Decay: Sibyll decay configuration:" << endl;
+    const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
+    const int absSibCode = abs(sibCode);
+    cout << vCode << " is ";
+    if (s_csydec_.idb[absSibCode - 1] <= 0)
+      cout << "stable" << endl;
+    else
+      cout << "unstable" << endl;
   }
 
   template <>
@@ -148,12 +135,16 @@ namespace corsika::process::sibyll {
     // remember position
     Point const decayPoint = vP.GetPosition();
     TimeType const t0 = vP.GetTime();
-    // set all particles/hadrons unstable
-    // setHadronsUnstable();
+    // remember if particles is unstable
+    // auto const priorIsUnstable = IsUnstable(pCode);
+    // switch on decay for this particle
     SetUnstable(pCode);
+    PrintDecayConfig(pCode);
+
     // call sibyll decay
     cout << "Decay: calling Sibyll decay routine.." << endl;
     decsib_();
+
     // reset to stable
     SetStable(pCode);
     // print output

Processes/Sibyll/Decay.h

@@ -13,6 +13,7 @@
 
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/process/DecayProcess.h>
+#include <corsika/process/SecondariesProcess.h>
 
 #include <vector>
 
@@ -21,27 +22,42 @@ namespace corsika::process {
   namespace sibyll {
 
     class Decay : public corsika::process::DecayProcess<Decay> {
-      std::vector<particles::Code> const fTrackedParticles;
       int fCount = 0;
 
     public:
-      Decay(std::vector<particles::Code>);
+      Decay();
       ~Decay();
       void Init();
 
-      void SetParticleListStable(const std::vector<particles::Code>);
+      void SetStable(const std::vector<particles::Code>);
+      void SetUnstable(const std::vector<particles::Code>);
       void SetUnstable(const corsika::particles::Code);
       void SetStable(const corsika::particles::Code);
+      void SetAllUnstable();
       void SetAllStable();
+      bool IsStable(const corsika::particles::Code);
+      bool IsUnstable(const corsika::particles::Code);
+      void SetDecay(const particles::Code, const bool);
+
+      void PrintDecayConfig(const corsika::particles::Code);
       void SetHadronsUnstable();
 
       template <typename TParticle>
       corsika::units::si::TimeType GetLifetime(TParticle const&) const;
 
+      /**
+       In this function SIBYLL is called to produce to decay the input particle.
+     */
+
       template <typename TProjectile>
       void DoDecay(TProjectile&);
+
+      template <typename TParticleView>
+      EProcessReturn DoSecondaries(TParticleView&);
     };
+
   } // namespace sibyll
+
 } // namespace corsika::process
 
 #endif

Processes/Sibyll/Interaction.cc

@@ -47,46 +47,38 @@ namespace corsika::process::sibyll {
     // initialize Sibyll
     if (!fInitialized) {
       sibyll_ini_();
-
-      // any decays in sibyll? if yes need to define which particles
-      if (fInternalDecays) {
-        // define which particles are passed to corsika, i.e. which particles make it into
-        // history even very shortlived particles like charm or pi0 are of interest here
-        const std::vector<particles::Code> hadronsWeWantTrackedByCorsika = {
-            particles::Code::PiPlus,     particles::Code::PiMinus,
-            particles::Code::Pi0,        particles::Code::KMinus,
-            particles::Code::KPlus,      particles::Code::K0Long,
-            particles::Code::K0Short,    particles::Code::SigmaPlus,
-            particles::Code::SigmaMinus, particles::Code::Lambda0,
-            particles::Code::Xi0,        particles::Code::XiMinus,
-            particles::Code::OmegaMinus, particles::Code::DPlus,
-            particles::Code::DMinus,     particles::Code::D0,
-            particles::Code::D0Bar};
-
-        Interaction::SetParticleListStable(hadronsWeWantTrackedByCorsika);
-      }
-
       fInitialized = true;
     }
   }
 
-  void Interaction::SetParticleListStable(
-      std::vector<particles::Code> const& particleList) {
-    for (auto p : particleList) Interaction::SetStable(p);
+  void Interaction::SetStable(std::vector<particles::Code> const& vParticleList) {
+    for (auto p : vParticleList) Interaction::SetStable(p);
+  }
+
+  void Interaction::SetUnstable(std::vector<particles::Code> const& vParticleList) {
+    for (auto p : vParticleList) Interaction::SetUnstable(p);
   }
 
-  void Interaction::SetUnstable(const particles::Code pCode) {
-    cout << "Sibyll::Interaction: setting " << pCode << " unstable.." << endl;
-    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+  void Interaction::SetUnstable(const particles::Code vCode) {
+    cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
+    const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
   }
 
-  void Interaction::SetStable(const particles::Code pCode) {
-    cout << "Sibyll::Interaction: setting " << pCode << " stable.." << endl;
-    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+  void Interaction::SetStable(const particles::Code vCode) {
+    cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
+    const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
     s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
   }
 
+  void Interaction::SetAllStable() {
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
+  }
+
+  void Interaction::SetAllUnstable() {
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+  }
+
   tuple<units::si::CrossSectionType, units::si::CrossSectionType>
   Interaction::GetCrossSection(const particles::Code BeamId,
                                const particles::Code TargetId,
@@ -325,8 +317,16 @@ namespace corsika::process::sibyll {
         const double sqs = Ecm / 1_GeV;
         // running sibyll, filling stack
         sibyll_(kBeam, targetSibCode, sqs);
-        // running decays
-        if (fInternalDecays) decsib_();
+        if (fInternalDecays) {
+          // particles that decay internally will never appear on the corsika stack
+          // switch on all decays except for the particles we want to take part in the
+          // tracking
+          SetAllUnstable();
+          SetStable(fTrackedParticles);
+          decsib_();
+          // reset
+          SetAllStable();
+        }
         // print final state
         int print_unit = 6;
         sib_list_(print_unit);

Processes/Sibyll/Interaction.h

@@ -31,9 +31,13 @@ namespace corsika::process::sibyll {
 
     void Init();
 
-    void SetParticleListStable(std::vector<particles::Code> const&);
+    void SetStable(std::vector<particles::Code> const&);
+    void SetUnstable(std::vector<particles::Code> const&);
+
     void SetUnstable(const corsika::particles::Code);
     void SetStable(const corsika::particles::Code);
+    void SetAllUnstable();
+    void SetAllStable();
 
     bool WasInitialized() { return fInitialized; }
     bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
@@ -65,7 +69,18 @@ namespace corsika::process::sibyll {
   private:
     corsika::random::RNG& fRNG =
         corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-
+    // FOR NOW keep trackedParticles private, could be configurable
+    std::vector<particles::Code> const fTrackedParticles = {
+        particles::Code::PiPlus,     particles::Code::PiMinus,
+        particles::Code::Pi0,        particles::Code::KMinus,
+        particles::Code::KPlus,      particles::Code::K0Long,
+        particles::Code::K0Short,    particles::Code::SigmaPlus,
+        particles::Code::SigmaMinus, particles::Code::Lambda0,
+        particles::Code::Xi0,        particles::Code::XiMinus,
+        particles::Code::OmegaMinus, particles::Code::DPlus,
+        particles::Code::DMinus,     particles::Code::D0,
+        particles::Code::MuMinus,    particles::Code::MuPlus,
+        particles::Code::D0Bar};
     const bool fInternalDecays = true;
     const corsika::units::si::HEPEnergyType fMinEnergyCoM =
         10. * 1e9 * corsika::units::si::electronvolt;

Processes/Sibyll/sibyll2.3c.f

@@ -478,7 +478,7 @@ C-----------------------------------------------------------------------
       CALL BLOCK_INI
       CALL NUC_GEOM_INI
       CALL SIG_AIR_INI
-      CALL DEC_INI
+c      CALL DEC_INI
 c...  charm frag. normalisation
       CALL ZNORMAL

Processes/Sibyll/sibyll_codes.dat

@@ -17,6 +17,8 @@ Pi0 6 1 2
 # rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
 Rho0 27 0 0
 K0Long 11 1 3
+K0 21 0 3
+K0Bar 22 0 3
 PiPlus 7 1 2
 PiMinus 8 1 2
 RhoPlus 25 0 0

Processes/Sibyll/testSibyll.cc

@@ -69,6 +69,7 @@ TEST_CASE("Sibyll", "[processes]") {
 #include <corsika/environment/Environment.h>
 #include <corsika/environment/HomogeneousMedium.h>
 #include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
 
 using namespace corsika::units::si;
 using namespace corsika::units;
@@ -157,18 +158,34 @@ TEST_CASE("SibyllInterface", "[processes]") {
     auto particle = stack.AddParticle(
         std::tuple<particles::Code, units::si::HEPEnergyType,
                    corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::Proton, E0, plab, pos, 0_ns});
+            particles::Code::Lambda0, E0, plab, pos, 0_ns});
     corsika::stack::SecondaryView view(particle);
     auto projectile = view.GetProjectile();
 
-    const std::vector<particles::Code> particleList = {
-        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-        particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-    Decay model(particleList);
+    Decay model;
 
     model.Init();
     /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
+    // run checks
     [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
   }
+
+  SECTION("DecayConfiguration") {
+
+    Decay model;
+
+    const std::vector<particles::Code> particleTestList = {
+        particles::Code::PiPlus,     particles::Code::PiMinus, particles::Code::KPlus,
+        particles::Code::Lambda0Bar, particles::Code::NuE,     particles::Code::D0Bar};
+
+    for (auto& pCode : particleTestList) {
+      model.SetUnstable(pCode);
+      // get state of sibyll internal config
+      REQUIRE(0 <= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
+
+      model.SetStable(pCode);
+      // get state of sibyll internal config
+      REQUIRE(0 >= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
+    }
+  }
 }