fixed bug in SetStable for antiparticles, repaired interference in decay...

fixed bug in SetStable for antiparticles, repaired interference in decay configuration between interaction and decay

Documentation/Examples/boundary_example.cc

@@ -123,9 +123,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);
 

Documentation/Examples/cascade_example.cc

@@ -100,12 +100,12 @@ int main() {
 
   const std::vector<particles::Code> trackedHadrons = {
       particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
+      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short,      particles::Code::MuMinus,particles::Code::MuPlus};
 
   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");

Documentation/Examples/vertical_EAS.cc

@@ -169,13 +169,9 @@ int main() {
   // setup processes, decays and interactions
   tracking_line::TrackingLine tracking;
 
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
   process::sibyll::Interaction sibyll;
   process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
+  process::sibyll::Decay decay;
 
   //~ process::pythia::Decay decay(trackedHadrons);
   process::particle_cut::ParticleCut cut(2_GeV);

Processes/Sibyll/Decay.cc

@@ -29,77 +29,62 @@ 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::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::SetStable(const vector<particles::Code> vParticleList) {
+    for (auto p : vParticleList) Decay::SetStable(p);
   }
 
-  void Decay::SetUnstable(const particles::Code pCode) {
-    int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+  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::SetDecay(const particles::Code vCode, const bool vMakeUnstable) {
+    vMakeUnstable ? SetUnstable(vCode) : SetStable(vCode);
+  }
+
+  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 <>
@@ -149,12 +134,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

@@ -22,23 +22,33 @@ 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&);
     };

Processes/Sibyll/Interaction.cc

@@ -48,46 +48,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,
@@ -326,8 +318,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

@@ -32,9 +32,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 {
@@ -66,7 +70,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/testSibyll.cc

@@ -72,6 +72,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;
@@ -164,14 +165,29 @@ TEST_CASE("SibyllInterface", "[processes]") {
     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);
     [[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]);
+    }
+  }
 }