diff --git a/PWGLF/Tasks/Resonances/chk892pp.cxx b/PWGLF/Tasks/Resonances/chk892pp.cxx index e888723a6b3..e0986eb7cbe 100644 --- a/PWGLF/Tasks/Resonances/chk892pp.cxx +++ b/PWGLF/Tasks/Resonances/chk892pp.cxx @@ -25,6 +25,7 @@ #include // #include // FIXME #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" #include "PWGLF/Utils/collisionCuts.h" #include "PWGLF/Utils/inelGt.h" @@ -59,6 +60,7 @@ #include "Math/Vector3D.h" #include "Math/Vector4D.h" #include "TF1.h" +#include "TParticlePDG.h" #include "TRandom3.h" #include "TVector2.h" #include @@ -107,9 +109,20 @@ struct Chk892pp { kTYEnd }; + enum EvtStep { + kAll = 0, + kZvtx, + kINELgt0, + kAssocReco, + kNSteps + }; + + const int nSteps = static_cast(EvtStep::kNSteps); + SliceCache cache; Preslice perCollision = aod::track::collisionId; Preslice perCollisionV0 = aod::v0data::collisionId; + Preslice perMCCollision = o2::aod::mcparticle::mcCollisionId; using EventCandidates = soa::Join; using TrackCandidates = soa::Join; @@ -164,7 +177,7 @@ struct Chk892pp { Configurable cfgincludeCentralityMC{"cfgincludeCentralityMC", false, "Include centrality in MC"}; Configurable cfgEvtCollInTimeRangeStandard{"cfgEvtCollInTimeRangeStandard", false, "Evt sel: apply NoCollInTimeRangeStandard"}; Configurable cfgEventCentralityMin{"cfgEventCentralityMin", 0.0f, "Event sel: minimum centrality"}; - Configurable cfgEventCentralityMax{"cfgEventCentralityMax", 80.0f, "Event sel: maximum centrality"}; + Configurable cfgEventCentralityMax{"cfgEventCentralityMax", 100.0f, "Event sel: maximum centrality"}; Configurable cfgEvtUseRCTFlagChecker{"cfgEvtUseRCTFlagChecker", false, "Evt sel: use RCT flag checker"}; Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; @@ -250,6 +263,9 @@ struct Chk892pp { Configurable cfgNrotBkg{"cfgNrotBkg", 4, "Number of rotated copies (background) per each original candidate"}; } BkgEstimationConfig; + Configurable cfgTruthUseInelGt0{"cfgTruthUseInelGt0", true, "Truth denominator: require INEL>0"}; + Configurable cfgTruthIncludeZvtx{"cfgTruthIncludeZvtx", true, "Truth denominator: also require |vtxz|(HIST("Correction/setSizes")); + hset->GetXaxis()->SetBinLabel(1, "refClassIds"); + hset->GetXaxis()->SetBinLabel(2, "allowedMcIds"); + hset->GetXaxis()->SetBinLabel(3, "intersection"); + hset->GetXaxis()->SetBinLabel(4, "allowed-only"); + + histos.add("Correction/hNEventsMCTruth", "hNEventsMCTruth", HistType::kTH1F, {AxisSpec{nSteps, 0.5, nSteps + 0.5, ""}}); + auto hstep = histos.get(HIST("Correction/hNEventsMCTruth")); + hstep->GetXaxis()->SetBinLabel(1, "All"); + hstep->GetXaxis()->SetBinLabel(2, "zvtx"); + hstep->GetXaxis()->SetBinLabel(3, "INEL>0"); + hstep->GetXaxis()->SetBinLabel(4, "Assoc with reco coll"); } ccdb->setURL(CCDBConfig.cfgURL); @@ -453,27 +491,6 @@ struct Chk892pp { } } - template - int getlDetId(DetNameType const& name) - { - LOGF(info, "GetDetID running"); - if (name.value == "FT0C") { - return 0; - } else if (name.value == "FT0A") { - return 1; - } else if (name.value == "FT0M") { - return 2; - } else if (name.value == "FV0A") { - return 3; - } else if (name.value == "TPCpos") { - return 4; - } else if (name.value == "TPCneg") { - return 5; - } else { - return false; - } - } - // Track selection template bool trackCut(TrackType const& track) @@ -689,62 +706,105 @@ struct Chk892pp { } std::unordered_set allowedMcIds; - std::unordered_map centByMcIds; + std::unordered_map centTruthByAllowed; + std::unordered_set refClassIds; + std::unordered_map refCentByMcId; - void buildAllowedMcIds(MCEventCandidates const& events) + template + void buildAllowedMcIds(RecoEventsT const& events) { allowedMcIds.clear(); - centByMcIds.clear(); + centTruthByAllowed.clear(); for (const auto& coll : events) { - lCentrality = getCentrality(coll); - histos.fill(HIST("QACent_woCut"), lCentrality); - histos.fill(HIST("QAvtxz_woCut"), coll.posZ()); + // lCentrality = getCentrality(coll); + + if (!coll.has_mcCollision()) + continue; + + const auto mcid = coll.mcCollisionId(); + const auto mccoll = coll.template mcCollision_as>(); + + const float lCentrality = mccoll.centFT0M(); + + if (doprocessMC) { + histos.fill(HIST("QACent_woCut"), lCentrality); + histos.fill(HIST("QAvtxz_woCut"), coll.posZ()); + } if (!colCuts.isSelected(coll)) continue; + if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(coll)) + continue; if (!coll.isInelGt0()) continue; - histos.fill(HIST("QACent_woCentCut"), lCentrality); - histos.fill(HIST("QAvtxz_wVtxzCut"), coll.posZ()); + if (doprocessMC) { + histos.fill(HIST("QACent_woCentCut"), lCentrality); + histos.fill(HIST("QAvtxz_wVtxzCut"), coll.posZ()); + } if (lCentrality < EventCuts.cfgEventCentralityMin || lCentrality > EventCuts.cfgEventCentralityMax) continue; - histos.fill(HIST("QACent_wCentCut"), lCentrality); - if (!coll.has_mcCollision()) - continue; - const auto mcid = coll.mcCollisionId(); - if (mcid < 0) - continue; + + if (doprocessMC) { + histos.fill(HIST("QACent_wCentCut"), lCentrality); + } allowedMcIds.insert(mcid); - centByMcIds[mcid] = lCentrality; + centTruthByAllowed.emplace(mcid, lCentrality); + } + } + + template + void buildReferenceMcIds(McCollsT const& mccolls, McPartsT const& mcparts) + { + refClassIds.clear(); + refCentByMcId.clear(); + + for (const auto& coll : mccolls) { + bool pass = true; + + if (cfgTruthIncludeZvtx && std::abs(coll.posZ()) >= EventCuts.cfgEvtZvtx) + pass = false; + + if (pass && cfgTruthUseInelGt0) { + auto partsThisMc = mcparts.sliceBy(perMCCollision, coll.globalIndex()); + if (!pwglf::isINELgtNmc(partsThisMc, 0, pdg)) + pass = false; + } + + if (!pass) + continue; + + const auto mcid = coll.globalIndex(); + refClassIds.insert(mcid); + const float lCentrality = coll.centFT0M(); + refCentByMcId.emplace(mcid, lCentrality); } } void effK0sProcessGen(MCTrueTrackCandidates const& mcparts) { for (const auto& part : mcparts) { - const auto mcid = part.mcCollisionId(); - if (!part.has_mcCollision()) continue; if (part.pdgCode() != kPDGK0s) continue; if (!part.isPhysicalPrimary()) continue; - if (allowedMcIds.find(mcid) == allowedMcIds.end()) + if (std::abs(part.y()) > SecondaryCuts.cfgSecondaryRapidityMax) continue; - auto iter = centByMcIds.find(mcid); - if (iter == centByMcIds.end()) + const auto mcid = part.mcCollisionId(); + if (allowedMcIds.count(mcid) == 0) continue; - const float lCentrality = iter->second; - - if (std::abs(part.y()) > SecondaryCuts.cfgSecondaryRapidityMax) + auto iter = centTruthByAllowed.find(mcid); + if (iter == centTruthByAllowed.end()) continue; + const float lCentrality = iter->second; + histos.fill(HIST("EffK0s/genK0s"), part.pt(), lCentrality); } } @@ -759,14 +819,11 @@ struct Chk892pp { const auto mcid = coll.mcCollisionId(); - if (allowedMcIds.find(mcid) == allowedMcIds.end()) - continue; - - auto iter = centByMcIds.find(mcid); - if (iter == centByMcIds.end()) + if (allowedMcIds.count(mcid) == 0) continue; - const float lCentrality = iter->second; + const auto mccoll = coll.template mcCollision_as>(); + const float lCentrality = mccoll.centFT0M(); const double ptreco = v0.pt(); const double yreco = v0.yK0Short(); @@ -845,23 +902,23 @@ struct Chk892pp { void effKstarProcessGen(MCTrueTrackCandidates const& mcparts) { for (const auto& part : mcparts) { - const auto mcid = part.mcCollisionId(); if (!part.has_mcCollision()) continue; if (std::abs(part.pdgCode()) != kKstarPlus) continue; - if (allowedMcIds.find(mcid) == allowedMcIds.end()) + if (std::abs(part.y()) > KstarCuts.cfgKstarMaxRap) continue; - auto iter = centByMcIds.find(mcid); - if (iter == centByMcIds.end()) + const auto mcid = part.mcCollisionId(); + if (allowedMcIds.count(mcid) == 0) continue; - const float lCentrality = iter->second; - - if (std::abs(part.y()) > KstarCuts.cfgKstarMaxRap) + auto iter = centTruthByAllowed.find(mcid); + if (iter == centTruthByAllowed.end()) continue; + const float lCentrality = iter->second; + histos.fill(HIST("EffKstar/genKstar"), part.pt(), lCentrality); } } // effKstarProcessGen @@ -877,14 +934,11 @@ struct Chk892pp { const auto mcid = coll.mcCollisionId(); - if (allowedMcIds.find(mcid) == allowedMcIds.end()) + if (allowedMcIds.count(mcid) == 0) continue; - auto iter = centByMcIds.find(mcid); - if (iter == centByMcIds.end()) - continue; - - const float lCentrality = iter->second; + const auto mccoll = coll.template mcCollision_as>(); + const float lCentrality = mccoll.centFT0M(); if (!SecondaryCuts.cfgByPassDauPIDSelection) { auto posDauTrack = v0.template posTrack_as(); @@ -901,7 +955,6 @@ struct Chk892pp { for (const auto& bTrack : trks) { if (bTrack.collisionId() != v0.collisionId()) continue; - if (!trackCut(bTrack)) continue; if (!selectionPIDPion(bTrack)) @@ -927,6 +980,54 @@ struct Chk892pp { } } // effKstarProcessReco + void fillSigLossNum(MCTrueTrackCandidates const& mcparts) + { + for (auto const& part : mcparts) { + if (!part.has_mcCollision()) + continue; + if (std::abs(part.pdgCode()) != kKstarPlus) + continue; + if (std::abs(part.y()) > KstarCuts.cfgKstarMaxRap) + continue; + + const auto mcid = part.mcCollisionId(); + if (allowedMcIds.count(mcid) == 0) + continue; + + auto iter = centTruthByAllowed.find(mcid); + if (iter == centTruthByAllowed.end()) + continue; + + const float lCentrality = iter->second; + + histos.fill(HIST("Correction/sigLoss_num"), part.pt(), lCentrality); + } + } // fillSigLossNum + + void fillSigLossDen(MCTrueTrackCandidates const& mcparts) + { + for (auto const& part : mcparts) { + if (!part.has_mcCollision()) + continue; + if (std::abs(part.pdgCode()) != kKstarPlus) + continue; + if (std::abs(part.y()) > KstarCuts.cfgKstarMaxRap) + continue; + + const auto mcid = part.mcCollisionId(); + if (refClassIds.count(mcid) == 0) + continue; + + auto iter = refCentByMcId.find(mcid); + if (iter == refCentByMcId.end()) + continue; + + const float lCentrality = iter->second; + + histos.fill(HIST("Correction/sigLoss_den"), part.pt(), lCentrality); + } + } // fillSigLossDen + int count = 0; template @@ -1157,7 +1258,6 @@ struct Chk892pp { } // fillHistograms - // process data void processData(EventCandidates::iterator const& collision, TrackCandidates const& tracks, V0Candidates const& v0s, @@ -1179,19 +1279,77 @@ struct Chk892pp { } PROCESS_SWITCH(Chk892pp, processData, "Process Event for data without Partitioning", false); - // process MC reconstructed level void processMC(MCTrueTrackCandidates const& mcpart, MCTrackCandidates const& tracks, MCV0Candidates const& v0s, - MCEventCandidates const& events) + MCEventCandidates const& events, + soa::Join const& mccolls) { buildAllowedMcIds(events); + buildReferenceMcIds(mccolls, mcpart); effK0sProcessGen(mcpart); effK0sProcessReco(v0s); effKstarProcessGen(mcpart); effKstarProcessReco(v0s, tracks); + fillSigLossNum(mcpart); + fillSigLossDen(mcpart); + + for (const auto& mcid : refClassIds) { + histos.fill(HIST("Correction/EF_den"), refCentByMcId[mcid]); + } + for (const auto& mcid : allowedMcIds) { + auto iter = centTruthByAllowed.find(mcid); + if (iter == centTruthByAllowed.end()) + continue; + + const float lCentrality = iter->second; + histos.fill(HIST("Correction/EF_num"), lCentrality); + } + + size_t nIntersect = 0; + for (const auto& mcid : allowedMcIds) + if (refClassIds.count(mcid)) + nIntersect++; + histos.fill(HIST("Correction/setSizes"), 0.0, refClassIds.size()); + histos.fill(HIST("Correction/setSizes"), 1.0, allowedMcIds.size()); + histos.fill(HIST("Correction/setSizes"), 2.0, nIntersect); + histos.fill(HIST("Correction/setSizes"), 3.0, allowedMcIds.size() - nIntersect); + + for (const auto& mcc : mccolls) { + histos.fill(HIST("Correction/MCTruthCent_all"), mcc.centFT0M()); + } + + for (const auto& mcid : refClassIds) { + auto iter = refCentByMcId.find(mcid); + if (iter == refCentByMcId.end()) + continue; + lCentrality = iter->second; + histos.fill(HIST("Correction/MCTruthCent_cut"), lCentrality); + } + + for (auto const& mcc : mccolls) { + const auto mcid = mcc.globalIndex(); + + histos.fill(HIST("Correction/hNEventsMCTruth"), 1.0); + + bool passZvtx = true; + if (cfgTruthIncludeZvtx && std::abs(mcc.posZ()) > EventCuts.cfgEvtZvtx) { + passZvtx = false; + } + if (passZvtx) { + histos.fill(HIST("Correction/hNEventsMCTruth"), 2.0); + + auto partsThisMc = mcpart.sliceBy(perMCCollision, mcid); + if (pwglf::isINELgtNmc(partsThisMc, 0, pdg)) { + histos.fill(HIST("Correction/hNEventsMCTruth"), 3.0); + } + } + if (allowedMcIds.count(mcid)) { + histos.fill(HIST("Correction/hNEventsMCTruth"), 4.0); + } + } } - PROCESS_SWITCH(Chk892pp, processMC, "Process Event for MC", false); + PROCESS_SWITCH(Chk892pp, processMC, "Process Event for MC", true); void processMCQA(MCEventCandidates::iterator const& collision, MCTrackCandidates const& tracks, @@ -1202,17 +1360,21 @@ struct Chk892pp { return; if (EventCuts.cfgEvtUseRCTFlagChecker && !rctChecker(collision)) return; + if (!collision.isInelGt0()) + return; - lCentrality = getCentrality(collision); - if (lCentrality < EventCuts.cfgEventCentralityMin || lCentrality > EventCuts.cfgEventCentralityMax) + if (!collision.has_mcCollision()) return; - if (!collision.isInelGt0()) + auto mccoll = collision.template mcCollision_as>(); + const float lCentrality = mccoll.centFT0M(); + + if (lCentrality < EventCuts.cfgEventCentralityMin || lCentrality > EventCuts.cfgEventCentralityMax) return; colCuts.fillQA(collision); fillHistograms(collision, tracks, v0s); } - PROCESS_SWITCH(Chk892pp, processMCQA, "Process Event for MC and fill QA plots", true); + PROCESS_SWITCH(Chk892pp, processMCQA, "Process Event for MC and fill QA plots", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)