Add configurations for MC production for SigmaC background studies.

MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C

@@ -33,22 +33,26 @@ public:
     mHadronPdgList = hadronPdgList;
     mPartPdgToReplaceList = partPdgToReplaceList;
     mFreqReplaceList = freqReplaceList;
-    // Ds1*(2700), Ds1*(2860), Ds3*(2860), Xic(3055)+, Xic(3080)+, Xic(3055)0, Xic(3080)0
-    mCustomPartPdgs = {30433, 40433, 437, 4315, 4316, 4325, 4326};
+    // Ds1*(2700), Ds1*(2860), Ds3*(2860), Xic(3055)+, Xic(3080)+, Xic(3055)0, Xic(3080)0, LambdaC(2625), LambdaC(2595)
+    mCustomPartPdgs = {30433, 40433, 437, 4315, 4316, 4325, 4326, 4124, 14122};
     mCustomPartMasses[30433] = 2.714f;
     mCustomPartMasses[40433] =  2.859f;
     mCustomPartMasses[437] = 2.860f;
     mCustomPartMasses[4315] = 3.0590f;
     mCustomPartMasses[4316] = 3.0799f;
     mCustomPartMasses[4325] = 3.0559f;
     mCustomPartMasses[4326] = 3.0772f;
+    mCustomPartMasses[4124] = 2.62810f;
+    mCustomPartMasses[14122] = 2.59225f;
     mCustomPartWidths[30433] = 0.122f;
     mCustomPartWidths[40433] =  0.160f;
     mCustomPartWidths[437] = 0.053f;
     mCustomPartWidths[4315] = 0.0064f;
     mCustomPartWidths[4316] = 0.0056f;
     mCustomPartWidths[4325] = 0.0078f;
     mCustomPartWidths[4326] = 0.0036f;
+    mCustomPartWidths[4124] = 0.00052f;
+    mCustomPartWidths[14122] = 0.0026f;
     Print();
   }    
 
@@ -436,4 +440,4 @@ FairGenerator *GeneratorPythia8GapHF(int inputTriggerRatio, float yQuarkMin = -1
   myGen->setHadronRapidity(yHadronMin, yHadronMax);
 
   return myGen;
-}
+}

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.ini

@@ -0,0 +1,10 @@
+
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5, -1.5, 1.5, {413}, {{413, 14122}, {413, 4124}}, {0.5, 0.5})
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkgSigmaC.cfg
+includePartonEvent=false
+### not needed for jet studies, hence no need to keep parton event

MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.C

@@ -0,0 +1,200 @@
+#include "TFile.h"
+#include "TTree.h"
+
+#include <string>
+#include <iostream>
+
+int External() {
+    std::string path{"/home/mattia/Documenti/cernbox/Documents/PostDoc/D2H/MC/corrBkgSigmaC/tf1/genevents_Kine.root"};
+
+    int checkPdgQuarkOne{4};
+    int checkPdgQuarkTwo{5};
+    float ratioTrigger = 1./5.; // one event triggered out of 5
+    std::array<std::array<int, 2>, 2> pdgReplParticles = {std::array{413, 14122}, std::array{413, 4124}};
+    std::array<std::array<int, 2>, 2> pdgReplPartCounters = {std::array{0, 0}, std::array{0, 0}};
+    std::array<float, 2> freqRepl = {0.5, 0.5};
+    std::map<int, int> sumOrigReplacedParticles = {{413, 0}};
+
+    std::array<int, 2> checkPdgHadron{14122, 4124};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted (!) pdg of daughters
+        //{14122, {{4222, -211}, {4112, 211}, {4122, 211, -211}}}, // Lc(2595)+
+        //{4124, {{4222, -211}, {4112, 211}, {4122, 211, -211}}} // Lc(2625)+
+        {14122, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}}, // Lc(2595)+
+        {4124, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}} // Lc(2625)+
+    };
+
+    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 nSignals{}, nSignalGoodDecay{};
+    auto nEvents = tree->GetEntries();
+
+    for (int i = 0; i < nEvents; i++) {
+
+        std::cout << std::endl;
+
+        tree->GetEntry(i);
+
+        // check subgenerator information
+        int subGeneratorId{-1};
+        if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
+            bool isValid = false;
+            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();
+            auto absPdg = std::abs(pdg);
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), absPdg) != checkPdgHadron.end()) { // found signal
+                nSignals++; // count signal PDG
+                std::cout << "==> signal " << absPdg << " found!" << std::endl;
+
+                if (subGeneratorId == checkPdgQuarkOne) { // replacement only for prompt  ---> BUT ALSO NON-PROMPT D* SEEM TO BE REPLACED
+                    for (int iRepl{0}; iRepl<2; ++iRepl) {
+                        if (absPdg == pdgReplParticles[iRepl][0]) {
+                            pdgReplPartCounters[iRepl][0]++;
+                            sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
+                        } else if (absPdg == pdgReplParticles[iRepl][1]) {
+                            pdgReplPartCounters[iRepl][1]++;
+                            sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
+                        }
+                    }
+                } else if (subGeneratorId == checkPdgQuarkTwo) {
+                    std::cout << "   NB: we have a " << absPdg << " also in event with quark " << checkPdgQuarkTwo << std::endl;
+                    std::cout << "          ### mother indices: ";
+                    int idFirstMother = track.getMotherTrackId();
+                    int idSecondMother = track.getSecondMotherTrackId();
+                    std::vector<int> motherIds = {};
+                    for(int i=idFirstMother; i<=idSecondMother; i++) {
+                        std::cout << i << " ";
+                        motherIds.push_back(i);
+                    }
+                    bool partonicEventOn = false;
+                    if(motherIds != std::vector<int>{-1, -1}) {
+                        std::cout << "The " << absPdg << " particle has mothers. This should mean that it comes directly from parton hadronization, and that the partonic event was kept in the MC production " << std::endl;
+                        partonicEventOn = true;
+                    }
+                    std::cout << "          ### mother PDG codes: ";
+                    std::vector<int> motherPdgCodes = {};
+                    if(partonicEventOn) {
+                        for(int i=idFirstMother; i<=idSecondMother; i++) {
+                            motherPdgCodes.push_back(tracks->at(i).GetPdgCode());
+                            std::cout << motherPdgCodes.back() << " ";
+                        }
+
+                        /// check that among the mothers there is a c/cbar quark
+                        /// This means that the charm hadron comes from the c-quark hadronization, where the c/cbar quark
+                        /// comes from a c-cbar pair present in the current event, tagged with a b-bbar (e.g. double-parton scattering)
+                        if(std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4) == motherPdgCodes.end()) {
+                            /// if the partinc event is not really saved and we arrive here, it means that  motherIds != {-1, -1} because 
+                            /// the hadron comes from the decay of a beauty hadron. This can happen if and only if this is not a replaced one (i.e. native from Lambdab0 decay)
+                            if (std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 5122) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -5122) == motherPdgCodes.end()) {
+                                std::cerr << "The particle " << absPdg << " does not originate neither from a c/c-bar quark (replaced) nor from a Lambda_b0 decay. There is something wrong, aborting..." << std::endl;
+                            return 1;
+                            }
+                        }
+                    }
+                    std::cout << std::endl;
+
+                    /// only if we arrive here it means that everything is ok, and we can safely update the counters for the final statistics
+                    for (int iRepl{0}; iRepl<2; ++iRepl) {
+                        if (absPdg == pdgReplParticles[iRepl][0]) {
+                            pdgReplPartCounters[iRepl][0]++;
+                            sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
+                        } else if (absPdg == pdgReplParticles[iRepl][1]) {
+                            pdgReplPartCounters[iRepl][1]++;
+                            sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
+                        }
+                    }
+                }
+
+
+                std::vector<int> pdgsDecay{};
+                std::vector<int> pdgsDecayAntiPart{};
+                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::cout << "   -- daughter " << j << ": " << pdgDau << std::endl;
+                        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++;
+                        std::cout << "     !!! GOOD DECAY FOUND !!!" << std::endl;
+                        break;
+                    }
+                }
+            }
+        } // end loop over tracks
+    }
+
+    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 <<"# 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;
+    }
+
+    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;
+    }
+
+    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]];
+            } 
+            std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0] << " into " << pdgReplParticles[iRepl][1] << " is " << fracMeas <<" (expected "<< freqRepl[iRepl] << ")\n";
+            return 1;    
+        }
+    }
+
+    return 0;
+}

MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkgSigmaC.cfg

@@ -0,0 +1,84 @@
+### authors: Mattia Faggin (mattia.faggin@cern.ch)
+### last update:  July 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
+
+# switch off charm-hadron decays
+4122:onMode = off
+14122:onMode = off  ### Lambda_c(2595)+ already present in PYTHIA
+4124:onMode = off   ### Lambda_c(2625)+ already present in PYTHIA
+
+## Λc decays
+### Λc+ -> p K- π+
+4122:oneChannel = 1 0.0350 0 2212 -321 211       ### Λc+ -> p K- π+ (non-resonant) 3.5%
+4122:addChannel = 1 0.0196 100 2212 -313         ### Λc+ -> p K*0(892) 1.96%
+4122:addChannel = 1 0.0108 100 2224 -321         ### Λc+ -> Delta++ K- 1.08%
+4122:addChannel = 1 0.0022 100 102134 211        ### Λc+ -> Lambda(1520) K- 2.20e-3
+4122:onIfMatch = 2212 321 211
+4122:onIfMatch = 2212 313
+4122:onIfMatch = 2224 321
+4122:onIfMatch = 102134 211
+### Λc(2595)+
+14122:oneChannel = 1 0.24 0 4222 -211        ### Λc(2595)+ -> Σc(2455)++ π- (24%)
+14122:addChannel = 1 0.24 0 4112 211         ### Λc(2595)+ -> Σc(2455)0 π+ (24%)
+14122:addChannel = 1 0.18 0 4122 211 -211    ### Λc(2595)+ -> Λc+ π+ π- (18%)
+14122:onIfMatch = 4222 211        ### Λc(2595)+ -> Σc(2455)++ π-
+14122:onIfMatch = 4112 211        ### Λc(2595)+ -> Σc(2455)0 π+ 
+14122:onIfMatch = 4122 211 211    ### Λc(2595)+ -> Λc+ π+ π-
+### Λc(2625)+
+4124:oneChannel = 1 0.24 0 4222 -211        ### Λc(2592)+ -> Σc(2455)++ π- (24%)
+4124:addChannel = 1 0.24 0 4112 211         ### Λc(2592)+ -> Σc(2455)0 π+ (24%)
+4124:addChannel = 1 0.18 0 4122 211 -211    ### Λc(2592)+ -> Λc+ π+ π- (18%)
+4124:onIfMatch = 4222 211        ### Λc(2625)+ -> Σc(2455)++ π-
+4124:onIfMatch = 4112 211        ### Λc(2625)+ -> Σc(2455)0 π+ 
+4124:onIfMatch = 4122 211 211    ### Λc(2625)+ -> Λc+ π+ π-
+
+# Correct decay lengths (wrong in PYTHIA8 decay table)
+# Lb
+5122:tau0 = 0.4390
+# Xic0
+4132:tau0 = 0.0455
+# OmegaC
+4332:tau0 = 0.0803
+
+# Correct Breit-Wigner width
+4124:mWidth = 0.00052 # <0.52 MeV
+
+# Force the decay of resonances in the decay chain
+### K*0(892) -> K- π+
+313:onMode = off
+313:onIfAll = 321 211
+### for Λc -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Λc -> Lambda(1520) K-
+102134:onMode = off
+102134:onIfAll = 2212 321