From ca18456bea30c26f1e6fa0746e0e20434237e246 Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Wed, 19 Nov 2025 15:22:41 +0100 Subject: [PATCH 1/9] Add configuration for GeneratorHF D2H Sigmac+/D* replaced by Omegac0/Xic0 --- .../ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini 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..ba5503948 --- /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{}, std::vector{}, std::vector>{{{423,4132}},{{4212,4332}}}, std::vector{1,1}) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg +includePartonEvent=true From b946646e7f8b6af010363c28979495b9e46621d9 Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Wed, 19 Nov 2025 15:31:58 +0100 Subject: [PATCH 2/9] Add Pythia8 configuration for charm hadronic decays This configuration file sets up the Pythia8 generator for charm hadronic processes, specifying decay modes and parameters for various charm and bottom hadrons. It includes detailed channels for D, Ds, Lambda_c, Xi_c, and Omega_c decays, as well as adjustments for resonance decays. --- ...eplaced_with_decays_Mode2_hardQCD_5TeV.cfg | 296 ++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100644 MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg new file mode 100644 index 000000000..90ae5b528 --- /dev/null +++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg @@ -0,0 +1,296 @@ +### 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 golden charm hadrons decay modes for D2H studies +# HF decays +### BR are set to yield 50% of signal in golden channel and uniform abundance of corr. bkg channels (weighted by BR from PDG) +### +### D0 decays +### D0 -> K- π+ (50%) +421:oneChannel = 1 0.50000 0 -321 211 ### D0 -> K- π+ 3.94% +### D0 -> K- π+ π0 (12.50%) +421:addChannel = 1 0.00625 0 -321 211 111 ### D0 -> K- π+ π0 (non-resonant) 1.15% (e.g. 115/(115+231+195+1120)*0.2) +421:addChannel = 1 0.08750 0 213 -321 ### D0 -> rho+ K- 11.2% +421:addChannel = 1 0.01250 0 -313 111 ### D0 -> antiK*0(892) π0 1.95% +421:addChannel = 1 0.01875 0 -323 211 ### D0 -> K*-(892) π+ 2.31% +### D0 -> π- π+ (12.50%) +421:addChannel = 1 0.12500 0 -211 211 ### D0 -> π- π+ (non-resonant) 1.0e-4 +### D0 -> π- π+ π0 (12.50%) +421:addChannel = 1 0.08750 0 213 -211 ### D0 -> rho+ π- 1.01% +421:addChannel = 1 0.03750 0 -211 211 111 ### D0 -> π- π+ π0 (non-resonant) 1.3e-4 +### D0 -> K- K+ (12.50%) +421:addChannel = 1 0.12500 0 -321 321 ### D0 -> K- K+ (non-resonant) 4.08e-3 + +### D+ decays +### D+ -> K- π+ π+ (50%) +411:oneChannel = 1 0.40189 0 -321 211 211 ### D+ -> K- π+ π+ 9.38% +411:addChannel = 1 0.05356 0 -10311 211 ### D+ -> antiK*0(1430) π+ 1.25% +411:addChannel = 1 0.04455 0 -313 211 ### D+ -> K*0(892) pi+ 1.04% +### D+ -> K- π+ π+ π0 (small, 3%) +411:addChannel = 1 0.03000 0 -321 211 211 111 ### D+ -> K- π+ π+ π0 6.25% +### D+ -> K- K+ π+ (36.00%, set 25% for D+ -> φ π+, 11% for the rest) +411:addChannel = 1 0.25000 0 333 211 ### D+ -> φ π+ 2.69e-3 !needed for signal +411:addChannel = 1 0.03929 0 -313 321 ### D+ -> K*0(892) K+ 2.49e-3 +411:addChannel = 1 0.02865 0 -10311 321 ### D+ -> antiK*0(1430) K+ 1.82e-3 +411:addChannel = 1 0.04206 0 -321 321 211 ### D+ -> K- K+ π+ (non-resonant) 2.68e-3 +### D+ -> π- π+ π+ (11.00%) +411:addChannel = 1 0.02911 0 113 211 ### D+ -> rho0 π+ 8.4e-4 +411:addChannel = 1 0.01618 0 225 211 ### D+ -> f2(1270) π+ 4.6e-4 +411:addChannel = 1 0.06471 0 -211 211 211 ### D+ -> π- π+ π+ (non-resonant) 1.0e-4 + +### Ds+ decays +### Ds+ -> K- K+ π+ (50%) +431:oneChannel = 1 0.50000 0 333 211 ### Ds+ -> φ(1020) π+ 2.21% +431:addChannel = 1 0.15000 0 -313 321 ### Ds+ -> antiK*(892) K+ 2.58% +### Ds+ -> K- K+ π+ π0 (small, 2%) +431:addChannel = 1 0.02000 0 333 213 ### Ds+ -> φ(1020) rho 5.50% +### Ds+ -> π- π+ π+ (11.00%) +431:addChannel = 1 0.00220 0 113 211 ### Ds+ -> rho π+ 1.1e-4 +431:addChannel = 1 0.00220 0 225 211 ### Ds+ -> f2(1270) π+ 1.4e-3 +431:addChannel = 1 0.10560 0 -211 211 211 ### Ds+ -> π- π+ π+ 9.12e-3 (s-wave) +### Ds+ -> π- K+ π+ (11.00%) +431:addChannel = 1 0.03080 0 313 211 ### Ds+ -> K*(892)0 π+ 1.67e-3 +431:addChannel = 1 0.02200 0 10221 321 ### Ds+ -> f0(1370) K+ 1.2e-3 +431:addChannel = 1 0.03960 0 113 321 ### Ds+ -> rho0 K+ 2.17e-3 +431:addChannel = 1 0.01760 0 -211 321 211 ### Ds+ -> π- K+ π+ (non-resonant) 1.16-3 +### Ds+ -> π+ π- π+ π0 (11.00%) +431:addChannel = 1 0.11000 0 221 211 ### Ds+ -> eta π+ -> π0 π+ π+ π- (affects D+ golden channel) + +## Λc decays +### Λ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 +### 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- π+ +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 +### 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 all decay channels +411:onMode = off +421:onMode = off +431:onMode = off +4122:onMode = off +4232:onMode = off +4132:onMode = off +443:onMode = off +4332: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* π +411:onIfMatch = 313 211 +### D+/- -> K* K +411:onIfMatch = 313 321 +### D+/- -> antiK* π +411:onIfMatch = 10311 211 +### D+/- -> antiK* K +411:onIfMatch = 10311 321 +### D+/- -> φ π +411:onIfMatch = 333 211 +### D+/- -> K K π +411:onIfMatch = 321 321 211 +### D+/- -> f2(1270) π +411:onIfMatch = 225 211 +### D+/- -> rho π +411:onIfMatch = 113 211 +### D+/- -> π π π +411:onIfMatch = 211 211 211 + +### Ds -> φ π +431:onIfMatch = 333 211 +### Ds -> K* K +431:onIfMatch = 313 321 +### Ds -> φ rho+ +431:onIfMatch = 333 213 +### Ds -> rho π +431:onIfMatch = 113 211 +### Ds -> f2 π +431:onIfMatch = 225 211 +### Ds -> π π π +431:onIfMatch = 211 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 -> K π π +431:onIfMatch = 321 211 211 +### Ds -> eta π+ +431:onIfMatch = 221 211 + +### Λ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 From 9b37b4bb80f1d532b2dc8e88a7d455624be35847 Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Wed, 19 Nov 2025 15:40:39 +0100 Subject: [PATCH 3/9] Add replacement config Updated the destructor to clean up the HF generator and modified the setupGeneratorEvHF method to include additional parameters for particle replacement. --- .../generator/generator_pythia8_embed_hf.C | 38 ++++++++++++------- 1 file changed, 25 insertions(+), 13 deletions(-) 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..95dbbdb8f 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,8 @@ using namespace Pythia8; +#include + namespace hf_generators { enum GenType : int { @@ -35,7 +37,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 +59,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 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 quarkPdgList = {}, std::vector hadronPdgList = {}, std::vector> partPdgToReplaceList = {}, std::vector 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(GeneratorPythia8GapTriggeredCharm(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList)); + mGeneratorEvHF = dynamic_cast(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(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList)); + mGeneratorEvHF = dynamic_cast(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(GeneratorPythia8GapTriggeredCharmAndBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList)); + mGeneratorEvHF = dynamic_cast(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(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList)); + mGeneratorEvHF = dynamic_cast(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList)); break; default: LOG(fatal) << "********** [GeneratorPythia8EmbedHF] bad configuration, fix it! **********"; @@ -111,7 +120,8 @@ 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); + // compute and round the number of signal events to an integer + mNumSigEvs = static_cast(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; }; @@ -142,6 +152,7 @@ bool isFromCharmOrBeauty(const int partId, std::vector const& particl std::vector initVec{partId}; arrayIds.push_back(initVec); // the first vector contains the index of the original particle int stage = 0; + while(arrayIds[-stage].size() > 0) { //LOG(info) << "### stage " << stage << ", arrayIds[-stage].size() = " << arrayIds[-stage].size(); @@ -149,6 +160,7 @@ bool isFromCharmOrBeauty(const int partId, std::vector const& particl std::vector arrayIdsStage{}; for (auto& iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage + const TParticle& partStage = particles.at(iPart); // check the first mother @@ -380,34 +392,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 hadronPdgList = {}) +FairGenerator * GeneratorPythia8EmbedHFCharm(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector quarkPdgList = {}, std::vector hadronPdgList = {}, std::vector> partPdgToReplaceList = {}, std::vector 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 hadronPdgList = {}) +FairGenerator * GeneratorPythia8EmbedHFBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector quarkPdgList = {}, std::vector hadronPdgList = {}, std::vector> partPdgToReplaceList = {}, std::vector 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 hadronPdgList = {}) +FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector quarkPdgList = {}, std::vector hadronPdgList = {}, std::vector> partPdgToReplaceList = {}, std::vector 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; } From 46c370b6d3cce29a7b65e1a56523cafdc1216e7e Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Wed, 19 Nov 2025 15:47:37 +0100 Subject: [PATCH 4/9] Update replaceParticle to include status code check Refactor replaceParticle function to check status code before replacing particles. --- .../generator/generator_pythia8_gaptriggered_hf.C | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) 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 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(); From 68a5f344a8f6f2b1f2958813735f4ab097467659 Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Fri, 21 Nov 2025 11:41:33 +0100 Subject: [PATCH 5/9] Update MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini Co-authored-by: Fabrizio --- .../PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index ba5503948..b8c9f3495 100644 --- 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 @@ -4,5 +4,5 @@ fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_py funcName=GeneratorPythia8EmbedHFCharmAndBeauty(false, -1.5, 1.5, -1.5, 1.5, std::vector{}, std::vector{}, std::vector>{{{423,4132}},{{4212,4332}}}, std::vector{1,1}) [GeneratorPythia8] -config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg includePartonEvent=true From 5a0bd9a24f9f581f58661c96b5bef3d11d9731bf Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Fri, 21 Nov 2025 11:44:43 +0100 Subject: [PATCH 6/9] Add pythia8 config for XicOmegaC decays --- ...pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename MC/config/PWGHF/pythia8/generator/{pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg => pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg} (100%) diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg similarity index 100% rename from MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_replaced_with_decays_Mode2_hardQCD_5TeV.cfg rename to MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg From 4a133cc7c27adfa1492a872ffb9a607d163d219e Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Fri, 21 Nov 2025 11:45:11 +0100 Subject: [PATCH 7/9] Update MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini Co-authored-by: Fabrizio --- .../PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index b8c9f3495..3cbed6171 100644 --- 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 @@ -1,7 +1,7 @@ ### 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{}, std::vector{}, std::vector>{{{423,4132}},{{4212,4332}}}, std::vector{1,1}) +funcName=GeneratorPythia8EmbedHFCharmAndBeauty(false, -1.5, 1.5, -1.5, 1.5, std::vector{}, std::vector{}, std::vector>{{423,4132},{423,4232},{4212,4332}}, std::vector{0.5,0.5,1}) [GeneratorPythia8] config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg From 7eb8a8cb2eabb8d96a61dc9fa1dd556bd5b0af51 Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Sun, 23 Nov 2025 14:29:07 +0100 Subject: [PATCH 8/9] Keep only charm baryon decay modes Updated decay modes for charm baryons and added new channels for Omega_c decays. --- ...ic_decays_Mode2_hardQCD_5TeV_XicOmegaC.cfg | 168 ++---------------- 1 file changed, 14 insertions(+), 154 deletions(-) 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 index 90ae5b528..7c6ae8f0e 100644 --- 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 @@ -45,62 +45,10 @@ BeamRemnants:saturation 5 # OmegaC 4332:tau0 = 0.0803 -### Force golden charm hadrons decay modes for D2H studies -# HF decays -### BR are set to yield 50% of signal in golden channel and uniform abundance of corr. bkg channels (weighted by BR from PDG) -### -### D0 decays -### D0 -> K- π+ (50%) -421:oneChannel = 1 0.50000 0 -321 211 ### D0 -> K- π+ 3.94% -### D0 -> K- π+ π0 (12.50%) -421:addChannel = 1 0.00625 0 -321 211 111 ### D0 -> K- π+ π0 (non-resonant) 1.15% (e.g. 115/(115+231+195+1120)*0.2) -421:addChannel = 1 0.08750 0 213 -321 ### D0 -> rho+ K- 11.2% -421:addChannel = 1 0.01250 0 -313 111 ### D0 -> antiK*0(892) π0 1.95% -421:addChannel = 1 0.01875 0 -323 211 ### D0 -> K*-(892) π+ 2.31% -### D0 -> π- π+ (12.50%) -421:addChannel = 1 0.12500 0 -211 211 ### D0 -> π- π+ (non-resonant) 1.0e-4 -### D0 -> π- π+ π0 (12.50%) -421:addChannel = 1 0.08750 0 213 -211 ### D0 -> rho+ π- 1.01% -421:addChannel = 1 0.03750 0 -211 211 111 ### D0 -> π- π+ π0 (non-resonant) 1.3e-4 -### D0 -> K- K+ (12.50%) -421:addChannel = 1 0.12500 0 -321 321 ### D0 -> K- K+ (non-resonant) 4.08e-3 +### Force charm baryon decay modes (Ξc and Ωc) +### D mesons use default PYTHIA8 branching ratios -### D+ decays -### D+ -> K- π+ π+ (50%) -411:oneChannel = 1 0.40189 0 -321 211 211 ### D+ -> K- π+ π+ 9.38% -411:addChannel = 1 0.05356 0 -10311 211 ### D+ -> antiK*0(1430) π+ 1.25% -411:addChannel = 1 0.04455 0 -313 211 ### D+ -> K*0(892) pi+ 1.04% -### D+ -> K- π+ π+ π0 (small, 3%) -411:addChannel = 1 0.03000 0 -321 211 211 111 ### D+ -> K- π+ π+ π0 6.25% -### D+ -> K- K+ π+ (36.00%, set 25% for D+ -> φ π+, 11% for the rest) -411:addChannel = 1 0.25000 0 333 211 ### D+ -> φ π+ 2.69e-3 !needed for signal -411:addChannel = 1 0.03929 0 -313 321 ### D+ -> K*0(892) K+ 2.49e-3 -411:addChannel = 1 0.02865 0 -10311 321 ### D+ -> antiK*0(1430) K+ 1.82e-3 -411:addChannel = 1 0.04206 0 -321 321 211 ### D+ -> K- K+ π+ (non-resonant) 2.68e-3 -### D+ -> π- π+ π+ (11.00%) -411:addChannel = 1 0.02911 0 113 211 ### D+ -> rho0 π+ 8.4e-4 -411:addChannel = 1 0.01618 0 225 211 ### D+ -> f2(1270) π+ 4.6e-4 -411:addChannel = 1 0.06471 0 -211 211 211 ### D+ -> π- π+ π+ (non-resonant) 1.0e-4 - -### Ds+ decays -### Ds+ -> K- K+ π+ (50%) -431:oneChannel = 1 0.50000 0 333 211 ### Ds+ -> φ(1020) π+ 2.21% -431:addChannel = 1 0.15000 0 -313 321 ### Ds+ -> antiK*(892) K+ 2.58% -### Ds+ -> K- K+ π+ π0 (small, 2%) -431:addChannel = 1 0.02000 0 333 213 ### Ds+ -> φ(1020) rho 5.50% -### Ds+ -> π- π+ π+ (11.00%) -431:addChannel = 1 0.00220 0 113 211 ### Ds+ -> rho π+ 1.1e-4 -431:addChannel = 1 0.00220 0 225 211 ### Ds+ -> f2(1270) π+ 1.4e-3 -431:addChannel = 1 0.10560 0 -211 211 211 ### Ds+ -> π- π+ π+ 9.12e-3 (s-wave) -### Ds+ -> π- K+ π+ (11.00%) -431:addChannel = 1 0.03080 0 313 211 ### Ds+ -> K*(892)0 π+ 1.67e-3 -431:addChannel = 1 0.02200 0 10221 321 ### Ds+ -> f0(1370) K+ 1.2e-3 -431:addChannel = 1 0.03960 0 113 321 ### Ds+ -> rho0 K+ 2.17e-3 -431:addChannel = 1 0.01760 0 -211 321 211 ### Ds+ -> π- K+ π+ (non-resonant) 1.16-3 -### Ds+ -> π+ π- π+ π0 (11.00%) -431:addChannel = 1 0.11000 0 221 211 ### Ds+ -> eta π+ -> π0 π+ π+ π- (affects D+ golden channel) - -## Λc decays +## Λ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% @@ -135,40 +83,19 @@ BeamRemnants:saturation 5 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 -### 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- π+ -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 +# 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 @@ -179,82 +106,13 @@ BeamRemnants:saturation 5 3334:onMode = off 3334:onIfAll = 3122 -321 -# Switch off all decay channels -411:onMode = off -421:onMode = off -431:onMode = off +# Switch off charm baryon decay channels (D mesons use PYTHIA8 defaults) 4122:onMode = off 4232:onMode = off 4132:onMode = off -443:onMode = off 4332: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* π -411:onIfMatch = 313 211 -### D+/- -> K* K -411:onIfMatch = 313 321 -### D+/- -> antiK* π -411:onIfMatch = 10311 211 -### D+/- -> antiK* K -411:onIfMatch = 10311 321 -### D+/- -> φ π -411:onIfMatch = 333 211 -### D+/- -> K K π -411:onIfMatch = 321 321 211 -### D+/- -> f2(1270) π -411:onIfMatch = 225 211 -### D+/- -> rho π -411:onIfMatch = 113 211 -### D+/- -> π π π -411:onIfMatch = 211 211 211 - -### Ds -> φ π -431:onIfMatch = 333 211 -### Ds -> K* K -431:onIfMatch = 313 321 -### Ds -> φ rho+ -431:onIfMatch = 333 213 -### Ds -> rho π -431:onIfMatch = 113 211 -### Ds -> f2 π -431:onIfMatch = 225 211 -### Ds -> π π π -431:onIfMatch = 211 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 -> K π π -431:onIfMatch = 321 211 211 -### Ds -> eta π+ -431:onIfMatch = 221 211 - +# Allow the decay of charm baryons ### Λc -> pK0s 4122:onIfMatch = 2212 311 ### Λc -> p K- π+ π0 @@ -294,3 +152,5 @@ BeamRemnants:saturation 5 4332:onIfMatch = 3334 211 ### Omega_c -> Xi pi 4332:onIfMatch = 3312 211 +### Omega_c -> Omega- Ka+ +4332:onIfMatch = 3334 321 From bc2281006134400529f26b284825594c3187a6a9 Mon Sep 17 00:00:00 2001 From: RuiqiYin <25110190115@m.fudan.edu.cn> Date: Sun, 23 Nov 2025 14:32:58 +0100 Subject: [PATCH 9/9] Add GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced teat Test if the replacements and the decay model are correct. --- ...atorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C | 215 ++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C 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..536825e92 --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_replaced.C @@ -0,0 +1,215 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + //std::string path{"tf1/sgn_Kine.root"}; + + // Particle replacement configuration: {original, replacement}, frequencies + std::array, 3> pdgReplParticles = { + std::array{423, 4132}, // D*0 -> Xic0 + std::array{423, 4232}, // D*0 -> Xic+ + std::array{4212, 4332} // Sigmac+ -> Omegac0 + }; + std::array, 3> pdgReplPartCounters = { + std::array{0, 0}, + std::array{0, 0}, + std::array{0, 0} + }; + std::array freqRepl = {0.5, 0.5, 1.0}; + std::map sumOrigReplacedParticles = {{423, 0}, {4212, 0}}; + + // Signal hadrons to check (only final charm baryons after replacement) + std::array checkPdgHadron{4122, 4132, 4232, 4332}; + + // Expected decay channels - will be sorted automatically + // Define both particle and antiparticle versions + std::map>> 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>> 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 *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nSignals{}, nSignalGoodDecay{}; + std::map failedDecayCount{{4122, 0}, {4132, 0}, {4232, 0}, {4332, 0}}; + std::map>> 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 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; + 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 < 3; ++iRepl) { + if (sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] == 0) { + continue; // Skip if no original particles found + } + + float expectedCount = freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]; + float sigma = std::sqrt(freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]); + + if (std::abs(pdgReplPartCounters[iRepl][1] - expectedCount) > 2 * sigma) { + float fracMeas = float(pdgReplPartCounters[iRepl][1]) / sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]; + std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0] + << " into " << pdgReplParticles[iRepl][1] + << " is " << fracMeas << " (expected " << freqRepl[iRepl] << ")\n"; + return 1; + } + } + + return 0; +}