diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
index 423c4509e..7a28d631a 100644
--- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
+++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
@@ -10,6 +10,7 @@
 
 using namespace Pythia8;
 
+#include <cmath>
 namespace hf_generators
 {
     enum GenType : int {
@@ -35,7 +36,13 @@ public:
     //}
 
     ///  Destructor
-    ~GeneratorPythia8EmbedHF() = default;
+    ~GeneratorPythia8EmbedHF() {
+        // Clean up the internally created HF generator if any
+        if (mGeneratorEvHF) {
+            delete mGeneratorEvHF;
+            mGeneratorEvHF = nullptr;
+        }
+    }
 
     /// Init
     bool Init() override
@@ -51,29 +58,30 @@ public:
     /// \param yHadronMin minimum hadron rapidity
     /// \param yHadronMax maximum hadron rapidity
     /// \param hadronPdgList list of PDG codes for hadrons to be used in trigger
-    void setupGeneratorEvHF(int genType, bool usePtHardBins, float yQuarkMin, float yQuarkMax, float yHadronMin, float yHadronMax, std::vector<int> hadronPdgList = {}) {
+    /// \param quarkPdgList list of PDG codes for quarks to be enriched in the trigger
+    void setupGeneratorEvHF(int genType, bool usePtHardBins, float yQuarkMin, float yQuarkMax, float yHadronMin, float yHadronMax, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int,2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {}) {
         mGeneratorEvHF = nullptr;
         switch (genType)
         {
         case hf_generators::GapTriggeredCharm:
             LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapTriggeredCharm **********";
             LOG(info) << "**********                           Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
-            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredCharm(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
+            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredCharm(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList, partPdgToReplaceList, freqReplaceList));
             break;
         case hf_generators::GapTriggeredBeauty:
             LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapTriggeredBeauty **********";
             LOG(info) << "**********                           Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
-            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
+            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList, partPdgToReplaceList, freqReplaceList));
             break;
         case hf_generators::GapTriggeredCharmAndBeauty:
             LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapTriggeredCharmAndBeauty **********";
             LOG(info) << "**********                           Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
-            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredCharmAndBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
+            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredCharmAndBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList, partPdgToReplaceList, freqReplaceList));
             break;
         case hf_generators::GapHF:
             LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapHF **********";
             LOG(info) << "**********                           Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
-            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
+            mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList, partPdgToReplaceList, freqReplaceList));
             break;
         default:
             LOG(fatal) << "********** [GeneratorPythia8EmbedHF] bad configuration, fix it! **********";
@@ -111,7 +119,7 @@ public:
         LOG(info) << "[notifyEmbedding] ----- Collision impact parameter: " << x;
 
         /// number of events to be embedded in a background event
-        mNumSigEvs = 5 + 0.886202881*std::pow(std::max(0.0f, 17.5f - x),1.7);
+        mNumSigEvs = static_cast<int>(std::lround(5.0 + 0.886202881 * std::pow(std::max(0.0f, 17.5f - x), 1.7)));
         LOG(info) << "[notifyEmbedding] ----- generating " << mNumSigEvs << " signal events " << std::endl;
     };
 
@@ -380,34 +388,34 @@ private:
 };
 
 // Charm enriched
-FairGenerator * GeneratorPythia8EmbedHFCharm(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
+FairGenerator * GeneratorPythia8EmbedHFCharm(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int,2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
 {
     auto myGen = new GeneratorPythia8EmbedHF();
 
     /// setup the internal generator for HF events
-    myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharm, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList);
+    myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharm, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList);
 
     return myGen;
 }
 
 // Beauty enriched
-FairGenerator * GeneratorPythia8EmbedHFBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
+FairGenerator * GeneratorPythia8EmbedHFBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int,2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
 {
     auto myGen = new GeneratorPythia8EmbedHF();
 
     /// setup the internal generator for HF events
-    myGen->setupGeneratorEvHF(hf_generators::GapTriggeredBeauty, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList);
+    myGen->setupGeneratorEvHF(hf_generators::GapTriggeredBeauty, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList);
 
     return myGen;
 }
 
 // Charm and beauty enriched (with same ratio)
-FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
+FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int,2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
 {
     auto myGen = new GeneratorPythia8EmbedHF();
 
     /// setup the internal generator for HF events
-    myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharmAndBeauty, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList);
+    myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharmAndBeauty, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList);
 
     return myGen;
 }
diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C b/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
index cdf472096..8d5060251 100644
--- a/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+++ b/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
@@ -258,17 +258,14 @@ protected:
 
   bool replaceParticle(int iPartToReplace, int pdgCodeNew) {
 
-    auto mothers = mPythia.event[iPartToReplace].motherList();
-
-    std::array<int, 25> pdgDiquarks = {1103, 2101, 2103, 2203, 3101, 3103, 3201, 3203, 3303, 4101, 4103, 4201, 4203, 4301, 4303, 4403, 5101, 5103, 5201, 5203, 5301, 5303, 5401, 5403, 5503};
-
-    for (auto const& mother: mothers) {
-      auto pdgMother = std::abs(mPythia.event[mother].id());
-      if (pdgMother > 100 && std::find(pdgDiquarks.begin(), pdgDiquarks.end(), pdgMother) == pdgDiquarks.end()) {
-        return false;
-      }
+    // Check status code: particles with status 91-99 come from decays and should not be replaced
+    int statusMother = std::abs(mPythia.event[iPartToReplace].status());
+    if (statusMother >= 91 && statusMother <= 99) {
+      return false;
     }
 
+    auto mothers = mPythia.event[iPartToReplace].motherList();
+
     int charge = mPythia.event[iPartToReplace].id() / std::abs(mPythia.event[iPartToReplace].id());
     float px = mPythia.event[iPartToReplace].px();
     float py = mPythia.event[iPartToReplace].py();
diff --git a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini
new file mode 100644
index 000000000..41fcc71cf
--- /dev/null
+++ b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini
@@ -0,0 +1,8 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
+funcName=GeneratorPythia8EmbedHFCharmAndBeauty(false, -1.5, 1.5, -1.5, 1.5, std::vector<int>{}, std::vector<int>{}, std::vector<std::array<int,2>>{{423,4132},{423,4232},{4212,4332}}, std::vector<float>{0.5,0.5,1})
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg
+includePartonEvent=true
\ No newline at end of file
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_Bforced.C b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_Bforced.C
index c3cb9bcfb..01fff33aa 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_Bforced.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_Bforced.C
@@ -90,8 +90,9 @@ int External()
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.75) // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.75) > uncFracForcedDecays) // we put some tolerance (e.g. due to oscillations which might change the final state)
     {
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_LambdaBToNuclei.C b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_LambdaBToNuclei.C
index 18a7dff79..af6b7e21e 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_LambdaBToNuclei.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_LambdaBToNuclei.C
@@ -45,8 +45,9 @@ int External()
     std::cout <<"# signal hadrons: " << nSignals << "\n";
     std::cout <<"# signal hadrons decaying into nuclei: " << nSignalGoodDecay << "\n";
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.8) // we put some tolerance (lambdaB in MB events do not coalesce)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.8) > uncFracForcedDecays) // we put some tolerance (lambdaB in MB events do not coalesce)
     {
         std::cerr << "Fraction of signals decaying into nuclei: " << fracForcedDecays << ", lower than expected\n";
         return 1;
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_bbbar.C b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_bbbar.C
index d67a9b48d..2a9c0f0d9 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_bbbar.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_bbbar.C
@@ -82,8 +82,9 @@ int External()
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) // we put some tolerance (e.g. due to oscillations which might change the final state)
     {
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_ccbar.C b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_ccbar.C
index 234b51bbe..73ee658f8 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_ccbar.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_ccbar.C
@@ -82,8 +82,9 @@ int External()
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.85) // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) // we put some tolerance (e.g. due to oscillations which might change the final state)
     {
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Bforced_gap5_Mode2.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Bforced_gap5_Mode2.C
index b5264c370..9f6a58d15 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Bforced_gap5_Mode2.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Bforced_gap5_Mode2.C
@@ -105,8 +105,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.85) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_BtoDK_gap5_Mode2.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_BtoDK_gap5_Mode2.C
index 36019e022..f3cdd9ea3 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_BtoDK_gap5_Mode2.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_BtoDK_gap5_Mode2.C
@@ -105,8 +105,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.85) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_OmegaC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_OmegaC.C
index e55bb2709..c0411d95e 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_OmegaC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_OmegaC.C
@@ -95,8 +95,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_XiC0_XiCplus.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_XiC0_XiCplus.C
index f2c79c2ad..6d83b21cc 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_XiC0_XiCplus.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_Mode2_XiC0_XiCplus.C
@@ -106,8 +106,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_gap5.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_gap5.C
index 7fb3eacb6..f92a829c5 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_gap5.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_bbbar_gap5.C
@@ -102,8 +102,9 @@ int External()
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.7) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.8) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_OmegaC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_OmegaC.C
index 081dfecae..c71340129 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_OmegaC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_OmegaC.C
@@ -95,8 +95,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_XiC0_XiCplus.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_XiC0_XiCplus.C
index e2f818fb2..17cddbb04 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_XiC0_XiCplus.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_Mode2_XiC0_XiCplus.C
@@ -106,8 +106,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_Mode2_ptHardBins.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_Mode2_ptHardBins.C
index 11cb92fc3..4e2ef1e28 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_Mode2_ptHardBins.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_Mode2_ptHardBins.C
@@ -123,13 +123,14 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
 
-    if (averagePt < 6.5) { // by testing locally it should be around 8.5 GeV/c with pthard bin 20-200 (contrary to 2-2.5 GeV/c of SoftQCD)
+    if (averagePt < 6) { // by testing locally it should be around 8.5 GeV/c with pthard bin 20-200 (contrary to 2-2.5 GeV/c of SoftQCD)
         std::cerr << "Average pT of charmed hadrons " << averagePt << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb.C
index d3f861f24..286e79eba 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb.C
@@ -181,8 +181,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.C
index 4437bf04f..49d95f084 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.C
@@ -172,8 +172,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C
index d3f861f24..286e79eba 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C
@@ -181,8 +181,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C
new file mode 100644
index 000000000..68cc27ce0
--- /dev/null
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C
@@ -0,0 +1,210 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+    //std::string path{"tf1/sgn_Kine.root"};
+
+    // Particle replacement configuration: {original, replacement}, frequencies
+    std::array<std::array<int, 2>, 3> pdgReplParticles = {
+        std::array{423, 4132},   // D*0 -> Xic0
+        std::array{423, 4232},   // D*0 -> Xic+
+        std::array{4212, 4332}   // Sigmac+ -> Omegac0
+    };
+    std::array<std::array<int, 2>, 3> pdgReplPartCounters = {
+        std::array{0, 0}, 
+        std::array{0, 0}, 
+        std::array{0, 0}
+    };
+    std::array<float, 3> freqRepl = {0.5, 0.5, 1.0};
+    std::map<int, int> sumOrigReplacedParticles = {{423, 0}, {4212, 0}};
+
+    // Signal hadrons to check (only final charm baryons after replacement)
+    std::array<int, 4> checkPdgHadron{4122, 4132, 4232, 4332};
+
+    // Expected decay channels - will be sorted automatically
+    // Define both particle and antiparticle versions
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecaysRaw{
+        // Λc+ decays (from cfg: 4122:addChannel + resonance decays)
+        {4122, {
+            {2212, 311}, {-2212, -311},                           // p K0s
+            {2212, -321, 211}, {-2212, 321, -211},                // p K- π+ (non-resonant)
+            {2212, 313}, {-2212, -313},                           // p K*0 (not decayed)
+            {2212, 321, 211}, {-2212, -321, -211},                // p K*0 -> p (K- π+) [K*0 decayed]
+            {2224, -321}, {-2224, 321},                           // Delta++ K- (not decayed)
+            {2212, 211, -321}, {-2212, -211, 321},                // Delta++ K- -> (p π+) K- [Delta decayed]
+            {102134, 211}, {-102134, -211},                       // Lambda(1520) π+ (not decayed)
+            {2212, 321, 211}, {-2212, -321, -211},                // Lambda(1520) π+ -> (p K-) π+ [Lambda* decayed]
+            {2212, -321, 211, 111}, {-2212, 321, -211, 111},      // p K- π+ π0
+            {2212, -211, 211}, {-2212, 211, -211},                // p π- π+ (cfg line 61: 2212 -211 211)
+            {2212, 333}, {-2212, 333},                            // p φ (not decayed)
+            {2212, 321, -321}, {-2212, -321, 321}                 // p φ -> p (K+ K-) [φ decayed]
+        }},
+        // Ξc0 decays (from cfg: 4132:onIfMatch)
+        {4132, {
+            {3312, 211}, {-3312, -211},                           // Ξ- π+
+            {3334, 321}, {-3334, -321}                            // Ω- K+
+        }},
+        // Ξc+ decays (from cfg: 4232:onIfMatch + resonance decays)
+        {4232, {
+            {2212, -321, 211}, {-2212, 321, -211},                // p K- π+
+            {2212, -313}, {-2212, 313},                           // p K̄*0 (not decayed)
+            {2212, -321, 211}, {-2212, 321, -211},                // p K̄*0 -> p (K+ π-) [K*0 decayed]
+            {2212, 333}, {-2212, 333},                            // p φ (not decayed)
+            {2212, 321, -321}, {-2212, -321, 321},                // p φ -> p (K+ K-) [φ decayed]
+            {3222, -211, 211}, {-3222, 211, -211},                // Σ+ π- π+
+            {3324, 211}, {-3324, -211},                           // Ξ*0 π+
+            {3312, 211, 211}, {-3312, -211, -211}                 // Ξ- π+ π+
+        }},
+        // Ωc0 decays (from cfg: 4332:onIfMatch)
+        {4332, {
+            {3334, 211}, {-3334, -211},                           // Ω- π+
+            {3312, 211}, {-3312, -211},                           // Ξ- π+
+            {3334, 321}, {-3334, -321}                            // Ω- K+
+        }}
+    };
+
+    // Sort all decay channels
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays;
+    for (auto &[pdg, decays] : checkHadronDecaysRaw) {
+        for (auto decay : decays) {
+            std::sort(decay.begin(), decay.end());
+            checkHadronDecays[pdg].push_back(decay);
+        }
+    }
+
+    TFile file(path.c_str(), "READ");
+    if (file.IsZombie())
+    {
+        std::cerr << "Cannot open ROOT file " << path << "\n";
+        return 1;
+    }
+
+    auto tree = (TTree *)file.Get("o2sim");
+    std::vector<o2::MCTrack> *tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+
+    int nSignals{}, nSignalGoodDecay{};
+    std::map<int, int> failedDecayCount{{4122, 0}, {4132, 0}, {4232, 0}, {4332, 0}};
+    std::map<int, std::set<std::vector<int>>> unknownDecays;
+    auto nEvents = tree->GetEntries();
+
+    for (int i = 0; i < nEvents; i++)
+    {
+        tree->GetEntry(i);
+        for (auto &track : *tracks)
+        {
+            auto pdg = track.GetPdgCode();
+            auto absPdg = std::abs(pdg);
+
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), absPdg) != checkPdgHadron.end())
+            {
+                nSignals++; // count signal PDG
+
+                // Count replacement particles (single-match per track)
+                int matchedIdx = -1;
+                bool matchedIsReplacement = false;
+                for (int iRepl{0}; iRepl < 3; ++iRepl) {
+                    if (absPdg == pdgReplParticles[iRepl][0]) {
+                        matchedIdx = iRepl;
+                        matchedIsReplacement = false;
+                        break;
+                    }
+                    if (absPdg == pdgReplParticles[iRepl][1]) {
+                        matchedIdx = iRepl;
+                        matchedIsReplacement = true;
+                        break;
+                    }
+                }
+                if (matchedIdx >= 0) {
+                    if (matchedIsReplacement) {
+                        pdgReplPartCounters[matchedIdx][1]++;
+                    } else {
+                        pdgReplPartCounters[matchedIdx][0]++;
+                    }
+                    // Count the original-particle population once for this matched group
+                    sumOrigReplacedParticles[pdgReplParticles[matchedIdx][0]]++;
+                }
+
+                // Collect decay products
+                std::vector<int> pdgsDecay{};
+                if (track.getFirstDaughterTrackId() >= 0 && track.getLastDaughterTrackId() >= 0) {
+                    for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                        auto pdgDau = tracks->at(j).GetPdgCode();
+                        pdgsDecay.push_back(pdgDau);
+                    }
+                }
+
+                std::sort(pdgsDecay.begin(), pdgsDecay.end());
+
+                // Check if decay matches expected channels
+                bool foundMatch = false;
+                for (auto &decay : checkHadronDecays[absPdg]) {
+                    if (pdgsDecay == decay) {
+                        nSignalGoodDecay++;
+                        foundMatch = true;
+                        break;
+                    }
+                }
+
+                // Record failed decays for debugging
+                if (!foundMatch && pdgsDecay.size() > 0) {
+                    failedDecayCount[absPdg]++;
+                    unknownDecays[absPdg].insert(pdgsDecay);
+                }
+            }
+        }
+    }
+
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    std::cout << "# signal charm baryons: " << nSignals << "\n";
+    std::cout << "# signal charm baryons decaying in the correct channel: " << nSignalGoodDecay << "\n";
+
+    // Print failed decay statistics
+    std::cout << "\nFailed decay counts:\n";
+    for (auto &[pdg, count] : failedDecayCount) {
+        if (count > 0) {
+            std::cout << "PDG " << pdg << ": " << count << " failed decays\n";
+            std::cout << "  Unknown decay channels (first 5):\n";
+            int printed = 0;
+            for (auto &decay : unknownDecays[pdg]) {
+                if (printed++ >= 5) break;
+                std::cout << "    [";
+                for (size_t i = 0; i < decay.size(); ++i) {
+                    std::cout << decay[i];
+                    if (i < decay.size()-1) std::cout << ", ";
+                }
+                std::cout << "]\n";
+            }
+        }
+    }
+    std::cout << "\n";
+
+    std::cout << "# D*0 (original): " << pdgReplPartCounters[0][0] << "\n";
+    std::cout << "# Xic0 (replaced from D*0): " << pdgReplPartCounters[0][1] << "\n";
+    std::cout << "# Xic+ (replaced from D*0): " << pdgReplPartCounters[1][1] << "\n";
+    std::cout << "# Sigmac+ (original): " << pdgReplPartCounters[2][0] << "\n";
+    std::cout << "# Omegac0 (replaced from Sigmac+): " << pdgReplPartCounters[2][1] << "\n";
+
+    // Check forced decay fraction
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    std::cout << "# fraction of signals decaying into the correct channel: " << fracForcedDecays
+              << " (" << fracForcedDecays * 100.0f << "%)\n";
+    if (fracForcedDecays < 0.9f) { // 90% threshold with tolerance
+        std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
+        return 1;
+    }
+
+    // Check particle replacement ratios (2-sigma statistical compatibility)
+    for (int iRepl{0}; iRepl<6; ++iRepl) {
+        float numPart = sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
+        float fracMeas = numPart ? float(pdgReplPartCounters[iRepl][1]) / numPart : 0.0f;
+        float fracMeasUnc = numPart ? std::sqrt(fracMeas * (1 - fracMeas) / numPart) : 1.0f;
+ 
+        if (std::abs(fracMeas - freqRepl[iRepl]) > fracMeasUnc) {
+            std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0] << " into " << pdgReplParticles[iRepl][1] << " is " << fracMeas <<" (expected "<< freqRepl[iRepl] << ")\n";
+            return 1;    
+        }
+    }
+
+    return 0;
+}
\ No newline at end of file
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_OO.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_OO.C
index 5cc0e18e6..186bfd58a 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_OO.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_OO.C
@@ -179,8 +179,9 @@ int External(std::string path = "o2sim_Kine.root") {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_pO.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_pO.C
index 5cc0e18e6..186bfd58a 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_pO.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_pO.C
@@ -179,8 +179,9 @@ int External(std::string path = "o2sim_Kine.root") {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3.C
index d29245b32..3369fcc6e 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3.C
@@ -118,8 +118,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2.C
index 5735c969d..e29821f63 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2.C
@@ -170,8 +170,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.85) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2_XiC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2_XiC.C
index c99426354..f6fd63c6b 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2_XiC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap3_Mode2_XiC.C
@@ -112,8 +112,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5.C
index e5c0bf08d..8c29ce209 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5.C
@@ -118,8 +118,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C
index 3c1cfd492..c88b50ce3 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C
@@ -130,8 +130,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DResoTrigger.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DResoTrigger.C
index f70d958d5..cd06ab74d 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DResoTrigger.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DResoTrigger.C
@@ -124,18 +124,19 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) : 0.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
 
     for (int iRepl{0}; iRepl<6; ++iRepl) {
-        if (std::abs(pdgReplPartCounters[iRepl][1] - freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]) > 2 * std::sqrt(freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]])) { // 2 sigma compatibility
-            float fracMeas = 0.;
-            if (sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] > 0.) {
-                fracMeas = float(pdgReplPartCounters[iRepl][1]) / sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
-            } 
+        float numPart = sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
+        float fracMeas = numPart ? float(pdgReplPartCounters[iRepl][1]) / numPart : 0.0f;
+        float fracMeasUnc = numPart ? std::sqrt(fracMeas * (1 - fracMeas) / numPart) : 1.0f;
+ 
+        if (std::abs(fracMeas - freqRepl[iRepl]) > fracMeasUnc) {
             std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0] << " into " << pdgReplParticles[iRepl][1] << " is " << fracMeas <<" (expected "<< freqRepl[iRepl] << ")\n";
             return 1;    
         }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2.C
index 73ca51134..f0f964432 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2.C
@@ -180,8 +180,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons.C
index e8a678315..edbb29cf9 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons.C
@@ -128,8 +128,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons_pp_ref.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons_pp_ref.C
index 9e8456fad..c8fde6c70 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons_pp_ref.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_CharmBaryons_pp_ref.C
@@ -179,8 +179,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC.C
index b9b398cfe..02dffdac1 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC.C
@@ -112,8 +112,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC_OmegaC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC_OmegaC.C
index 909ba4ccb..ff6ad16da 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC_OmegaC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_XiC_OmegaC.C
@@ -113,8 +113,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.C
index eeb37acfc..79485d1fe 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.C
@@ -170,8 +170,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.8) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.C
index 31aba24ad..2a5a5ae64 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.C
@@ -199,23 +199,22 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.5) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) : 1.0f;
+    if (std::abs(fracForcedDecays - 0.5) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
 
-    for (int iRepl{0}; iRepl<2; ++iRepl) {
-
-        std::cout << " --- pdgReplPartCounters[" << iRepl << "][1] = " << pdgReplPartCounters[iRepl][1] << ", freqRepl[" << iRepl <<"] = " << freqRepl[iRepl] << ", sumOrigReplacedParticles[pdgReplParticles[" << iRepl << "][0]] = " << sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] << std::endl; 
-        
-        if (std::abs(pdgReplPartCounters[iRepl][1] - freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]) > 2 * std::sqrt(freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]])) { // 2 sigma compatibility
-            float fracMeas = 0.;
-            if (sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] > 0.) {
-                fracMeas = float(pdgReplPartCounters[iRepl][1]) / sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
-            } 
+    // Check particle replacement ratios
+    for (int iRepl{0}; iRepl < 3; ++iRepl) {
+        float numPart = sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
+        float fracMeas = numPart ? float(pdgReplPartCounters[iRepl][1]) / numPart : 0.0f;
+        float fracMeasUnc = numPart ? std::sqrt(fracMeas * (1 - fracMeas) / numPart) : 1.0f;
+ 
+        if (std::abs(fracMeas - freqRepl[iRepl]) > fracMeasUnc) {
             std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0] << " into " << pdgReplParticles[iRepl][1] << " is " << fracMeas <<" (expected "<< freqRepl[iRepl] << ")\n";
-            return 1;    
+            return 1;
         }
     }
 
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_pp_ref.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_pp_ref.C
index 6f13ae05e..5886cbe8f 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_pp_ref.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_pp_ref.C
@@ -180,8 +180,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8.C
index b6828e205..1906c59cf 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8.C
@@ -118,8 +118,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8_Mode2_XiC.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8_Mode2_XiC.C
index 596639d12..ab818b8f2 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8_Mode2_XiC.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap8_Mode2_XiC.C
@@ -112,8 +112,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_gap5.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_gap5.C
index 99a039e6f..e0f07b6af 100644
--- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_gap5.C
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_gap5.C
@@ -101,8 +101,9 @@ int External() {
         return 1;
     }
 
-    float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
-    if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
+    float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
+    float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) / nSignals : 1.0f;
+    if (std::abs(fracForcedDecays - 0.85) > uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
         std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
         return 1;
     }
diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg
new file mode 100644
index 000000000..84ac21536
--- /dev/null
+++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg
@@ -0,0 +1,156 @@
+### authors: Fabrizio Grosa (fabrizio.grosa@cern.ch)
+###          Ruiqi Yin (ruiqi.yin@cern.ch)
+### last update:  October 2025
+
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 5360. 		# GeV
+
+### processes: c-cbar and b-bbar processes
+HardQCD:hardccbar on
+HardQCD:hardbbbar on
+HardQCD:gg2ccbar on
+HardQCD:qqbar2ccbar on
+HardQCD:gg2bbbar on
+HardQCD:qqbar2bbbar on
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.
+
+### switching on Pythia Mode2
+ColourReconnection:mode 1
+ColourReconnection:allowDoubleJunRem off
+ColourReconnection:m0 0.3
+ColourReconnection:allowJunctions on
+ColourReconnection:junctionCorrection 1.20
+ColourReconnection:timeDilationMode 2
+ColourReconnection:timeDilationPar 0.18
+StringPT:sigma 0.335
+StringZ:aLund 0.36
+StringZ:bLund 0.56
+StringFlav:probQQtoQ 0.078
+StringFlav:ProbStoUD 0.2
+StringFlav:probQQ1toQQ0join 0.0275,0.0275,0.0275,0.0275
+MultiPartonInteractions:pT0Ref 2.15
+BeamRemnants:remnantMode 1
+BeamRemnants:saturation 5
+
+# Correct decay lengths (wrong in PYTHIA8 decay table)
+# Lb
+5122:tau0 = 0.4390
+# Xic0
+4132:tau0 = 0.0455
+# OmegaC
+4332:tau0 = 0.0803
+
+### Force charm baryon decay modes (Ξc and Ωc)
+### D mesons use default PYTHIA8 branching ratios
+
+## Λc decays (optional, can be kept for reference)
+### Λc+ -> p K- π+ (36%)
+4122:oneChannel = 1 0.14400 0 2212 -321 211       ### Λc+ -> p K- π+ (non-resonant) 3.5%
+4122:addChannel = 1 0.08100 100 2212 -313         ### Λc+ -> p K*0(892) 1.96%
+4122:addChannel = 1 0.04500 100 2224 -321         ### Λc+ -> Delta++ K- 1.08%
+4122:addChannel = 1 0.09000 100 102134 211        ### Λc+ -> Lambda(1520) K- 2.20e-3
+### Λc+ -> p K0S (36%)
+4122:addChannel = 1 0.36000 0 2212 311            ### Λc+ -> p K0S 1.59%
+### Λc+ -> p K- π+ π0 (small, 3%)
+4122:addChannel = 1 0.03000 0 2212 -321 211 111   ### Λc+ -> p K- π+ π0 (non-resonant) 4.6%
+### Λc+ -> p π- π+ (12.50%)
+4122:addChannel = 1 0.12500 0 2212 -211 211       ### Λc+ -> p π+ π+ 4.59%
+### Λc+ -> p K- K+ (12.50%)
+4122:addChannel = 1 0.12500 0 2212 333            ### Λc+ -> p phi 1.06%
+
+## Xic decays
+### Ξc+ -> p K- π+ (35%)
+4232:oneChannel = 1 0.17500 0 2212 -321 211        ### Ξc+ -> p K- π+ 6.18e-3
+4232:addChannel = 1 0.17500 0 2212 -313            ### Ξc+ -> p antiK*0(892)
+### Ξc+ -> Ξ- π+ π+ (35%) (set the same as Ξc+ -> p K- π+)
+4232:addChannel = 1 0.35000 0 3312 211 211        ### Ξc+ -> Ξ- π+ π+ 2.86%
+### Ξc+ -> p φ (10%)
+4232:addChannel = 1 0.10000 0 2212 333             ### Ξc+ -> p φ
+### Ξc+ -> sigma+ π+ π- (10%)
+4232:addChannel = 1 0.12500 0 3222 -211 211        ### Ξc+ -> sigma+ π- π+ 1.37%
+### Ξc+ -> Ξ*0 π+ (10%)
+4232:addChannel = 1 0.12500 0 3324 211
+
+### add Xic0 decays absent in PYTHIA8 decay table
+4132:oneChannel = 1 0.0143 0 3312 211
+
+### add OmegaC decays absent in PYTHIA8 decay table
+4332:oneChannel = 1 0.5 0 3334 211
+4332:addChannel = 1 0.5 0 3312 211
+
+# Allow the decay of resonances in the decay chain (only for charm baryons)
+### for Λc -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Λc -> Lambda(1520) K-
+102134:onMode = off
+102134:onIfAll = 2212 321
+### for Λc -> p K*0(892)
+313:onMode = off
+313:onIfAll = 321 211
+### for Λc -> p phi
+333:onMode = off
+333:onIfAll = 321 321
+### for Xic0 -> pi Xi -> pi pi Lambda -> pi pi pi p 
+### and Omega_c -> pi Xi -> pi pi Lambda -> pi pi pi p
+3312:onMode = off
+3312:onIfAll = 3122 -211
+3122:onMode = off
+3122:onIfAll = 2212 -211
+### for Omega_c -> pi Omega -> pi K Lambda -> pi K pi p
+3334:onMode = off
+3334:onIfAll = 3122 -321
+
+# Switch off charm baryon decay channels (D mesons use PYTHIA8 defaults)
+4122:onMode = off
+4232:onMode = off
+4132:onMode = off
+4332:onMode = off
+
+# Allow the decay of charm baryons
+### Λc -> pK0s
+4122:onIfMatch = 2212 311 
+### Λc -> p K- π+ π0
+4122:onIfMatch = 2212 321 211
+### Λc -> p K*
+4122:onIfMatch = 2212 313
+### Λc -> Delta++ K
+4122:onIfMatch = 2224 321
+### Λc -> Lambda(1520) π
+4122:onIfMatch = 102134 211
+### Λc -> p K- π+ π0
+4122:onIfMatch = 2212 321 211 111
+### Λc -> p π π
+4122:onIfMatch = 2212 211 211
+### Λc -> p K K
+4122:onIfMatch = 2212 333
+
+### Ξc+ -> p K- π+
+4232:onIfMatch = 2212 321 211
+### Ξc+ -> p antiK*0(892)
+4232:onIfMatch = 2212 313
+### Ξc+ -> p φ
+4232:onIfMatch = 2212 333
+### Ξc+ -> sigma- π+ π+
+4232:onIfMatch = 3222 211 211
+### Ξc+ -> Ξ*0 π+, Ξ*0 -> Ξ- π+
+4232:onIfMatch = 3324 211
+### Ξc+ -> Ξ- π+ π+
+4232:onIfMatch = 3312 211 211
+
+### Xic0 -> Xi- pi+
+4132:onIfMatch = 3312 211
+### Xic0 -> Omega- Ka+
+4132:onIfMatch = 3334 321
+
+### Omega_c -> Omega pi
+4332:onIfMatch = 3334 211
+### Omega_c -> Xi pi
+4332:onIfMatch = 3312 211
+### Omega_c -> Omega- Ka+
+4332:onIfMatch = 3334 321
\ No newline at end of file