PWGHF: add ini for D resonances MCs for pp collisions

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.ini

@@ -0,0 +1,8 @@
+### The external generator derives from GeneratorPythia8.
+[GeneratorExternal]
+fileName=${O2DPG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
+funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5)
+
+[GeneratorPythia8]
+config=${O2DPG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg
+includePartonEvent=true

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

@@ -0,0 +1,140 @@
+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, 415, 421, 425, 431, 435, 511, 521, 531, 4122, 10411, 10421, 10433, 20423, 20433};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {411, {{-321, 211, 211}, {-313, 211}, {211, 311}, {211, 333}}}, // D+
+        {415, {{211, 421}}}, // D2*(2460)+
+        {421, {{-321, 211}, {-321, 111, 211}}}, // D0
+        {425, {{-211, 413},{-211, 411}}}, // D2*(2460)0
+        {431, {{211, 333}, {-313, 321}}}, // Ds+
+        {425, {{311, 413}, {311, 411}, {321,421}}}, // Ds2*(2573)
+        {511, {{-415, -11, 12}, {-10411, -11, 12}, {-415, -13, 14}, {-10411, -13, 14}, {-415, -15, 16}, {-10411, -15, 16},
+               {-10411, 211}, {-10421, 211}, {-415, 433}, { -415, 431}, {-415, 211}, {-415, 213} }}, // B0
+        {521, {{-20423, -11, 12}, {-425, -11, 12}, {-10421, -11, 12}, {-20423, -13, 14}, {-425, -13, 14}, {-10421, -13, 14},
+               {-20423, -15, 16}, {-425, -15, 16}, {-10421, -15, 16}, {-20423, 211}, {-20423, 213}, {-20423, 431}, {-20423, 433}, 
+               {-425, 211}, {-425, 213}, {-425, 431}, {-425, 433}}}, // B+
+        {531, {{-435, -11, 12}, {-10433, -11, 12}, {-435, -13, 14}, {-10433, -13, 14}, {-435, -15, 16}, {-10433, -15, 16},
+               {-435, 211}, {-20433, 211}, {-20433, 213}}},// Bs0
+        {4122, {{-313, 2212}, {-321, 2224}, {211, 3124}, {-321, 211, 2212}, {311, 2212}}}, // Lc+
+        {10411, {{211, 421}}}, // D0*+
+        {10421, {{-211, 411}}}, // D0*0
+        {10433, {{311, 413}}}, // Ds1(2536)
+        {20423, {{-211, 413}}}, // D1(2430)0
+        {20433, {{22, 431}, {-211, 211, 431} }} // Ds1 (2460)
+    };
+
+    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;
+}

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

@@ -0,0 +1,237 @@
+### author: Stefano Politano, Luca Aglietta ([email protected], [email protected])
+### since:  July 2024
+
+### 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.
+
+# parameters to boost resonances production
+StringFlav:mesonCL1S0J1= 3 
+StringFlav:mesonCL1S1J2= 3.2 
+StringFlav:mesonCL1S1J0= 0.75 
+StringFlav:mesonCL1S1J1= 1. 
+
+
+### turn off all charm resonances decays
+10433:onMode = off  # Ds1(2536)
+435:onMode = off    # Ds2*(2573)
+425:onMode = off    # D2*(2460)0
+415:onMode = off    # D2*(2460)+
+10411:onMode = off  # D0*+
+10421:onMode = off  # D0*0
+20433:onMode = off  # Ds1 (2460)
+20423:onMode = off  # D1(2430)0
+
+### turn off all beauty hadron decays
+531:onMode = off
+511:onMode = off
+521:onMode = off
+
+### add D1(2430)0
+20423:oneChannel = 1 1 0  413 -211
+20423:onIfMatch = 413 -211
+
+### add Ds1(2536)
+10433:oneChannel = 1 1 0 413 311
+10433:onIfMatch = 413 311
+
+### Ds1 (2460)
+20433:oneChannel = 1 0.5 0 431 22
+20433:addChannel = 1 0.5 0 431 211 -211
+20433:onIfMatch = 431 22
+20433:onIfMatch = 431 211 211
+
+### add Ds2*(2573)
+435:oneChannel = 1 0.0500000 0 413 311
+435:addChannel = 1 0.4500000 0 411 311
+435:addChannel = 1 0.4500000 0 421 321
+435:onIfMatch = 413 311
+435:onIfMatch = 411 311
+435:onIfMatch = 421 321
+
+### add D2*(2460)0
+425:oneChannel = 1 0.5 0 413 -211
+425:addChannel = 1 0.5 0 411 -211
+425:onIfMatch = 413 211
+425:onIfMatch = 411 211
+
+### add D2*(2460)+
+415:oneChannel = 1 1 0 421 211
+415:onIfMatch = 421 211
+
+
+### add D0*+
+10411:oneChannel = 1 1 0 421 211
+10411:onIfMatch = 421 211
+
+### add D0*0
+10421:oneChannel = 1 1 0 411 -211
+10421:onIfMatch = 411 211
+
+
+### add Bs0
+531:oneChannel = 1 0.0070000 0 12 -11 -435 
+531:addChannel = 1 0.0070000 0 12 -11 -10433
+531:addChannel = 1 0.0070000 0 14 -13 -435
+531:addChannel = 1 0.0070000 0 14 -13 -10433
+531:addChannel = 1 0.0040000 0 14 -13 -20433
+531:addChannel = 1 0.0160000 0 16 -15 -433
+531:addChannel = 1 0.0028000 0 16 -15 -435
+531:addChannel = 1 0.0028000 0 16 -15 -10433
+531:addChannel = 1 0.0013000 0 -435 211
+531:addChannel = 1 0.0008000 0 -20433 211
+531:addChannel = 1 0.0021000 0 -20433 213
+
+531:onIfMatch = 12 11 435 
+531:onIfMatch = 12 11 10433
+531:onIfMatch = 14 13 435
+531:onIfMatch = 14 13 10433
+531:onIfMatch = 14 13 20433
+531:onIfMatch = 16 15 433
+531:onIfMatch = 16 15 435
+531:onIfMatch = 16 15 10433
+531:onIfMatch = 435 211
+531:onIfMatch = 20433 211
+531:onIfMatch = 20433 213
+
+### add B0
+511:oneChannel = 1 0.0023000 0 12 -11 -415
+511:addChannel = 1 0.0045000 0 12 -11 -10411
+511:addChannel = 1 0.0023000 0 14 -13 -415
+511:addChannel = 1 0.0045000 0 14 -13 -10411
+511:addChannel = 1 0.0020000 0 16 -15 -415
+511:addChannel = 1 0.0013000 0 16 -15 -10411
+511:addChannel = 1 0.0002000 0 -10411 211
+511:addChannel = 1 0.0009100 0 -10421 211
+511:addChannel = 1 0.0040000 0 433 -415
+511:addChannel = 1 0.0042000 0 431 -415
+511:addChannel = 1 0.0009000 0 -415 211
+511:addChannel = 1 0.0022000 0 -415 213
+
+511:onIfMatch = 12 11 415
+511:onIfMatch = 12 11 10411
+511:onIfMatch = 14 13 415
+511:onIfMatch = 14 13 10411
+511:onIfMatch = 16 15 415
+511:onIfMatch = 16 15 10411
+511:onIfMatch = 10411 211
+511:onIfMatch = 10421 211
+511:onIfMatch = 433 415
+511:onIfMatch = 431 415
+511:onIfMatch = 415 211
+511:onIfMatch = 415 213
+
+### add B+
+521:oneChannel = 1 0.0090000 0 12 -11 -20423
+521:addChannel = 1 0.0090000 0 14 -13 -20423
+521:addChannel = 1 0.0020000 0 16 -15 -20423
+521:addChannel = 1 0.0030000 0 12 -11 -425
+521:addChannel = 1 0.0030000 0 14 -13 -425
+521:addChannel = 1 0.0020000 0 16 -15 -425
+521:addChannel = 1 0.0049000 0 12 -11 -10421
+521:addChannel = 1 0.0049000 0 14 -13 -10421
+521:addChannel = 1 0.0013000 0 16 -15 -10421
+521:addChannel = 1 0.0007500 0 -20423 211
+521:addChannel = 1 0.0022000 0 -20423 213
+521:addChannel = 1 0.0006000 0 -20423 431
+521:addChannel = 1 0.0012000 0 -20423 433
+521:addChannel = 1 0.0008000 0 -425 211
+521:addChannel = 1 0.0038000 0 -425 213
+521:addChannel = 1 0.0042000 0 431 -425
+521:addChannel = 1 0.0040000 0 433 -425
+
+521:onIfMatch = 12 11 20423
+521:onIfMatch = 14 13 20423
+521:onIfMatch = 16 15 20423
+521:onIfMatch = 12 11 425
+521:onIfMatch = 14 13 425
+521:onIfMatch = 16 15 425
+521:onIfMatch = 12 11 10421
+521:onIfMatch = 14 13 10421
+521:onIfMatch = 16 15 10421
+521:onIfMatch = 20423 211
+521:onIfMatch = 20423 213
+521:onIfMatch = 20423 431
+521:onIfMatch = 20423 433
+521:onIfMatch = 425 211
+521:onIfMatch = 425 213
+521:onIfMatch = 431 425
+521:onIfMatch = 433 425
+
+# 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
+### add D+ decays absent in PYTHIA8 decay table and set BRs from PDG for other
+411:oneChannel = 1 0.0752 0 -321 211 211
+411:addChannel = 1 0.0104 0 -313 211
+411:addChannel = 1 0.0156 0 311 211
+411:addChannel = 1 0.0752 0 333 211 # to have the same amount of D+->KKpi and D+->Kpipi
+## add Lc decays absent in PYTHIA8 decay table and set BRs from PDG for other
+4122:oneChannel = 1 0.0196 100 2212 -313
+4122:addChannel = 1 0.0108 100 2224 -321
+4122:addChannel = 1 0.022 100 102134 211
+4122:addChannel = 1 0.035 0 2212 -321 211
+4122:addChannel = 1 0.0159 0 2212 311
+
+### K* -> K pi
+313:onMode = off
+313:onIfAll = 321 211
+### for Ds -> Phi pi+
+333:onMode = off
+333:onIfAll = 321 321
+### for D0 -> rho0 pi+ k-
+113:onMode = off
+113:onIfAll = 211 211
+### for Lambda_c -> Delta++ K-
+2224:onMode = off
+2224:onIfAll = 2212 211
+### for Lambda_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
+
+### D0 -> K pi
+421:onIfMatch = 321 211
+
+### D+/- -> K pi pi 
+411:onIfMatch = 321 211 211
+### D+/- -> K* pi
+411:onIfMatch = 313 211
+### D+/- -> phi pi
+411:onIfMatch = 333 211
+
+### D_s -> K K*
+431:onIfMatch = 321 313
+### D_s -> Phi pi
+431:onIfMatch = 333 211
+
+### Lambda_c -> p K*
+4122:onIfMatch = 2212 313
+### Lambda_c -> Delta K
+4122:onIfMatch = 2224 321
+### Lambda_c -> Lambda(1520) pi
+4122:onIfMatch = 102134 211
+### Lambda_c -> p K pi
+4122:onIfMatch = 2212 321 211
+### Lambda_c -> pK0s
+4122:onIfMatch = 2212 311 
+