[PWGHF] Add configurations for D2H XiC0, XiC+ and OmegaC MCs (gen onl…

…y, no gap trigger)

MC/config/PWGHF/ini/GeneratorHF_D2H_bbbar_Mode2_OmegaC.ini

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

MC/config/PWGHF/ini/GeneratorHF_D2H_bbbar_Mode2_XiC0_XiCplus.ini

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

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_Mode2_OmegaC.ini

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

MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_Mode2_XiC0_XiCplus.ini

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

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

@@ -0,0 +1,105 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuark{5};
+    float ratioTrigger = 1.; // only enriched events
+
+    std::vector<int> checkPdgHadron{4332};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {4332, {{211, 3334}, {211, 3312}}}, // Omegac0
+    };
+
+    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{}, nEventsInj{};
+    int nQuarks{}, 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 == checkPdgQuark) {
+                nEventsInj++;
+            }
+        }
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::abs(pdg) == checkPdgQuark) {
+                nQuarks++;
+                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: ", checkPdgQuark) << nEventsInj << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuark) << nQuarks << "\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 (nEventsInj < nEvents * ratioTrigger * 0.95 || nEventsInj > nEvents * ratioTrigger * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuark << " different than expected\n";
+        return 1;
+    }
+
+    if (nQuarks < 2 * 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 " << checkPdgQuark << " 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/ini/tests/GeneratorHF_D2H_bbbar_Mode2_XiC0_XiCplus.C

@@ -0,0 +1,106 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuark{5};
+    float ratioTrigger = 1.; // only enriched events
+
+    std::vector<int> checkPdgHadron{4132, 4232};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {4132, {{211, 3312}}}, // Xic0
+        {4232, {{-313, 2212}, {-321, 3324}, {211, 211, 3312}, {-321, 211, 2212}}} // Xic+
+    };
+
+    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{}, nEventsInj{};
+    int nQuarks{}, 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 == checkPdgQuark) {
+                nEventsInj++;
+            }
+        }
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::abs(pdg) == checkPdgQuark) {
+                nQuarks++;
+                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: ", checkPdgQuark) << nEventsInj << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuark) << nQuarks << "\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 (nEventsInj < nEvents * ratioTrigger * 0.95 || nEventsInj > nEvents * ratioTrigger * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuark << " different than expected\n";
+        return 1;
+    }
+
+    if (nQuarks < 2 * 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 " << checkPdgQuark << " 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/ini/tests/GeneratorHF_D2H_ccbar_Mode2_OmegaC.C

@@ -0,0 +1,105 @@
+int External() {
+    std::string path{"o2sim_Kine.root"};
+
+    int checkPdgQuark{4};
+    float ratioTrigger = 1.; // only enriched events
+
+    std::vector<int> checkPdgHadron{4332};
+    std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
+        {4332, {{211, 3334}, {211, 3312}}}, // Omegac0
+    };
+
+    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{}, nEventsInj{};
+    int nQuarks{}, 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 == checkPdgQuark) {
+                nEventsInj++;
+            }
+        }
+
+        for (auto &track : *tracks) {
+            auto pdg = track.GetPdgCode();
+            if (std::abs(pdg) == checkPdgQuark) {
+                nQuarks++;
+                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: ", checkPdgQuark) << nEventsInj << "\n";
+    std::cout << Form("# %d (anti)quarks: ", checkPdgQuark) << nQuarks << "\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 (nEventsInj < nEvents * ratioTrigger * 0.95 || nEventsInj > nEvents * ratioTrigger * 1.05) {
+        std::cerr << "Number of generated events injected with " << checkPdgQuark << " different than expected\n";
+        return 1;
+    }
+
+    if (nQuarks < 2 * 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 " << checkPdgQuark << " 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;
+}