Configs for Lc resonance prod.

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

@@ -33,8 +33,8 @@ public:
     mHadronPdgList = hadronPdgList;
     mPartPdgToReplaceList = partPdgToReplaceList;
     mFreqReplaceList = freqReplaceList;
-    // 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};
+    // Ds1*(2700), Ds1*(2860), Ds3*(2860), Xic(3055)+, Xic(3080)+, Xic(3055)0, Xic(3080)0, LambdaC(2625), LambdaC(2595), LambdaC(2860), LambdaC(2880), LambdaC(2940), ThetaC(3100)
+    mCustomPartPdgs = {30433, 40433, 437, 4315, 4316, 4325, 4326, 4124, 14122, 34122, 44122, 54122, 9422111};
     mCustomPartMasses[30433] = 2.714f;
     mCustomPartMasses[40433] =  2.859f;
     mCustomPartMasses[437] = 2.860f;
@@ -44,6 +44,10 @@ public:
     mCustomPartMasses[4326] = 3.0772f;
     mCustomPartMasses[4124] = 2.62810f;
     mCustomPartMasses[14122] = 2.59225f;
+    mCustomPartMasses[34122] = 2.8561f;
+    mCustomPartMasses[44122] = 2.8816;
+    mCustomPartMasses[54124] = 2.9396f;
+    mCustomPartMasses[9422111] = 3.099f;
     mCustomPartWidths[30433] = 0.122f;
     mCustomPartWidths[40433] =  0.160f;
     mCustomPartWidths[437] = 0.053f;
@@ -53,6 +57,10 @@ public:
     mCustomPartWidths[4326] = 0.0036f;
     mCustomPartWidths[4124] = 0.00052f;
     mCustomPartWidths[14122] = 0.0026f;
+    mCustomPartWidths[34122] = 0.0676f;
+    mCustomPartWidths[44122] = 0.0056f;
+    mCustomPartWidths[54122] = 0.017f;
+    mCustomPartWidths[9422111] = 0.0000083f;
     Print();
   }    
 

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_LcResoTrigger.ini

@@ -0,0 +1,10 @@
+#NEV_TEST> 20
+### 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, {34122, 44122, 54122, 9422111}, {{4122, 34122}, {4122, 44122}, {4122, 54122}, {4122, 9422111}}, {0.1, 0.1, 0.1, 0.1})
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_LcResoTrigger.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_LcResoTrigger.C

@@ -0,0 +1,131 @@
+int ExternalLc() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuarkOne{4};  // c quark injection
+    int checkPdgQuarkTwo{5};  // b quark injection
+    float ratioTrigger = 1./5.; // one event triggered out of 5
+
+    // PDG replacements: Λc(4122) -> resonances
+    std::array<std::array<int, 2>, 4> pdgReplParticles = {
+        std::array{4122, 34122},   // Λc -> Λc(2860)
+        std::array{4122, 44122},   // Λc -> Λc(2880)
+        std::array{4122, 54122},   // Λc -> Λc(2940)
+        std::array{4122, 9422111}  // Λc -> Tc(3100)
+    };
+
+    // Counters for replacements
+    std::array<std::array<int, 2>, 4> pdgReplPartCounters = { 
+        std::array{0,0}, std::array{0,0}, std::array{0,0}, std::array{0,0} 
+    };
+
+    std::array<float, 4> freqRepl = {0.2, 0.2, 0.2, 0.2};
+    std::map<int,int> sumOrigReplacedParticles = {{4122,0}};
+
+    // Hadrons to check
+    std::array<int, 5> checkPdgHadron{34122, 44122, 54122, 9422111, 5122};
+
+    // Correct decays
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{
+        {34122, {{421, 2212}}},    
+        {44122, {{421, 2212}}},    
+        {54122, {{421, 2212}}},    
+        {9422111, {{413, 2212}}},  
+        {5122, {{421, 2212, -211}, {41221, -211}, {34122, -211}, {44122, -211}, {54122, -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 nSignals{}, nSignalGoodDecay{};
+    auto nEvents = tree->GetEntries();
+
+    std::map<int,int> signalHadronsPerType;
+    std::map<int,int> signalGoodDecayPerType;
+    for (auto pdg : checkPdgHadron) { signalHadronsPerType[pdg] = 0; signalGoodDecayPerType[pdg] = 0; }
+
+    for (int i=0; i<nEvents; i++) {
+        tree->GetEntry(i);
+
+        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) {
+            int pdg = track.GetPdgCode();
+            int absPdg = std::abs(pdg);
+
+            if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), absPdg) != checkPdgHadron.end()) {
+                nSignals++;
+                signalHadronsPerType[absPdg]++;
+
+                if (subGeneratorId==checkPdgQuarkOne) {
+                    for (int iRepl=0; iRepl<3; ++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> pdgsDecayAnti;
+                if (track.getFirstDaughterTrackId()>=0) {
+                    for (int j=track.getFirstDaughterTrackId(); j<=track.getLastDaughterTrackId(); ++j) {
+                        int pdgDau = tracks->at(j).GetPdgCode();
+                        pdgsDecay.push_back(pdgDau);
+                        pdgsDecayAnti.push_back((pdgDau==111||pdgDau==333)? pdgDau : -pdgDau);
+                    }
+                }
+
+                std::sort(pdgsDecay.begin(), pdgsDecay.end());
+                std::sort(pdgsDecayAnti.begin(), pdgsDecayAnti.end());
+
+                bool matchedDecay = false;
+                for (auto& decay : checkHadronDecays[absPdg]) {
+                    std::vector<int> decayCopy = decay;
+                    std::sort(decayCopy.begin(), decayCopy.end());
+                    if (pdgsDecay == decayCopy || pdgsDecayAnti == decayCopy) {
+                        matchedDecay = true;
+                        nSignalGoodDecay++;
+                        signalGoodDecayPerType[absPdg]++;
+                        break;
+                    }
+                }
+
+                // Print daughters
+                std::cout << "Particle " << absPdg << " daughters: ";
+                for (auto d : pdgsDecay) std::cout << d << " ";
+                if (matchedDecay) std::cout << "(matches expected decay)";
+                else std::cout << "(does NOT match expected decay)";
+                std::cout << "\n";
+            }
+        }
+    }
+
+    std::cout << "--------------------------------\n";
+    std::cout << "# Events: " << nEvents << "\n";
+    std::cout << "# MB events: " << nEventsMB << "\n";
+    std::cout << "# events injected with " << checkPdgQuarkOne << ": " << nEventsInjOne << "\n";
+    std::cout << "# events injected with " << checkPdgQuarkTwo << ": " << nEventsInjTwo << "\n";
+    std::cout << "# signal hadrons: " << nSignals << "\n";
+    std::cout << "# signal hadrons decaying in correct channels: " << nSignalGoodDecay << "\n";
+
+    for (auto& [pdg, count] : signalHadronsPerType) {
+        int good = signalGoodDecayPerType[pdg];
+        float frac = count>0 ? float(good)/count : 0.;
+        std::cout << "Particle " << pdg << ": " << count << " signals, " << good << " good decays, fraction: " << frac << "\n";
+    }
+
+    // Optional: sanity checks...
+    return 0;
+}

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

@@ -0,0 +1,118 @@
+### author: Antonio Palasciano (antonio.palasciano@cern.ch)
+### since:  July 2024
+
+### beams
+Beams:idA 2212			# proton
+Beams:idB 2212 			# proton
+Beams:eCM 13600. 		# GeV
+
+### processes
+SoftQCD:inelastic on		# all inelastic processes
+
+### 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
+
+### decays
+ParticleDecays:limitTau0 on	
+ParticleDecays:tau0Max 10.
+
+### add resonances not present in PYTHIA
+### id:all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0
+
+### LambdaC(2860)
+34122:all = Lambda_c(2860)+ Lambda_c(2860)- 3 1 0 2.8561 0.0676 2.7 10. 0.
+34122:mayDecay = on
+
+### LambdaC(2880)
+44122:all = Lambda_c(2880)+ Lambda_c(2880)- 5 1 0 2.8816 0.0056 2.85 10. 0.
+44122:mayDecay = on
+
+### LambdaC(2940)
+54122:all = Lambda_c(2940)+ Lambda_c(2940)- 3 1 0 2.9396 0.020 2.9 10. 0.
+54122:mayDecay = on
+
+### TC(3100)
+9422111:all = Theta_c(3100) 3 0 0 3.099 0.0000083 2.95 10. 0.
+9422111:mayDecay = on
+
+### turn off all charm resonances decays
+# 34122:onMode = off    # LambdaC(2860)
+# 44122:onMode = off    # LambdaC(2880)
+# 54122:onMode = off    # LambdaC(2940)
+# 9422111:onMode = off  # Theta_c(3100)
+
+### turn off all beauty hadron decays
+531:onMode = off
+521:onMode = off
+5122:onMode = off
+
+### add LambdaC(2860)
+34122:oneChannel = 1 1 0 421 2212
+34122:onIfMatch = 421 2212
+
+### add LambdaC(2880)
+44122:oneChannel = 1 1 0 421 2212
+44122:onIfMatch = 421 2212
+
+### add LambdaC(2940)
+54122:oneChannel = 1 1 0 421 2212
+54122:onIfMatch = 421 2212
+
+### add ThetaC(3100)
+9422111:oneChannel = 1 1 0 413 2212
+9422111:onIfMatch = 413 2212
+
+### add Lb
+5122:addChannel = 1 0.2 0 421 2212 -211        # Lb -> D0 p pi- non resonant
+5122:addChannel = 1 0.4 0 34122 -211           # Lb -> LambdaC(2860)+ pi-
+5122:addChannel = 1 0.25 0 44122 -211          # Lb -> LambdaC(2880)+ pi-
+5122:addChannel = 1 0.15 0 54122 -211          # Lb -> LambdaC(2940)+ pi-
+
+5122:onIfMatch = 421 2212 -211    # D0 p pi-
+5122:onIfMatch = 34122 -211       # LambdaC(2860)+ pi-
+5122:onIfMatch = 44122 -211       # LambdaC(2880)+ pi-
+5122:onIfMatch = 54122 -211       # LambdaC(2940)+ pi-
+
+# 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
+### set D0 BRs
+421:oneChannel = 1 0.9 0 -321 211
+421:addChannel = 1 0.1 0 -321 211 111
+
+### K* -> K pi
+313:onMode = off
+313:onIfAll = 321 211
+
+### switch off all decay channels
+411:onMode = off
+413:onMode = off
+421:onMode = off
+
+### D0 -> K pi
+421:onIfMatch = 321 211
+### D0 -> K pi pi0
+421:onIfMatch = 321 211 111
+### D*+0 -> D0 pi
+413:onIfMatch = 421 211