From e379563c5a6a6d608e9f4f67d10846cd4f1568f3 Mon Sep 17 00:00:00 2001 From: Alicja Plachta Date: Thu, 18 Dec 2025 14:16:38 +0100 Subject: [PATCH] Add pair fractions calculation for track-track and extend it for track-V0 and V0-V0 --- ...iversePairTaskTrackTrackMultKtExtended.cxx | 348 ++++++++++++++---- .../femtoUniversePairTaskTrackV0Extended.cxx | 30 +- 2 files changed, 314 insertions(+), 64 deletions(-) diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx index 5ff117461ca..7eee98d1f62 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackMultKtExtended.cxx @@ -25,20 +25,15 @@ #include "PWGCF/FemtoUniverse/Core/FemtoUniversePairWithCentMultKt.h" #include "PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h" #include "PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h" -#include "PWGCF/FemtoUniverse/Core/femtoUtils.h" #include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h" #include "Framework/ASoAHelpers.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/O2DatabasePDGPlugin.h" -#include "Framework/RunningWorkflowInfo.h" -#include "Framework/StepTHn.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/PID.h" -#include "TDatabasePDG.h" - #include #include @@ -87,10 +82,16 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { SliceCache cache; Preslice perCol = aod::femtouniverseparticle::fdCollisionId; + using FemtoTruthParticles = soa::Filtered>; + Preslice perColMCTruth = aod::femtouniverseparticle::fdCollisionId; + + using FemtoRecoParticles = soa::Filtered>; + Preslice perColMC = aod::femtouniverseparticle::fdCollisionId; + /// Particle 1 struct : o2::framework::ConfigurableGroup { Configurable confPDGCodePartOne{"confPDGCodePartOne", 211, "Particle 1 -- PDG code"}; - // Configurable ConfCutPartOne{"ConfCutPartOne", 5542474, "Particle 1 -- Selection bit from cutCulator"}; + // Configurable confCutPartOne{"confCutPartOne", 5542474, "Particle 1 -- Selection bit from cutCulator"}; Configurable confPIDPartOne{"confPIDPartOne", 2, "Particle 1 -- Read from cutCulator"}; // we also need the possibility to specify whether the bit is true/false ->std>>vector>int>> Configurable confpLowPart1{"confpLowPart1", 0.14, "Lower limit for Pt for the first particle"}; Configurable confPtHighPart1{"confPtHighPart1", 1.5, "Higher limit for Pt for the first particle"}; @@ -100,7 +101,10 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { /// Partition for particle 1 Partition partsOne = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && aod::femtouniverseparticle::sign == as(trackonefilter.confChargePart1) && aod::femtouniverseparticle::pt < trackonefilter.confPtHighPart1 && aod::femtouniverseparticle::pt > trackonefilter.confpLowPart1; - Partition> partsOneMC = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && aod::femtouniverseparticle::sign == as(trackonefilter.confChargePart1) && aod::femtouniverseparticle::pt < trackonefilter.confPtHighPart1 && aod::femtouniverseparticle::pt > trackonefilter.confpLowPart1; + Partition partsOneMC = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && aod::femtouniverseparticle::sign == as(trackonefilter.confChargePart1) && aod::femtouniverseparticle::pt < trackonefilter.confPtHighPart1 && aod::femtouniverseparticle::pt > trackonefilter.confpLowPart1; + + Partition partsOneMCTruth = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack) && aod::femtouniverseparticle::pt < trackonefilter.confPtHighPart1 && aod::femtouniverseparticle::pt > trackonefilter.confpLowPart1; + ; /// Histogramming for particle 1 FemtoUniverseParticleHisto trackHistoPartOne; @@ -108,7 +112,7 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { /// Particle 2 struct : o2::framework::ConfigurableGroup { Configurable confPDGCodePartTwo{"confPDGCodePartTwo", 211, "Particle 2 -- PDG code"}; - // Configurable ConfCutPartTwo{"ConfCutPartTwo", 5542474, "Particle 2 -- Selection bit"}; + // Configurable confCutPartTwo{"confCutPartTwo", 5542474, "Particle 2 -- Selection bit"}; Configurable confPIDPartTwo{"confPIDPartTwo", 2, "Particle 2 -- Read from cutCulator"}; // we also need the possibility to specify whether the bit is true/false ->std>>vector> Configurable confpLowPart2{"confpLowPart2", 0.14, "Lower limit for Pt for the second particle"}; Configurable confPtHighPart2{"confPtHighPart2", 1.5, "Higher limit for Pt for the second particle"}; @@ -118,7 +122,9 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { /// Partition for particle 2 Partition partsTwo = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && (aod::femtouniverseparticle::sign == as(tracktwofilter.confChargePart2)) && aod::femtouniverseparticle::pt < tracktwofilter.confPtHighPart2 && aod::femtouniverseparticle::pt > tracktwofilter.confpLowPart2; - Partition> partsTwoMC = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack) && (aod::femtouniverseparticle::sign == as(tracktwofilter.confChargePart2)) && aod::femtouniverseparticle::pt < tracktwofilter.confPtHighPart2 && aod::femtouniverseparticle::pt > tracktwofilter.confpLowPart2; + Partition partsTwoMC = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack) && (aod::femtouniverseparticle::sign == as(tracktwofilter.confChargePart2)) && aod::femtouniverseparticle::pt < tracktwofilter.confPtHighPart2 && aod::femtouniverseparticle::pt > tracktwofilter.confpLowPart2; + + Partition partsTwoMCTruth = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack) && aod::femtouniverseparticle::pt < tracktwofilter.confPtHighPart2 && aod::femtouniverseparticle::pt > tracktwofilter.confpLowPart2; /// Histogramming for particle 2 FemtoUniverseParticleHisto trackHistoPartTwo; @@ -131,13 +137,15 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { std::vector kNsigma; /// Event part - Configurable confV0MLow{"confV0MLow", 0.0, "Lower limit for V0M multiplicity"}; - Configurable confV0MHigh{"confV0MHigh", 25000.0, "Upper limit for V0M multiplicity"}; - Configurable confSphericityCutMin{"confSphericityCutMin", 0, "Min. sphericity"}; - Configurable confSphericityCutMax{"confSphericityCutMax", 3, "Max. sphericity"}; - - Filter collV0Mfilter = ((o2::aod::femtouniversecollision::multV0M > confV0MLow) && (o2::aod::femtouniversecollision::multV0M < confV0MHigh)); - Filter colSpherfilter = ((o2::aod::femtouniversecollision::sphericity > confSphericityCutMin) && (o2::aod::femtouniversecollision::sphericity < confSphericityCutMax)); + struct : o2::framework::ConfigurableGroup { + Configurable confV0MLow{"confV0MLow", 0.0, "Lower limit for V0M multiplicity"}; + Configurable confV0MHigh{"confV0MHigh", 25000.0, "Upper limit for V0M multiplicity"}; + Configurable confSphericityCutMin{"confSphericityCutMin", 0, "Min. sphericity"}; + Configurable confSphericityCutMax{"confSphericityCutMax", 3, "Max. sphericity"}; + } eventSel; + + Filter collV0Mfilter = ((o2::aod::femtouniversecollision::multV0M > eventSel.confV0MLow) && (o2::aod::femtouniversecollision::multV0M < eventSel.confV0MHigh)); + Filter colSpherfilter = ((o2::aod::femtouniversecollision::sphericity > eventSel.confSphericityCutMin) && (o2::aod::femtouniversecollision::sphericity < eventSel.confSphericityCutMax)); // Filter trackAdditionalfilter = (nabs(aod::femtouniverseparticle::eta) < twotracksconfigs.confEtaMax); // example filtering on configurable /// Particle part @@ -168,11 +176,14 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { Configurable confCPRChosenRadii{"confCPRChosenRadii", 0.80, "Delta Eta cut for Close Pair Rejection"}; Configurable isPairIdentical{"isPairIdentical", true, "'true' for identical particles, 'false' for non-identical particles"}; - Configurable cfgProcessPM{"cfgProcessPM", true, "Process differently charged particles (plus-minus)"}; - Configurable cfgProcessPP{"cfgProcessPP", true, "Process positively charged particles (plus-plus)"}; - Configurable cfgProcessMM{"cfgProcessMM", true, "Process negatively charged particles (minus-minus)"}; - Configurable cfgProcessMultBins{"cfgProcessMultBins", true, "Process kstar histograms (in multiplicity bins)"}; - Configurable cfgProcessKtBins{"cfgProcessKtBins", true, "Process kstar histograms in kT bins (if 'cfgProcessMultBins' is false, it will not be processed regardless of 'cfgProcessKtBins' state)"}; + struct : o2::framework::ConfigurableGroup { + Configurable cfgProcessPM{"cfgProcessPM", true, "Process particles of the opposite charge"}; + Configurable cfgProcessPP{"cfgProcessPP", false, "Process particles of the same, positice charge"}; + Configurable cfgProcessMM{"cfgProcessMM", false, "Process particles of the same, positice charge"}; + } processPair; + Configurable confFillDebug{"confFillDebug", false, "Fill additional histograms: event mixing, sphericity"}; + Configurable cfgProcessMultBins{"cfgProcessMultBins", false, "Process kstar histograms (in multiplicity bins)"}; + Configurable cfgProcessKtBins{"cfgProcessKtBins", false, "Process kstar histograms in kT bins (if 'cfgProcessMultBins' is false, it will not be processed regardless of 'cfgProcessKtBins' state)"}; Configurable cfgProcessKtMt3DCF{"cfgProcessKtMt3DCF", false, "Process 3D histograms in kT and MultBins"}; FemtoUniverseFemtoContainer sameEventCont; @@ -240,15 +251,19 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { /// TPC Kaon Sigma selection (stricter cuts for K+ and K-) -- based on Run2 results bool isKaonNsigma(float mom, float nsigmaTPCK, float nsigmaTOFK) { + double nSigmaTPCRanges[3] = {1, 2, 3}; + double nSigmaTOFRanges[2] = {1.5, 2}; + double momRanges[4] = {0.4, 0.45, 0.8, 1.5}; + if (twotracksconfigs.isKaonNsigma == true) { - if (mom < 0.4) { - return std::abs(nsigmaTPCK) < 2; - } else if (mom > 0.4 && mom < 0.45) { - return std::abs(nsigmaTPCK) < 1; - } else if (mom > 0.45 && mom < 0.8) { - return (std::abs(nsigmaTPCK) < 3 && std::abs(nsigmaTOFK) < 2); - } else if (mom > 0.8 && mom < 1.5) { - return (std::abs(nsigmaTPCK) < 3 && std::abs(nsigmaTOFK) < 1.5); + if (mom < momRanges[0]) { + return std::abs(nsigmaTPCK) < nSigmaTPCRanges[1]; + } else if (mom > momRanges[0] && mom < momRanges[1]) { + return std::abs(nsigmaTPCK) < nSigmaTPCRanges[0]; + } else if (mom > momRanges[1] && mom < momRanges[2]) { + return (std::abs(nsigmaTPCK) < nSigmaTPCRanges[2] && std::abs(nsigmaTOFK) < nSigmaTOFRanges[1]); + } else if (mom > momRanges[2] && mom < momRanges[3]) { + return (std::abs(nsigmaTPCK) < nSigmaTPCRanges[2] && std::abs(nsigmaTOFK) < nSigmaTOFRanges[0]); } else { return false; } @@ -259,7 +274,9 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { bool isParticleNSigma(int8_t particle_number, float mom, float nsigmaTPCPr, float nsigmaTOFPr, float nsigmaTPCPi, float nsigmaTOFPi, float nsigmaTPCK, float nsigmaTOFK) { - if (particle_number == 1) { + enum ParticleNumber { partOne = 1, + partTwo = 2 }; + if (particle_number == partOne) { switch (trackonefilter.confPDGCodePartOne) { case 2212: // Proton+ case -2212: // Proton- @@ -282,7 +299,7 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { return false; } return false; - } else if (particle_number == 2) { + } else if (particle_number == partTwo) { switch (tracktwofilter.confPDGCodePartTwo) { case 2212: // Proton+ case -2212: // Proton- @@ -313,19 +330,21 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { void init(InitContext&) { eventHisto.init(&qaRegistry); - sphericityRegistry.add("sphericity", ";Sphericity;Entries", kTH1F, {{150, 0.0, 3, "Sphericity"}}); trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true); trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true); - mixQaRegistry.add("MixingQA/hSECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); - mixQaRegistry.add("MixingQA/hMECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); + if (confFillDebug) { + sphericityRegistry.add("sphericity", ";Sphericity;Entries", kTH1F, {{150, 0.0, 3, "Sphericity"}}); + mixQaRegistry.add("MixingQA/hSECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); + mixQaRegistry.add("MixingQA/hMECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); + } mass1 = pdg->Mass(trackonefilter.confPDGCodePartOne); mass2 = pdg->Mass(tracktwofilter.confPDGCodePartTwo); - if (cfgProcessPM) { + if (processPair.cfgProcessPM) { sameEventCont.init(&resultRegistryPM, confkstarBins, confMultBins, confkTBins, confmTBins, confmultBins3D, confmTBins3D, twotracksconfigs.confIsMC, twotracksconfigs.confUse3D); mixedEventCont.init(&resultRegistryPM, confkstarBins, confMultBins, confkTBins, confmTBins, confmultBins3D, confmTBins3D, twotracksconfigs.confIsMC, twotracksconfigs.confUse3D); @@ -336,9 +355,23 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { sameEventMultCont.init(&sameMultRegistryPM, confkstarBins, confMultKstarBins, confKtKstarBins, cfgProcessKtBins, cfgProcessKtMt3DCF); mixedEventMultCont.init(&mixedMultRegistryPM, confkstarBins, confMultKstarBins, confKtKstarBins, cfgProcessKtBins, cfgProcessKtMt3DCF); } + if (doprocessFractions) { + mixedMultRegistryPM.add("MCreco/motherParticle", "pair fractions; part1 mother PDG; part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + mixedMultRegistryPM.add("MCreco/motherParticlePDGCheck", "pair fractions;part1 mother PDG;part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + } + if (doprocessFractionsMCTruth) { + mixedMultRegistryPM.add("MCtruth/motherParticle", "pair fractions; part1 mother PDG; part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + } } - if (cfgProcessPP) { + if (processPair.cfgProcessPP) { + if (doprocessFractions) { + mixedMultRegistryPP.add("MCreco/motherParticle", "pair fractions; part1 mother PDG; part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + mixedMultRegistryPP.add("MCreco/motherParticlePDGCheck", "pair fractions;part1 mother PDG;part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + } + if (doprocessFractionsMCTruth) { + mixedMultRegistryPP.add("MCtruth/motherParticle", "pair fractions; part1 mother PDG; part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + } sameEventContPP.init(&resultRegistryPP, confkstarBins, confMultBins, confkTBins, confmTBins, confmultBins3D, confmTBins3D, twotracksconfigs.confIsMC, twotracksconfigs.confUse3D); mixedEventContPP.init(&resultRegistryPP, confkstarBins, confMultBins, confkTBins, confmTBins, confmultBins3D, confmTBins3D, twotracksconfigs.confIsMC, twotracksconfigs.confUse3D); sameEventContPP.setPDGCodes(trackonefilter.confPDGCodePartOne, tracktwofilter.confPDGCodePartTwo); @@ -350,7 +383,14 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { } } - if (cfgProcessMM) { + if (processPair.cfgProcessMM) { + if (doprocessFractions) { + mixedMultRegistryMM.add("MCreco/motherParticle", "pair fractions; part1 mother PDG; part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + mixedMultRegistryMM.add("MCreco/motherParticlePDGCheck", "pair fractions;part1 mother PDG;part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + } + if (doprocessFractionsMCTruth) { + mixedMultRegistryMM.add("MCtruth/motherParticle", "pair fractions; part1 mother PDG; part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + } sameEventContMM.init(&resultRegistryMM, confkstarBins, confMultBins, confkTBins, confmTBins, confmultBins3D, confmTBins3D, twotracksconfigs.confIsMC, twotracksconfigs.confUse3D); mixedEventContMM.init(&resultRegistryMM, confkstarBins, confMultBins, confkTBins, confmTBins, confmultBins3D, confmTBins3D, twotracksconfigs.confIsMC, twotracksconfigs.confUse3D); sameEventContMM.setPDGCodes(trackonefilter.confPDGCodePartOne, tracktwofilter.confPDGCodePartTwo); @@ -375,7 +415,9 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { template void fillCollision(CollisionType col) { - mixQaRegistry.fill(HIST("MixingQA/hSECollisionBins"), colBinning.getBin({col.posZ(), col.multV0M()})); + if (confFillDebug) { + mixQaRegistry.fill(HIST("MixingQA/hSECollisionBins"), colBinning.getBin({col.posZ(), col.multV0M()})); + } eventHisto.fillQA(col); } @@ -395,8 +437,11 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { void doSameEvent(PartitionType groupPartsOne, PartitionType groupPartsTwo, PartType parts, float magFieldTesla, int multCol, int pairType, bool fillQA) { + enum PairTypes { type1 = 1, + type2 = 2, + type3 = 3 }; /// Histogramming same event - if ((pairType == 1 || pairType == 2) && fillQA) { + if ((pairType == type1 || pairType == type2) && fillQA) { for (const auto& part : groupPartsOne) { if (!isParticleNSigma((int8_t)1, part.p(), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon))) { continue; @@ -405,7 +450,7 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { } } - if ((pairType == 1 || pairType == 3) && fillQA) { + if ((pairType == type1 || pairType == type3) && fillQA) { for (const auto& part : groupPartsTwo) { if (!isParticleNSigma((int8_t)2, part.p(), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon))) { continue; @@ -414,7 +459,7 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { } } - if (pairType == 1) { + if (pairType == type1) { /// Now build the combinations for non-identical particle pairs for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) { @@ -522,30 +567,32 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { FilteredFemtoFullParticles const& parts) { fillCollision(col); - sphericityRegistry.fill(HIST("sphericity"), col.sphericity()); + if (confFillDebug) { + sphericityRegistry.fill(HIST("sphericity"), col.sphericity()); + } auto thegroupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); auto thegroupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); bool fillQA = true; - if (cfgProcessPM) { + if (processPair.cfgProcessPM) { doSameEvent(thegroupPartsOne, thegroupPartsTwo, parts, col.magField(), col.multV0M(), 1, fillQA); fillQA = false; } - if (cfgProcessPP) + if (processPair.cfgProcessPP) doSameEvent(thegroupPartsOne, thegroupPartsOne, parts, col.magField(), col.multV0M(), 2, fillQA); - if (cfgProcessMM) + if (processPair.cfgProcessMM) doSameEvent(thegroupPartsTwo, thegroupPartsTwo, parts, col.magField(), col.multV0M(), 3, fillQA); } PROCESS_SWITCH(FemtoUniversePairTaskTrackTrackMultKtExtended, processSameEvent, "Enable processing same event", true); /// process function for to call doSameEvent with Monte Carlo /// \param col subscribe to the collision table (Monte Carlo Reconstructed reconstructed) - /// \param parts subscribe to joined table FemtoUniverseParticles and FemtoUniverseMCLables to access Monte Carlo truth + /// \param parts subscribe to joined table FemtoUniverseParticles and FemtoUniverseMCLabels to access Monte Carlo truth /// \param FemtoUniverseMCParticles subscribe to the Monte Carlo truth table void processSameEventMC(o2::aod::FdCollision const& col, - soa::Join const& parts, + FemtoRecoParticles const& parts, o2::aod::FdMCParticles const&) { fillCollision(col); @@ -554,13 +601,13 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { auto thegroupPartsTwo = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); bool fillQA = true; - if (cfgProcessPM) { - doSameEvent(thegroupPartsOne, thegroupPartsTwo, parts, col.magField(), col.multV0M(), 1, fillQA); + if (processPair.cfgProcessPM) { + doSameEvent(thegroupPartsOne, thegroupPartsTwo, parts, col.magField(), col.multV0M(), 1, fillQA); fillQA = false; } - if (cfgProcessPP) + if (processPair.cfgProcessPP) doSameEvent(thegroupPartsOne, thegroupPartsOne, parts, col.magField(), col.multV0M(), 2, fillQA); - if (cfgProcessMM) + if (processPair.cfgProcessMM) doSameEvent(thegroupPartsTwo, thegroupPartsTwo, parts, col.magField(), col.multV0M(), 3, fillQA); } PROCESS_SWITCH(FemtoUniversePairTaskTrackTrackMultKtExtended, processSameEventMC, "Enable processing same event for Monte Carlo", false); @@ -661,7 +708,9 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { for (const auto& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, cols, cols)) { const int multiplicityCol = collision1.multV0M(); - mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + if (confFillDebug) { + mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + } const auto& magFieldTesla1 = collision1.magField(); const auto& magFieldTesla2 = collision2.magField(); @@ -670,17 +719,17 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { continue; } - if (cfgProcessPM) { + if (processPair.cfgProcessPM) { auto groupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); doMixedEvent(groupPartsOne, groupPartsTwo, parts, magFieldTesla1, multiplicityCol, 1); } - if (cfgProcessPP) { + if (processPair.cfgProcessPP) { auto groupPartsOne = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsOne->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); doMixedEvent(groupPartsOne, groupPartsTwo, parts, magFieldTesla1, multiplicityCol, 2); } - if (cfgProcessMM) { + if (processPair.cfgProcessMM) { auto groupPartsOne = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsTwo->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); doMixedEvent(groupPartsOne, groupPartsTwo, parts, magFieldTesla1, multiplicityCol, 3); @@ -694,14 +743,15 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { /// \param parts subscribe to joined table FemtoUniverseParticles and FemtoUniverseMCLables to access Monte Carlo truth /// \param FemtoUniverseMCParticles subscribe to the Monte Carlo truth table void processMixedEventMC(o2::aod::FdCollisions const& cols, - soa::Join const& parts, + FemtoRecoParticles const& parts, o2::aod::FdMCParticles const&) { for (const auto& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, cols, cols)) { const int multiplicityCol = collision1.multV0M(); - mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); - + if (confFillDebug) { + mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + } const auto& magFieldTesla1 = collision1.magField(); const auto& magFieldTesla2 = collision2.magField(); @@ -711,24 +761,198 @@ struct FemtoUniversePairTaskTrackTrackMultKtExtended { /// \todo before mixing we should check whether both collisions contain a pair of particles! // if (partsOne.size() == 0 || nPart2Evt1 == 0 || nPart1Evt2 == 0 || partsTwo.size() == 0 ) continue; - if (cfgProcessPM) { + if (processPair.cfgProcessPM) { auto groupPartsOne = partsOneMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); doMixedEvent(groupPartsOne, groupPartsTwo, parts, magFieldTesla1, multiplicityCol, 1); } - if (cfgProcessPP) { + if (processPair.cfgProcessPP) { auto groupPartsOne = partsOneMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsOneMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); doMixedEvent(groupPartsOne, groupPartsTwo, parts, magFieldTesla1, multiplicityCol, 2); } - if (cfgProcessMM) { + if (processPair.cfgProcessMM) { auto groupPartsOne = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); auto groupPartsTwo = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); doMixedEvent(groupPartsOne, groupPartsTwo, parts, magFieldTesla1, multiplicityCol, 3); } } } - PROCESS_SWITCH(FemtoUniversePairTaskTrackTrackMultKtExtended, processMixedEventMC, "Enable processing mixed events MC", false); + PROCESS_SWITCH(FemtoUniversePairTaskTrackTrackMultKtExtended, processMixedEventMC, "Enable processing mixed events MC reco", false); + + /// brief process function for to call process fractions with Monte Carlo reco + /// \param cols subscribe to the collisions table (Monte Carlo Reconstructed reconstructed) + /// \param parts subscribe to joined table FemtoUniverseParticles and FemtoUniverseMCLables to access Monte Carlo truth + /// \param FemtoUniverseMCParticles subscribe to the Monte Carlo truth table + void processFractions(o2::aod::FdCollisions const& cols, + FemtoRecoParticles const& parts, + o2::aod::FdMCParticles const& mparts) + { + + auto doFractions = [&](auto& p1, auto& p2, auto& magFieldTesla, int partType) -> void { + if (!isParticleNSigma((int8_t)2, p1.p(), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon))) { + return; + } + + if (!isParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon))) { + return; + } + + if (confIsCPR.value) { + if (pairCloseRejection.isClosePair(p1, p2, parts, magFieldTesla, femto_universe_container::EventType::mixed)) { + return; + } + } + switch (partType) { + case 1: + mixedMultRegistryPM.fill(HIST("MCreco/motherParticle"), p1.motherPDG(), p2.motherPDG()); + break; + case 2: + mixedMultRegistryPP.fill(HIST("MCreco/motherParticle"), p1.motherPDG(), p2.motherPDG()); + break; + case 3: + mixedMultRegistryMM.fill(HIST("MCreco/motherParticle"), p1.motherPDG(), p2.motherPDG()); + break; + } + auto mcPartId1 = p1.fdMCParticleId(); + if (mcPartId1 == -1) + return; + auto mcPartId2 = p2.fdMCParticleId(); + if (mcPartId2 == -1) + return; + const auto& mcParticle1 = mparts.iteratorAt(mcPartId1); + const auto& mcParticle2 = mparts.iteratorAt(mcPartId2); + switch (partType) { + case 1: + if ((trackonefilter.confPDGCodePartOne != mcParticle1.pdgMCTruth()) || (tracktwofilter.confPDGCodePartTwo != mcParticle2.pdgMCTruth())) + return; + mixedMultRegistryPM.fill(HIST("MCreco/motherParticlePDGCheck"), p1.motherPDG(), p2.motherPDG()); + break; + case 2: + if ((trackonefilter.confPDGCodePartOne != mcParticle1.pdgMCTruth()) || (trackonefilter.confPDGCodePartOne != mcParticle2.pdgMCTruth())) + return; + mixedMultRegistryPP.fill(HIST("MCreco/motherParticlePDGCheck"), p1.motherPDG(), p2.motherPDG()); + break; + case 3: + if ((tracktwofilter.confPDGCodePartTwo != mcParticle1.pdgMCTruth()) || (tracktwofilter.confPDGCodePartTwo != mcParticle2.pdgMCTruth())) + return; + mixedMultRegistryMM.fill(HIST("MCreco/motherParticlePDGCheck"), p1.motherPDG(), p2.motherPDG()); + break; + } + }; + + for (const auto& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, cols, cols)) { + const int multiplicityCol = collision1.multV0M(); + if (confFillDebug) { + mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + } + const auto& magFieldTesla1 = collision1.magField(); + const auto& magFieldTesla2 = collision2.magField(); + + if (magFieldTesla1 != magFieldTesla2) { + continue; + } + + auto groupPartsOne = partsOneMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); + auto groupPartsTwo = partsTwoMC->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); + + if (processPair.cfgProcessPM) { + for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) { + doFractions(p1, p2, magFieldTesla1, 1); + } + } + if (processPair.cfgProcessPP) { + for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsOne))) { + doFractions(p1, p2, magFieldTesla1, 2); + } + } + if (processPair.cfgProcessMM) { + for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsTwo, groupPartsTwo))) { + doFractions(p1, p2, magFieldTesla1, 3); + } + } + } + } + PROCESS_SWITCH(FemtoUniversePairTaskTrackTrackMultKtExtended, processFractions, "Enable processing fractions for MC reco", false); + + /// This function processes fills histograms for fractions analysis + /// \todo the trivial loops over the collisions and tracks should be factored out since they will be common to all combinations of T-T, T-V0, V0-V0, ... + /// \tparam PartitionType + /// \tparam PartType + /// \tparam isMC: enables Monte Carlo truth specific histograms + /// \param groupPartsOne partition for the first particle passed by the process function + /// \param groupPartsTwo partition for the second particle passed by the process function + template + void doFractionsMCTruth(PartitionType groupPartsOne, PartitionType groupPartsTwo, int ContType) + { + + for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) { + + int pdgCodePartOne = static_cast(p1.pidCut()); + const auto& pdgParticleOne = pdg->GetParticle(pdgCodePartOne); + int pdgCodePartTwo = static_cast(p2.pidCut()); + const auto& pdgParticleTwo = pdg->GetParticle(pdgCodePartTwo); + switch (ContType) { + case 1: { + if (pdgParticleOne && pdgParticleTwo && (pdgCodePartOne == trackonefilter.confPDGCodePartOne) && (pdgCodePartTwo == tracktwofilter.confPDGCodePartTwo)) { + continue; + } + mixedMultRegistryPM.fill(HIST("MCtruth/motherParticle"), p1.tempFitVar(), p2.tempFitVar()); + break; + } + case 2: { + if (pdgParticleOne && pdgParticleTwo && (pdgCodePartOne == trackonefilter.confPDGCodePartOne) && (pdgCodePartTwo == trackonefilter.confPDGCodePartOne)) { + continue; + } + mixedMultRegistryPP.fill(HIST("MCtruth/motherParticle"), p1.tempFitVar(), p2.tempFitVar()); + break; + } + case 3: { + if (pdgParticleOne && pdgParticleTwo && (pdgCodePartOne == tracktwofilter.confPDGCodePartTwo) && (pdgCodePartTwo == tracktwofilter.confPDGCodePartTwo)) { + continue; + } + mixedMultRegistryMM.fill(HIST("MCtruth/motherParticle"), p1.tempFitVar(), p2.tempFitVar()); + } + default: + break; + } + } + } + + /// brief process function for to call doFractionsMCTruth with Monte Carlo + /// \param cols subscribe to the collisions table (Monte Carlo Reconstructed reconstructed) + /// \param FemtoTruthParticles subscribe to the Monte Carlo truth particle table + void processFractionsMCTruth(o2::aod::FdCollisions const& cols, + FemtoTruthParticles const&) + { + for (const auto& [collision1, collision2] : soa::selfCombinations(colBinning, 5, -1, cols, cols)) { + + const int multiplicityCol = collision1.multV0M(); + if (confFillDebug) { + mixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinning.getBin({collision1.posZ(), multiplicityCol})); + } + const auto& magFieldTesla1 = collision1.magField(); + const auto& magFieldTesla2 = collision2.magField(); + + if (magFieldTesla1 != magFieldTesla2) { + continue; + } + /// \todo before mixing we should check whether both collisions contain a pair of particles! + // if (partsOne.size() == 0 || nPart2Evt1 == 0 || nPart1Evt2 == 0 || partsTwo.size() == 0 ) continue; + auto groupPartsOne = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache); + auto groupPartsTwo = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache); + if (processPair.cfgProcessPM) { + doFractionsMCTruth(groupPartsOne, groupPartsTwo, 1); + } + if (processPair.cfgProcessPP) { + doFractionsMCTruth(groupPartsOne, groupPartsOne, 2); + } + if (processPair.cfgProcessMM) { + doFractionsMCTruth(groupPartsTwo, groupPartsTwo, 3); + } + } + } + PROCESS_SWITCH(FemtoUniversePairTaskTrackTrackMultKtExtended, processFractionsMCTruth, "Enable processing residual correlations fractions for MC truth", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Extended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Extended.cxx index f158b23b726..34639df3919 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Extended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Extended.cxx @@ -282,6 +282,7 @@ struct FemtoUniversePairTaskTrackV0Extended { registryMCreco.add("minus/MCrecoPrPt", "MC reco protons;#it{p}_{T} (GeV/c)", {HistType::kTH1F, {{500, 0, 5}}}); registryMCreco.add("mothersReco/motherParticle", "pair fractions;part1 mother PDG;part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); + registryMCreco.add("mothersReco/motherParticlePDGCheck", "pair fractions;part1 mother PDG;part2 mother PDG", {HistType::kTH2F, {{8001, -4000, 4000}, {8001, -4000, 4000}}}); sameEventCont.init(&resultRegistry, confkstarBins, confMultBins, confkTBins, confmTBins, confMultBins3D, confmTBins3D, confEtaBins, confPhiBins, confIsMC, confUse3D); sameEventCont.setPDGCodes(confTrkPDGCodePartOne, ConfV0Selection.confV0PDGCodePartTwo); @@ -923,7 +924,7 @@ struct FemtoUniversePairTaskTrackV0Extended { PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Extended, processMCTruth, "Process MC truth data", false); - void processPairFractions(FilteredFDCollisions const& cols, FemtoRecoParticles const& parts) + void processPairFractions(FilteredFDCollisions const& cols, FemtoRecoParticles const& parts, aod::FdMCParticles const& mcparts) { ColumnBinningPolicy colBinningMult{{confVtxBins, confMultBins}, true}; ColumnBinningPolicy colBinningCent{{confVtxBins, confMultBins}, true}; @@ -963,6 +964,18 @@ struct FemtoUniversePairTaskTrackV0Extended { } } registryMCreco.fill(HIST("mothersReco/motherParticle"), p1.motherPDG(), p2.motherPDG()); + + auto mcPartId1 = p1.fdMCParticleId(); + if (mcPartId1 == -1) + continue; + auto mcPartId2 = p2.fdMCParticleId(); + if (mcPartId2 == -1) + continue; + const auto& mcParticle1 = mcparts.iteratorAt(mcPartId1); + const auto& mcParticle2 = mcparts.iteratorAt(mcPartId2); + if (mcParticle1.pdgMCTruth() == confTrkPDGCodePartOne && mcParticle2.pdgMCTruth() == ConfV0Selection.confV0PDGCodePartTwo) { + registryMCreco.fill(HIST("mothersReco/motherParticlePDGCheck"), p1.motherPDG(), p2.motherPDG()); + } } }; @@ -972,7 +985,7 @@ struct FemtoUniversePairTaskTrackV0Extended { } PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Extended, processPairFractions, "Process MC data to obtain pair fractions", false); - void processPairFractionsV0(FilteredFDCollisions const& cols, FemtoRecoParticles const& parts) + void processPairFractionsV0(FilteredFDCollisions const& cols, FemtoRecoParticles const& parts, aod::FdMCParticles const& mcparts) { ColumnBinningPolicy colBinningMult{{confVtxBins, confMultBins}, true}; ColumnBinningPolicy colBinningCent{{confVtxBins, confMultBins}, true}; @@ -1021,6 +1034,19 @@ struct FemtoUniversePairTaskTrackV0Extended { } registryMCreco.fill(HIST("mothersReco/motherParticle"), p1.motherPDG(), p2.motherPDG()); + auto mcPartId1 = p1.fdMCParticleId(); + if (mcPartId1 == -1) + continue; + auto mcPartId2 = p2.fdMCParticleId(); + if (mcPartId2 == -1) + continue; + const auto& mcParticle1 = mcparts.iteratorAt(mcPartId1); + const auto& mcParticle2 = mcparts.iteratorAt(mcPartId2); + if ((ConfV0Selection.confV0Type1 == 0 && mcParticle1.pdgMCTruth() != kLambda0) || (ConfV0Selection.confV0Type1 == 1 && mcParticle1.pdgMCTruth() != kLambda0Bar)) + continue; + if ((ConfV0Selection.confV0Type2 == 0 && mcParticle2.pdgMCTruth() != kLambda0) || (ConfV0Selection.confV0Type2 == 1 && mcParticle2.pdgMCTruth() != kLambda0Bar)) + continue; + registryMCreco.fill(HIST("mothersReco/motherParticlePDGCheck"), p1.motherPDG(), p2.motherPDG()); } };