From 04d537f5be1e20eea189235b4539c3038cbf8392 Mon Sep 17 00:00:00 2001 From: Anna-Mariia Andrushko <01175703@pw.edu.pl> Date: Wed, 3 Sep 2025 12:55:22 +0200 Subject: [PATCH] Add new MCTruth analysis --- .../femtoUniverseProducerTask.cxx | 212 +++++++++++++----- .../femtoUniversePairTaskTrackV0Helicity.cxx | 176 ++++++++------- 2 files changed, 246 insertions(+), 142 deletions(-) diff --git a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx index 37a16abc46b..23d14062fee 100644 --- a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx +++ b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx @@ -74,24 +74,16 @@ using namespace o2::constants::physics; namespace o2::aod { -using FemtoFullCollision = - soa::Join::iterator; -using FemtoFullCollisionCentPP = - soa::Join::iterator; -using FemtoFullCollisionCentRun2 = - soa::Join::iterator; -using FemtoFullCollisionCentRun3 = - soa::Join::iterator; +using FemtoFullCollision = soa::Join::iterator; +using FemtoFullCollisionCentPP = soa::Join::iterator; +using FemtoFullCollisionCentRun2 = soa::Join::iterator; +using FemtoFullCollisionCentRun3 = soa::Join::iterator; using FemtoFullCollisionMC = soa::Join::iterator; -using FemtoFullCollisionCentRun3MCs = - soa::Join; -using FemtoFullCollisionCentRun3MC = - soa::Join::iterator; -using FemtoFullTracks = - soa::Join; +using FemtoFullCollisionCentRun3MCs = soa::Join; +using FemtoFullCollisionCentRun3MC = soa::Join::iterator; +using FemtoFullTracks = soa::Join; // using FilteredFullV0s = soa::Filtered; /// predefined Join // table for o2::aod::V0s = soa::Join @@ -232,6 +224,8 @@ struct FemtoUniverseProducerTask { Configurable confV0RejectKaons{"confV0RejectKaons", false, "Switch to reject kaons"}; Configurable confV0InvKaonMassLowLimit{"confV0InvKaonMassLowLimit", 0.48, "Lower limit of the V0 invariant mass for Kaon rejection"}; Configurable confV0InvKaonMassUpLimit{"confV0InvKaonMassUpLimit", 0.515, "Upper limit of the V0 invariant mass for Kaon rejection"}; + + Configurable> confV0PDGMCTruth{"confV0PDGMCTruth", std::vector{2212, -211, 3122}, "PDG codes of V0 daughters and mother, the order must be as follows -- positive daughter, negative daughter, mother"}; } ConfV0Selection; struct : o2::framework::ConfigurableGroup { @@ -875,6 +869,37 @@ struct FemtoUniverseProducerTask { } } + template + void fillMCTruthParticle(MCParticleType const& mcparticle, o2::aod::femtouniverseparticle::ParticleType fdparttype) + { + auto pdgCode = mcparticle.pdgCode(); + int particleOrigin = 99; + + // Determine particle origin + if (std::abs(pdgCode) == std::abs(confPDGCodePartOne.value) || std::abs(pdgCode) == std::abs(confPDGCodePartTwo.value)) { + if (mcparticle.isPhysicalPrimary()) { + particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kPrimary; + } else { + auto mothers = mcparticle.template mothers_as(); + if (!mothers.empty()) { + auto mother = mothers.front(); + if (mother.producedByGenerator()) { + particleOrigin = checkDaughterType(fdparttype, mother.pdgCode()); + } else { + particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kMaterial; + } + } + } + } else { + particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kFake; + } + + // Fill MC companion tables + outputPartsMC(particleOrigin, pdgCode, mcparticle.pt(), mcparticle.eta(), mcparticle.phi()); + fillDebugParticleMC(mcparticle); + outputPartsMCLabels(outputPartsMC.lastIndex()); + } + template void fillMCParticlePhi(ParticleType const& kaon1, ParticleType const& kaon2) { @@ -1244,6 +1269,77 @@ struct FemtoUniverseProducerTask { } } + template + void fillV0MCTruth(MCParticlesType const& mcParticles) + { + for (const auto& mc : mcParticles) { // Loop over all MC Truth particles + if (mc.pdgCode() != ConfV0Selection.confV0PDGMCTruth.value[2]) + continue; // Artificially single out V0s + + auto daughters = mc.template daughters_as(); // Access daughters (no differentiation of signs, it needs to be checked separately) + + bool foundPos = false, foundNeg = false; + + std::vector childIDs = {0, 0}; + int rowPos = 0; + int rowNeg = 0; + + for (auto const& d : daughters) { // Loop over daughters + if (d.pdgCode() == ConfV0Selection.confV0PDGMCTruth.value[0]) { // Check for a positive child + foundPos = true; + + outputParts(outputCollision.lastIndex(), + d.pt(), + d.eta(), + d.phi(), + aod::femtouniverseparticle::ParticleType::kMCTruthTrack, + 0, + 0, + d.pdgCode(), + childIDs, // {0, 0} + 0, + 0); + rowPos = outputParts.lastIndex(); + fillMCTruthParticle(d, aod::femtouniverseparticle::ParticleType::kMCTruthTrack); + } else if (d.pdgCode() == ConfV0Selection.confV0PDGMCTruth.value[1]) { // Check for a negative child + foundNeg = true; + + outputParts(outputCollision.lastIndex(), + d.pt(), + d.eta(), + d.phi(), + aod::femtouniverseparticle::ParticleType::kMCTruthTrack, + 0, + 0, + d.pdgCode(), + childIDs, // {0, 0} + 0, + 0); + rowNeg = outputParts.lastIndex(); + fillMCTruthParticle(d, aod::femtouniverseparticle::ParticleType::kMCTruthTrack); + } + } + + if (!foundPos || !foundNeg) + continue; + + childIDs[0] = rowPos; + childIDs[1] = rowNeg; + outputParts(outputCollision.lastIndex(), + mc.pt(), + mc.eta(), + mc.phi(), + aod::femtouniverseparticle::ParticleType::kMCTruthTrack, + 0, + 0, + mc.pdgCode(), + childIDs, + 0, + 0); + fillMCTruthParticle(mc, aod::femtouniverseparticle::ParticleType::kMCTruthTrack); + } + } + template void fillCascade(CollisionType const& col, CascadeType const& fullCascades, TrackType const&) { @@ -2429,47 +2525,6 @@ struct FemtoUniverseProducerTask { } PROCESS_SWITCH(FemtoUniverseProducerTask, processTruthAndFullMCV0, "Provide both MC truth and reco for tracks and V0s", false); - void processTruthAndFullMCCentRun3V0( - aod::McCollisions const& mccols, - aod::McParticles const& mcParticles, - aod::FemtoFullCollisionCentRun3MCs const& collisions, - soa::Filtered> const& tracks, - soa::Join const& fullV0s, - aod::BCsWithTimestamps const&) - { - - // recos - std::set recoMcIds; - for (const auto& col : collisions) { - auto groupedTracks = tracks.sliceBy(perCollisionTracks, col.globalIndex()); - auto groupedV0Parts = fullV0s.sliceBy(perCollisionV0s, col.globalIndex()); - getMagneticFieldTesla(col.bc_as()); - const auto colcheck = fillCollisionsCentRun3(col); - if (colcheck) { - fillTracks(groupedTracks); - fillV0(col, groupedV0Parts, groupedTracks); - } - for (const auto& track : groupedTracks) { - if (trackCuts.isSelectedMinimal(track)) - recoMcIds.insert(track.mcParticleId()); - } - } - - // truth - for (const auto& mccol : mccols) { - auto groupedCollisions = collisions.sliceBy(recoCollsPerMCCollCentPbPb, mccol.globalIndex()); - for (const auto& col : groupedCollisions) { - const auto colcheck = fillMCTruthCollisionsCentRun3(col); // fills the reco collisions for mc collision - if (colcheck) { - auto groupedMCParticles = mcParticles.sliceBy(perMCCollision, mccol.globalIndex()); - outputCollExtra(1.0, 1.0); - fillParticles(groupedMCParticles, recoMcIds); // fills mc particles - } - } - } - } - PROCESS_SWITCH(FemtoUniverseProducerTask, processTruthAndFullMCCentRun3V0, "Provide both MC truth and reco for tracks and V0s with centrality", false); - Preslice> perCollisionCascs = aod::track::collisionId; void processTruthAndFullMCCasc( aod::McCollisions const& mccols, @@ -2673,6 +2728,47 @@ struct FemtoUniverseProducerTask { } PROCESS_SWITCH(FemtoUniverseProducerTask, processV0CentRun3Data, "Provide experimental data for Run 3 with centrality for track track", false); + void processTruthAndFullMCCentRun3V0( + aod::McCollisions const& mccols, + aod::McParticles const& mcParticles, + aod::FemtoFullCollisionCentRun3MCs const& collisions, + soa::Filtered> const& tracks, + soa::Join const& fullV0s, + aod::BCsWithTimestamps const&) + { + + // MCReco + std::set recoMcIds; + for (const auto& col : collisions) { // loop over collisions + auto groupedTracks = tracks.sliceBy(perCollisionTracks, col.globalIndex()); // slicing for tracks + auto groupedV0Parts = fullV0s.sliceBy(perCollisionV0s, col.globalIndex()); // slicing for V0 + getMagneticFieldTesla(col.bc_as()); + const auto colcheck = fillCollisionsCentRun3(col); + if (colcheck) { + fillTracks(groupedTracks); + fillV0(col, groupedV0Parts, groupedTracks); + } + for (const auto& track : groupedTracks) { + if (trackCuts.isSelectedMinimal(track)) + recoMcIds.insert(track.mcParticleId()); + } + } + + // MCTruth + for (const auto& mccol : mccols) { + auto groupedCollisions = collisions.sliceBy(recoCollsPerMCCollCentPbPb, mccol.globalIndex()); // slicing for MC collisions + auto groupedMCParticles = mcParticles.sliceBy(perMCCollision, mccol.globalIndex()); // slicing for MC particles + for (const auto& col : groupedCollisions) { + const auto colcheck = fillMCTruthCollisionsCentRun3(col); + if (colcheck) { + outputCollExtra(1.0, 1.0); + fillV0MCTruth(groupedMCParticles); // fills MC V0s and its daughters + } + } + } + } + PROCESS_SWITCH(FemtoUniverseProducerTask, processTruthAndFullMCCentRun3V0, "Provide both MC truth and reco for tracks and V0s with centrality", false); + void processCascadeCentRun3Data(aod::FemtoFullCollisionCentRun3 const& col, aod::BCsWithTimestamps const&, soa::Filtered const& tracks, diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Helicity.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Helicity.cxx index bd66de7dc37..0801d404a4a 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Helicity.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackV0Helicity.cxx @@ -53,12 +53,16 @@ struct FemtoUniversePairTaskTrackV0Helicity { Service pdgMC; SliceCache cache; + using FemtoFullParticles = soa::Join; Preslice perCol = aod::femtouniverseparticle::fdCollisionId; using FemtoRecoParticles = soa::Join; Preslice perColMC = aod::femtouniverseparticle::fdCollisionId; + using FemtoTruthParticles = soa::Join; + Preslice perColMCTruth = aod::femtouniverseparticle::fdCollisionId; + /// To apply narrow cut Configurable confZVertexCut{"confZVertexCut", 10.f, "Event sel: Maximum z-Vertex (cm)"}; Configurable confEta{"confEta", 0.8, "Eta cut for the global track"}; @@ -89,25 +93,28 @@ struct FemtoUniversePairTaskTrackV0Helicity { FemtoUniverseParticleHisto trackHistoPartOneNeg; /// Particle 2 (V0) - Configurable confV0PDGCodePartTwo{"confV0PDGCodePartTwo", 3122, "Particle 2 (V0) - PDG code"}; - ConfigurableAxis confV0TempFitVarBins{"confV0TempFitVarBins", {300, 0.95, 1.}, "V0: binning of the TempFitVar in the pT vs. TempFitVar plot"}; - ConfigurableAxis confV0TempFitVarpTBins{"confV0TempFitVarpTBins", {20, 0.5, 4.05}, "V0: pT binning of the pT vs. TempFitVar plot"}; - Configurable confV0Type1{"confV0Type1", 0, "select one of the V0s (lambda = 0, anti-lambda = 1, k0 = 2) for v0-v0 and Track-v0 combination"}; - Configurable confV0Type2{"confV0Type2", 0, "select one of the V0s (lambda = 0, anti-lambda = 1, k0 = 2) for v0-v0 combination"}; - Configurable confV0InvMassLowLimit{"confV0InvMassLowLimit", 1.10, "Lower limit of the V0 invariant mass"}; - Configurable confV0InvMassUpLimit{"confV0InvMassUpLimit", 1.13, "Upper limit of the V0 invariant mass"}; - ConfigurableAxis confChildTempFitVarBins{"confChildTempFitVarBins", {300, -0.15, 0.15}, "V0 child: binning of the TempFitVar in the pT vs. TempFitVar plot"}; - ConfigurableAxis confChildTempFitVarpTBins{"confChildTempFitVarpTBins", {20, 0.5, 4.05}, "V0 child: pT binning of the pT vs. TempFitVar plot"}; - Configurable confHPtPart2{"confHPtPart2", 4.0f, "higher limit for pt of particle 2"}; - Configurable confLPtPart2{"confLPtPart2", 0.3f, "lower limit for pt of particle 2"}; - Configurable confPDGCodeV0{"confPDGCodeV0", 3122, "V0 -- PDG code"}; - Configurable confPDGCodePosChild{"confPDGCodePosChild", 2212, "Positive Child -- PDG code"}; - Configurable confPDGCodeNegChild{"confPDGCodeNegChild", 211, "Negative Child -- PDG code"}; + struct : o2::framework::ConfigurableGroup { + Configurable confV0PDGCodePartTwo{"confV0PDGCodePartTwo", 3122, "Particle 2 (V0) - PDG code"}; + ConfigurableAxis confV0TempFitVarBins{"confV0TempFitVarBins", {300, 0.95, 1.}, "V0: binning of the TempFitVar in the pT vs. TempFitVar plot"}; + ConfigurableAxis confV0TempFitVarpTBins{"confV0TempFitVarpTBins", {20, 0.5, 4.05}, "V0: pT binning of the pT vs. TempFitVar plot"}; + Configurable confV0Type1{"confV0Type1", 0, "select one of the V0s (lambda = 0, anti-lambda = 1, k0 = 2) for v0-v0 and Track-v0 combination"}; + Configurable confV0Type2{"confV0Type2", 0, "select one of the V0s (lambda = 0, anti-lambda = 1, k0 = 2) for v0-v0 combination"}; + Configurable confV0InvMassLowLimit{"confV0InvMassLowLimit", 1.10, "Lower limit of the V0 invariant mass"}; + Configurable confV0InvMassUpLimit{"confV0InvMassUpLimit", 1.13, "Upper limit of the V0 invariant mass"}; + ConfigurableAxis confChildTempFitVarBins{"confChildTempFitVarBins", {300, -0.15, 0.15}, "V0 child: binning of the TempFitVar in the pT vs. TempFitVar plot"}; + ConfigurableAxis confChildTempFitVarpTBins{"confChildTempFitVarpTBins", {20, 0.5, 4.05}, "V0 child: pT binning of the pT vs. TempFitVar plot"}; + Configurable confHPtPart2{"confHPtPart2", 4.0f, "higher limit for pt of particle 2"}; + Configurable confLPtPart2{"confLPtPart2", 0.3f, "lower limit for pt of particle 2"}; + Configurable confPDGCodeV0{"confPDGCodeV0", 3122, "V0 -- PDG code"}; + Configurable confPDGCodePosChild{"confPDGCodePosChild", 2212, "Positive Child -- PDG code"}; + Configurable confPDGCodeNegChild{"confPDGCodeNegChild", 211, "Negative Child -- PDG code"}; + } V0configs; /// Partition for particle 2 - Partition partsTwo = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kV0)) && (aod::femtouniverseparticle::pt < confHPtPart2) && (aod::femtouniverseparticle::pt > confLPtPart2); - Partition partsTwoMC = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)) && (aod::femtouniverseparticle::pt < confHPtPart2) && (aod::femtouniverseparticle::pt > confLPtPart2); - Partition partsTwoMCReco = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kV0)) && (aod::femtouniverseparticle::pt < confHPtPart2) && (aod::femtouniverseparticle::pt > confLPtPart2); + Partition partsTwo = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kV0)) && (aod::femtouniverseparticle::pt < V0configs.confHPtPart2) && (aod::femtouniverseparticle::pt > V0configs.confLPtPart2); + Partition partsTwoMC = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)) && (aod::femtouniverseparticle::pt < V0configs.confHPtPart2) && (aod::femtouniverseparticle::pt > V0configs.confLPtPart2); + Partition partsTwoMCReco = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kV0)) && (aod::femtouniverseparticle::pt < V0configs.confHPtPart2) && (aod::femtouniverseparticle::pt > V0configs.confLPtPart2); + Partition partsTwoMCTruth = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)) && (aod::femtouniverseparticle::pt < V0configs.confHPtPart2) && (aod::femtouniverseparticle::pt > V0configs.confLPtPart2); /// Histogramming for particle 2 FemtoUniverseParticleHisto trackHistoPartTwo; @@ -182,7 +189,7 @@ struct FemtoUniversePairTaskTrackV0Helicity { bool invMLambda(float invMassLambda, float invMassAntiLambda) { - if ((invMassLambda < confV0InvMassLowLimit || invMassLambda > confV0InvMassUpLimit) && (invMassAntiLambda < confV0InvMassLowLimit || invMassAntiLambda > confV0InvMassUpLimit)) { + if ((invMassLambda < V0configs.confV0InvMassLowLimit || invMassLambda > V0configs.confV0InvMassUpLimit) && (invMassAntiLambda < V0configs.confV0InvMassLowLimit || invMassAntiLambda > V0configs.confV0InvMassUpLimit)) { return false; } return true; @@ -224,16 +231,16 @@ struct FemtoUniversePairTaskTrackV0Helicity { qaRegistry.add("Tracks_neg/nSigmaTOF", "; #it{p} (GeV/#it{c}); n#sigma_{TOF}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); trackHistoPartOnePos.init(&qaRegistry, confTrkTempFitVarpTBins, confTrkTempFitVarBins, confIsMC, confTrkPDGCodePartOne); trackHistoPartOneNeg.init(&qaRegistry, confTrkTempFitVarpTBins, confTrkTempFitVarBins, confIsMC, confTrkPDGCodePartOne); - trackHistoPartTwo.init(&qaRegistry, confV0TempFitVarpTBins, confV0TempFitVarBins, confIsMC, confV0PDGCodePartTwo, true); - posChildHistos.init(&qaRegistry, confChildTempFitVarpTBins, confChildTempFitVarBins, false, 0, true); - negChildHistos.init(&qaRegistry, confChildTempFitVarpTBins, confChildTempFitVarBins, false, 0, true); + trackHistoPartTwo.init(&qaRegistry, V0configs.confV0TempFitVarpTBins, V0configs.confV0TempFitVarBins, confIsMC, V0configs.confV0PDGCodePartTwo, true); + posChildHistos.init(&qaRegistry, V0configs.confChildTempFitVarpTBins, V0configs.confChildTempFitVarBins, false, 0, true); + negChildHistos.init(&qaRegistry, V0configs.confChildTempFitVarpTBins, V0configs.confChildTempFitVarBins, false, 0, true); - trackHistoV0Type1.init(&qaRegistry, confV0TempFitVarpTBins, confV0TempFitVarBins, confIsMC, confV0PDGCodePartTwo, true, "V0Type1"); - posChildV0Type1.init(&qaRegistry, confChildTempFitVarpTBins, confChildTempFitVarBins, false, 0, true, "posChildV0Type1"); - negChildV0Type1.init(&qaRegistry, confChildTempFitVarpTBins, confChildTempFitVarBins, false, 0, true, "negChildV0Type1"); - trackHistoV0Type2.init(&qaRegistry, confV0TempFitVarpTBins, confV0TempFitVarBins, confIsMC, confV0PDGCodePartTwo, true, "V0Type2"); - posChildV0Type2.init(&qaRegistry, confChildTempFitVarpTBins, confChildTempFitVarBins, false, 0, true, "posChildV0Type2"); - negChildV0Type2.init(&qaRegistry, confChildTempFitVarpTBins, confChildTempFitVarBins, false, 0, true, "negChildV0Type2"); + trackHistoV0Type1.init(&qaRegistry, V0configs.confV0TempFitVarpTBins, V0configs.confV0TempFitVarBins, confIsMC, V0configs.confV0PDGCodePartTwo, true, "V0Type1"); + posChildV0Type1.init(&qaRegistry, V0configs.confChildTempFitVarpTBins, V0configs.confChildTempFitVarBins, false, 0, true, "posChildV0Type1"); + negChildV0Type1.init(&qaRegistry, V0configs.confChildTempFitVarpTBins, V0configs.confChildTempFitVarBins, false, 0, true, "negChildV0Type1"); + trackHistoV0Type2.init(&qaRegistry, V0configs.confV0TempFitVarpTBins, V0configs.confV0TempFitVarBins, confIsMC, V0configs.confV0PDGCodePartTwo, true, "V0Type2"); + posChildV0Type2.init(&qaRegistry, V0configs.confChildTempFitVarpTBins, V0configs.confChildTempFitVarBins, false, 0, true, "posChildV0Type2"); + negChildV0Type2.init(&qaRegistry, V0configs.confChildTempFitVarpTBins, V0configs.confChildTempFitVarBins, false, 0, true, "negChildV0Type2"); // Helicity angle thetaRegistry.add("Theta/hTheta", " ; p (GeV/#it{c}); cos(#theta)", kTH2F, {{100, 0, 10}, {110, -1.1, 1.1}}); @@ -310,9 +317,9 @@ struct FemtoUniversePairTaskTrackV0Helicity { registryMCreco.add("ThetaMCReco/NegativeChild/hThetaPhi", " ; #phi; cos(#theta)", kTH2F, {{100, -1, 7}, {110, -1.1, 1.1}}); sameEventCont.init(&resultRegistry, confkstarBins, confMultBins, confkTBins, confmTBins, confMultBins3D, confmTBins3D, confEtaBins, confPhiBins, confIsMC, confUse3D); - sameEventCont.setPDGCodes(confTrkPDGCodePartOne, confV0PDGCodePartTwo); + sameEventCont.setPDGCodes(confTrkPDGCodePartOne, V0configs.confV0PDGCodePartTwo); mixedEventCont.init(&resultRegistry, confkstarBins, confMultBins, confkTBins, confmTBins, confMultBins3D, confmTBins3D, confEtaBins, confPhiBins, confIsMC, confUse3D); - mixedEventCont.setPDGCodes(confTrkPDGCodePartOne, confV0PDGCodePartTwo); + mixedEventCont.setPDGCodes(confTrkPDGCodePartOne, V0configs.confV0PDGCodePartTwo); pairCleaner.init(&qaRegistry); pairCleanerV0.init(&qaRegistry); @@ -327,11 +334,11 @@ struct FemtoUniversePairTaskTrackV0Helicity { LOGF(fatal, "Could not load efficiency histogram from %s", confLocalEfficiency.value.c_str()); if (doprocessSameEvent || doprocessMixedEvent) { plocalEffp1 = (confChargePart1 > 0) ? std::unique_ptr(plocalEffFile.get()->Get("PrPlus")) : std::unique_ptr(plocalEffFile.get()->Get("PrMinus")); // note: works only for protons for now - plocalEffp2 = (confV0Type1 == 0) ? std::unique_ptr(plocalEffFile.get()->Get("Lambda")) : std::unique_ptr(plocalEffFile.get()->Get("AntiLambda")); + plocalEffp2 = (V0configs.confV0Type1 == 0) ? std::unique_ptr(plocalEffFile.get()->Get("Lambda")) : std::unique_ptr(plocalEffFile.get()->Get("AntiLambda")); LOGF(info, "Loaded efficiency histograms for track-V0."); } else if (doprocessSameEventV0 || doprocessMixedEventV0) { - plocalEffp1 = (confV0Type1 == 0) ? std::unique_ptr(plocalEffFile.get()->Get("Lambda")) : std::unique_ptr(plocalEffFile.get()->Get("AntiLambda")); - plocalEffp2 = (confV0Type2 == 0) ? std::unique_ptr(plocalEffFile.get()->Get("Lambda")) : std::unique_ptr(plocalEffFile.get()->Get("AntiLambda")); + plocalEffp1 = (V0configs.confV0Type1 == 0) ? std::unique_ptr(plocalEffFile.get()->Get("Lambda")) : std::unique_ptr(plocalEffFile.get()->Get("AntiLambda")); + plocalEffp2 = (V0configs.confV0Type2 == 0) ? std::unique_ptr(plocalEffFile.get()->Get("Lambda")) : std::unique_ptr(plocalEffFile.get()->Get("AntiLambda")); LOGF(info, "Loaded efficiency histograms for V0-V0."); } } @@ -353,11 +360,11 @@ struct FemtoUniversePairTaskTrackV0Helicity { const auto& posChild = parts.iteratorAt(part.index() - 2); const auto& negChild = parts.iteratorAt(part.index() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild, V0ChildTable[confV0Type1][0]) || !isParticleTPC(negChild, V0ChildTable[confV0Type1][1])) + if (!isParticleTPC(posChild, V0ChildTable[V0configs.confV0Type1][0]) || !isParticleTPC(negChild, V0ChildTable[V0configs.confV0Type1][1])) continue; - auto posChildMass = pdg->Mass(confPDGCodePosChild); - auto negChildMass = pdg->Mass(confPDGCodeNegChild); + auto posChildMass = pdg->Mass(V0configs.confPDGCodePosChild); + auto negChildMass = pdg->Mass(V0configs.confPDGCodeNegChild); auto posChildBoosted = FemtoUniverseMath::boostPRF(posChild, posChildMass, negChild, negChildMass); auto cosineTheta = (posChildBoosted.Px() * part.px() + posChildBoosted.Py() * part.py() + posChildBoosted.Pz() * part.pz()) / (posChildBoosted.P() * part.p()); @@ -415,7 +422,7 @@ struct FemtoUniversePairTaskTrackV0Helicity { const auto& negChild = parts.iteratorAt(p2.index() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild, V0ChildTable[confV0Type1][0]) || !isParticleTPC(negChild, V0ChildTable[confV0Type1][1])) + if (!isParticleTPC(posChild, V0ChildTable[V0configs.confV0Type1][0]) || !isParticleTPC(negChild, V0ChildTable[V0configs.confV0Type1][1])) continue; float weight = 1.0f; @@ -462,13 +469,13 @@ struct FemtoUniversePairTaskTrackV0Helicity { const auto& negChild = parts.iteratorAt(part.index() - 1); /// Check daughters of first V0 particle - if (isParticleTPC(posChild, V0ChildTable[confV0Type1][0]) && isParticleTPC(negChild, V0ChildTable[confV0Type1][1])) { + if (isParticleTPC(posChild, V0ChildTable[V0configs.confV0Type1][0]) && isParticleTPC(negChild, V0ChildTable[V0configs.confV0Type1][1])) { trackHistoV0Type1.fillQABase(part, HIST("V0Type1")); posChildV0Type1.fillQABase(posChild, HIST("posChildV0Type1")); negChildV0Type1.fillQABase(negChild, HIST("negChildV0Type1")); } /// Check daughters of second V0 particle - if (isParticleTPC(posChild, V0ChildTable[confV0Type2][0]) && isParticleTPC(negChild, V0ChildTable[confV0Type2][1])) { + if (isParticleTPC(posChild, V0ChildTable[V0configs.confV0Type2][0]) && isParticleTPC(negChild, V0ChildTable[V0configs.confV0Type2][1])) { trackHistoV0Type2.fillQABase(part, HIST("V0Type2")); posChildV0Type2.fillQABase(posChild, HIST("posChildV0Type2")); negChildV0Type2.fillQABase(negChild, HIST("negChildV0Type2")); @@ -494,18 +501,18 @@ struct FemtoUniversePairTaskTrackV0Helicity { const auto& posChild1 = parts.iteratorAt(p1.index() - 2); const auto& negChild1 = parts.iteratorAt(p1.index() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild1, V0ChildTable[confV0Type1][0]) || !isParticleTPC(negChild1, V0ChildTable[confV0Type1][1])) + if (!isParticleTPC(posChild1, V0ChildTable[V0configs.confV0Type1][0]) || !isParticleTPC(negChild1, V0ChildTable[V0configs.confV0Type1][1])) return; const auto& posChild2 = parts.iteratorAt(p2.index() - 2); const auto& negChild2 = parts.iteratorAt(p2.index() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild2, V0ChildTable[confV0Type2][0]) || !isParticleTPC(negChild2, V0ChildTable[confV0Type2][1])) + if (!isParticleTPC(posChild2, V0ChildTable[V0configs.confV0Type2][0]) || !isParticleTPC(negChild2, V0ChildTable[V0configs.confV0Type2][1])) return; sameEventCont.setPair(p1, p2, multCol, confUse3D); }; - if (confV0Type1 == confV0Type2) { + if (V0configs.confV0Type1 == V0configs.confV0Type2) { /// Now build the combinations for identical V0s for (const auto& [p1, p2] : combinations(CombinationsStrictlyUpperIndexPolicy(groupPartsTwo, groupPartsTwo))) { pairProcessFunc(p1, p2); @@ -534,7 +541,7 @@ struct FemtoUniversePairTaskTrackV0Helicity { /// Histogramming same event for (const auto& part : groupPartsTwo) { int pdgCode = static_cast(part.pidCut()); - if ((confV0Type1 == 0 && pdgCode != confPDGCodeV0) || (confV0Type1 == 1 && pdgCode != -confPDGCodeV0)) + if ((V0configs.confV0Type1 == 0 && pdgCode != V0configs.confPDGCodeV0) || (V0configs.confV0Type1 == 1 && pdgCode != -V0configs.confPDGCodeV0)) continue; trackHistoPartTwo.fillQA(part); } @@ -560,7 +567,7 @@ struct FemtoUniversePairTaskTrackV0Helicity { if (static_cast(p1.pidCut()) != confTrkPDGCodePartOne) continue; int pdgCode2 = static_cast(p2.pidCut()); - if ((confV0Type1 == 0 && pdgCode2 != confPDGCodeV0) || (confV0Type1 == 1 && pdgCode2 != -confPDGCodeV0)) + if ((V0configs.confV0Type1 == 0 && pdgCode2 != V0configs.confPDGCodeV0) || (V0configs.confV0Type1 == 1 && pdgCode2 != -V0configs.confPDGCodeV0)) continue; // track cleaning if (confIsCPR.value) { @@ -585,22 +592,22 @@ struct FemtoUniversePairTaskTrackV0Helicity { /// Histogramming same event for (const auto& part : groupPartsTwo) { int pdgCode = static_cast(part.pidCut()); - if ((confV0Type1 == 0 && pdgCode != confPDGCodeV0) || (confV0Type1 == 1 && pdgCode != -confPDGCodeV0)) + if ((V0configs.confV0Type1 == 0 && pdgCode != V0configs.confPDGCodeV0) || (V0configs.confV0Type1 == 1 && pdgCode != -V0configs.confPDGCodeV0)) continue; trackHistoPartTwo.fillQA(part); } auto pairProcessFunc = [&](auto& p1, auto& p2) -> void { int pdgCode1 = static_cast(p1.pidCut()); - if ((confV0Type1 == 0 && pdgCode1 != confPDGCodeV0) || (confV0Type1 == 1 && pdgCode1 != -confPDGCodeV0)) + if ((V0configs.confV0Type1 == 0 && pdgCode1 != V0configs.confPDGCodeV0) || (V0configs.confV0Type1 == 1 && pdgCode1 != -V0configs.confPDGCodeV0)) return; int pdgCode2 = static_cast(p2.pidCut()); - if ((confV0Type2 == 0 && pdgCode2 != confPDGCodeV0) || (confV0Type2 == 1 && pdgCode2 != -confPDGCodeV0)) + if ((V0configs.confV0Type2 == 0 && pdgCode2 != V0configs.confPDGCodeV0) || (V0configs.confV0Type2 == 1 && pdgCode2 != -V0configs.confPDGCodeV0)) return; sameEventCont.setPair(p1, p2, multCol, confUse3D); }; /// Now build the combinations - if (confV0Type1 == confV0Type2) { + if (V0configs.confV0Type1 == V0configs.confV0Type2) { /// Now build the combinations for identical V0s for (const auto& [p1, p2] : combinations(CombinationsStrictlyUpperIndexPolicy(groupPartsTwo, groupPartsTwo))) { pairProcessFunc(p1, p2); @@ -646,7 +653,7 @@ struct FemtoUniversePairTaskTrackV0Helicity { const auto& posChild = parts.iteratorAt(p2.globalIndex() - 2); const auto& negChild = parts.iteratorAt(p2.globalIndex() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild, V0ChildTable[confV0Type1][0]) || !isParticleTPC(negChild, V0ChildTable[confV0Type1][1])) + if (!isParticleTPC(posChild, V0ChildTable[V0configs.confV0Type1][0]) || !isParticleTPC(negChild, V0ChildTable[V0configs.confV0Type1][1])) continue; // track cleaning @@ -726,13 +733,13 @@ struct FemtoUniversePairTaskTrackV0Helicity { const auto& posChild1 = parts.iteratorAt(p1.globalIndex() - 2); const auto& negChild1 = parts.iteratorAt(p1.globalIndex() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild1, V0ChildTable[confV0Type1][0]) || !isParticleTPC(negChild1, V0ChildTable[confV0Type1][1])) + if (!isParticleTPC(posChild1, V0ChildTable[V0configs.confV0Type1][0]) || !isParticleTPC(negChild1, V0ChildTable[V0configs.confV0Type1][1])) continue; const auto& posChild2 = parts.iteratorAt(p2.globalIndex() - 2); const auto& negChild2 = parts.iteratorAt(p2.globalIndex() - 1); /// Daughters that do not pass this condition are not selected - if (!isParticleTPC(posChild2, V0ChildTable[confV0Type2][0]) || !isParticleTPC(negChild2, V0ChildTable[confV0Type2][1])) + if (!isParticleTPC(posChild2, V0ChildTable[V0configs.confV0Type2][0]) || !isParticleTPC(negChild2, V0ChildTable[V0configs.confV0Type2][1])) continue; // track cleaning @@ -784,7 +791,7 @@ struct FemtoUniversePairTaskTrackV0Helicity { if (static_cast(p1.pidCut()) != confTrkPDGCodePartOne) continue; int pdgCode2 = static_cast(p2.pidCut()); - if ((confV0Type1 == 0 && pdgCode2 != confPDGCodeV0) || (confV0Type1 == 1 && pdgCode2 != -confPDGCodeV0)) + if ((V0configs.confV0Type1 == 0 && pdgCode2 != V0configs.confPDGCodeV0) || (V0configs.confV0Type1 == 1 && pdgCode2 != -V0configs.confPDGCodeV0)) continue; if (confIsCPR.value) { if (pairCloseRejection.isClosePair(p1, p2, parts, magFieldTesla1, femto_universe_container::EventType::mixed)) { @@ -824,10 +831,10 @@ struct FemtoUniversePairTaskTrackV0Helicity { for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) { int pdgCode1 = static_cast(p1.pidCut()); - if ((confV0Type1 == 0 && pdgCode1 != confPDGCodeV0) || (confV0Type1 == 1 && pdgCode1 != -confPDGCodeV0)) + if ((V0configs.confV0Type1 == 0 && pdgCode1 != V0configs.confPDGCodeV0) || (V0configs.confV0Type1 == 1 && pdgCode1 != -V0configs.confPDGCodeV0)) continue; int pdgCode2 = static_cast(p2.pidCut()); - if ((confV0Type2 == 0 && pdgCode2 != confPDGCodeV0) || (confV0Type2 == 1 && pdgCode2 != -confPDGCodeV0)) + if ((V0configs.confV0Type2 == 0 && pdgCode2 != V0configs.confPDGCodeV0) || (V0configs.confV0Type2 == 1 && pdgCode2 != -V0configs.confPDGCodeV0)) continue; mixedEventCont.setPair(p1, p2, multCol, confUse3D); } @@ -847,23 +854,48 @@ struct FemtoUniversePairTaskTrackV0Helicity { } PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Helicity, processMCMixedEventV0, "Enable processing mixed events for MC truth V0 - V0", false); - ///--------------------------------------------MC-------------------------------------------------/// - /// This function fills MC truth particles from derived MC table - void processMCTruth(aod::FDParticles const& parts) + /// --------------------------------------------------------------- MC --------------------------------------------------------------- /// + + /// This function fills MC Truth particles from derived MC table + void processMCTruth(o2::aod::FdCollision const& col, FemtoTruthParticles const& parts) { - for (const auto& part : parts) { + eventHisto.fillQA(col); + auto groupPartsTwo = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); + for (const auto& part : groupPartsTwo) { if (part.partType() != uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)) continue; - int pdgCode = static_cast(part.pidCut()); + int pdgCode = static_cast(part.tempFitVar()); const auto& pdgParticle = pdgMC->GetParticle(pdgCode); - if (!pdgParticle) { + if (!pdgParticle) continue; - } if (pdgCode == 3122) { registryMCtruth.fill(HIST("plus/MCtruthLambda"), part.pt(), part.eta()); + + // Helicity angle + const auto& posChild = parts.iteratorAt(part.index() - 2); + const auto& negChild = parts.iteratorAt(part.index() - 1); + + registryMCtruth.fill(HIST("PosChildMCTruth/hPt"), posChild.pt()); + registryMCtruth.fill(HIST("NegChildMCTruth/hPt"), negChild.pt()); + + auto posChildMass = pdg->Mass(V0configs.confPDGCodePosChild); + auto negChildMass = pdg->Mass(V0configs.confPDGCodeNegChild); + auto posChildBoosted = FemtoUniverseMath::boostPRF(posChild, posChildMass, negChild, negChildMass); + auto cosineTheta = (posChildBoosted.Px() * part.px() + posChildBoosted.Py() * part.py() + posChildBoosted.Pz() * part.pz()) / (posChildBoosted.P() * part.p()); + + registryMCtruth.fill(HIST("ThetaMCTruth/hTheta"), part.p(), cosineTheta); + registryMCtruth.fill(HIST("ThetaMCTruth/hTheta3D_PosChild"), part.p(), cosineTheta, posChild.pt()); + registryMCtruth.fill(HIST("ThetaMCTruth/hTheta3D_NegChild"), part.p(), cosineTheta, negChild.pt()); + registryMCtruth.fill(HIST("ThetaMCTruth/PositiveChild/hThetaPt"), posChild.pt(), cosineTheta); + registryMCtruth.fill(HIST("ThetaMCTruth/PositiveChild/hThetaEta"), posChild.eta(), cosineTheta); + registryMCtruth.fill(HIST("ThetaMCTruth/PositiveChild/hThetaPhi"), posChild.phi(), cosineTheta); + registryMCtruth.fill(HIST("ThetaMCTruth/NegativeChild/hThetaPt"), negChild.pt(), cosineTheta); + registryMCtruth.fill(HIST("ThetaMCTruth/NegativeChild/hThetaEta"), negChild.eta(), cosineTheta); + registryMCtruth.fill(HIST("ThetaMCTruth/NegativeChild/hThetaPhi"), negChild.phi(), cosineTheta); + continue; } else if (pdgCode == -3122) { registryMCtruth.fill(HIST("minus/MCtruthLambda"), part.pt(), part.eta()); @@ -893,31 +925,8 @@ struct FemtoUniversePairTaskTrackV0Helicity { registryMCtruth.fill(HIST("minus/MCtruthPr"), part.pt(), part.eta()); registryMCtruth.fill(HIST("minus/MCtruthPrPt"), part.pt()); } - - // Helicity angle - const auto& posChild = parts.iteratorAt(part.index() - 2); - const auto& negChild = parts.iteratorAt(part.index() - 1); - - registryMCtruth.fill(HIST("PosChildMCTruth/hPt"), posChild.pt()); - registryMCtruth.fill(HIST("NegChildMCTruth/hPt"), negChild.pt()); - - auto posChildMass = pdg->Mass(confPDGCodePosChild); - auto negChildMass = pdg->Mass(confPDGCodeNegChild); - auto posChildBoosted = FemtoUniverseMath::boostPRF(posChild, posChildMass, negChild, negChildMass); - auto cosineTheta = (posChildBoosted.Px() * part.px() + posChildBoosted.Py() * part.py() + posChildBoosted.Pz() * part.pz()) / (posChildBoosted.P() * part.p()); - - registryMCtruth.fill(HIST("ThetaMCTruth/hTheta"), part.p(), cosineTheta); - registryMCtruth.fill(HIST("ThetaMCTruth/hTheta3D_PosChild"), part.p(), cosineTheta, posChild.p()); - registryMCtruth.fill(HIST("ThetaMCTruth/hTheta3D_NegChild"), part.p(), cosineTheta, negChild.p()); - registryMCtruth.fill(HIST("ThetaMCTruth/PositiveChild/hThetaPt"), posChild.pt(), cosineTheta); - registryMCtruth.fill(HIST("ThetaMCTruth/PositiveChild/hThetaEta"), posChild.eta(), cosineTheta); - registryMCtruth.fill(HIST("ThetaMCTruth/PositiveChild/hThetaPhi"), posChild.phi(), cosineTheta); - registryMCtruth.fill(HIST("ThetaMCTruth/NegativeChild/hThetaPt"), negChild.pt(), cosineTheta); - registryMCtruth.fill(HIST("ThetaMCTruth/NegativeChild/hThetaEta"), negChild.eta(), cosineTheta); - registryMCtruth.fill(HIST("ThetaMCTruth/NegativeChild/hThetaPhi"), negChild.phi(), cosineTheta); } } - PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Helicity, processMCTruth, "Process MC truth data", false); void processMCReco(FemtoRecoParticles const& parts, aod::FdMCParticles const& mcparts) @@ -985,7 +994,6 @@ struct FemtoUniversePairTaskTrackV0Helicity { } // partType } } - PROCESS_SWITCH(FemtoUniversePairTaskTrackV0Helicity, processMCReco, "Process MC reco data", false); };