diff --git a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.ini b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.ini
new file mode 100644
index 000000000..143f9589f
--- /dev/null
+++ b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.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()
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV_CorrBkg.cfg
+includePartonEvent=true
diff --git a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.ini b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.ini
index 49102a6b2..a88c48b63 100644
--- a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.ini
+++ b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.ini
@@ -4,5 +4,5 @@ fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_py
 funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5)
 
 [GeneratorPythia8]
-config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_CorrBkg.cfg
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkg.cfg
 includePartonEvent=true
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
new file mode 100644
index 000000000..ef17d7a24
--- /dev/null
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.C
@@ -0,0 +1,172 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+    //std::string path{"tf1/sgn_1_Kine.root"};
+
+    int checkPdgQuarkOne{4};
+    int checkPdgQuarkTwo{5};
+    float ratioTrigger = 1.; // one event triggered out of 1
+
+    std::vector<int> checkPdgHadron{411, 421, 431, 4122, 4232};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {411, {
+            {-321, 211, 211},               // K- π+ π+ (non-resonant)
+            {-313, 321},                    // K*0(892) K+
+            {-10311, 321},                  // K*0(1430) K+
+            {211, 333},                     // φ π+
+            {-321, 321, 211},              // K- K+ π+ (non-resonant)
+            {113, 211},                     // ρ0 π+
+            {225, 211},                     // f2(1270) π+
+            {-211, 211, 211}               // π- π+ π+ (non-resonant)
+        }},
+        {421, {
+            {-321, 211},                   // K- π+ (non-resonant)
+            {-321, 111, 211},              // K- π+ π0
+            {213, -321},                   // ρ+ K-
+            {-313, 111},                   // antiK*0(892) π0
+            {-323, 211},                   // K*-(892) π+
+            {-211, 211},                   // π- π+
+            {213, -211},                   // ρ+ π-
+            {-211, 211, 111},              // π- π+ π0
+            {-321, 321}                   // K- K+
+        }},
+        {431, {
+            {211, 333},                    // φ π+
+            {-313, 321},                   // antiK*(892) K+
+            {333, 213},                    // φ ρ
+            {113, 211},                    // ρ π+
+            {225, 211},                    // f2(1270) π+
+            {-211, 211, 211},              // π- π+ π+ (s-wave)
+            {313, 211},                    // K*(892)0 π+
+            {10221, 321},                  // f0(1370) K+
+            {113, 321},                    // ρ0 K+
+            {-211, 321, 211},              // π- K+ π+ (non-resonant)
+            {221, 211}                    // η π+
+        }},
+        {4122, {
+            {2212, -321, 211},             // p K- π+ (non-resonant)
+            {2212, -313},                  // p K*0(892)
+            {2224, -321},                  // Δ++ K-
+            {102134, 211},                 // Λ(1520) K-
+            {2212, -321, 211, 111},        // p K- π+ π0
+            {2212, -211, 211},             // p π- π+
+            {2212, 333}                   // p φ
+        }},
+        {4232, {
+            {-313, 2212},                  // antiK*0(892) p
+            {2212, -321, 211},             // p K- π+
+            {2212, 333},                   // p φ
+            {3222, -211, 211}             // Σ+ π- π+
+        }},
+    };
+    
+
+    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);
+    o2::dataformats::MCEventHeader *eventHeader = nullptr;
+    tree->SetBranchAddress("MCEventHeader.", &eventHeader);
+
+    int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
+    int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{};
+    auto nEvents = tree->GetEntries();
+
+    for (int i = 0; i < nEvents; i++) {
+        tree->GetEntry(i);
+
+        // check subgenerator information
+        //if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
+        //    bool isValid = false;
+        //    int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
+        //    if (subGeneratorId == 0) {
+        //        nEventsMB++;
+        //    } else if (subGeneratorId == checkPdgQuarkOne) {
+        //        nEventsInjOne++;
+        //    } else if (subGeneratorId == checkPdgQuarkTwo) {
+        //        nEventsInjTwo++;
+        //    }
+        //}
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::abs(pdg) == checkPdgQuarkOne) {
+                nQuarksOne++;
+                continue;
+            }
+            if (std::abs(pdg) == checkPdgQuarkTwo) {
+                nQuarksTwo++;
+                continue;
+            }
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
+                nSignals++; // count signal PDG
+
+                std::vector<int> pdgsDecay{};
+                std::vector<int> pdgsDecayAntiPart{};
+                for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                    auto pdgDau = tracks->at(j).GetPdgCode();
+                    pdgsDecay.push_back(pdgDau);
+                    if (pdgDau != 333) { // phi is antiparticle of itself
+                        pdgsDecayAntiPart.push_back(-pdgDau);
+                    } else {
+                        pdgsDecayAntiPart.push_back(pdgDau);
+                    }
+                }
+
+                std::sort(pdgsDecay.begin(), pdgsDecay.end());
+                std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
+
+                for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
+                    if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
+                        nSignalGoodDecay++;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    //std::cout << "# MB events: " << nEventsMB << "\n";
+    //std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
+    //std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n";
+    std::cout <<"# signal hadrons: " << nSignals << "\n";
+    std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
+
+    //if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.05) { // we put some tolerance since the number of generated events is small
+    //    std::cerr << "Number of generated MB events different than expected\n";
+    //    return 1;
+    //}
+    //if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
+    //    std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
+    //    return 1;
+    //}
+    //if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
+    //    std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
+    //    return 1;
+    //}
+
+    if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
+        std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n";
+        return 1;
+    }
+    if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
+        std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n";
+        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)
+        std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
+        return 1;
+    }
+
+    return 0;
+}
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
new file mode 100644
index 000000000..688939eae
--- /dev/null
+++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.C
@@ -0,0 +1,170 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuarkOne{4};
+    int checkPdgQuarkTwo{5};
+    float ratioTrigger = 1./5; // one event triggered out of 5
+
+    std::vector<int> checkPdgHadron{411, 421, 431, 4122, 4232};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {411, {
+            {-321, 211, 211},               // K- π+ π+ (non-resonant)
+            {-313, 321},                    // K*0(892) K+
+            {-10311, 321},                  // K*0(1430) K+
+            {211, 333},                     // φ π+
+            {-321, 321, 211},              // K- K+ π+ (non-resonant)
+            {113, 211},                     // ρ0 π+
+            {225, 211},                     // f2(1270) π+
+            {-211, 211, 211}               // π- π+ π+ (non-resonant)
+        }},
+        {421, {
+            {-321, 211},                   // K- π+ (non-resonant)
+            {-321, 111, 211},              // K- π+ π0
+            {213, -321},                   // ρ+ K-
+            {-313, 111},                   // antiK*0(892) π0
+            {-323, 211},                   // K*-(892) π+
+            {-211, 211},                   // π- π+
+            {213, -211},                   // ρ+ π-
+            {-211, 211, 111},              // π- π+ π0
+            {-321, 321}                   // K- K+
+        }},
+        {431, {
+            {211, 333},                    // φ π+
+            {-313, 321},                   // antiK*(892) K+
+            {333, 213},                    // φ ρ
+            {113, 211},                    // ρ π+
+            {225, 211},                    // f2(1270) π+
+            {-211, 211, 211},              // π- π+ π+ (s-wave)
+            {313, 211},                    // K*(892)0 π+
+            {10221, 321},                  // f0(1370) K+
+            {113, 321},                    // ρ0 K+
+            {-211, 321, 211},              // π- K+ π+ (non-resonant)
+            {221, 211}                    // η π+
+        }},
+        {4122, {
+            {2212, -321, 211},             // p K- π+ (non-resonant)
+            {2212, -313},                  // p K*0(892)
+            {2224, -321},                  // Δ++ K-
+            {102134, 211},                 // Λ(1520) K-
+            {2212, -321, 211, 111},        // p K- π+ π0
+            {2212, -211, 211},             // p π- π+
+            {2212, 333}                   // p φ
+        }},
+        {4232, {
+            {-313, 2212},                  // antiK*0(892) p
+            {2212, -321, 211},             // p K- π+
+            {2212, 333},                   // p φ
+            {3222, -211, 211}             // Σ+ π- π+
+        }},
+    };
+
+    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);
+    o2::dataformats::MCEventHeader *eventHeader = nullptr;
+    tree->SetBranchAddress("MCEventHeader.", &eventHeader);
+
+    int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
+    int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{};
+    auto nEvents = tree->GetEntries();
+
+    for (int i = 0; i < nEvents; i++) {
+        tree->GetEntry(i);
+
+        // check subgenerator information
+        if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
+            bool isValid = false;
+            int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
+            if (subGeneratorId == 0) {
+                nEventsMB++;
+            } else if (subGeneratorId == checkPdgQuarkOne) {
+                nEventsInjOne++;
+            } else if (subGeneratorId == checkPdgQuarkTwo) {
+                nEventsInjTwo++;
+            }
+        }
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::abs(pdg) == checkPdgQuarkOne) {
+                nQuarksOne++;
+                continue;
+            }
+            if (std::abs(pdg) == checkPdgQuarkTwo) {
+                nQuarksTwo++;
+                continue;
+            }
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
+                nSignals++; // count signal PDG
+
+                std::vector<int> pdgsDecay{};
+                std::vector<int> pdgsDecayAntiPart{};
+                for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
+                    auto pdgDau = tracks->at(j).GetPdgCode();
+                    pdgsDecay.push_back(pdgDau);
+                    if (pdgDau != 333) { // phi is antiparticle of itself
+                        pdgsDecayAntiPart.push_back(-pdgDau);
+                    } else {
+                        pdgsDecayAntiPart.push_back(pdgDau);
+                    }
+                }
+
+                std::sort(pdgsDecay.begin(), pdgsDecay.end());
+                std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
+
+                for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
+                    if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
+                        nSignalGoodDecay++;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    std::cout << "# MB events: " << nEventsMB << "\n";
+    std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
+    std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n";
+    std::cout <<"# signal hadrons: " << nSignals << "\n";
+    std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
+
+    if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.05) { // we put some tolerance since the number of generated events is small
+        std::cerr << "Number of generated MB events different than expected\n";
+        return 1;
+    }
+    if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
+        return 1;
+    }
+    if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
+        return 1;
+    }
+
+    if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
+        std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n";
+        return 1;
+    }
+    if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
+        std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n";
+        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)
+        std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
+        return 1;
+    }
+
+    return 0;
+}
diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_CorrBkg.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_CorrBkg.cfg
deleted file mode 100644
index 4e39ac46f..000000000
--- a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_CorrBkg.cfg
+++ /dev/null
@@ -1,158 +0,0 @@
-### authors: Stefano Politano (stefano.politano@cern.ch)
-### last update:  April 2025
-
-### beams
-Beams:idA 2212			# proton
-Beams:idB 2212 			# proton
-Beams:eCM 13600. 		# GeV
-
-### processes
-SoftQCD:inelastic on		# all inelastic processes
-
-### 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
-
-### HF decays
-### BR are set to yield the same amount of counts for each final state
-###
-### D0 decays
-421:oneChannel = 1 0.170 0 -321 211             ### D0 -> K- π+
-421:addChannel = 1 0.166 0 -321 211 111         ### D0 -> K- π+ π0
-421:addChannel = 1 0.166 0 -211 211             ### D0 -> π- π+
-421:addChannel = 1 0.166 0 -211 211 111         ### D0 -> π- π+ π0
-421:addChannel = 1 0.166 0 -321 321             ### D0 -> K- K+
-421:addChannel = 1 0.166 0 -321 321 111         ### D0 -> K- K+ π0
-
-### D+ decays
-411:oneChannel = 1 0.17 0 -321 211 211          ### D+ -> K- π+ π+
-411:addChannel = 1 0.33 0 -211 211 211          ### D+ -> π- π+ π+
-411:addChannel = 1 0.17 0  313 211              ### D+ -> K* π+ -> K- π+ π+
-411:addChannel = 1 0.33 0  333 211              ### D+ -> φ π+ -> K- K+ π+
-
-### Ds+ decays
-431:oneChannel = 1 0.125 0  333 211             ### Ds+ -> φ π+ -> K- K+ π+
-431:addChannel = 1 0.125 0  313 321             ### Ds+ -> K* K+ -> K- K+ π+
-431:addChannel = 1 0.250 0 -211 211 211         ### Ds+ -> π- π+ π+
-431:addChannel = 1 0.250 0 -211 321 211         ### Ds+ -> π- K+ π+
-431:addChannel = 1 0.250 0 -321 321 321         ### Ds+ -> K- K+ K+
-
-## Lc decays
-4122:oneChannel = 1 0.0825 0 2212 -321 211      ### Λc+ -> p K- π+
-4122:addChannel = 1 0.0825 100 2212 -313        ### Λc+ -> p K* -> p K- π+
-4122:addChannel = 1 0.0825 100 2224 -321        ### Λc+ -> Delta++ K- -> p K- π+
-4122:addChannel = 1 0.0825 100 102134 211       ### Λc+ -> Lambda(1520) K- -> p K- π+
-4122:addChannel = 1 0.1650 0 2212 211 211       ### Λc+ -> p π+ π+
-4122:addChannel = 1 0.1650 0 2212 311           ### Λc+ -> p K0s -> p π+ π-
-4122:addChannel = 1 0.3300 0 2212 -321 321      ### Λc+ -> p K- K+
-
-### Ξc+ decays
-4232:oneChannel = 1 0.34 0 2212 -321 211        ### Ξc+ -> p K- π+
-4232:addChannel = 1 0.33 0 -3312 211 211        ### Ξc+ -> Ξ- π+ π+
-4232:addChannel = 1 0.33 0 2212 333             ### Ξc+ -> p φ -> p K- K+
-
-# Allow the decay of resonances in the decay chain
-### K* -> K π
-313:onMode = off
-313:onIfAll = 321 211
-### for Ds -> φ π+
-333:onMode = off
-333:onIfAll = 321 321
-### for D0 -> rho0 π+ k-
-113:onMode = off
-113:onIfAll = 211 211
-### for Λc -> Delta++ K-
-2224:onMode = off
-2224:onIfAll = 2212 211
-### for Λc -> Lambda(1520) K-
-102134:onMode = off
-102134:onIfAll = 2212 321
-
-### Switch off all decay channels
-411:onMode = off
-421:onMode = off
-431:onMode = off
-4122:onMode = off
-4232:onMode = off
-
-# Allow the decay of HF
-### D0 -> K π
-421:onIfMatch = 321 211
-### D0 -> K π π0
-421:onIfMatch = 321 211 111
-### D0 -> π π
-421:onIfMatch = 211 211
-### D0 -> π π π0
-421:onIfMatch = 211 211 111
-### D0 -> K K
-421:onIfMatch = 321 321
-### D0 -> K K π0
-421:onIfMatch = 321 321 111
-
-### D+/- -> K π π 
-411:onIfMatch = 321 211 211
-### D+/- -> K* π
-411:onIfMatch = 313 211
-### D+/- -> φ π
-411:onIfMatch = 333 211
-### D+/- -> π π π
-411:onIfMatch = 211 211 211
-
-### Ds -> φ π
-431:onIfMatch = 333 211
-### Ds -> K* K
-431:onIfMatch = 321 313
-### Ds -> π π π
-431:onIfMatch = 211 211 211
-### Ds -> K π π
-431:onIfMatch = 321 211 211
-### Ds -> K K K
-431:onIfMatch = 321 321 321
-
-### Λc -> p K π
-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 π π
-4122:onIfMatch = 2212 211 211
-### Λc -> pK0s
-4122:onIfMatch = 2212 311 
-### Λc -> p K K
-4122:onIfMatch = 2212 321 321
-
-### Ξc+ -> p K- π+
-4232:onIfMatch = 2212 321 211
-### Ξc+ -> Ξ- π+ π+
-4232:onIfMatch = 3312 211 211
-### Ξc+ -> p φ
-4232:onIfMatch = 2212 333
diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkg.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkg.cfg
new file mode 100644
index 000000000..778aaddba
--- /dev/null
+++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkg.cfg
@@ -0,0 +1,240 @@
+### authors: Stefano Politano (stefano.politano@cern.ch), Mattia Faggin (mattia.faggin@cern.ch), Fabrizio Grosa (fabrizio.grosa@cern.ch), Marcello Di Costanzo (marcello.di.costanzo@cern.ch)
+### last update:  April 2025
+
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 13600. 		# GeV
+
+### processes
+SoftQCD:inelastic on		# all inelastic processes
+
+### 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
+
+### HF decays
+### BR are set to yield 20% of signal in golden channel and uniform abundance of corr. bkg channels (weighted by BR from PDG) 
+###
+### D0 decays
+### D0 -> K- π+ (20%)
+421:oneChannel = 1 0.20 0 -321 211              ### D0 -> K- π+ 3.94% 
+### D0 -> K- π+ π0 (20%)
+421:addChannel = 1 0.01 0 -321 211 111          ### D0 -> K- π+ π0 (non-resonant) 1.15% (e.g. 115/(115+231+195+1120)*0.2)
+421:addChannel = 1 0.14 0  213 -321             ### D0 -> rho+ K- 11.2%
+421:addChannel = 1 0.02 0 -313 111              ### D0 -> antiK*0(892) π0 1.95%
+421:addChannel = 1 0.03 0 -323 211              ### D0 -> K*-(892) π+ 2.31%
+### D0 -> π- π+ (20%)
+421:addChannel = 1 0.20 0 -211 211              ### D0 -> π- π+ (non-resonant) 1.0e-4
+### D0 -> π- π+ π0 (20%)
+421:addChannel = 1 0.14 0  213 -211             ### D0 -> rho+ π- 1.01%
+421:addChannel = 1 0.06 0 -211 211 111          ### D0 -> π- π+ π0 (non-resonant) 1.3e-4
+### D0 -> K- K+ (20%)
+421:addChannel = 1 0.20 0 -321 321              ### D0 -> K- K+ (non-resonant) 4.08e-3
+
+### D+ decays
+### D+ -> K- π+ π+ (28%)
+411:oneChannel = 1 0.28 0 -321 211 211          ### D+ -> K- π+ π+ 9.38%
+### D+ -> K- π+ π+ π0 (small, 5%)
+411:addChannel = 1 0.05 0 -321 211 211 111      ### D+ -> K- π+ π+ π0 6.25%
+### D+ -> K- K+ π+ (33%)
+411:addChannel = 1 0.085 0 -313 321              ### D+ -> K*0(892) K+ 2.49e-3
+411:addChannel = 1 0.062 0 -10311 321            ### D+ -> antiK*0(1430) K+ 1.82e-3
+411:addChannel = 1 0.092 0  333 211              ### D+ -> φ π+ 2.69e-3
+411:addChannel = 1 0.091 0  -321 321 211         ### D+ -> K- K+ π+ (non-resonant) 2.68e-3 
+### D+ -> π- π+ π+ (34%)
+411:addChannel = 1 0.09 0  113 211               ### D+ -> rho0 π+ 8.4e-4
+411:addChannel = 1 0.05 0  225 211               ### D+ -> f2(1270) π+ 4.6e-4
+411:addChannel = 1 0.20 0 -211 211 211           ### D+ -> π- π+ π+ (non-resonant) 1.0e-4
+
+### Ds+ decays
+### Ds+ -> K- K+ π+ (20%)
+431:oneChannel = 1 0.100 0  333 211             ### Ds+ -> φ(1020) π+ 2.21%
+431:addChannel = 1 0.100 0 -313 321             ### Ds+ -> antiK*(892) K+ 2.58%
+### Ds+ -> K- K+ π+ π0 (small, 5%)
+431:addChannel = 1 0.050 0 333 213              ### Ds+ -> φ(1020) rho 5.50%
+### Ds+ -> π- π+ π+ (25%)
+431:addChannel = 1 0.005 0  113 211             ### Ds+ -> rho π+ 1.1e-4
+431:addChannel = 1 0.005 0  225 211             ### Ds+ -> f2(1270) π+ 1.4e-3
+431:addChannel = 1 0.240 0 -211 211 211         ### Ds+ -> π- π+ π+ 9.12e-3 (s-wave)
+### Ds+ -> π- K+ π+ (25%)
+431:addChannel = 1 0.07 0  313 211             ### Ds+ -> K*(892)0 π+ 1.67e-3
+431:addChannel = 1 0.05 0  10221 321           ### Ds+ -> f0(1370) K+ 1.2e-3
+431:addChannel = 1 0.09 0  113 321             ### Ds+ -> rho0 K+ 2.17e-3
+431:addChannel = 1 0.04 0 -211 321 211         ### Ds+ -> π- K+ π+ (non-resonant) 1.16-3
+### Ds+ -> π+ π- π+ π0 (25%)
+431:addChannel = 1 0.250 0  221 211             ### Ds+ -> eta π+ -> π0 π+ π+ π- (affects D+ golden channel) 
+
+## Λc decays
+### Λc+ -> p K- π+ (20%, BR set up to 20%)
+4122:oneChannel = 1 0.080 0 2212 -321 211       ### Λc+ -> p K- π+ (non-resonant) 3.5%
+4122:addChannel = 1 0.045 100 2212 -313         ### Λc+ -> p K*0(892) 1.96%
+4122:addChannel = 1 0.025 100 2224 -321         ### Λc+ -> Delta++ K- 1.08%
+4122:addChannel = 1 0.050 100 102134 211        ### Λc+ -> Lambda(1520) K- 2.20e-3
+### Λc+ -> p K- π+ π0 (small, 5%)
+4122:oneChannel = 1 0.050 0 2212 -321 211 111   ### Λc+ -> p K- π+ π0 (non-resonant) 4.6%
+### Λc+ -> p π- π+ (25%)
+4122:addChannel = 1 0.250 0 2212 -211 211       ### Λc+ -> p π+ π+ 4.59%
+### Λc+ -> p K- K+ (25%)
+4122:addChannel = 1 0.250 0 2212 333            ### Λc+ -> p phi 1.06%
+
+### Ξc+ decays (unfiormly distributed)
+### Ξc+ -> p K- π+ (34%)
+4232:oneChannel = 1 0.17 0 2212 -321 211        ### Ξc+ -> p K- π+
+4232:oneChannel = 1 0.17 0 2212 -313            ### Ξc+ -> p antiK*0(892)
+### Ξc+ -> p K- K+ (33%)
+4232:addChannel = 1 0.33 0 2212 333             ### Ξc+ -> p φ
+### Ξc+ -> sigma+ π+ π- (33%)
+4232:addChannel = 1 0.33 0 3222 -211 211        ### Ξc+ -> sigma+ π- π+
+
+# Allow the decay of resonances in the decay chain
+### K*0(892) -> K- π+
+313:onMode = off
+313:onIfAll = 321 211
+### K*(892)+ -> K- π0
+323:onMode = off
+323:onIfAll = 321 111
+### K*(1430)0 -> K- π0
+10311:onMode = off
+10311:onIfAll = 321 211
+### rho+ -> π+ π0
+213:onMode = off
+213:onIfAll = 211 111
+### φ -> K+ K-
+333:onMode = off
+333:onIfAll = 321 321
+### rho0 -> π+ π-
+113:onMode = off
+113:onIfAll = 211 211
+### f2(1270) -> π+ π-
+225:onMode = off
+225:onIfAll = 211 211
+### f0(1370) -> π+ π-
+10221:onMode = off
+10221:onIfAll = 211 211
+### eta -> π+ π-
+221:onMode = off
+221:onIfAll =  111 211 211
+### for Λc -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Λc -> Lambda(1520) K-
+102134:onMode = off
+102134:onIfAll = 2212 321
+
+### Switch off all decay channels
+411:onMode = off
+421:onMode = off
+431:onMode = off
+4122:onMode = off
+4232:onMode = off
+
+# Allow the decay of HF
+### D0 -> K π
+421:onIfMatch = 321 211
+### D0 -> K π π0
+421:onIfMatch = 321 211 111
+### D0 -> rho K
+421:onIfMatch = 213 321
+### D0 -> antiK*0 π0
+421:onIfMatch = 313 111
+### D0 -> K*- π+
+421:onIfMatch = 323 211 
+### D0 -> π π
+421:onIfMatch = 211 211
+### D0 -> rho+ π-
+421:onIfMatch = 213 -211
+### D0 -> π π π0
+421:onIfMatch = 211 211 111
+### D0 -> K K
+421:onIfMatch = 321 321
+
+### D+/- -> K π π 
+411:onIfMatch = 321 211 211
+### D+/- -> K π π 
+411:onIfMatch = 321 211 211 111
+### D+/- -> K* K
+411:onIfMatch = 313 321
+### D+/- -> antiK* π
+411:onIfMatch = 10311 321
+### D+/- -> φ π
+411:onIfMatch = 333 211
+### D+/- -> π π π
+411:onIfMatch = 211 211 211
+### D+/- -> f2(1270) π
+411:onIfMatch = 225 211
+### D+/- -> rho π
+411:onIfMatch = 113 211
+
+### Ds -> φ π
+431:onIfMatch = 333 211
+### Ds -> K* K
+431:onIfMatch = 313 321
+### Ds -> φ rho+
+431:onIfMatch = 333 213
+### Ds -> π π π
+431:onIfMatch = 211 211 211
+### Ds -> f2 π
+431:onIfMatch = 225 211
+### Ds -> rho π
+431:onIfMatch = 113 211
+### Ds -> K π π
+431:onIfMatch = 321 211 211
+### Ds -> K*(892)0 π+
+431:onIfMatch = 313 211
+### Ds -> f0(1370) K+
+431:onIfMatch = 10221 321
+### Ds -> rho0 K+
+431:onIfMatch = 113 321
+### Ds -> eta π+
+431:onIfMatch = 221 211
+
+### Λc -> p K π
+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+ -> sigma- π+ π+
+4232:onIfMatch = 3222 211 211
+### Ξc+ -> p φ
+4232:onIfMatch = 2212 333
\ No newline at end of file
diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV_CorrBkg.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV_CorrBkg.cfg
new file mode 100644
index 000000000..592238787
--- /dev/null
+++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV_CorrBkg.cfg
@@ -0,0 +1,248 @@
+### authors: Fabrizio Grosa (fabrizio.grosa@cern.ch)
+###          Mattia Faggin (mattia.faggin@cern.ch)
+###          Stefano Politano (stefano.politano@cern.ch)
+###50        Marcello Di Costanzo (marcello.di.costanzo@cern.ch)
+### last update:  June 2024
+
+### 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
+
+### HF decays
+### BR are set to yield 20% of signal in golden channel and uniform abundance of corr. bkg channels (weighted by BR from PDG) 
+###
+### D0 decays
+### D0 -> K- π+ (20%)
+421:oneChannel = 1 0.20 0 -321 211              ### D0 -> K- π+ 3.94% 
+### D0 -> K- π+ π0 (20%)
+421:addChannel = 1 0.01 0 -321 211 111          ### D0 -> K- π+ π0 (non-resonant) 1.15% (e.g. 115/(115+231+195+1120)*0.2)
+421:addChannel = 1 0.14 0  213 -321             ### D0 -> rho+ K- 11.2%
+421:addChannel = 1 0.02 0 -313 111              ### D0 -> antiK*0(892) π0 1.95%
+421:addChannel = 1 0.03 0 -323 211              ### D0 -> K*-(892) π+ 2.31%
+### D0 -> π- π+ (20%)
+421:addChannel = 1 0.20 0 -211 211              ### D0 -> π- π+ (non-resonant) 1.0e-4
+### D0 -> π- π+ π0 (20%)
+421:addChannel = 1 0.14 0  213 -211             ### D0 -> rho+ π- 1.01%
+421:addChannel = 1 0.06 0 -211 211 111          ### D0 -> π- π+ π0 (non-resonant) 1.3e-4
+### D0 -> K- K+ (20%)
+421:addChannel = 1 0.20 0 -321 321              ### D0 -> K- K+ (non-resonant) 4.08e-3
+
+### D+ decays
+### D+ -> K- π+ π+ (28%)
+411:oneChannel = 1 0.28 0 -321 211 211          ### D+ -> K- π+ π+ 9.38%
+### D+ -> K- π+ π+ π0 (small, 5%)
+411:addChannel = 1 0.05 0 -321 211 211 111      ### D+ -> K- π+ π+ π0 6.25%
+### D+ -> K- K+ π+ (33%)
+411:addChannel = 1 0.085 0 -313 321              ### D+ -> K*0(892) K+ 2.49e-3
+411:addChannel = 1 0.062 0 -10311 321            ### D+ -> antiK*0(1430) K+ 1.82e-3
+411:addChannel = 1 0.092 0  333 211              ### D+ -> φ π+ 2.69e-3
+411:addChannel = 1 0.091 0  -321 321 211         ### D+ -> K- K+ π+ (non-resonant) 2.68e-3 
+### D+ -> π- π+ π+ (34%)
+411:addChannel = 1 0.09 0  113 211               ### D+ -> rho0 π+ 8.4e-4
+411:addChannel = 1 0.05 0  225 211               ### D+ -> f2(1270) π+ 4.6e-4
+411:addChannel = 1 0.20 0 -211 211 211           ### D+ -> π- π+ π+ (non-resonant) 1.0e-4
+
+### Ds+ decays
+### Ds+ -> K- K+ π+ (20%)
+431:oneChannel = 1 0.100 0  333 211             ### Ds+ -> φ(1020) π+ 2.21%
+431:addChannel = 1 0.100 0 -313 321             ### Ds+ -> antiK*(892) K+ 2.58%
+### Ds+ -> K- K+ π+ π0 (small, 5%)
+431:addChannel = 1 0.050 0 333 213              ### Ds+ -> φ(1020) rho 5.50%
+### Ds+ -> π- π+ π+ (25%)
+431:addChannel = 1 0.005 0  113 211             ### Ds+ -> rho π+ 1.1e-4
+431:addChannel = 1 0.005 0  225 211             ### Ds+ -> f2(1270) π+ 1.4e-3
+431:addChannel = 1 0.240 0 -211 211 211         ### Ds+ -> π- π+ π+ 9.12e-3 (s-wave)
+### Ds+ -> π- K+ π+ (25%)
+431:addChannel = 1 0.07 0  313 211             ### Ds+ -> K*(892)0 π+ 1.67e-3
+431:addChannel = 1 0.05 0  10221 321           ### Ds+ -> f0(1370) K+ 1.2e-3
+431:addChannel = 1 0.09 0  113 321             ### Ds+ -> rho0 K+ 2.17e-3
+431:addChannel = 1 0.04 0 -211 321 211         ### Ds+ -> π- K+ π+ (non-resonant) 1.16-3
+### Ds+ -> π+ π- π+ π0 (25%)
+431:addChannel = 1 0.250 0  221 211             ### Ds+ -> eta π+ -> π0 π+ π+ π- (affects D+ golden channel) 
+
+## Λc decays
+### Λc+ -> p K- π+ (20%, BR set up to 20%)
+4122:oneChannel = 1 0.080 0 2212 -321 211       ### Λc+ -> p K- π+ (non-resonant) 3.5%
+4122:addChannel = 1 0.045 100 2212 -313         ### Λc+ -> p K*0(892) 1.96%
+4122:addChannel = 1 0.025 100 2224 -321         ### Λc+ -> Delta++ K- 1.08%
+4122:addChannel = 1 0.050 100 102134 211        ### Λc+ -> Lambda(1520) K- 2.20e-3
+### Λc+ -> p K- π+ π0 (small, 5%)
+4122:oneChannel = 1 0.050 0 2212 -321 211 111   ### Λc+ -> p K- π+ π0 (non-resonant) 4.6%
+### Λc+ -> p π- π+ (25%)
+4122:addChannel = 1 0.250 0 2212 -211 211       ### Λc+ -> p π+ π+ 4.59%
+### Λc+ -> p K- K+ (25%)
+4122:addChannel = 1 0.250 0 2212 333            ### Λc+ -> p phi 1.06%
+
+### Ξc+ decays (unfiormly distributed)
+### Ξc+ -> p K- π+ (34%)
+4232:oneChannel = 1 0.17 0 2212 -321 211        ### Ξc+ -> p K- π+
+4232:oneChannel = 1 0.17 0 2212 -313            ### Ξc+ -> p antiK*0(892)
+### Ξc+ -> p K- K+ (33%)
+4232:addChannel = 1 0.33 0 2212 333             ### Ξc+ -> p φ
+### Ξc+ -> sigma+ π+ π- (33%)
+4232:addChannel = 1 0.33 0 3222 -211 211        ### Ξc+ -> sigma+ π- π+
+
+# Allow the decay of resonances in the decay chain
+### K*0(892) -> K- π+
+313:onMode = off
+313:onIfAll = 321 211
+### K*(892)+ -> K- π0
+323:onMode = off
+323:onIfAll = 321 111
+### K*(1430)0 -> K- π0
+10311:onMode = off
+10311:onIfAll = 321 211
+### rho+ -> π+ π0
+213:onMode = off
+213:onIfAll = 211 111
+### φ -> K+ K-
+333:onMode = off
+333:onIfAll = 321 321
+### rho0 -> π+ π-
+113:onMode = off
+113:onIfAll = 211 211
+### f2(1270) -> π+ π-
+225:onMode = off
+225:onIfAll = 211 211
+### f0(1370) -> π+ π-
+10221:onMode = off
+10221:onIfAll = 211 211
+### eta -> π+ π-
+221:onMode = off
+221:onIfAll =  111 211 211
+### for Λc -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Λc -> Lambda(1520) K-
+102134:onMode = off
+102134:onIfAll = 2212 321
+
+### Switch off all decay channels
+411:onMode = off
+421:onMode = off
+431:onMode = off
+4122:onMode = off
+4232:onMode = off
+
+# Allow the decay of HF
+### D0 -> K π
+421:onIfMatch = 321 211
+### D0 -> K π π0
+421:onIfMatch = 321 211 111
+### D0 -> rho K
+421:onIfMatch = 213 321
+### D0 -> antiK*0 π0
+421:onIfMatch = 313 111
+### D0 -> K*- π+
+421:onIfMatch = 323 211 
+### D0 -> π π
+421:onIfMatch = 211 211
+### D0 -> rho+ π-
+421:onIfMatch = 213 -211
+### D0 -> π π π0
+421:onIfMatch = 211 211 111
+### D0 -> K K
+421:onIfMatch = 321 321
+
+### D+/- -> K π π 
+411:onIfMatch = 321 211 211
+### D+/- -> K π π 
+411:onIfMatch = 321 211 211 111
+### D+/- -> K* K
+411:onIfMatch = 313 321
+### D+/- -> antiK* π
+411:onIfMatch = 10311 321
+### D+/- -> φ π
+411:onIfMatch = 333 211
+### D+/- -> π π π
+411:onIfMatch = 211 211 211
+### D+/- -> f2(1270) π
+411:onIfMatch = 225 211
+### D+/- -> rho π
+411:onIfMatch = 113 211
+
+### Ds -> φ π
+431:onIfMatch = 333 211
+### Ds -> K* K
+431:onIfMatch = 313 321
+### Ds -> φ rho+
+431:onIfMatch = 333 213
+### Ds -> π π π
+431:onIfMatch = 211 211 211
+### Ds -> f2 π
+431:onIfMatch = 225 211
+### Ds -> rho π
+431:onIfMatch = 113 211
+### Ds -> K π π
+431:onIfMatch = 321 211 211
+### Ds -> K*(892)0 π+
+431:onIfMatch = 313 211
+### Ds -> f0(1370) K+
+431:onIfMatch = 10221 321
+### Ds -> rho0 K+
+431:onIfMatch = 113 321
+### Ds -> eta π+
+431:onIfMatch = 221 211
+
+### Λc -> p K π
+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+ -> sigma- π+ π+
+4232:onIfMatch = 3222 211 211
+### Ξc+ -> p φ
+4232:onIfMatch = 2212 333
\ No newline at end of file