From f69a3526db6810b5cca33674971dc72758977072 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Mon, 28 Sep 2026 23:55:48 +0200 Subject: [PATCH 01/13] Fix builder --- .../Strangeness/sigma0builder.cxx | 116 +++++++++--------- 1 file changed, 57 insertions(+), 59 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index 2808e65d893..bdac95866ef 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -21,17 +21,16 @@ #include "PWGLF/DataModel/LFStrangenessMLTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "Common/Tools/StandardCCDBLoader.h" +#include #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" -#include "Common/Tools/StandardCCDBLoader.h" #include #include #include -#include -#include #include #include #include @@ -43,7 +42,6 @@ #include #include #include -#include #include // IWYU pragma: keep (do not replace with Math/Vector3Dfwd.h) #include @@ -115,7 +113,7 @@ struct sigma0builder { //__________________________________________________ // Xi0 specific Produces xi0cores; // xi0 candidates info for analysis - Produces xi0collRefs; // references to straCollision + Produces xi0collRefs; // references to straCollision Produces xi0Indices; // references to V0Cores Produces xi0mccores; // Reco xi0 MC properties Produces xi0mccollRefs; // references to straMCCollision @@ -894,12 +892,12 @@ struct sigma0builder { // ______________________________________________________ // Struct to store V0Pair properties struct V0PairTopoInfo { - std::array position{-999.f, -999.f, -999.f}; + std::array position {-999.f, -999.f, -999.f}; float dcaPi0ToPV = 999.f; float daughterDCA = 999.f; float cosPA = -1.f; - float pi0Mass = 999.f; - float pi0Y = 999.f; + float pi0Mass= 999.f; + float pi0Y= 999.f; int v01Index = 0; // index to de-reference V0Cores table int v02Index = 0; // index to de-reference V0Cores table @@ -950,8 +948,6 @@ struct sigma0builder { bool IsSigma0 = false; bool IsAntiSigma0 = false; bool IsKStar = false; - bool IsLambdaStar = false; - bool IsAntiLambdaStar = false; bool IsProducedByGenerator = false; bool IsSterile = false; int MCProcess = -1; @@ -968,31 +964,31 @@ struct sigma0builder { // ______________________________________________________ // Struct to store V0Pair properties struct Xi0Info { - std::array gamma1Momentum{999.f, 999.f, 999.f}; - std::array gamma2Momentum{999.f, 999.f, 999.f}; - std::array lambdaMomentum{999.f, 999.f, 999.f}; - std::array pi0Position{0., 0., 0.}; - std::array cascadePosition{0., 0., 0.}; + std::array gamma1Momentum {999.f, 999.f, 999.f}; + std::array gamma2Momentum {999.f, 999.f, 999.f}; + std::array lambdaMomentum {999.f, 999.f, 999.f}; + std::array pi0Position {0., 0., 0.}; + std::array cascadePosition {0., 0., 0.}; float xi0Mass = 999.f; float dcaPi0ToPV = 999.f; float pi0DaughterDCA = 999.f; float cascadeDaughterDCA = 999.f; float cascadeDCAxy = 999.f; float cascadeDCAz = 999.f; - + int collisionId = 0; int gamma1Index = 0; // index to de-reference V0Cores table int gamma2Index = 0; // index to de-reference V0Cores table int lambdaIndex = 0; // index to de-reference V0Cores table - std::array covariance{0.}; + std::array covariance {0.}; }; struct Xi0MCInfo { int gamma1Index = 0; // index to de-reference V0MCCores table int gamma2Index = 0; // index to de-reference V0MCCores table int lambdaIndex = 0; // index to de-reference V0MCCores table - + // Basic kinematic info float gamma1MCpx = -999.f; float gamma1MCpy = -999.f; @@ -1098,9 +1094,8 @@ struct sigma0builder { auto v02MC = v02.template v0MCCore_as>(); // Sanity check: Is V0Pair <-> Mother assignment correct? - int expectedPairPDG = doLambdaStar ? 3124 : PDG_t::kSigma0; bool fIsSigma0 = false; - if ((v01MC.pdgCode() == PDG_t::kGamma) && (v01MC.pdgCodeMother() == expectedPairPDG) && (v02MC.pdgCode() == PDG_t::kLambda0) && (v02MC.pdgCodeMother() == expectedPairPDG) && (v01.motherMCPartId() == v02.motherMCPartId())) + if ((v01MC.pdgCode() == PDG_t::kGamma) && (v01MC.pdgCodeMother() == PDG_t::kSigma0) && (v02MC.pdgCode() == PDG_t::kLambda0) && (v02MC.pdgCodeMother() == PDG_t::kSigma0) && (v01.motherMCPartId() == v02.motherMCPartId())) fIsSigma0 = true; bool fIsKStar = false; @@ -1251,11 +1246,11 @@ struct sigma0builder { } } // Check association correctness - if (fIsSigma0 && (MCinfo.V0PairPDGCode == expectedPairPDG)) + if (fIsSigma0 && (MCinfo.V0PairPDGCode == PDG_t::kSigma0)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 1); // match - if (fIsSigma0 && !(MCinfo.V0PairPDGCode == expectedPairPDG)) + if (fIsSigma0 && !(MCinfo.V0PairPDGCode == PDG_t::kSigma0)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 2); // mismatch - if (!fIsSigma0 && (MCinfo.V0PairPDGCode == expectedPairPDG)) + if (!fIsSigma0 && (MCinfo.V0PairPDGCode == PDG_t::kSigma0)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 3); // mismatch // Check association correctness @@ -1431,8 +1426,8 @@ struct sigma0builder { return MCinfo; } - template - Xi0MCInfo getXi0MCInfo(TV0 const& gamma1, TV0 const& gamma2, TV0 const& lambda, TMCParticles const& mcparticles) + template + Xi0MCInfo getXi0MCInfo(TV0 const& gamma1, TV0 const& gamma2, TV0 const& lambda, TCollision const& collision, TMCParticles const& mcparticles) { Xi0MCInfo MCinfo; @@ -1477,8 +1472,8 @@ struct sigma0builder { auto MCParticle_gamma2 = mcparticles.rawIteratorAt(gamma2MC.particleIdMC()); auto MCParticle_lambda = mcparticles.rawIteratorAt(lambdaMC.particleIdMC()); - if ((std::abs(MCParticle_gamma1.pdgCode()) == PDG_t::kGamma) && - (std::abs(MCParticle_gamma2.pdgCode()) == PDG_t::kGamma) && + if ((std::abs(MCParticle_gamma1.pdgCode()) == PDG_t::kGamma) && + (std::abs(MCParticle_gamma2.pdgCode()) == PDG_t::kGamma) && (std::abs(MCParticle_lambda.pdgCode()) == PDG_t::kLambda0)) { // Get MC Mothers @@ -1487,9 +1482,9 @@ struct sigma0builder { auto const& MCMothersList_lambda = MCParticle_lambda.template mothers_as(); if (!MCMothersList_gamma1.empty() && !MCMothersList_gamma2.empty() && !MCMothersList_lambda.empty()) { // Are there mothers? - auto const& MCMother_gamma1 = MCMothersList_gamma1.front(); // First mother - auto const& MCMother_gamma2 = MCMothersList_gamma2.front(); // First mother - auto const& MCMother_lambda = MCMothersList_lambda.front(); // First mother + auto const& MCMother_gamma1 = MCMothersList_gamma1.front(); // First mother + auto const& MCMother_gamma2 = MCMothersList_gamma2.front(); // First mother + auto const& MCMother_lambda = MCMothersList_lambda.front(); // First mother MCinfo.gamma1PDGCodeMother = MCMother_gamma1.pdgCode(); MCinfo.gamma2PDGCodeMother = MCMother_gamma2.pdgCode(); @@ -1505,13 +1500,13 @@ struct sigma0builder { MCinfo.pi0Z = MCMother_gamma1.vz(); auto const& MCMothersList_pi0 = MCMother_gamma1.template mothers_as(); // get pi0 mother list - if (!MCMothersList_pi0.empty()) { // Are there mothers? - auto const& MCMother_pi0 = MCMothersList_pi0.front(); // get pi0 mother - + if (!MCMothersList_pi0.empty()) { // Are there mothers? + auto const& MCMother_pi0 = MCMothersList_pi0.front(); // get pi0 mother + if (MCMother_pi0.globalIndex() == MCMother_lambda.globalIndex()) { // check that lambda and pi0 have the same mother. MCinfo.xi0PDGCode = MCMother_pi0.pdgCode(); MCinfo.xi0MCProcess = MCMother_pi0.getProcess(); - MCinfo.pi0IsPhysicalPrimary = MCMother_pi0.isPhysicalPrimary(); + MCinfo.xi0IsPhysicalPrimary = MCMother_pi0.isPhysicalPrimary(); MCinfo.xi0X = MCMother_pi0.vx(); MCinfo.xi0Y = MCMother_pi0.vy(); MCinfo.xi0Z = MCMother_pi0.vz(); @@ -1910,8 +1905,6 @@ struct sigma0builder { GenInfo.IsPi0 = mcParticle.pdgCode() == PDG_t::kPi0; // 111; GenInfo.IsSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0; // PDG_t::kSigma0 GenInfo.IsAntiSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0Bar; //-3212 - GenInfo.IsLambdaStar = mcParticle.pdgCode() == 3124; // 102134 (PYTHIA8) - GenInfo.IsAntiLambdaStar = mcParticle.pdgCode() == -3124; // -102134 GenInfo.IsKStar = std::abs(mcParticle.pdgCode()) == o2::constants::physics::Pdg::kK0Star892; // 313; GenInfo.IsProducedByGenerator = mcParticle.producedByGenerator(); GenInfo.MCProcess = mcParticle.getProcess(); @@ -1923,7 +1916,7 @@ struct sigma0builder { GenInfo.MCCollId = mcParticle.mcCollisionId(); // save this reference, please // Checking decay mode if sigma0 or pi0 (it is easier here) - if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0 || GenInfo.IsPi0 || GenInfo.IsKStar || GenInfo.IsLambdaStar || GenInfo.IsAntiLambdaStar) { + if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0 || GenInfo.IsPi0 || GenInfo.IsKStar) { // This is a costly operation, so we do it only for pi0s and sigma0s auto const& daughters = mcParticle.template daughters_as(); @@ -1933,7 +1926,7 @@ struct sigma0builder { auto const& GenMothersList = mcParticle.template mothers_as(); GenInfo.PDGCodeMother = (!GenMothersList.empty()) ? GenMothersList.front().pdgCode() : 0; - if ((doLambdaStar ? (GenInfo.IsLambdaStar || GenInfo.IsAntiLambdaStar) : (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0)) && genSelections.doQA) { + if ((GenInfo.IsSigma0 || GenInfo.IsAntiSigma0) && genSelections.doQA) { histos.fill(HIST("GenQA/h2dSigma0MCSourceVsPDGMother"), GenInfo.IsProducedByGenerator, GenInfo.PDGCodeMother); // Checking decay modes and getting daughter pTs @@ -1996,7 +1989,7 @@ struct sigma0builder { histos.fill(HIST("GenQA/hGenSpeciesKStar"), 0); // Checking decay mode - if (doLambdaStar ? (GenInfo.IsLambdaStar || GenInfo.IsAntiLambdaStar) : (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0)) { + if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0) { histos.fill(HIST("GenQA/hSigma0NDau"), GenInfo.NDaughters); histos.fill(HIST("GenQA/h2dSigma0NDauVsProcess"), GenInfo.NDaughters, GenInfo.MCProcess); @@ -2016,7 +2009,7 @@ struct sigma0builder { histos.fill(HIST("GenQA/h2DGenSigma0TypeVsProducedByGen"), typeIndex, genIndex); // Fill histograms - if (doLambdaStar ? GenInfo.IsLambdaStar : GenInfo.IsSigma0) { + if (GenInfo.IsSigma0) { histos.fill(HIST("GenQA/hGenSpecies"), 2); histos.fill(HIST("GenQA/hGenSigma0"), GenInfo.MCPt); histos.fill(HIST("GenQA/h3dGenSigma0_pTMap"), GenInfo.MCPt, GenInfo.MCDau1Pt, GenInfo.MCDau2Pt); @@ -2025,7 +2018,7 @@ struct sigma0builder { if (GenInfo.IsPrimary) histos.fill(HIST("GenQA/hPrimarySigma0s"), 1); } - if (doLambdaStar ? GenInfo.IsAntiLambdaStar : GenInfo.IsAntiSigma0) { + if (GenInfo.IsAntiSigma0) { histos.fill(HIST("GenQA/hGenSpecies"), 3); histos.fill(HIST("GenQA/hGenAntiSigma0"), GenInfo.MCPt); histos.fill(HIST("GenQA/h3dGenASigma0_pTMap"), GenInfo.MCPt, GenInfo.MCDau1Pt, GenInfo.MCDau2Pt); @@ -2102,10 +2095,9 @@ struct sigma0builder { pi0GenCollRefs(MCGenInfo.MCCollId); // link to stramccollision table } - // Sigma0/ASigma0 (Lambda(1520)/ALambda(1520)) - bool fIsGenSigma0Like = doLambdaStar ? (MCGenInfo.IsLambdaStar || MCGenInfo.IsAntiLambdaStar) : (MCGenInfo.IsSigma0 || MCGenInfo.IsAntiSigma0); - if (fillSigma0Tables && fIsGenSigma0Like) { - sigma0Gens(doLambdaStar ? MCGenInfo.IsLambdaStar : MCGenInfo.IsSigma0, MCGenInfo.IsProducedByGenerator, MCGenInfo.MCPt, mcParticle.y()); + // Sigma0/ASigma0 + if (fillSigma0Tables && (MCGenInfo.IsSigma0 || MCGenInfo.IsAntiSigma0)) { + sigma0Gens(MCGenInfo.IsSigma0, MCGenInfo.IsProducedByGenerator, MCGenInfo.MCPt, mcParticle.y()); sigma0GenCollRefs(MCGenInfo.MCCollId); // link to stramccollision table } @@ -2609,7 +2601,7 @@ struct sigma0builder { //_______________________________________________ // Build pi0 candidate for QA template - bool buildPi0ForXi0(TV0Object const& gamma1, TV0Object const& gamma2, TCollision const& collision, V0PairTopoInfo& info) + bool buildPi0ForXi0(TV0Object const& gamma1, TV0Object const& gamma2, TCollision const& collision, V0PairTopoInfo &info) { //_______________________________________________ // Check if both V0s are made of the same tracks @@ -2637,7 +2629,7 @@ struct sigma0builder { info.v01Index = gamma1.globalIndex(); info.v02Index = gamma2.globalIndex(); - + //_______________________________________________ // Pi0-specific selections: if (std::abs(info.pi0Y) > cascadeSelections.maxPi0Rapidity) @@ -2683,7 +2675,7 @@ struct sigma0builder { if (std::hypot(info.position[0], info.position[1]) < cascadeSelections.radiusPi0) { return false; } - + float x = vtxPi0[0]; float y = vtxPi0[1]; float z = vtxPi0[2]; @@ -2703,7 +2695,7 @@ struct sigma0builder { //_______________________________________________ // Build Xi0 candidate for analysis template - bool buildXi0(TV0Object const& gamma1, TV0Object const& gamma2, TV0Object const& lambda, TCollision const& collision, TMCParticles const& mcparticles) + bool buildXi0(TV0Object const& gamma1, TV0Object const& gamma2, TV0Object const& lambda, V0PairTopoInfo &pi0Info, TCollision const& collision, TMCParticles const& mcparticles) { Xi0Info cascade; //_______________________________________________ @@ -2764,6 +2756,8 @@ struct sigma0builder { return false; } + cascade.pi0DaughterDCA = pi0Info.daughterDCA; + lGamma1Track.getPxPyPzGlo(cascade.gamma1Momentum); lGamma2Track.getPxPyPzGlo(cascade.gamma2Momentum); lLambdaTrack.getPxPyPzGlo(cascade.lambdaMomentum); @@ -2776,6 +2770,10 @@ struct sigma0builder { return false; } + cascade.pi0Position[0] = pi0Info.position[0]; + cascade.pi0Position[1] = pi0Info.position[1]; + cascade.pi0Position[2] = pi0Info.position[2]; + double cosPA = RecoDecay::cpa( std::array{collision.posX(), collision.posY(), collision.posZ()}, std::array{cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2]}, @@ -2788,9 +2786,9 @@ struct sigma0builder { // Calculate DCAxy of the cascade (with bending) auto lCascadeTrack = fitter3Prongs.createParentTrackParCov(); - lCascadeTrack.setAbsCharge(0); // to be sure + lCascadeTrack.setAbsCharge(0); // to be sure lCascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas - std::array dcaInfo{999.f, 999.f}; + std::array dcaInfo {999.f, 999.f}; o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, fitter3Prongs.getMatCorrType(), &dcaInfo); cascade.cascadeDCAxy = dcaInfo[0]; @@ -2803,10 +2801,10 @@ struct sigma0builder { cascade.gamma2Index = gamma2.globalIndex(); cascade.lambdaIndex = lambda.globalIndex(); - std::array momPi0{gamma1.px() + gamma2.px(), gamma1.py() + gamma2.py(), gamma1.pz() + gamma2.pz()}; + std::array momPi0{gamma1.px()+gamma2.px(), gamma1.py()+gamma2.py(), gamma1.pz()+gamma2.pz()}; auto arrMom = std::array{momPi0, momLambda}; - cascade.xi0Mass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassLambda0}); - + cascade.xi0Mass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassXi0}); + // Calculate position covariance matrix auto covVtxV = fitter3Prongs.calcPCACovMatrix(0); // std::array positionCovariance; @@ -2835,8 +2833,8 @@ struct sigma0builder { } // Check if MC data and populate corresponding table - if constexpr (requires { collision.straMCCollisionId(); gamma1.motherMCPartId(); gamma2.motherMCPartId(); lambda.motherMCPartId(); }) { - auto xi0MCInfo = getXi0MCInfo(gamma1, gamma2, lambda, mcparticles); + if constexpr (requires { collision.straMCCollisionId(); gamma1.motherMCPartId(); gamma2.motherMCPartId(); lambda.motherMCPartId();}) { + auto xi0MCInfo = getXi0MCInfo(gamma1, gamma2, lambda, collision, mcparticles); xi0mccores(xi0MCInfo.gamma1MCpx, xi0MCInfo.gamma1MCpy, xi0MCInfo.gamma1MCpz, xi0MCInfo.gamma1IsPhysicalPrimary, xi0MCInfo.gamma1PDGCodePos, xi0MCInfo.gamma1PDGCodeNeg, xi0MCInfo.gamma1PDGCode, xi0MCInfo.gamma1PDGCodeMother, @@ -2846,7 +2844,7 @@ struct sigma0builder { xi0MCInfo.lambdaIsPhysicalPrimary, xi0MCInfo.lambdaPDGCodePos, xi0MCInfo.lambdaPDGCodeNeg, xi0MCInfo.lambdaPDGCode, xi0MCInfo.lambdaPDGCodeMother, xi0MCInfo.pi0X, xi0MCInfo.pi0Y, xi0MCInfo.pi0Z, xi0MCInfo.pi0IsPhysicalPrimary, xi0MCInfo.pi0MCProcess, xi0MCInfo.pi0PDGCode, xi0MCInfo.xi0X, xi0MCInfo.xi0Y, xi0MCInfo.xi0Z, xi0MCInfo.xi0IsPhysicalPrimary, xi0MCInfo.xi0MCProcess, xi0MCInfo.xi0PDGCode, xi0MCInfo.xi0PDGCodeMother); - + xi0mcIndices(xi0MCInfo.gamma1Index, xi0MCInfo.gamma2Index, xi0MCInfo.lambdaIndex); int mcCollisionIndex = -1; @@ -3410,7 +3408,7 @@ struct sigma0builder { pi0Candidates.push_back(pi0Info); } } - + // second build Xi0 candidates for (size_t i = 0; i < bestLambdasArray.size(); ++i) { // loop over Lambda auto lambda = fullV0s.rawIteratorAt(bestLambdasArray[i]); @@ -3418,7 +3416,7 @@ struct sigma0builder { auto gamma1 = fullV0s.rawIteratorAt(pi0Candidates[j].v01Index); auto gamma2 = fullV0s.rawIteratorAt(pi0Candidates[j].v02Index); // Building pi0 candidate & filling tables - if (!buildXi0(gamma1, gamma2, lambda, coll, mcparticles)) + if (!buildXi0(gamma1, gamma2, lambda, pi0Candidates[j], coll, mcparticles)) continue; } } From 66e561fce2a2f0c20cd7ec9fbb1266f475316b0e Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Mon, 28 Sep 2026 23:56:28 +0200 Subject: [PATCH 02/13] Fix data model for Xi0 --- PWGLF/DataModel/LFSigmaTables.h | 93 ++++++++++++--------------------- 1 file changed, 32 insertions(+), 61 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index de5d0f46acd..cfb9c5668e0 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -107,22 +107,6 @@ DECLARE_SOA_DYNAMIC_COLUMN(OPAngle, opAngle, return v1.Angle(v2); }); -// Armenteros-Podolanski variables (photon = positive daughter, lambda = negative daughter) -DECLARE_SOA_DYNAMIC_COLUMN(LStarAlpha, lStarAlpha, //! Armenteros Alpha - [](float photonPx, float photonPy, float photonPz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { - float momTot = RecoDecay::p(photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz); - float lQlNeg = RecoDecay::dotProd(std::array{lambdaPx, lambdaPy, lambdaPz}, std::array{photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz}) / momTot; - float lQlPos = RecoDecay::dotProd(std::array{photonPx, photonPy, photonPz}, std::array{photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz}) / momTot; - return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); - }); - -DECLARE_SOA_DYNAMIC_COLUMN(LStarQtArm, lStarQtarm, //! Armenteros Qt - [](float photonPx, float photonPy, float photonPz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { - float momTot = RecoDecay::p2(photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz); - float dp = RecoDecay::dotProd(std::array{lambdaPx, lambdaPy, lambdaPz}, std::array{photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz}); - return std::sqrt(RecoDecay::p2(lambdaPx, lambdaPy, lambdaPz) - dp * dp / momTot); // qtarm - }); - // Photon DECLARE_SOA_DYNAMIC_COLUMN(PhotonPt, photonPt, //! Transverse momentum in GeV/c [](float photonPx, float photonPy) -> float { @@ -195,8 +179,6 @@ DECLARE_SOA_TABLE(Sigma0Cores, "AOD", "SIGMA0CORES", sigma0Core::Eta, sigma0Core::Radius, sigma0Core::OPAngle, - sigma0Core::LStarAlpha, - sigma0Core::LStarQtArm, sigma0Core::PhotonPt, sigma0Core::PhotonP, @@ -674,12 +656,6 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsSigma0, isSigma0, DECLARE_SOA_DYNAMIC_COLUMN(IsAntiSigma0, isAntiSigma0, //! IsASigma0 [](int pdgCode) -> bool { return pdgCode == PDG_t::kSigma0Bar; }); //-3212 -DECLARE_SOA_DYNAMIC_COLUMN(IsLambdaStar, isLambdaStar, //! IsLambdaStar - [](int pdgCode) -> bool { return pdgCode == 3124; }); // PYTHIA8 code for Lambda(1520) - -DECLARE_SOA_DYNAMIC_COLUMN(IsAntiLambdaStar, isAntiLambdaStar, //! IsAntiLambdaStar - [](int pdgCode) -> bool { return pdgCode == -3124; }); // PYTHIA8 code for AntiLambda(1520) - DECLARE_SOA_DYNAMIC_COLUMN(MCPx, mcpx, //! Sigma0 px [](float photonMCPx, float lambdaMCPx) -> float { return photonMCPx + lambdaMCPx; }); DECLARE_SOA_DYNAMIC_COLUMN(MCPy, mcpy, //! Sigma0 py @@ -710,11 +686,6 @@ DECLARE_SOA_DYNAMIC_COLUMN(Sigma0MCY, sigma0MCY, return RecoDecay::y(std::array{photonMCPx + lambdaMCPx, photonMCPy + lambdaMCPy, photonMCPz + lambdaMCPz}, o2::constants::physics::MassSigma0); }); -DECLARE_SOA_DYNAMIC_COLUMN(LambdaStarMCY, lambdaStarMCY, - [](float photonMCPx, float photonMCPy, float photonMCPz, float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { - return RecoDecay::y(std::array{photonMCPx + lambdaMCPx, photonMCPy + lambdaMCPy, photonMCPz + lambdaMCPz}, o2::constants::physics::MassLambda1520); - }); - DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcphi, //! Phi in the range [0, 2pi) [](float photonMCPx, float photonMCPy, float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(photonMCPx + lambdaMCPx, photonMCPy + lambdaMCPy); }); @@ -796,8 +767,6 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", // Dynamic columns sigma0MCCore::IsSigma0, sigma0MCCore::IsAntiSigma0, - sigma0MCCore::IsLambdaStar, - sigma0MCCore::IsAntiLambdaStar, sigma0MCCore::MCPx, sigma0MCCore::MCPy, @@ -806,7 +775,6 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", sigma0MCCore::MCP, sigma0MCCore::Sigma0MCMass, sigma0MCCore::Sigma0MCY, - sigma0MCCore::LambdaStarMCY, sigma0MCCore::MCPhi, sigma0MCCore::MCEta, sigma0MCCore::MCOPAngle, @@ -1395,6 +1363,7 @@ DECLARE_SOA_TABLE(Pi0Gens, "AOD", "PI0GENS", DECLARE_SOA_TABLE(Pi0GenCollRef, "AOD", "PI0GENCOLLREF", //! optional table to refer back to a collision o2::soa::Index<>, v0data::StraMCCollisionId); + // ___________________________________________________________________________ // Xi0 namespace Xi0Core @@ -1449,11 +1418,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Total momentum in GeV/c DECLARE_SOA_DYNAMIC_COLUMN(Xi0Mass, xi0Mass, [](float photon1Px, float photon1Py, float photon1Pz, float photon2Px, float photon2Py, float photon2Pz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { - std::array pVecPi0{photon1Px + photon2Px, photon1Py + photon2Py, photon1Pz + photon2Pz}; + std::array pVecPi0{photon1Px+photon2Px, photon1Py+photon2Py, photon1Pz+photon2Pz}; std::array pVecLambda{lambdaPx, lambdaPy, lambdaPz}; auto arrMom = std::array{pVecPi0, pVecLambda}; - return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassXi0}); - }); + return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassLambda0}); + }); DECLARE_SOA_DYNAMIC_COLUMN(PxPi0, pxPi0, //! Pi0 px [](float photon1Px, float photon2Px) -> float { return photon1Px + photon2Px; }); @@ -1487,11 +1456,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(Radius, radius, //! Xi0 decay radius (2D, centered at [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(RadiusPi0, radiusPi0, //! Pi0 decay radius (2D, centered at zero) - [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(RadiusLambda, radiusLambda, //! Pi0 decay radius (2D, centered at zero) - [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); - + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + // CosPAs DECLARE_SOA_DYNAMIC_COLUMN(Pi0CosPA, pi0CosPA, //! [](float xPi0, float yPi0, float zPi0, float pxGamma1, float pyGamma1, float pzGamma1, float pxGamma2, float pyGamma2, float pzGamma2, float pvX, float pvY, float pvZ) -> float { @@ -1510,14 +1479,14 @@ DECLARE_SOA_DYNAMIC_COLUMN(DCALambdaToPV, dcaLambdaToPV, //! [](float X, float Y, float Z, float px, float py, float pz, float pvX, float pvY, float pvZ) -> float { return std::sqrt((std::pow((pvY - Y) * pz - (pvZ - Z) * py, 2) + std::pow((pvX - X) * pz - (pvZ - Z) * px, 2) + std::pow((pvX - X) * py - (pvY - Y) * px, 2)) / (px * px + py * py + pz * pz)); }); - + DECLARE_SOA_DYNAMIC_COLUMN(DCAPi0ToPV, dcaPi0ToPV, //! [](float X, float Y, float Z, float pxGamma1, float pyGamma1, float pzGamma1, float pxGamma2, float pyGamma2, float pzGamma2, float pvX, float pvY, float pvZ) -> float { float px = pxGamma1 + pxGamma2; float py = pyGamma1 + pyGamma2; float pz = pzGamma1 + pzGamma2; return std::sqrt((std::pow((pvY - Y) * pz - (pvZ - Z) * py, 2) + std::pow((pvX - X) * pz - (pvZ - Z) * px, 2) + std::pow((pvX - X) * py - (pvY - Y) * px, 2)) / (px * px + py * py + pz * pz)); - }); + }); DECLARE_SOA_DYNAMIC_COLUMN(Photon1Pt, photon1Pt, //! Transverse momentum in GeV/c [](float photon1Px, float photon1Py) -> float { @@ -1586,12 +1555,12 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaY, lambdaY, //! Rapidity }); DECLARE_SOA_DYNAMIC_COLUMN(LambdaPhi, lambdaPhi, //! Phi in the range [0, 2pi) - [](float lambdaPx, float lambdaPy) -> float { return RecoDecay::phi(lambdaPx, lambdaPy); }); + [](float lambdaPx, float lambdaPy) -> float { return RecoDecay::phi(lambdaPx, lambdaPy); }); -} // namespace Xi0Core +} // namespace Pi0Core DECLARE_SOA_TABLE(Xi0Cores, "AOD", "XI0CORES", - Xi0Core::X, Xi0Core::Y, Xi0Core::Z, + Xi0Core::X, Xi0Core::Y, Xi0Core::Z, Xi0Core::XPi0, Xi0Core::YPi0, Xi0Core::ZPi0, Xi0Core::XLambda, Xi0Core::YLambda, Xi0Core::ZLambda, Xi0Core::DCADaughters, Xi0Core::DCADaughtersPi0, Xi0Core::DCAXYCascToPV, Xi0Core::DCAZCascToPV, @@ -1607,11 +1576,12 @@ DECLARE_SOA_TABLE(Xi0Cores, "AOD", "XI0CORES", Xi0Core::Pz, Xi0Core::Pt, Xi0Core::P, + Xi0Core::Xi0Mass, Xi0Core::Rapidity, Xi0Core::Phi, Xi0Core::Eta, Xi0Core::Radius, - Xi0Core::CascCosPA, + Xi0Core::CascCosPA, // Dynamic columns for pi0 Xi0Core::PxPi0, @@ -1620,7 +1590,7 @@ DECLARE_SOA_TABLE(Xi0Cores, "AOD", "XI0CORES", Xi0Core::Pi0Mass, Xi0Core::RadiusPi0, Xi0Core::Pi0CosPA, - Xi0Core::DCAPi0ToPV, + Xi0Core::DCAPi0ToPV, // Dynamic columns for Lambda Xi0Core::RadiusLambda, @@ -1650,14 +1620,15 @@ DECLARE_SOA_TABLE(Xi0CollRefs, "AOD", "XI0COLLREFS", //! optional table to refer o2::soa::Index<>, v0data::StraCollisionId); DECLARE_SOA_TABLE(Xi0Indices, "AOD", "XI0INDICES", //! optional table to refer back to V0Cores - o2::soa::Index<>, Xi0Core::Photon1Index, Xi0Core::Photon2Index, Xi0Core::LambdaIndex); - + o2::soa::Index<>, Xi0Core::Photon1Index, Xi0Core::Photon2Index, Xi0Core::LambdaIndex); + + // for MC namespace Xi0MCCore { DECLARE_SOA_COLUMN(Photon1Index, photon1Index, int); DECLARE_SOA_COLUMN(Photon2Index, photon2Index, int); -DECLARE_SOA_COLUMN(LambdaIndex, lambdaIndex, int); +DECLARE_SOA_COLUMN(LambdaIndex, lambdaIndex, int); DECLARE_SOA_COLUMN(Photon1MCPx, photon1mcpx, float); DECLARE_SOA_COLUMN(Photon1MCPy, photon1mcpy, float); @@ -1730,13 +1701,13 @@ DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcPhi, //! Phi in the range [0, 2pi) DECLARE_SOA_DYNAMIC_COLUMN(MCEta, mcEta, //! Pseudorapidity [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz, float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { return RecoDecay::eta(std::array{photon1MCPx + photon2MCPx + lambdaMCPx, photon1MCPy + photon2MCPy + lambdaMCPy, photon1MCPz + photon2MCPz + lambdaMCPz}); - }); - + }); + DECLARE_SOA_DYNAMIC_COLUMN(MCRadius, mcRadius, //! Xi0 decay radius (2D, centered at zero) [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(MCRadiusPi0, mcRadiusPi0, //! Pi0 decay radius (2D, centered at zero) - [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPx, pi0mcpx, //! Pi0 MC px [](float photon1MCPx, float photon2MCPx) -> float { return photon1MCPx + photon2MCPx; }); @@ -1850,9 +1821,9 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCY, lambdaMCY, //! Rapidity }); DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCPhi, lambdaMCPhi, //! Phi in the range [0, 2pi) - [](float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(lambdaMCPx, lambdaMCPy); }); + [](float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(lambdaMCPx, lambdaMCPy); }); -} // namespace Xi0MCCore +} // namespace Xi0CoreMC DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", // Basic properties @@ -1869,7 +1840,7 @@ DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", Xi0MCCore::IsPi0Primary, Xi0MCCore::Pi0MCprocess, Xi0MCCore::Pi0PDGCode, Xi0MCCore::MCx, Xi0MCCore::MCy, Xi0MCCore::MCz, - Xi0MCCore::IsPhysicalPrimary, Xi0MCCore::MCprocess, Xi0MCCore::PDGCode, Xi0MCCore::PDGCodeMother, + Xi0MCCore::IsPhysicalPrimary, Xi0MCCore::MCprocess, Xi0MCCore::PDGCode, Xi0MCCore::PDGCodeMother, // Dynamic columns Xi0MCCore::MCPx, @@ -1884,11 +1855,11 @@ DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", Xi0MCCore::MCRadius, Xi0MCCore::MCRadiusPi0, - Xi0MCCore::MCPx, - Xi0MCCore::MCPy, - Xi0MCCore::MCPz, - Xi0MCCore::MCPt, - Xi0MCCore::MCP, + Xi0MCCore::Pi0MCPx, + Xi0MCCore::Pi0MCPy, + Xi0MCCore::Pi0MCPz, + Xi0MCCore::Pi0MCPt, + Xi0MCCore::Pi0MCP, Xi0MCCore::Pi0MCMass, Xi0MCCore::Pi0RapidityMC, Xi0MCCore::Pi0MCPhi, @@ -1906,7 +1877,7 @@ DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", Xi0MCCore::Photon2MCEta, Xi0MCCore::Photon2MCY, Xi0MCCore::Photon2MCPhi, - + Xi0MCCore::LambdaMCPt, Xi0MCCore::LambdaMCP, Xi0MCCore::LambdaMCEta, @@ -1917,7 +1888,7 @@ DECLARE_SOA_TABLE(Xi0MCCollRefs, "AOD", "XI0MCCOLLREFS", //! optional table to r o2::soa::Index<>, v0data::StraMCCollisionId); DECLARE_SOA_TABLE(Xi0MCIndices, "AOD", "XI0MCINDICES", //! optional table to refer back to V0MCCores - o2::soa::Index<>, Xi0MCCore::Photon1Index, Xi0MCCore::Photon2Index, Xi0MCCore::LambdaIndex); + o2::soa::Index<>, Xi0MCCore::Photon1Index, Xi0MCCore::Photon2Index, Xi0MCCore::LambdaIndex); } // namespace o2::aod From 10c47571dd0b401c48505720cf084b3c1be72a8a Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Mon, 28 Sep 2026 23:57:23 +0200 Subject: [PATCH 03/13] Add derivedXi0analysis in CMakeLists.txt --- PWGLF/Tasks/Strangeness/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index 9b08a271008..9ddda410938 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -69,6 +69,11 @@ o2physics_add_dpl_workflow(derivedcascadeanalysis PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore O2Physics::AnalysisCCDB COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(derivedxi0analysis + SOURCES derivedXi0Analysis.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(cascpostprocessing SOURCES cascpostprocessing.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore From 235be03f6581c285858101f0d37e66b9a678d47f Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 00:02:30 +0200 Subject: [PATCH 04/13] Add derivedXi0analysis task --- .../Tasks/Strangeness/derivedXi0analysis.cxx | 1274 +++++++++++++++++ 1 file changed, 1274 insertions(+) create mode 100644 PWGLF/Tasks/Strangeness/derivedXi0analysis.cxx diff --git a/PWGLF/Tasks/Strangeness/derivedXi0analysis.cxx b/PWGLF/Tasks/Strangeness/derivedXi0analysis.cxx new file mode 100644 index 00000000000..cfc54fd89fa --- /dev/null +++ b/PWGLF/Tasks/Strangeness/derivedXi0analysis.cxx @@ -0,0 +1,1274 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +/// \file derivedXi0Analysis.cxx +/// \brief Xi0 (--> Lambda pi0 --> Lambda gamma gamma) analysis task using strangeness derived data (produced with the Xi0 builder) +/// +/// \author Romain Schotter , Austrian Academy of Sciences & MBI +// +// Xi0 analysis task +// ================= +// +// This code loops over the Xi0Cores table produced by the sigma0builder +// (PWGLF/TableProducer/Strangeness/sigma0builder.cxx, fillXi0Tables = true) +// and produces some standard analysis output. It is meant to be run over +// strangeness derived data. +// +// The Xi0 candidate is built out of three V0s: two photon conversions +// (forming the pi0) and one Lambda. The daughter V0s are reached by +// de-referencing the V0Cores table through the Xi0Indices table. +// +// Three process functions are provided: +// - processRealData : real data +// - processMonteCarlo : reconstructed information in MC +// - processGenerated : pure generated information in MC (from CascMCCores) +// +// N.B.: when running over reconstructed MC information, the rapidity and the +// pT used both in the selections and in the histograms are the generated +// ones, so that the numerator and the denominator of the +// acceptance x efficiency share exactly the same definition. +// +// Comments, questions, complaints, suggestions? +// Please write to: +// romain.schotter@cern.ch +// + +#include "PWGLF/DataModel/LFSigmaTables.h" +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/CCDB/ctpRateFetcher.h" +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +// simple checkers, but ensure 64 bit integers +#define BITSET(var, nbit) ((var) |= (static_cast(1) << static_cast(nbit))) +#define BITCHECK(var, nbit) ((var) & (static_cast(1) << static_cast(nbit))) + +// Xi0 candidates: cores + collision reference + indices to the daughter V0s. +// N.B.: aod::Xi0Indices and aod::Xi0MCIndices must never be joined together, +// since they share the very same column names (photon1Index, ...). +using Xi0Candidates = soa::Join; +using Xi0McCandidates = soa::Join; + +using V0Candidates = soa::Join; +using DauTracks = soa::Join; + +enum CentEstimator { + kCentFT0C = 0, + kCentFT0M, + kCentFT0CVariant1, + kCentMFT, + kCentNGlobal, + kCentFV0A +}; + +struct DerivedXi0Analysis { + Service ccdb; + ctpRateFetcher rateFetcher; + + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // Xi0 proper decay length, in cm (PDG: tau = 2.90e-10 s) + static constexpr float CtauXi0 = 8.71f; + // number of QA histogram sets: before and after the candidate selections + static constexpr int NSelectionStages = 2; + + //__________________________________________________ + // Selection bits: one per selection criterion, so that the full selection is + // a single mask operation and the bookkeeping fits in a single histogram + enum Selectionbits : int { // o2-linter: disable=name/enum (bit names follow derivedlambdakzeroanalysis) + // photon selections (both photons of the pi0 have to pass) + selPhotonV0Type = 0, + selPhotonMass, + selPhotonRapidity, + selPhotonDauEta, + selPhotonDCADauToPV, + selPhotonDCADau, + selPhotonRadius, + selPhotonCosPA, + selPhotonArmenteros, + selPhotonPt, + selPhotonTPCCrossedRows, + selPhotonTPCPID, + // Lambda selections, common to Lambda and anti-Lambda + selLambdaV0Type, + selLambdaRapidity, + selLambdaDauEta, + selLambdaDCAPosToPV, + selLambdaDCANegToPV, + selLambdaDCADau, + selLambdaRadius, + selLambdaCosPA, + selLambdaLifetime, + selLambdaTPCCrossedRows, + selLambdaITSClusters, + // Lambda selections, specific to one of the two mass hypotheses + selLambdaMass, + selAntiLambdaMass, + selLambdaArmenteros, // Armenteros alpha > 0 --> the positive daughter is the baryon + selAntiLambdaArmenteros, // Armenteros alpha < 0 --> the negative daughter is the baryon + selTPCPIDPositiveProton, + selTPCPIDNegativePion, + selTPCPIDNegativeProton, + selTPCPIDPositivePion, + // pi0 and Xi0 (cascade) selections + selPi0Mass, + selPi0Radius, + selPi0CosPA, + selPi0DCADau, + selXi0Radius, + selXi0CosPA, + selXi0DCADau, + selXi0DCAxyToPV, + selXi0DCAzToPV, + selXi0Lifetime, + selXi0Rapidity, + selXi0Pt, + // MC tagging + selConsiderXi0, + selConsiderAntiXi0, + selPhysPrimXi0, + selPhysPrimAntiXi0, + }; + + uint64_t maskTopological = 0; + uint64_t maskXi0Specific = 0; + uint64_t maskAntiXi0Specific = 0; + uint64_t maskSelectionXi0 = 0; + uint64_t maskSelectionAntiXi0 = 0; + + //__________________________________________________ + // Event level + Configurable centralityEstimator{"centralityEstimator", kCentFT0M, "Run 3 centrality estimator (0:CentFT0C, 1:CentFT0M, 2:CentFT0CVariant1, 3:CentMFT, 4:CentNGlobal, 5:CentFV0A)"}; + Configurable fGetIR{"fGetIR", false, "Flag to retrieve the IR info."}; + Configurable fIRCrashOnNull{"fIRCrashOnNull", false, "Flag to avoid CTP RateFetcher crash."}; + Configurable irSource{"irSource", "T0VTX", "Estimator of the interaction rate (Recommended: pp --> T0VTX, Pb-Pb --> ZNC hadronic)"}; + + struct : ConfigurableGroup { + std::string prefix = "eventSelections"; // JSON group name + Configurable requireSel8{"requireSel8", true, "require sel8 event selection"}; + Configurable requireTriggerTVX{"requireTriggerTVX", true, "require FT0 vertex (acceptable FT0C-FT0A time difference) at trigger level"}; + Configurable rejectITSROFBorder{"rejectITSROFBorder", true, "reject events at ITS ROF border"}; + Configurable rejectTFBorder{"rejectTFBorder", true, "reject events at TF border"}; + Configurable requireIsVertexITSTPC{"requireIsVertexITSTPC", false, "require events with at least one ITS-TPC track"}; + Configurable requireIsGoodZvtxFT0VsPV{"requireIsGoodZvtxFT0VsPV", true, "require events with PV position along z consistent (within 1 cm) between PV reconstructed using tracks and PV using FT0 A-C time difference"}; + Configurable requireIsVertexTOFmatched{"requireIsVertexTOFmatched", false, "require events with at least one of vertex contributors matched to TOF"}; + Configurable requireIsVertexTRDmatched{"requireIsVertexTRDmatched", false, "require events with at least one of vertex contributors matched to TRD"}; + Configurable rejectSameBunchPileup{"rejectSameBunchPileup", true, "reject collisions in case of pileup with another collision in the same foundBC"}; + Configurable requireNoCollInTimeRangeStd{"requireNoCollInTimeRangeStd", false, "reject collisions corrupted by the cannibalism, with other collisions within +/- 2 microseconds or mult above a certain threshold in -4 - -2 microseconds"}; + Configurable requireNoCollInTimeRangeStrict{"requireNoCollInTimeRangeStrict", false, "reject collisions corrupted by the cannibalism, with other collisions within +/- 10 microseconds"}; + Configurable requireNoCollInTimeRangeNarrow{"requireNoCollInTimeRangeNarrow", false, "reject collisions corrupted by the cannibalism, with other collisions within +/- 2 microseconds"}; + Configurable requireNoCollInROFStd{"requireNoCollInROFStd", false, "reject collisions corrupted by the cannibalism, with other collisions within the same ITS ROF with mult. above a certain threshold"}; + Configurable requireNoCollInROFStrict{"requireNoCollInROFStrict", false, "reject collisions corrupted by the cannibalism, with other collisions within the same ITS ROF"}; + Configurable requireINEL0{"requireINEL0", true, "require INEL>0 event selection"}; + Configurable requireINEL1{"requireINEL1", false, "require INEL>1 event selection"}; + Configurable maxZVtxPosition{"maxZVtxPosition", 10., "max Z vtx position"}; + Configurable useEvtSelInDenomEff{"useEvtSelInDenomEff", false, "Consider event selections in the recoed <-> gen collision association for the denominator (or numerator) of the acc. x eff. (or signal loss)?"}; + Configurable applyZVtxSelOnMCPV{"applyZVtxSelOnMCPV", false, "Apply Z-vtx cut on the PV of the generated collision?"}; + Configurable useFT0CbasedOccupancy{"useFT0CbasedOccupancy", false, "Use sum of FT0-C amplitudes for estimating occupancy? (if not, use track-based definition)"}; + Configurable minOccupancy{"minOccupancy", -1, "minimum occupancy from neighbouring collisions"}; + Configurable maxOccupancy{"maxOccupancy", -1, "maximum occupancy from neighbouring collisions"}; + Configurable minIR{"minIR", -1, "minimum IR collisions"}; + Configurable maxIR{"maxIR", -1, "maximum IR collisions"}; + } eventSelections; + + //__________________________________________________ + // Photon (from the pi0) selections, applied on the de-referenced V0Cores + struct : ConfigurableGroup { + std::string prefix = "photonSelections"; // JSON group name + Configurable v0TypeSelection{"v0TypeSelection", 7, "select on a certain V0 type (leave negative if no selection desired)"}; + Configurable maxMass{"maxMass", 0.1, "Max photon mass (GeV/c^2)"}; + Configurable minRapidity{"minRapidity", -0.8, "Min photon rapidity"}; + Configurable maxRapidity{"maxRapidity", 0.8, "Max photon rapidity"}; + Configurable maxDauEta{"maxDauEta", 0.8, "Max |eta| of the daughter tracks"}; + Configurable minDCADauToPV{"minDCADauToPV", 0.0, "Min DCA of the daughter tracks to PV (cm)"}; + Configurable maxDCADau{"maxDCADau", 3.5, "Max DCA between the V0 daughters (cm)"}; + Configurable minRadius{"minRadius", 0.0, "Min photon conversion radius (cm)"}; + Configurable maxRadius{"maxRadius", 240.0, "Max photon conversion radius (cm)"}; + Configurable minCosPA{"minCosPA", 0.8, "Min photon cosine of pointing angle"}; + Configurable maxQt{"maxQt", 0.05, "Max qt (Armenteros-Podolanski) (GeV/c)"}; + Configurable maxAlpha{"maxAlpha", 0.95, "Max |alpha| (Armenteros-Podolanski)"}; + Configurable minPt{"minPt", 0.0, "Min photon pT (GeV/c)"}; + Configurable maxPt{"maxPt", 50.0, "Max photon pT (GeV/c)"}; + Configurable minTPCCrossedRows{"minTPCCrossedRows", 30, "Min number of TPC crossed rows of the daughter tracks"}; + Configurable minTPCNSigmaEl{"minTPCNSigmaEl", -7, "Min TPC NSigma (electron) of the daughter tracks"}; + Configurable maxTPCNSigmaEl{"maxTPCNSigmaEl", +7, "Max TPC NSigma (electron) of the daughter tracks"}; + } photonSelections; + + //__________________________________________________ + // Lambda selections, applied on the de-referenced V0Cores + struct : ConfigurableGroup { + std::string prefix = "lambdaSelections"; // JSON group name + Configurable v0TypeSelection{"v0TypeSelection", 1, "select on a certain V0 type (leave negative if no selection desired)"}; + Configurable massWindow{"massWindow", 0.015, "Lambda mass window around the PDG value (GeV/c^2)"}; + Configurable minRapidity{"minRapidity", -0.8, "Min Lambda rapidity"}; + Configurable maxRapidity{"maxRapidity", 0.8, "Max Lambda rapidity"}; + Configurable maxDauEta{"maxDauEta", 0.8, "Max |eta| of the daughter tracks"}; + Configurable minDCAPosToPV{"minDCAPosToPV", 0.05, "Min DCA of the positive daughter to PV (cm)"}; + Configurable minDCANegToPV{"minDCANegToPV", 0.05, "Min DCA of the negative daughter to PV (cm)"}; + Configurable maxDCADau{"maxDCADau", 1.0, "Max DCA between the V0 daughters (cm)"}; + Configurable minRadius{"minRadius", 0.5, "Min Lambda decay radius (cm)"}; + Configurable maxRadius{"maxRadius", 200.0, "Max Lambda decay radius (cm)"}; + Configurable minCosPA{"minCosPA", 0.95, "Min Lambda cosine of pointing angle (w.r.t. the PV)"}; + Configurable maxLifetime{"maxLifetime", 30.0, "Max Lambda proper lifetime, m*L/p (cm)"}; + Configurable minTPCCrossedRows{"minTPCCrossedRows", 70, "Min number of TPC crossed rows of the daughter tracks"}; + Configurable minITSclusters{"minITSclusters", -1, "Min number of ITS clusters of the daughter tracks (leave negative if no selection desired)"}; + Configurable maxTPCNSigmaPr{"maxTPCNSigmaPr", 5, "Max |TPC NSigma| (proton) of the baryon daughter"}; + Configurable maxTPCNSigmaPi{"maxTPCNSigmaPi", 5, "Max |TPC NSigma| (pion) of the meson daughter"}; + // Lambda / anti-Lambda are told apart with the Armenteros-Podolanski alpha: + // the baryon takes most of the momentum, so alpha > 0 for Lambda (positive + // daughter = proton) and alpha < 0 for anti-Lambda (negative daughter = anti-proton) + Configurable minAbsAlpha{"minAbsAlpha", 0.25, "Min |alpha| (Armenteros-Podolanski) to tag the Lambda charge"}; + Configurable maxAbsAlpha{"maxAbsAlpha", 1.0, "Max |alpha| (Armenteros-Podolanski)"}; + } lambdaSelections; + + //__________________________________________________ + // Xi0 selections, applied on the cascade itself + struct : ConfigurableGroup { + std::string prefix = "xi0Selections"; // JSON group name + Configurable minRapidity{"minRapidity", -0.5, "Min Xi0 rapidity"}; + Configurable maxRapidity{"maxRapidity", 0.5, "Max Xi0 rapidity"}; + Configurable minPt{"minPt", 0.0, "Min Xi0 pT (GeV/c)"}; + Configurable maxPt{"maxPt", 50.0, "Max Xi0 pT (GeV/c)"}; + Configurable pi0MassWindow{"pi0MassWindow", 0.035, "pi0 mass window around the PDG value (GeV/c^2)"}; + Configurable minPi0Radius{"minPi0Radius", -1, "Min pi0 decay radius (cm), leave negative if no selection desired"}; + Configurable maxPi0Radius{"maxPi0Radius", -1, "Max pi0 decay radius (cm), leave negative if no selection desired"}; + Configurable minPi0CosPA{"minPi0CosPA", -2, "Min pi0 cosine of pointing angle, leave below -1 if no selection desired"}; + Configurable maxDCAPi0Daughters{"maxDCAPi0Daughters", -1, "Max DCA between the two photons (cm), leave negative if no selection desired"}; + Configurable minCascRadius{"minCascRadius", 0.5, "Min Xi0 decay radius (cm)"}; + Configurable maxCascRadius{"maxCascRadius", 200.0, "Max Xi0 decay radius (cm)"}; + Configurable minCascCosPA{"minCascCosPA", 0.98, "Min Xi0 cosine of pointing angle"}; + Configurable maxDCACascDaughters{"maxDCACascDaughters", 1.0, "Max DCA between the Xi0 daughters (cm)"}; + Configurable maxDCAxyCascToPV{"maxDCAxyCascToPV", 1.0, "Max |DCAxy| of the Xi0 to the PV (cm), leave negative if no selection desired"}; + Configurable maxDCAzCascToPV{"maxDCAzCascToPV", -1, "Max |DCAz| of the Xi0 to the PV (cm), leave negative if no selection desired"}; + Configurable maxLifetime{"maxLifetime", 3.0, "Max Xi0 proper lifetime, in units of c*tau"}; + } xi0Selections; + + //__________________________________________________ + // MC-specific selections + struct : ConfigurableGroup { + std::string prefix = "mcSelections"; // JSON group name + Configurable doMCAssociation{"doMCAssociation", false, "Keep only candidates whose full decay chain is matched to a true Xi0/anti-Xi0"}; + Configurable requirePhysicalPrimary{"requirePhysicalPrimary", true, "Keep only physical primary Xi0, at reconstructed and generated level"}; + } mcSelections; + + //__________________________________________________ + // QA switches + Configurable doEventQA{"doEventQA", true, "do event QA histograms"}; + Configurable doCandidateQA{"doCandidateQA", true, "do candidate-level QA histograms"}; + Configurable doMCQA{"doMCQA", true, "do MC-specific QA histograms (pT resolution, background composition)"}; + + //__________________________________________________ + // Axes + ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0.0f, 1.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f}, "Centrality (%)"}; + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.5f, 1.0f, 1.5f, 2.0f, 2.5f, 3.0f, 3.5f, 4.0f, 4.5f, 5.0f, 6.0f, 7.0f, 8.0f, 10.0f, 12.0f, 15.0f, 20.0f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis axisXi0Mass{"axisXi0Mass", {200, 1.2f, 1.5f}, "#it{M}_{#Lambda#pi^{0}} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisPi0Mass{"axisPi0Mass", {200, 0.0f, 0.4f}, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.07f, 1.17f}, "#it{M}_{p#pi} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisRadius{"axisRadius", {200, 0.0f, 100.0f}, "Decay radius (cm)"}; + ConfigurableAxis axisCosPA{"axisCosPA", {200, 0.9f, 1.0f}, "cos(#theta_{PA})"}; + ConfigurableAxis axisDCADau{"axisDCADau", {200, 0.0f, 5.0f}, "DCA between daughters (cm)"}; + ConfigurableAxis axisDCAToPV{"axisDCAToPV", {200, -5.0f, 5.0f}, "DCA to PV (cm)"}; + ConfigurableAxis axisLifetime{"axisLifetime", {200, 0.0f, 50.0f}, "Proper lifetime (cm)"}; + ConfigurableAxis axisRapidity{"axisRapidity", {200, -1.0f, 1.0f}, "#it{y}"}; + ConfigurableAxis axisArmAlpha{"axisArmAlpha", {200, -1.0f, 1.0f}, "#alpha (Armenteros-Podolanski)"}; + ConfigurableAxis axisArmQt{"axisArmQt", {200, 0.0f, 0.5f}, "#it{q}_{T} (GeV/#it{c})"}; + ConfigurableAxis axisNch{"axisNch", {300, 0.0f, 3000.0f}, "#it{N}_{ch}"}; + ConfigurableAxis axisPtResolution{"axisPtResolution", {200, -1.0f, 1.0f}, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"}; + ConfigurableAxis axisVtxZ{"axisVtxZ", {40, -20.0f, 20.0f}, "PV #it{z} (cm)"}; + ConfigurableAxis axisPDGCode{"axisPDGCode", {10001, -5000.5f, 5000.5f}, "PDG code"}; + ConfigurableAxis axisNCandidates{"axisNCandidates", {10, -0.5f, 9.5f}, "Number of candidates"}; + + void init(InitContext const&) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + //_______________________________________________ + // Assemble the selection masks once and for all + maskTopological = 0; + BITSET(maskTopological, selPhotonV0Type); + BITSET(maskTopological, selPhotonMass); + BITSET(maskTopological, selPhotonRapidity); + BITSET(maskTopological, selPhotonDauEta); + BITSET(maskTopological, selPhotonDCADauToPV); + BITSET(maskTopological, selPhotonDCADau); + BITSET(maskTopological, selPhotonRadius); + BITSET(maskTopological, selPhotonCosPA); + BITSET(maskTopological, selPhotonArmenteros); + BITSET(maskTopological, selPhotonPt); + BITSET(maskTopological, selPhotonTPCCrossedRows); + BITSET(maskTopological, selPhotonTPCPID); + BITSET(maskTopological, selLambdaV0Type); + BITSET(maskTopological, selLambdaRapidity); + BITSET(maskTopological, selLambdaDauEta); + BITSET(maskTopological, selLambdaDCAPosToPV); + BITSET(maskTopological, selLambdaDCANegToPV); + BITSET(maskTopological, selLambdaDCADau); + BITSET(maskTopological, selLambdaRadius); + BITSET(maskTopological, selLambdaCosPA); + BITSET(maskTopological, selLambdaLifetime); + BITSET(maskTopological, selLambdaTPCCrossedRows); + BITSET(maskTopological, selLambdaITSClusters); + BITSET(maskTopological, selPi0Mass); + BITSET(maskTopological, selPi0Radius); + BITSET(maskTopological, selPi0CosPA); + BITSET(maskTopological, selPi0DCADau); + BITSET(maskTopological, selXi0Radius); + BITSET(maskTopological, selXi0CosPA); + BITSET(maskTopological, selXi0DCADau); + BITSET(maskTopological, selXi0DCAxyToPV); + BITSET(maskTopological, selXi0DCAzToPV); + BITSET(maskTopological, selXi0Lifetime); + BITSET(maskTopological, selXi0Rapidity); + BITSET(maskTopological, selXi0Pt); + + maskXi0Specific = 0; + BITSET(maskXi0Specific, selLambdaMass); + BITSET(maskXi0Specific, selLambdaArmenteros); + BITSET(maskXi0Specific, selTPCPIDPositiveProton); + BITSET(maskXi0Specific, selTPCPIDNegativePion); + BITSET(maskXi0Specific, selConsiderXi0); + + maskAntiXi0Specific = 0; + BITSET(maskAntiXi0Specific, selAntiLambdaMass); + BITSET(maskAntiXi0Specific, selAntiLambdaArmenteros); + BITSET(maskAntiXi0Specific, selTPCPIDNegativeProton); + BITSET(maskAntiXi0Specific, selTPCPIDPositivePion); + BITSET(maskAntiXi0Specific, selConsiderAntiXi0); + + if (mcSelections.requirePhysicalPrimary) { + BITSET(maskXi0Specific, selPhysPrimXi0); + BITSET(maskAntiXi0Specific, selPhysPrimAntiXi0); + } + + maskSelectionXi0 = maskTopological | maskXi0Specific; + maskSelectionAntiXi0 = maskTopological | maskAntiXi0Specific; + + //_______________________________________________ + // Event selection bookkeeping + auto hEventSelection = histos.add("hEventSelection", "hEventSelection", kTH1D, {{23, -0.5f, +22.5f}}); + hEventSelection->GetXaxis()->SetBinLabel(1, "All collisions"); + hEventSelection->GetXaxis()->SetBinLabel(2, "sel8 cut"); + hEventSelection->GetXaxis()->SetBinLabel(3, "kIsTriggerTVX"); + hEventSelection->GetXaxis()->SetBinLabel(4, "kNoITSROFrameBorder"); + hEventSelection->GetXaxis()->SetBinLabel(5, "kNoTimeFrameBorder"); + hEventSelection->GetXaxis()->SetBinLabel(6, "posZ cut"); + hEventSelection->GetXaxis()->SetBinLabel(7, "kIsVertexITSTPC"); + hEventSelection->GetXaxis()->SetBinLabel(8, "kIsGoodZvtxFT0vsPV"); + hEventSelection->GetXaxis()->SetBinLabel(9, "kIsVertexTOFmatched"); + hEventSelection->GetXaxis()->SetBinLabel(10, "kIsVertexTRDmatched"); + hEventSelection->GetXaxis()->SetBinLabel(11, "kNoSameBunchPileup"); + hEventSelection->GetXaxis()->SetBinLabel(12, "kNoCollInTimeRangeStd"); + hEventSelection->GetXaxis()->SetBinLabel(13, "kNoCollInTimeRangeStrict"); + hEventSelection->GetXaxis()->SetBinLabel(14, "kNoCollInTimeRangeNarrow"); + hEventSelection->GetXaxis()->SetBinLabel(15, "kNoCollInRofStd"); + hEventSelection->GetXaxis()->SetBinLabel(16, "kNoCollInRofStrict"); + hEventSelection->GetXaxis()->SetBinLabel(17, "INEL>0"); + hEventSelection->GetXaxis()->SetBinLabel(18, "INEL>1"); + hEventSelection->GetXaxis()->SetBinLabel(19, "Below min occup."); + hEventSelection->GetXaxis()->SetBinLabel(20, "Above max occup."); + hEventSelection->GetXaxis()->SetBinLabel(21, "Below min IR"); + hEventSelection->GetXaxis()->SetBinLabel(22, "Above max IR"); + hEventSelection->GetXaxis()->SetBinLabel(23, "Selected collisions"); + + histos.add("hEventCentrality", "hEventCentrality", kTH1D, {axisCentrality}); + histos.add("hEventPVz", "hEventPVz", kTH1D, {axisVtxZ}); + if (doEventQA) { + histos.add("EventQA/hCentralityVsNch", "hCentralityVsNch", kTH2D, {axisCentrality, axisNch}); + histos.add("EventQA/hCentralityVsPVz", "hCentralityVsPVz", kTH2D, {axisCentrality, axisVtxZ}); + } + + //_______________________________________________ + // Single candidate-selection bookkeeping histogram + auto hSelections = histos.add("hSelections", "hSelections", kTH1D, {{static_cast(selPhysPrimAntiXi0) + 3, -0.5f, static_cast(selPhysPrimAntiXi0) + 2.5f}}); + hSelections->GetXaxis()->SetBinLabel(1, "All"); + hSelections->GetXaxis()->SetBinLabel(selPhotonV0Type + 2, "#gamma V0 type"); + hSelections->GetXaxis()->SetBinLabel(selPhotonMass + 2, "#gamma mass"); + hSelections->GetXaxis()->SetBinLabel(selPhotonRapidity + 2, "#gamma rapidity"); + hSelections->GetXaxis()->SetBinLabel(selPhotonDauEta + 2, "#gamma dau. #eta"); + hSelections->GetXaxis()->SetBinLabel(selPhotonDCADauToPV + 2, "#gamma DCA dau. to PV"); + hSelections->GetXaxis()->SetBinLabel(selPhotonDCADau + 2, "#gamma DCA dau."); + hSelections->GetXaxis()->SetBinLabel(selPhotonRadius + 2, "#gamma radius"); + hSelections->GetXaxis()->SetBinLabel(selPhotonCosPA + 2, "#gamma cosPA"); + hSelections->GetXaxis()->SetBinLabel(selPhotonArmenteros + 2, "#gamma Arm. pod."); + hSelections->GetXaxis()->SetBinLabel(selPhotonPt + 2, "#gamma #it{p}_{T}"); + hSelections->GetXaxis()->SetBinLabel(selPhotonTPCCrossedRows + 2, "#gamma TPC rows"); + hSelections->GetXaxis()->SetBinLabel(selPhotonTPCPID + 2, "#gamma TPC PID"); + hSelections->GetXaxis()->SetBinLabel(selLambdaV0Type + 2, "#Lambda V0 type"); + hSelections->GetXaxis()->SetBinLabel(selLambdaRapidity + 2, "#Lambda rapidity"); + hSelections->GetXaxis()->SetBinLabel(selLambdaDauEta + 2, "#Lambda dau. #eta"); + hSelections->GetXaxis()->SetBinLabel(selLambdaDCAPosToPV + 2, "#Lambda DCA pos. to PV"); + hSelections->GetXaxis()->SetBinLabel(selLambdaDCANegToPV + 2, "#Lambda DCA neg. to PV"); + hSelections->GetXaxis()->SetBinLabel(selLambdaDCADau + 2, "#Lambda DCA dau."); + hSelections->GetXaxis()->SetBinLabel(selLambdaRadius + 2, "#Lambda radius"); + hSelections->GetXaxis()->SetBinLabel(selLambdaCosPA + 2, "#Lambda cosPA"); + hSelections->GetXaxis()->SetBinLabel(selLambdaLifetime + 2, "#Lambda lifetime"); + hSelections->GetXaxis()->SetBinLabel(selLambdaTPCCrossedRows + 2, "#Lambda TPC rows"); + hSelections->GetXaxis()->SetBinLabel(selLambdaITSClusters + 2, "#Lambda ITS clusters"); + hSelections->GetXaxis()->SetBinLabel(selLambdaMass + 2, "#Lambda mass"); + hSelections->GetXaxis()->SetBinLabel(selAntiLambdaMass + 2, "#bar{#Lambda} mass"); + hSelections->GetXaxis()->SetBinLabel(selLambdaArmenteros + 2, "#Lambda Arm. pod."); + hSelections->GetXaxis()->SetBinLabel(selAntiLambdaArmenteros + 2, "#bar{#Lambda} Arm. pod."); + hSelections->GetXaxis()->SetBinLabel(selTPCPIDPositiveProton + 2, "TPC PID p"); + hSelections->GetXaxis()->SetBinLabel(selTPCPIDNegativePion + 2, "TPC PID #pi^{-}"); + hSelections->GetXaxis()->SetBinLabel(selTPCPIDNegativeProton + 2, "TPC PID #bar{p}"); + hSelections->GetXaxis()->SetBinLabel(selTPCPIDPositivePion + 2, "TPC PID #pi^{+}"); + hSelections->GetXaxis()->SetBinLabel(selPi0Mass + 2, "#pi^{0} mass"); + hSelections->GetXaxis()->SetBinLabel(selPi0Radius + 2, "#pi^{0} radius"); + hSelections->GetXaxis()->SetBinLabel(selPi0CosPA + 2, "#pi^{0} cosPA"); + hSelections->GetXaxis()->SetBinLabel(selPi0DCADau + 2, "#pi^{0} DCA dau."); + hSelections->GetXaxis()->SetBinLabel(selXi0Radius + 2, "#Xi^{0} radius"); + hSelections->GetXaxis()->SetBinLabel(selXi0CosPA + 2, "#Xi^{0} cosPA"); + hSelections->GetXaxis()->SetBinLabel(selXi0DCADau + 2, "#Xi^{0} DCA dau."); + hSelections->GetXaxis()->SetBinLabel(selXi0DCAxyToPV + 2, "#Xi^{0} DCA_{xy} to PV"); + hSelections->GetXaxis()->SetBinLabel(selXi0DCAzToPV + 2, "#Xi^{0} DCA_{z} to PV"); + hSelections->GetXaxis()->SetBinLabel(selXi0Lifetime + 2, "#Xi^{0} lifetime"); + hSelections->GetXaxis()->SetBinLabel(selXi0Rapidity + 2, "#Xi^{0} rapidity"); + hSelections->GetXaxis()->SetBinLabel(selXi0Pt + 2, "#Xi^{0} #it{p}_{T}"); + hSelections->GetXaxis()->SetBinLabel(selConsiderXi0 + 2, "True #Xi^{0}"); + hSelections->GetXaxis()->SetBinLabel(selConsiderAntiXi0 + 2, "True #bar{#Xi^{0}}"); + hSelections->GetXaxis()->SetBinLabel(selPhysPrimXi0 + 2, "Phys. prim. #Xi^{0}"); + hSelections->GetXaxis()->SetBinLabel(selPhysPrimAntiXi0 + 2, "Phys. prim. #bar{#Xi^{0}}"); + hSelections->GetXaxis()->SetBinLabel(selPhysPrimAntiXi0 + 3, "Cand. selected"); + + //_______________________________________________ + // Main analysis output + histos.add("h3dMassXi0", "h3dMassXi0", kTH3D, {axisCentrality, axisPt, axisXi0Mass}); + histos.add("h3dMassAntiXi0", "h3dMassAntiXi0", kTH3D, {axisCentrality, axisPt, axisXi0Mass}); + histos.add("h2dNbrOfXi0VsCentrality", "h2dNbrOfXi0VsCentrality", kTH2D, {axisCentrality, axisNCandidates}); + histos.add("h2dNbrOfAntiXi0VsCentrality", "h2dNbrOfAntiXi0VsCentrality", kTH2D, {axisCentrality, axisNCandidates}); + + //_______________________________________________ + // Candidate QA, before and after the selections + if (doCandidateQA) { + for (int mode = 0; mode < NSelectionStages; mode++) { + const std::string dir = (mode == 0) ? "QA/BeforeSel/" : "QA/AfterSel/"; + histos.add(dir + "h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisXi0Mass}); + histos.add(dir + "hMassPi0", "hMassPi0", kTH1D, {axisPi0Mass}); + histos.add(dir + "hMassLambda", "hMassLambda", kTH1D, {axisLambdaMass}); + histos.add(dir + "hPt", "hPt", kTH1D, {axisPt}); + histos.add(dir + "hRapidity", "hRapidity", kTH1D, {axisRapidity}); + histos.add(dir + "hCascRadius", "hCascRadius", kTH1D, {axisRadius}); + histos.add(dir + "hCascCosPA", "hCascCosPA", kTH1D, {axisCosPA}); + histos.add(dir + "hDCACascDaughters", "hDCACascDaughters", kTH1D, {axisDCADau}); + histos.add(dir + "hDCAxyCascToPV", "hDCAxyCascToPV", kTH1D, {axisDCAToPV}); + histos.add(dir + "hDCAzCascToPV", "hDCAzCascToPV", kTH1D, {axisDCAToPV}); + histos.add(dir + "hLifetime", "hLifetime", kTH1D, {axisLifetime}); + histos.add(dir + "hLambdaRadius", "hLambdaRadius", kTH1D, {axisRadius}); + histos.add(dir + "hLambdaCosPA", "hLambdaCosPA", kTH1D, {axisCosPA}); + histos.add(dir + "hPi0Radius", "hPi0Radius", kTH1D, {axisRadius}); + histos.add(dir + "hPi0CosPA", "hPi0CosPA", kTH1D, {axisCosPA}); + histos.add(dir + "hDCAPi0Daughters", "hDCAPi0Daughters", kTH1D, {axisDCADau}); + histos.add(dir + "h2dArmenterosLambda", "h2dArmenterosLambda", kTH2D, {axisArmAlpha, axisArmQt}); + } + } + + //_______________________________________________ + // MC-specific histograms + if (doprocessMonteCarlo && doMCQA) { + histos.add("MCQA/h2dPtResolution", "h2dPtResolution", kTH2D, {axisPt, axisPtResolution}); + histos.add("MCQA/h2dPtRecoVsPtGen", "h2dPtRecoVsPtGen", kTH2D, {axisPt, axisPt}); + histos.add("MCQA/h3dMassXi0VsPtReco", "h3dMassXi0VsPtReco", kTH3D, {axisCentrality, axisPt, axisXi0Mass}); + histos.add("MCQA/h3dMassAntiXi0VsPtReco", "h3dMassAntiXi0VsPtReco", kTH3D, {axisCentrality, axisPt, axisXi0Mass}); + histos.add("MCQA/hXi0PDGCode", "hXi0PDGCode", kTH1D, {axisPDGCode}); + histos.add("MCQA/hXi0PDGCodeMother", "hXi0PDGCodeMother", kTH1D, {axisPDGCode}); + histos.add("MCQA/h2dPhoton1VsPhoton2PDGCode", "h2dPhoton1VsPhoton2PDGCode", kTH2D, {axisPDGCode, axisPDGCode}); + histos.add("MCQA/hLambdaPDGCode", "hLambdaPDGCode", kTH1D, {axisPDGCode}); + } + + //_______________________________________________ + // Generated-level histograms + if (doprocessGenerated) { + histos.add("Gen/hGenEvents", "hGenEvents", kTH2D, {axisNch, {2, -0.5f, +1.5f}}); + histos.add("Gen/hGenEventCentrality", "hGenEventCentrality", kTH1D, {axisCentrality}); + histos.add("Gen/hCentralityVsNcoll_beforeEvSel", "hCentralityVsNcoll_beforeEvSel", kTH2D, {axisCentrality, {50, -0.5f, +49.5f}}); + histos.add("Gen/hCentralityVsNcoll_afterEvSel", "hCentralityVsNcoll_afterEvSel", kTH2D, {axisCentrality, {50, -0.5f, +49.5f}}); + histos.add("Gen/hCentralityVsMultMC", "hCentralityVsMultMC", kTH2D, {axisCentrality, axisNch}); + + histos.add("Gen/h2dGenXi0", "h2dGenXi0", kTH2D, {axisCentrality, axisPt}); + histos.add("Gen/h2dGenAntiXi0", "h2dGenAntiXi0", kTH2D, {axisCentrality, axisPt}); + histos.add("Gen/h2dGenXi0VsMultMC", "h2dGenXi0VsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("Gen/h2dGenAntiXi0VsMultMC", "h2dGenAntiXi0VsMultMC", kTH2D, {axisNch, axisPt}); + histos.add("Gen/h2dGenXi0_RecoedEvt", "h2dGenXi0_RecoedEvt", kTH2D, {axisCentrality, axisPt}); + histos.add("Gen/h2dGenAntiXi0_RecoedEvt", "h2dGenAntiXi0_RecoedEvt", kTH2D, {axisCentrality, axisPt}); + histos.add("Gen/hGenXi0Rapidity", "hGenXi0Rapidity", kTH1D, {axisRapidity}); + } + } + + //_______________________________________________ + // Check that all the bits of a mask are set in a bitmap + bool verifyMask(uint64_t bitmap, uint64_t mask) + { + return (bitmap & mask) == mask; + } + + //_______________________________________________ + // Centrality getter (Run 3) + template + float getCentralityRun3(TCollision const& collision) + { + if (centralityEstimator == kCentFT0C) + return collision.centFT0C(); + else if (centralityEstimator == kCentFT0M) + return collision.centFT0M(); + else if (centralityEstimator == kCentFT0CVariant1) + return collision.centFT0CVariant1(); + else if (centralityEstimator == kCentMFT) + return collision.centMFT(); + else if (centralityEstimator == kCentNGlobal) + return collision.centNGlobal(); + else if (centralityEstimator == kCentFV0A) + return collision.centFV0A(); + + return -1.f; + } + + //_______________________________________________ + // Check whether the collision passes our collision selections + template + bool isEventAccepted(TCollision const& collision, bool fillHists) + { + if (fillHists) + histos.fill(HIST("hEventSelection"), 0. /* all collisions */); + if (eventSelections.requireSel8 && !collision.sel8()) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 1 /* sel8 collisions */); + if (eventSelections.requireTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 2 /* FT0 vertex (acceptable FT0C-FT0A time difference) collisions */); + if (eventSelections.rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 3 /* Not at ITS ROF border */); + if (eventSelections.rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 4 /* Not at TF border */); + if (std::abs(collision.posZ()) > eventSelections.maxZVtxPosition) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 5 /* vertex-Z selected */); + if (eventSelections.requireIsVertexITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 6 /* Contains at least one ITS-TPC track */); + if (eventSelections.requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 7 /* PV position consistency check */); + if (eventSelections.requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 8 /* PV with at least one contributor matched with TOF */); + if (eventSelections.requireIsVertexTRDmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTRDmatched)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 9 /* PV with at least one contributor matched with TRD */); + if (eventSelections.rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 10 /* Not at same bunch pile-up */); + if (eventSelections.requireNoCollInTimeRangeStd && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 11 /* No other collision within +/- 2 microseconds or mult above a certain threshold in -4 - -2 microseconds */); + if (eventSelections.requireNoCollInTimeRangeStrict && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStrict)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 12 /* No other collision within +/- 10 microseconds */); + if (eventSelections.requireNoCollInTimeRangeNarrow && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 13 /* No other collision within +/- 2 microseconds */); + if (eventSelections.requireNoCollInROFStd && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 14 /* No other collision within the same ITS ROF with mult. above a certain threshold */); + if (eventSelections.requireNoCollInROFStrict && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStrict)) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 15 /* No other collision within the same ITS ROF */); + if (eventSelections.requireINEL0 && collision.multNTracksPVeta1() < 1) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 16 /* INEL > 0 */); + if (eventSelections.requireINEL1 && collision.multNTracksPVeta1() < 2) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 17 /* INEL > 1 */); + + float collisionOccupancy = eventSelections.useFT0CbasedOccupancy ? collision.ft0cOccupancyInTimeRange() : collision.trackOccupancyInTimeRange(); + if (eventSelections.minOccupancy >= 0 && collisionOccupancy < eventSelections.minOccupancy) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 18 /* Below min occupancy */); + if (eventSelections.maxOccupancy >= 0 && collisionOccupancy > eventSelections.maxOccupancy) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 19 /* Above max occupancy */); + + // Fetch interaction rate only if required (in order to limit ccdb calls) + float interactionRate = (fGetIR) ? rateFetcher.fetch(ccdb.service, collision.timestamp(), collision.runNumber(), irSource, fIRCrashOnNull) * 1.e-3 : -1; + if (eventSelections.minIR >= 0 && interactionRate < eventSelections.minIR) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 20 /* Below min IR */); + if (eventSelections.maxIR >= 0 && interactionRate > eventSelections.maxIR) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 21 /* Above max IR */); + + if (fillHists) { + histos.fill(HIST("hEventSelection"), 22 /* selected collisions */); + float centrality = getCentralityRun3(collision); + histos.fill(HIST("hEventCentrality"), centrality); + histos.fill(HIST("hEventPVz"), collision.posZ()); + if (doEventQA) { + histos.fill(HIST("EventQA/hCentralityVsNch"), centrality, collision.multNTracksPVeta1()); + histos.fill(HIST("EventQA/hCentralityVsPVz"), centrality, collision.posZ()); + } + } + + return true; + } + + //_______________________________________________ + // Photon (pi0 daughter) selection bitmap, computed on the de-referenced V0Cores entry. + // The bitmaps of the two photons are AND-ed, so that a bit is set only if + // both photons of the pi0 satisfy the corresponding criterion. + template + uint64_t computePhotonBitmap(TV0Object const& gamma) + { + uint64_t bitMap = 0; + + if (photonSelections.v0TypeSelection < 0 || gamma.v0Type() == photonSelections.v0TypeSelection) + BITSET(bitMap, selPhotonV0Type); + + if (gamma.mGamma() > 0 && gamma.mGamma() < photonSelections.maxMass) + BITSET(bitMap, selPhotonMass); + + const float photonY = RecoDecay::y(std::array{gamma.px(), gamma.py(), gamma.pz()}, o2::constants::physics::MassGamma); + if (photonY > photonSelections.minRapidity && photonY < photonSelections.maxRapidity) + BITSET(bitMap, selPhotonRapidity); + + if (std::abs(gamma.positiveeta()) < photonSelections.maxDauEta && std::abs(gamma.negativeeta()) < photonSelections.maxDauEta) + BITSET(bitMap, selPhotonDauEta); + + if (std::abs(gamma.dcapostopv()) > photonSelections.minDCADauToPV && std::abs(gamma.dcanegtopv()) > photonSelections.minDCADauToPV) + BITSET(bitMap, selPhotonDCADauToPV); + + if (std::abs(gamma.dcaV0daughters()) < photonSelections.maxDCADau) + BITSET(bitMap, selPhotonDCADau); + + if (gamma.v0radius() > photonSelections.minRadius && gamma.v0radius() < photonSelections.maxRadius) + BITSET(bitMap, selPhotonRadius); + + if (gamma.v0cosPA() > photonSelections.minCosPA) + BITSET(bitMap, selPhotonCosPA); + + if (gamma.qtarm() < photonSelections.maxQt && std::abs(gamma.alpha()) < photonSelections.maxAlpha) + BITSET(bitMap, selPhotonArmenteros); + + if (gamma.pt() > photonSelections.minPt && gamma.pt() < photonSelections.maxPt) + BITSET(bitMap, selPhotonPt); + + auto posTrack = gamma.template posTrackExtra_as(); + auto negTrack = gamma.template negTrackExtra_as(); + + if (posTrack.tpcCrossedRows() >= photonSelections.minTPCCrossedRows && negTrack.tpcCrossedRows() >= photonSelections.minTPCCrossedRows) + BITSET(bitMap, selPhotonTPCCrossedRows); + + if (posTrack.tpcNSigmaEl() > photonSelections.minTPCNSigmaEl && posTrack.tpcNSigmaEl() < photonSelections.maxTPCNSigmaEl && + negTrack.tpcNSigmaEl() > photonSelections.minTPCNSigmaEl && negTrack.tpcNSigmaEl() < photonSelections.maxTPCNSigmaEl) + BITSET(bitMap, selPhotonTPCPID); + + return bitMap; + } + + //_______________________________________________ + // Lambda selection bitmap, computed on the de-referenced V0Cores entry + template + uint64_t computeLambdaBitmap(TV0Object const& lambda, TCollision const& collision) + { + uint64_t bitMap = 0; + + if (lambdaSelections.v0TypeSelection < 0 || lambda.v0Type() == lambdaSelections.v0TypeSelection) + BITSET(bitMap, selLambdaV0Type); + + if (lambda.yLambda() > lambdaSelections.minRapidity && lambda.yLambda() < lambdaSelections.maxRapidity) + BITSET(bitMap, selLambdaRapidity); + + if (std::abs(lambda.positiveeta()) < lambdaSelections.maxDauEta && std::abs(lambda.negativeeta()) < lambdaSelections.maxDauEta) + BITSET(bitMap, selLambdaDauEta); + + if (std::abs(lambda.dcapostopv()) > lambdaSelections.minDCAPosToPV) + BITSET(bitMap, selLambdaDCAPosToPV); + + if (std::abs(lambda.dcanegtopv()) > lambdaSelections.minDCANegToPV) + BITSET(bitMap, selLambdaDCANegToPV); + + if (std::abs(lambda.dcaV0daughters()) < lambdaSelections.maxDCADau) + BITSET(bitMap, selLambdaDCADau); + + if (lambda.v0radius() > lambdaSelections.minRadius && lambda.v0radius() < lambdaSelections.maxRadius) + BITSET(bitMap, selLambdaRadius); + + if (lambda.v0cosPA() > lambdaSelections.minCosPA) + BITSET(bitMap, selLambdaCosPA); + + const float lambdaLifetime = lambda.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassLambda0; + if (lambdaSelections.maxLifetime < 0 || lambdaLifetime < lambdaSelections.maxLifetime) + BITSET(bitMap, selLambdaLifetime); + + // Lambda / anti-Lambda separation via the Armenteros-Podolanski alpha: + // in the Lambda decay the baryon carries most of the longitudinal momentum + if (lambda.alpha() > lambdaSelections.minAbsAlpha && lambda.alpha() < lambdaSelections.maxAbsAlpha) + BITSET(bitMap, selLambdaArmenteros); + if (lambda.alpha() < -lambdaSelections.minAbsAlpha && lambda.alpha() > -lambdaSelections.maxAbsAlpha) + BITSET(bitMap, selAntiLambdaArmenteros); + + if (std::abs(lambda.mLambda() - o2::constants::physics::MassLambda0) < lambdaSelections.massWindow) + BITSET(bitMap, selLambdaMass); + if (std::abs(lambda.mAntiLambda() - o2::constants::physics::MassLambda0) < lambdaSelections.massWindow) + BITSET(bitMap, selAntiLambdaMass); + + auto posTrack = lambda.template posTrackExtra_as(); + auto negTrack = lambda.template negTrackExtra_as(); + + if (posTrack.tpcCrossedRows() >= lambdaSelections.minTPCCrossedRows && negTrack.tpcCrossedRows() >= lambdaSelections.minTPCCrossedRows) + BITSET(bitMap, selLambdaTPCCrossedRows); + + if (lambdaSelections.minITSclusters < 0 || + (posTrack.itsNCls() >= lambdaSelections.minITSclusters && negTrack.itsNCls() >= lambdaSelections.minITSclusters)) + BITSET(bitMap, selLambdaITSClusters); + + // TPC PID, for both mass hypotheses + if (std::abs(posTrack.tpcNSigmaPr()) < lambdaSelections.maxTPCNSigmaPr) + BITSET(bitMap, selTPCPIDPositiveProton); + if (std::abs(negTrack.tpcNSigmaPi()) < lambdaSelections.maxTPCNSigmaPi) + BITSET(bitMap, selTPCPIDNegativePion); + if (std::abs(negTrack.tpcNSigmaPr()) < lambdaSelections.maxTPCNSigmaPr) + BITSET(bitMap, selTPCPIDNegativeProton); + if (std::abs(posTrack.tpcNSigmaPi()) < lambdaSelections.maxTPCNSigmaPi) + BITSET(bitMap, selTPCPIDPositivePion); + + return bitMap; + } + + //_______________________________________________ + // pi0 and Xi0 (cascade) selection bitmap. + // rapidity and pt are the generated ones when running over MC. + template + uint64_t computeCascadeBitmap(TXi0Object const& xi0, TCollision const& collision, float rapidity, float pt) + { + uint64_t bitMap = 0; + + if (std::abs(xi0.pi0Mass() - o2::constants::physics::MassPi0) < xi0Selections.pi0MassWindow) + BITSET(bitMap, selPi0Mass); + + if ((xi0Selections.minPi0Radius < 0 || xi0.radiusPi0() > xi0Selections.minPi0Radius) && + (xi0Selections.maxPi0Radius < 0 || xi0.radiusPi0() < xi0Selections.maxPi0Radius)) + BITSET(bitMap, selPi0Radius); + + if (xi0Selections.minPi0CosPA < -1 || xi0.pi0CosPA(collision.posX(), collision.posY(), collision.posZ()) > xi0Selections.minPi0CosPA) + BITSET(bitMap, selPi0CosPA); + + if (xi0Selections.maxDCAPi0Daughters < 0 || xi0.dcadaughtersPi0() < xi0Selections.maxDCAPi0Daughters) + BITSET(bitMap, selPi0DCADau); + + if (xi0.radius() > xi0Selections.minCascRadius && xi0.radius() < xi0Selections.maxCascRadius) + BITSET(bitMap, selXi0Radius); + + if (xi0.cascCosPA(collision.posX(), collision.posY(), collision.posZ()) > xi0Selections.minCascCosPA) + BITSET(bitMap, selXi0CosPA); + + if (xi0.dcadaughters() < xi0Selections.maxDCACascDaughters) + BITSET(bitMap, selXi0DCADau); + + if (xi0Selections.maxDCAxyCascToPV < 0 || std::abs(xi0.dcaXYCascToPV()) < xi0Selections.maxDCAxyCascToPV) + BITSET(bitMap, selXi0DCAxyToPV); + + if (xi0Selections.maxDCAzCascToPV < 0 || std::abs(xi0.dcaZCascToPV()) < xi0Selections.maxDCAzCascToPV) + BITSET(bitMap, selXi0DCAzToPV); + + if (xi0Selections.maxLifetime < 0 || getProperLifetime(xi0, collision) < xi0Selections.maxLifetime * CtauXi0) + BITSET(bitMap, selXi0Lifetime); + + if (rapidity > xi0Selections.minRapidity && rapidity < xi0Selections.maxRapidity) + BITSET(bitMap, selXi0Rapidity); + + if (pt > xi0Selections.minPt && pt < xi0Selections.maxPt) + BITSET(bitMap, selXi0Pt); + + return bitMap; + } + + //_______________________________________________ + // MC association bitmap. + // The builder sets the Xi0 PDG code only if the full decay chain is matched, + // i.e. two true photons from the same true pi0 and a true Lambda sharing the + // same true Xi0 mother, so a single PDG check is enough here. + template + uint64_t computeMCAssociation(TXi0Object const& xi0) + { + uint64_t bitMap = 0; + + if (xi0.pdgCode() == o2::constants::physics::Pdg::kXi0) { + BITSET(bitMap, selConsiderXi0); + if (xi0.isPhysicalPrimary()) + BITSET(bitMap, selPhysPrimXi0); + } + if (xi0.pdgCode() == -o2::constants::physics::Pdg::kXi0) { + BITSET(bitMap, selConsiderAntiXi0); + if (xi0.isPhysicalPrimary()) + BITSET(bitMap, selPhysPrimAntiXi0); + } + + return bitMap; + } + + //_______________________________________________ + // Xi0 proper lifetime, m*L/p, in cm + template + float getProperLifetime(TXi0Object const& xi0, TCollision const& collision) + { + const float decayLength = RecoDecay::distance(std::array{collision.posX(), collision.posY(), collision.posZ()}, + std::array{xi0.x(), xi0.y(), xi0.z()}); + const float totalMomentum = xi0.p(); + return (totalMomentum > 0) ? o2::constants::physics::MassXi0 * decayLength / totalMomentum : 1e6f; + } + + //_______________________________________________ + // Fill the candidate QA histograms. mode = 0: before selections, 1: after selections + template + void fillCandidateQA(TXi0Object const& xi0, TV0Object const& lambda, TCollision const& collision, float pt, float rapidity, float centrality, float massLambda) + { + static constexpr std::string_view MainDir[] = {"QA/BeforeSel", "QA/AfterSel"}; + + histos.fill(HIST(MainDir[mode]) + HIST("/h3dMass"), centrality, pt, xi0.xi0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/hMassPi0"), xi0.pi0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/hMassLambda"), massLambda); + histos.fill(HIST(MainDir[mode]) + HIST("/hPt"), pt); + histos.fill(HIST(MainDir[mode]) + HIST("/hRapidity"), rapidity); + histos.fill(HIST(MainDir[mode]) + HIST("/hCascRadius"), xi0.radius()); + histos.fill(HIST(MainDir[mode]) + HIST("/hCascCosPA"), xi0.cascCosPA(collision.posX(), collision.posY(), collision.posZ())); + histos.fill(HIST(MainDir[mode]) + HIST("/hDCACascDaughters"), xi0.dcadaughters()); + histos.fill(HIST(MainDir[mode]) + HIST("/hDCAxyCascToPV"), xi0.dcaXYCascToPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/hDCAzCascToPV"), xi0.dcaZCascToPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/hLifetime"), getProperLifetime(xi0, collision)); + histos.fill(HIST(MainDir[mode]) + HIST("/hLambdaRadius"), xi0.radiusLambda()); + histos.fill(HIST(MainDir[mode]) + HIST("/hLambdaCosPA"), xi0.lambdaCosPA(collision.posX(), collision.posY(), collision.posZ())); + histos.fill(HIST(MainDir[mode]) + HIST("/hPi0Radius"), xi0.radiusPi0()); + histos.fill(HIST(MainDir[mode]) + HIST("/hPi0CosPA"), xi0.pi0CosPA(collision.posX(), collision.posY(), collision.posZ())); + histos.fill(HIST(MainDir[mode]) + HIST("/hDCAPi0Daughters"), xi0.dcadaughtersPi0()); + histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenterosLambda"), lambda.alpha(), lambda.qtarm()); + } + + //_______________________________________________ + // Main reconstructed-level analysis function, common to data and MC + template + void analyseRecoedXi0s(TCollisions const& collisions, TXi0s const& fullXi0s, TV0s const& fullV0s) + { + // Custom grouping: the Xi0 candidates are not sorted per collision + std::vector> xi0grouped(collisions.size()); + for (const auto& xi0 : fullXi0s) { + xi0grouped[xi0.straCollisionId()].push_back(xi0.globalIndex()); + } + + const int64_t nV0s = fullV0s.size(); + + for (const auto& coll : collisions) { + // Event selection + if (!isEventAccepted(coll, true)) { + continue; + } + const float centrality = getCentralityRun3(coll); + + int nXi0s = 0; + int nAntiXi0s = 0; + + // Xi0 candidates loop + for (size_t i = 0; i < xi0grouped[coll.globalIndex()].size(); i++) { + auto xi0 = fullXi0s.rawIteratorAt(xi0grouped[coll.globalIndex()][i]); + + histos.fill(HIST("hSelections"), 0 /* all candidates */); + + //_______________________________________________ + // De-reference the three daughter V0s + if (xi0.photon1Index() < 0 || xi0.photon1Index() >= nV0s || + xi0.photon2Index() < 0 || xi0.photon2Index() >= nV0s || + xi0.lambdaIndex() < 0 || xi0.lambdaIndex() >= nV0s) { + continue; + } + auto photon1 = fullV0s.rawIteratorAt(xi0.photon1Index()); + auto photon2 = fullV0s.rawIteratorAt(xi0.photon2Index()); + auto lambda = fullV0s.rawIteratorAt(xi0.lambdaIndex()); + + //_______________________________________________ + // Kinematics entering both the selections and the histograms. + // Over MC, the generated rapidity and pT are used instead of the + // reconstructed ones, so that the acceptance x efficiency numerator + // and denominator are defined identically. + const float ptRecoed = xi0.pt(); + float pt = ptRecoed; + float rapidity = xi0.rapidity(); + bool hasMCInfo = false; + if constexpr (requires { xi0.pdgCode(); }) { + // the MC information is available only if the three daughter V0s + // could be associated to a V0MCCore entry by the builder + hasMCInfo = (xi0.photon1PDGCode() != 0 && xi0.photon2PDGCode() != 0 && xi0.lambdaPDGCode() != 0); + if (hasMCInfo) { + pt = xi0.mcpt(); + rapidity = xi0.rapidityMC(); + } + } + + //_______________________________________________ + // Assemble the selection bitmap + uint64_t selMap = computePhotonBitmap(photon1) & computePhotonBitmap(photon2); + selMap |= computeLambdaBitmap(lambda, coll); + selMap |= computeCascadeBitmap(xi0, coll, rapidity, pt); + + if constexpr (requires { xi0.pdgCode(); }) { + selMap |= computeMCAssociation(xi0); + + if (doMCQA) { + histos.fill(HIST("MCQA/hXi0PDGCode"), xi0.pdgCode()); + histos.fill(HIST("MCQA/hXi0PDGCodeMother"), xi0.pdgCodeMother()); + histos.fill(HIST("MCQA/h2dPhoton1VsPhoton2PDGCode"), xi0.photon1PDGCode(), xi0.photon2PDGCode()); + histos.fill(HIST("MCQA/hLambdaPDGCode"), xi0.lambdaPDGCode()); + } + + // disregard the MC association if not explicitly asked for + if (!mcSelections.doMCAssociation) { + BITSET(selMap, selConsiderXi0); + BITSET(selMap, selConsiderAntiXi0); + BITSET(selMap, selPhysPrimXi0); + BITSET(selMap, selPhysPrimAntiXi0); + } + } else { + // real data: every candidate is considered for both hypotheses, + // the Armenteros bits alone tell Xi0 and anti-Xi0 apart + BITSET(selMap, selConsiderXi0); + BITSET(selMap, selConsiderAntiXi0); + BITSET(selMap, selPhysPrimXi0); + BITSET(selMap, selPhysPrimAntiXi0); + } + + // Selection bookkeeping: one bin per criterion + for (int bit = 0; bit <= static_cast(selPhysPrimAntiXi0); bit++) { + if (BITCHECK(selMap, bit)) { + histos.fill(HIST("hSelections"), bit + 1); + } + } + + const bool passXi0 = verifyMask(selMap, maskSelectionXi0); + const bool passAntiXi0 = verifyMask(selMap, maskSelectionAntiXi0); + // the sign of the Armenteros alpha already tells the two hypotheses apart + const float massLambda = (lambda.alpha() < 0) ? lambda.mAntiLambda() : lambda.mLambda(); + + if (doCandidateQA) { + fillCandidateQA<0>(xi0, lambda, coll, pt, rapidity, centrality, massLambda); + } + + if (!passXi0 && !passAntiXi0) { + continue; + } + histos.fill(HIST("hSelections"), static_cast(selPhysPrimAntiXi0) + 2 /* candidate selected */); + + if (doCandidateQA) { + fillCandidateQA<1>(xi0, lambda, coll, pt, rapidity, centrality, massLambda); + } + + if (passXi0) { + nXi0s++; + histos.fill(HIST("h3dMassXi0"), centrality, pt, xi0.xi0Mass()); + } + if (passAntiXi0) { + nAntiXi0s++; + histos.fill(HIST("h3dMassAntiXi0"), centrality, pt, xi0.xi0Mass()); + } + + //_______________________________________________ + // MC-specific QA + if constexpr (requires { xi0.pdgCode(); }) { + if (doMCQA) { + if (passXi0) { + histos.fill(HIST("MCQA/h3dMassXi0VsPtReco"), centrality, ptRecoed, xi0.xi0Mass()); + } + if (passAntiXi0) { + histos.fill(HIST("MCQA/h3dMassAntiXi0VsPtReco"), centrality, ptRecoed, xi0.xi0Mass()); + } + if (hasMCInfo && pt > 0) { + histos.fill(HIST("MCQA/h2dPtResolution"), pt, (ptRecoed - pt) / pt); + histos.fill(HIST("MCQA/h2dPtRecoVsPtGen"), pt, ptRecoed); + } + } + } + } + + histos.fill(HIST("h2dNbrOfXi0VsCentrality"), centrality, nXi0s); + histos.fill(HIST("h2dNbrOfAntiXi0VsCentrality"), centrality, nAntiXi0s); + } + } + + //_______________________________________________ + // Generated-level processing + // Return the list of indices of the reconstructed collision associated to each MC collision + template + std::vector getListOfRecoCollIndices(TMCCollisions const& mcCollisions, TCollisions const& collisions) + { + std::vector listBestCollisionIdx(mcCollisions.size(), -1); + std::vector> groupedCollisions(mcCollisions.size()); + + for (const auto& coll : collisions) { + if (coll.straMCCollisionId() < 0) { + continue; + } + groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); + } + + for (auto const& mcCollision : mcCollisions) { + int biggestNContribs = -1; + int bestCollisionIndex = -1; + for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { + auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); + + if (eventSelections.useEvtSelInDenomEff && !isEventAccepted(collision, false)) { + continue; + } + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + bestCollisionIndex = collision.globalIndex(); + } + } + listBestCollisionIdx[mcCollision.globalIndex()] = bestCollisionIndex; + } + return listBestCollisionIdx; + } + + //_______________________________________________ + // Event selection applied at generated level + template + bool isGeneratedEventAccepted(TMCCollision const& mcCollision) + { + if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) { + return false; + } + if (eventSelections.requireINEL0 && mcCollision.multMCNParticlesEta10() < 1) { + return false; + } + if (eventSelections.requireINEL1 && mcCollision.multMCNParticlesEta10() < 2) { + return false; + } + return true; + } + + //_______________________________________________ + // Generated-level processing: generated event properties (event loss/splitting) + template + void fillGeneratedEventProperties(TMCCollisions const& mcCollisions, TCollisions const& collisions) + { + std::vector> groupedCollisions(mcCollisions.size()); + for (const auto& coll : collisions) { + if (coll.straMCCollisionId() < 0) { + continue; + } + groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); + } + + for (auto const& mcCollision : mcCollisions) { + if (!isGeneratedEventAccepted(mcCollision)) { + continue; + } + + histos.fill(HIST("Gen/hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events */); + + bool atLeastOne = false; + int biggestNContribs = -1; + float centrality = 100.5f; + int nCollisions = 0; + for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { + auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); + + if (!isEventAccepted(collision, false)) { + continue; + } + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + centrality = getCentralityRun3(collision); + } + nCollisions++; + atLeastOne = true; + } + + histos.fill(HIST("Gen/hCentralityVsNcoll_beforeEvSel"), centrality, static_cast(groupedCollisions[mcCollision.globalIndex()].size())); + histos.fill(HIST("Gen/hCentralityVsNcoll_afterEvSel"), centrality, nCollisions); + histos.fill(HIST("Gen/hCentralityVsMultMC"), centrality, mcCollision.multMCNParticlesEta05()); + + if (atLeastOne) { + histos.fill(HIST("Gen/hGenEvents"), mcCollision.multMCNParticlesEta05(), 1 /* at least 1 rec. event */); + histos.fill(HIST("Gen/hGenEventCentrality"), centrality); + } + } + } + + //_______________________________________________ + // Main generated-level analysis function. + // The generated Xi0 are read from the CascMCCores table. + template + void analyseGeneratedXi0s(TMCCollisions const& mcCollisions, TCollisions const& collisions, TCascMCs const& cascMCCores) + { + fillGeneratedEventProperties(mcCollisions, collisions); + std::vector listBestCollisionIdx = getListOfRecoCollIndices(mcCollisions, collisions); + + for (auto const& cascMC : cascMCCores) { + if (std::abs(cascMC.pdgCode()) != o2::constants::physics::Pdg::kXi0) { + continue; + } + if (!cascMC.has_straMCCollision()) { + continue; + } + if (mcSelections.requirePhysicalPrimary && !cascMC.isPhysicalPrimary()) { + continue; + } + + // N.B.: cascdata::RapidityMC only knows about the Xi- and Omega- masses, + // so the Xi0 rapidity is computed here explicitly + const float ptMC = cascMC.ptMC(); + const float rapidityMC = RecoDecay::y(std::array{cascMC.pxMC(), cascMC.pyMC(), cascMC.pzMC()}, o2::constants::physics::MassXi0); + + histos.fill(HIST("Gen/hGenXi0Rapidity"), rapidityMC); + + // same rapidity window as at reconstructed level + if (rapidityMC < xi0Selections.minRapidity || rapidityMC > xi0Selections.maxRapidity) { + continue; + } + + auto mcCollision = cascMC.template straMCCollision_as(); + if (!isGeneratedEventAccepted(mcCollision)) { + continue; + } + + float centrality = 100.5f; + const int bestCollisionIdx = listBestCollisionIdx[mcCollision.globalIndex()]; + if (bestCollisionIdx > -1) { + auto collision = collisions.rawIteratorAt(bestCollisionIdx); + centrality = getCentralityRun3(collision); + + if (cascMC.pdgCode() > 0) { + histos.fill(HIST("Gen/h2dGenXi0_RecoedEvt"), centrality, ptMC); + } else { + histos.fill(HIST("Gen/h2dGenAntiXi0_RecoedEvt"), centrality, ptMC); + } + } + + if (cascMC.pdgCode() > 0) { + histos.fill(HIST("Gen/h2dGenXi0"), centrality, ptMC); + histos.fill(HIST("Gen/h2dGenXi0VsMultMC"), mcCollision.multMCNParticlesEta05(), ptMC); + } else { + histos.fill(HIST("Gen/h2dGenAntiXi0"), centrality, ptMC); + histos.fill(HIST("Gen/h2dGenAntiXi0VsMultMC"), mcCollision.multMCNParticlesEta05(), ptMC); + } + } + } + + //_______________________________________________ + // Process functions + void processRealData(soa::Join const& collisions, + Xi0Candidates const& fullXi0s, + V0Candidates const& fullV0s, + DauTracks const&) + { + analyseRecoedXi0s(collisions, fullXi0s, fullV0s); + } + + void processMonteCarlo(soa::Join const& collisions, + Xi0McCandidates const& fullXi0s, + V0Candidates const& fullV0s, + DauTracks const&) + { + analyseRecoedXi0s(collisions, fullXi0s, fullV0s); + } + + void processGenerated(soa::Join const& mcCollisions, + soa::Join const& collisions, + soa::Join const& cascMCCores) + { + analyseGeneratedXi0s(mcCollisions, collisions, cascMCCores); + } + + PROCESS_SWITCH(DerivedXi0Analysis, processRealData, "process as if real data", true); + PROCESS_SWITCH(DerivedXi0Analysis, processMonteCarlo, "process reconstructed information in MC", false); + PROCESS_SWITCH(DerivedXi0Analysis, processGenerated, "process pure generated information in MC", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} From 740c23eb127d814838b66b2887f2037af5e4f053 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 00:08:10 +0200 Subject: [PATCH 05/13] Fix datamodel --- PWGLF/DataModel/LFSigmaTables.h | 36 +++++++++++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index cfb9c5668e0..3806acfcf3a 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -107,6 +107,22 @@ DECLARE_SOA_DYNAMIC_COLUMN(OPAngle, opAngle, return v1.Angle(v2); }); +// Armenteros-Podolanski variables (photon = positive daughter, lambda = negative daughter) +DECLARE_SOA_DYNAMIC_COLUMN(LStarAlpha, lStarAlpha, //! Armenteros Alpha + [](float photonPx, float photonPy, float photonPz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { + float momTot = RecoDecay::p(photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz); + float lQlNeg = RecoDecay::dotProd(std::array{lambdaPx, lambdaPy, lambdaPz}, std::array{photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz}) / momTot; + float lQlPos = RecoDecay::dotProd(std::array{photonPx, photonPy, photonPz}, std::array{photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz}) / momTot; + return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LStarQtArm, lStarQtarm, //! Armenteros Qt + [](float photonPx, float photonPy, float photonPz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { + float momTot = RecoDecay::p2(photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz); + float dp = RecoDecay::dotProd(std::array{lambdaPx, lambdaPy, lambdaPz}, std::array{photonPx + lambdaPx, photonPy + lambdaPy, photonPz + lambdaPz}); + return std::sqrt(RecoDecay::p2(lambdaPx, lambdaPy, lambdaPz) - dp * dp / momTot); // qtarm + }); + // Photon DECLARE_SOA_DYNAMIC_COLUMN(PhotonPt, photonPt, //! Transverse momentum in GeV/c [](float photonPx, float photonPy) -> float { @@ -163,6 +179,8 @@ DECLARE_SOA_TABLE(Sigma0Cores, "AOD", "SIGMA0CORES", // Basic properties sigma0Core::X, sigma0Core::Y, sigma0Core::Z, sigma0Core::DCADaughters, + sigma0Core::LStarAlpha, + sigma0Core::LStarQtArm, sigma0Core::PhotonPx, sigma0Core::PhotonPy, sigma0Core::PhotonPz, sigma0Core::PhotonMass, sigma0Core::LambdaPx, sigma0Core::LambdaPy, sigma0Core::LambdaPz, sigma0Core::LambdaMass, sigma0Core::AntiLambdaMass, @@ -656,6 +674,12 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsSigma0, isSigma0, DECLARE_SOA_DYNAMIC_COLUMN(IsAntiSigma0, isAntiSigma0, //! IsASigma0 [](int pdgCode) -> bool { return pdgCode == PDG_t::kSigma0Bar; }); //-3212 +DECLARE_SOA_DYNAMIC_COLUMN(IsLambdaStar, isLambdaStar, //! IsLambdaStar + [](int pdgCode) -> bool { return pdgCode == 3124; }); // PYTHIA8 code for Lambda(1520) + +DECLARE_SOA_DYNAMIC_COLUMN(IsAntiLambdaStar, isAntiLambdaStar, //! IsAntiLambdaStar + [](int pdgCode) -> bool { return pdgCode == -3124; }); // PYTHIA8 code for AntiLambda(1520) + DECLARE_SOA_DYNAMIC_COLUMN(MCPx, mcpx, //! Sigma0 px [](float photonMCPx, float lambdaMCPx) -> float { return photonMCPx + lambdaMCPx; }); DECLARE_SOA_DYNAMIC_COLUMN(MCPy, mcpy, //! Sigma0 py @@ -686,6 +710,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(Sigma0MCY, sigma0MCY, return RecoDecay::y(std::array{photonMCPx + lambdaMCPx, photonMCPy + lambdaMCPy, photonMCPz + lambdaMCPz}, o2::constants::physics::MassSigma0); }); +DECLARE_SOA_DYNAMIC_COLUMN(LambdaStarMCY, lambdaStarMCY, + [](float photonMCPx, float photonMCPy, float photonMCPz, float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { + return RecoDecay::y(std::array{photonMCPx + lambdaMCPx, photonMCPy + lambdaMCPy, photonMCPz + lambdaMCPz}, o2::constants::physics::MassLambda1520); + }); + DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcphi, //! Phi in the range [0, 2pi) [](float photonMCPx, float photonMCPy, float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(photonMCPx + lambdaMCPx, photonMCPy + lambdaMCPy); }); @@ -767,6 +796,8 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", // Dynamic columns sigma0MCCore::IsSigma0, sigma0MCCore::IsAntiSigma0, + sigma0MCCore::IsLambdaStar, + sigma0MCCore::IsAntiLambdaStar, sigma0MCCore::MCPx, sigma0MCCore::MCPy, @@ -775,6 +806,7 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", sigma0MCCore::MCP, sigma0MCCore::Sigma0MCMass, sigma0MCCore::Sigma0MCY, + sigma0MCCore::LambdaStarMCY, sigma0MCCore::MCPhi, sigma0MCCore::MCEta, sigma0MCCore::MCOPAngle, @@ -1557,7 +1589,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaY, lambdaY, //! Rapidity DECLARE_SOA_DYNAMIC_COLUMN(LambdaPhi, lambdaPhi, //! Phi in the range [0, 2pi) [](float lambdaPx, float lambdaPy) -> float { return RecoDecay::phi(lambdaPx, lambdaPy); }); -} // namespace Pi0Core +} // namespace Xi0Core DECLARE_SOA_TABLE(Xi0Cores, "AOD", "XI0CORES", Xi0Core::X, Xi0Core::Y, Xi0Core::Z, @@ -1823,7 +1855,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCY, lambdaMCY, //! Rapidity DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCPhi, lambdaMCPhi, //! Phi in the range [0, 2pi) [](float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(lambdaMCPx, lambdaMCPy); }); -} // namespace Xi0CoreMC +} // namespace Xi0MCCore DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", // Basic properties From bda492668b45f87fe07eb4b4d6a2cbcd3bf788be Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 28 Sep 2026 22:08:44 +0000 Subject: [PATCH 06/13] Please consider the following formatting changes --- PWGLF/DataModel/LFSigmaTables.h | 42 ++++++------ .../Strangeness/sigma0builder.cxx | 66 +++++++++---------- 2 files changed, 53 insertions(+), 55 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index 3806acfcf3a..f36c89c38c2 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -1395,7 +1395,6 @@ DECLARE_SOA_TABLE(Pi0Gens, "AOD", "PI0GENS", DECLARE_SOA_TABLE(Pi0GenCollRef, "AOD", "PI0GENCOLLREF", //! optional table to refer back to a collision o2::soa::Index<>, v0data::StraMCCollisionId); - // ___________________________________________________________________________ // Xi0 namespace Xi0Core @@ -1450,11 +1449,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Total momentum in GeV/c DECLARE_SOA_DYNAMIC_COLUMN(Xi0Mass, xi0Mass, [](float photon1Px, float photon1Py, float photon1Pz, float photon2Px, float photon2Py, float photon2Pz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { - std::array pVecPi0{photon1Px+photon2Px, photon1Py+photon2Py, photon1Pz+photon2Pz}; + std::array pVecPi0{photon1Px + photon2Px, photon1Py + photon2Py, photon1Pz + photon2Pz}; std::array pVecLambda{lambdaPx, lambdaPy, lambdaPz}; auto arrMom = std::array{pVecPi0, pVecLambda}; return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassLambda0}); - }); + }); DECLARE_SOA_DYNAMIC_COLUMN(PxPi0, pxPi0, //! Pi0 px [](float photon1Px, float photon2Px) -> float { return photon1Px + photon2Px; }); @@ -1488,11 +1487,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(Radius, radius, //! Xi0 decay radius (2D, centered at [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(RadiusPi0, radiusPi0, //! Pi0 decay radius (2D, centered at zero) - [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(RadiusLambda, radiusLambda, //! Pi0 decay radius (2D, centered at zero) - [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); - + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + // CosPAs DECLARE_SOA_DYNAMIC_COLUMN(Pi0CosPA, pi0CosPA, //! [](float xPi0, float yPi0, float zPi0, float pxGamma1, float pyGamma1, float pzGamma1, float pxGamma2, float pyGamma2, float pzGamma2, float pvX, float pvY, float pvZ) -> float { @@ -1511,14 +1510,14 @@ DECLARE_SOA_DYNAMIC_COLUMN(DCALambdaToPV, dcaLambdaToPV, //! [](float X, float Y, float Z, float px, float py, float pz, float pvX, float pvY, float pvZ) -> float { return std::sqrt((std::pow((pvY - Y) * pz - (pvZ - Z) * py, 2) + std::pow((pvX - X) * pz - (pvZ - Z) * px, 2) + std::pow((pvX - X) * py - (pvY - Y) * px, 2)) / (px * px + py * py + pz * pz)); }); - + DECLARE_SOA_DYNAMIC_COLUMN(DCAPi0ToPV, dcaPi0ToPV, //! [](float X, float Y, float Z, float pxGamma1, float pyGamma1, float pzGamma1, float pxGamma2, float pyGamma2, float pzGamma2, float pvX, float pvY, float pvZ) -> float { float px = pxGamma1 + pxGamma2; float py = pyGamma1 + pyGamma2; float pz = pzGamma1 + pzGamma2; return std::sqrt((std::pow((pvY - Y) * pz - (pvZ - Z) * py, 2) + std::pow((pvX - X) * pz - (pvZ - Z) * px, 2) + std::pow((pvX - X) * py - (pvY - Y) * px, 2)) / (px * px + py * py + pz * pz)); - }); + }); DECLARE_SOA_DYNAMIC_COLUMN(Photon1Pt, photon1Pt, //! Transverse momentum in GeV/c [](float photon1Px, float photon1Py) -> float { @@ -1587,12 +1586,12 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaY, lambdaY, //! Rapidity }); DECLARE_SOA_DYNAMIC_COLUMN(LambdaPhi, lambdaPhi, //! Phi in the range [0, 2pi) - [](float lambdaPx, float lambdaPy) -> float { return RecoDecay::phi(lambdaPx, lambdaPy); }); + [](float lambdaPx, float lambdaPy) -> float { return RecoDecay::phi(lambdaPx, lambdaPy); }); } // namespace Xi0Core DECLARE_SOA_TABLE(Xi0Cores, "AOD", "XI0CORES", - Xi0Core::X, Xi0Core::Y, Xi0Core::Z, + Xi0Core::X, Xi0Core::Y, Xi0Core::Z, Xi0Core::XPi0, Xi0Core::YPi0, Xi0Core::ZPi0, Xi0Core::XLambda, Xi0Core::YLambda, Xi0Core::ZLambda, Xi0Core::DCADaughters, Xi0Core::DCADaughtersPi0, Xi0Core::DCAXYCascToPV, Xi0Core::DCAZCascToPV, @@ -1613,7 +1612,7 @@ DECLARE_SOA_TABLE(Xi0Cores, "AOD", "XI0CORES", Xi0Core::Phi, Xi0Core::Eta, Xi0Core::Radius, - Xi0Core::CascCosPA, + Xi0Core::CascCosPA, // Dynamic columns for pi0 Xi0Core::PxPi0, @@ -1652,15 +1651,14 @@ DECLARE_SOA_TABLE(Xi0CollRefs, "AOD", "XI0COLLREFS", //! optional table to refer o2::soa::Index<>, v0data::StraCollisionId); DECLARE_SOA_TABLE(Xi0Indices, "AOD", "XI0INDICES", //! optional table to refer back to V0Cores - o2::soa::Index<>, Xi0Core::Photon1Index, Xi0Core::Photon2Index, Xi0Core::LambdaIndex); - - + o2::soa::Index<>, Xi0Core::Photon1Index, Xi0Core::Photon2Index, Xi0Core::LambdaIndex); + // for MC namespace Xi0MCCore { DECLARE_SOA_COLUMN(Photon1Index, photon1Index, int); DECLARE_SOA_COLUMN(Photon2Index, photon2Index, int); -DECLARE_SOA_COLUMN(LambdaIndex, lambdaIndex, int); +DECLARE_SOA_COLUMN(LambdaIndex, lambdaIndex, int); DECLARE_SOA_COLUMN(Photon1MCPx, photon1mcpx, float); DECLARE_SOA_COLUMN(Photon1MCPy, photon1mcpy, float); @@ -1733,13 +1731,13 @@ DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcPhi, //! Phi in the range [0, 2pi) DECLARE_SOA_DYNAMIC_COLUMN(MCEta, mcEta, //! Pseudorapidity [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz, float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { return RecoDecay::eta(std::array{photon1MCPx + photon2MCPx + lambdaMCPx, photon1MCPy + photon2MCPy + lambdaMCPy, photon1MCPz + photon2MCPz + lambdaMCPz}); - }); - + }); + DECLARE_SOA_DYNAMIC_COLUMN(MCRadius, mcRadius, //! Xi0 decay radius (2D, centered at zero) [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(MCRadiusPi0, mcRadiusPi0, //! Pi0 decay radius (2D, centered at zero) - [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPx, pi0mcpx, //! Pi0 MC px [](float photon1MCPx, float photon2MCPx) -> float { return photon1MCPx + photon2MCPx; }); @@ -1853,7 +1851,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCY, lambdaMCY, //! Rapidity }); DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCPhi, lambdaMCPhi, //! Phi in the range [0, 2pi) - [](float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(lambdaMCPx, lambdaMCPy); }); + [](float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(lambdaMCPx, lambdaMCPy); }); } // namespace Xi0MCCore @@ -1872,7 +1870,7 @@ DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", Xi0MCCore::IsPi0Primary, Xi0MCCore::Pi0MCprocess, Xi0MCCore::Pi0PDGCode, Xi0MCCore::MCx, Xi0MCCore::MCy, Xi0MCCore::MCz, - Xi0MCCore::IsPhysicalPrimary, Xi0MCCore::MCprocess, Xi0MCCore::PDGCode, Xi0MCCore::PDGCodeMother, + Xi0MCCore::IsPhysicalPrimary, Xi0MCCore::MCprocess, Xi0MCCore::PDGCode, Xi0MCCore::PDGCodeMother, // Dynamic columns Xi0MCCore::MCPx, @@ -1909,7 +1907,7 @@ DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", Xi0MCCore::Photon2MCEta, Xi0MCCore::Photon2MCY, Xi0MCCore::Photon2MCPhi, - + Xi0MCCore::LambdaMCPt, Xi0MCCore::LambdaMCP, Xi0MCCore::LambdaMCEta, @@ -1920,7 +1918,7 @@ DECLARE_SOA_TABLE(Xi0MCCollRefs, "AOD", "XI0MCCOLLREFS", //! optional table to r o2::soa::Index<>, v0data::StraMCCollisionId); DECLARE_SOA_TABLE(Xi0MCIndices, "AOD", "XI0MCINDICES", //! optional table to refer back to V0MCCores - o2::soa::Index<>, Xi0MCCore::Photon1Index, Xi0MCCore::Photon2Index, Xi0MCCore::LambdaIndex); + o2::soa::Index<>, Xi0MCCore::Photon1Index, Xi0MCCore::Photon2Index, Xi0MCCore::LambdaIndex); } // namespace o2::aod diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index bdac95866ef..c5c65b528e7 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -21,16 +21,16 @@ #include "PWGLF/DataModel/LFStrangenessMLTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" -#include "Common/Tools/StandardCCDBLoader.h" -#include #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" +#include "Common/Tools/StandardCCDBLoader.h" #include #include #include +#include #include #include #include @@ -113,7 +113,7 @@ struct sigma0builder { //__________________________________________________ // Xi0 specific Produces xi0cores; // xi0 candidates info for analysis - Produces xi0collRefs; // references to straCollision + Produces xi0collRefs; // references to straCollision Produces xi0Indices; // references to V0Cores Produces xi0mccores; // Reco xi0 MC properties Produces xi0mccollRefs; // references to straMCCollision @@ -892,12 +892,12 @@ struct sigma0builder { // ______________________________________________________ // Struct to store V0Pair properties struct V0PairTopoInfo { - std::array position {-999.f, -999.f, -999.f}; + std::array position{-999.f, -999.f, -999.f}; float dcaPi0ToPV = 999.f; float daughterDCA = 999.f; float cosPA = -1.f; - float pi0Mass= 999.f; - float pi0Y= 999.f; + float pi0Mass = 999.f; + float pi0Y = 999.f; int v01Index = 0; // index to de-reference V0Cores table int v02Index = 0; // index to de-reference V0Cores table @@ -964,31 +964,31 @@ struct sigma0builder { // ______________________________________________________ // Struct to store V0Pair properties struct Xi0Info { - std::array gamma1Momentum {999.f, 999.f, 999.f}; - std::array gamma2Momentum {999.f, 999.f, 999.f}; - std::array lambdaMomentum {999.f, 999.f, 999.f}; - std::array pi0Position {0., 0., 0.}; - std::array cascadePosition {0., 0., 0.}; + std::array gamma1Momentum{999.f, 999.f, 999.f}; + std::array gamma2Momentum{999.f, 999.f, 999.f}; + std::array lambdaMomentum{999.f, 999.f, 999.f}; + std::array pi0Position{0., 0., 0.}; + std::array cascadePosition{0., 0., 0.}; float xi0Mass = 999.f; float dcaPi0ToPV = 999.f; float pi0DaughterDCA = 999.f; float cascadeDaughterDCA = 999.f; float cascadeDCAxy = 999.f; float cascadeDCAz = 999.f; - + int collisionId = 0; int gamma1Index = 0; // index to de-reference V0Cores table int gamma2Index = 0; // index to de-reference V0Cores table int lambdaIndex = 0; // index to de-reference V0Cores table - std::array covariance {0.}; + std::array covariance{0.}; }; struct Xi0MCInfo { int gamma1Index = 0; // index to de-reference V0MCCores table int gamma2Index = 0; // index to de-reference V0MCCores table int lambdaIndex = 0; // index to de-reference V0MCCores table - + // Basic kinematic info float gamma1MCpx = -999.f; float gamma1MCpy = -999.f; @@ -1472,8 +1472,8 @@ struct sigma0builder { auto MCParticle_gamma2 = mcparticles.rawIteratorAt(gamma2MC.particleIdMC()); auto MCParticle_lambda = mcparticles.rawIteratorAt(lambdaMC.particleIdMC()); - if ((std::abs(MCParticle_gamma1.pdgCode()) == PDG_t::kGamma) && - (std::abs(MCParticle_gamma2.pdgCode()) == PDG_t::kGamma) && + if ((std::abs(MCParticle_gamma1.pdgCode()) == PDG_t::kGamma) && + (std::abs(MCParticle_gamma2.pdgCode()) == PDG_t::kGamma) && (std::abs(MCParticle_lambda.pdgCode()) == PDG_t::kLambda0)) { // Get MC Mothers @@ -1482,9 +1482,9 @@ struct sigma0builder { auto const& MCMothersList_lambda = MCParticle_lambda.template mothers_as(); if (!MCMothersList_gamma1.empty() && !MCMothersList_gamma2.empty() && !MCMothersList_lambda.empty()) { // Are there mothers? - auto const& MCMother_gamma1 = MCMothersList_gamma1.front(); // First mother - auto const& MCMother_gamma2 = MCMothersList_gamma2.front(); // First mother - auto const& MCMother_lambda = MCMothersList_lambda.front(); // First mother + auto const& MCMother_gamma1 = MCMothersList_gamma1.front(); // First mother + auto const& MCMother_gamma2 = MCMothersList_gamma2.front(); // First mother + auto const& MCMother_lambda = MCMothersList_lambda.front(); // First mother MCinfo.gamma1PDGCodeMother = MCMother_gamma1.pdgCode(); MCinfo.gamma2PDGCodeMother = MCMother_gamma2.pdgCode(); @@ -1500,9 +1500,9 @@ struct sigma0builder { MCinfo.pi0Z = MCMother_gamma1.vz(); auto const& MCMothersList_pi0 = MCMother_gamma1.template mothers_as(); // get pi0 mother list - if (!MCMothersList_pi0.empty()) { // Are there mothers? - auto const& MCMother_pi0 = MCMothersList_pi0.front(); // get pi0 mother - + if (!MCMothersList_pi0.empty()) { // Are there mothers? + auto const& MCMother_pi0 = MCMothersList_pi0.front(); // get pi0 mother + if (MCMother_pi0.globalIndex() == MCMother_lambda.globalIndex()) { // check that lambda and pi0 have the same mother. MCinfo.xi0PDGCode = MCMother_pi0.pdgCode(); MCinfo.xi0MCProcess = MCMother_pi0.getProcess(); @@ -2601,7 +2601,7 @@ struct sigma0builder { //_______________________________________________ // Build pi0 candidate for QA template - bool buildPi0ForXi0(TV0Object const& gamma1, TV0Object const& gamma2, TCollision const& collision, V0PairTopoInfo &info) + bool buildPi0ForXi0(TV0Object const& gamma1, TV0Object const& gamma2, TCollision const& collision, V0PairTopoInfo& info) { //_______________________________________________ // Check if both V0s are made of the same tracks @@ -2629,7 +2629,7 @@ struct sigma0builder { info.v01Index = gamma1.globalIndex(); info.v02Index = gamma2.globalIndex(); - + //_______________________________________________ // Pi0-specific selections: if (std::abs(info.pi0Y) > cascadeSelections.maxPi0Rapidity) @@ -2675,7 +2675,7 @@ struct sigma0builder { if (std::hypot(info.position[0], info.position[1]) < cascadeSelections.radiusPi0) { return false; } - + float x = vtxPi0[0]; float y = vtxPi0[1]; float z = vtxPi0[2]; @@ -2695,7 +2695,7 @@ struct sigma0builder { //_______________________________________________ // Build Xi0 candidate for analysis template - bool buildXi0(TV0Object const& gamma1, TV0Object const& gamma2, TV0Object const& lambda, V0PairTopoInfo &pi0Info, TCollision const& collision, TMCParticles const& mcparticles) + bool buildXi0(TV0Object const& gamma1, TV0Object const& gamma2, TV0Object const& lambda, V0PairTopoInfo& pi0Info, TCollision const& collision, TMCParticles const& mcparticles) { Xi0Info cascade; //_______________________________________________ @@ -2786,9 +2786,9 @@ struct sigma0builder { // Calculate DCAxy of the cascade (with bending) auto lCascadeTrack = fitter3Prongs.createParentTrackParCov(); - lCascadeTrack.setAbsCharge(0); // to be sure + lCascadeTrack.setAbsCharge(0); // to be sure lCascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas - std::array dcaInfo {999.f, 999.f}; + std::array dcaInfo{999.f, 999.f}; o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, fitter3Prongs.getMatCorrType(), &dcaInfo); cascade.cascadeDCAxy = dcaInfo[0]; @@ -2801,10 +2801,10 @@ struct sigma0builder { cascade.gamma2Index = gamma2.globalIndex(); cascade.lambdaIndex = lambda.globalIndex(); - std::array momPi0{gamma1.px()+gamma2.px(), gamma1.py()+gamma2.py(), gamma1.pz()+gamma2.pz()}; + std::array momPi0{gamma1.px() + gamma2.px(), gamma1.py() + gamma2.py(), gamma1.pz() + gamma2.pz()}; auto arrMom = std::array{momPi0, momLambda}; cascade.xi0Mass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassXi0}); - + // Calculate position covariance matrix auto covVtxV = fitter3Prongs.calcPCACovMatrix(0); // std::array positionCovariance; @@ -2833,7 +2833,7 @@ struct sigma0builder { } // Check if MC data and populate corresponding table - if constexpr (requires { collision.straMCCollisionId(); gamma1.motherMCPartId(); gamma2.motherMCPartId(); lambda.motherMCPartId();}) { + if constexpr (requires { collision.straMCCollisionId(); gamma1.motherMCPartId(); gamma2.motherMCPartId(); lambda.motherMCPartId(); }) { auto xi0MCInfo = getXi0MCInfo(gamma1, gamma2, lambda, collision, mcparticles); xi0mccores(xi0MCInfo.gamma1MCpx, xi0MCInfo.gamma1MCpy, xi0MCInfo.gamma1MCpz, @@ -2844,7 +2844,7 @@ struct sigma0builder { xi0MCInfo.lambdaIsPhysicalPrimary, xi0MCInfo.lambdaPDGCodePos, xi0MCInfo.lambdaPDGCodeNeg, xi0MCInfo.lambdaPDGCode, xi0MCInfo.lambdaPDGCodeMother, xi0MCInfo.pi0X, xi0MCInfo.pi0Y, xi0MCInfo.pi0Z, xi0MCInfo.pi0IsPhysicalPrimary, xi0MCInfo.pi0MCProcess, xi0MCInfo.pi0PDGCode, xi0MCInfo.xi0X, xi0MCInfo.xi0Y, xi0MCInfo.xi0Z, xi0MCInfo.xi0IsPhysicalPrimary, xi0MCInfo.xi0MCProcess, xi0MCInfo.xi0PDGCode, xi0MCInfo.xi0PDGCodeMother); - + xi0mcIndices(xi0MCInfo.gamma1Index, xi0MCInfo.gamma2Index, xi0MCInfo.lambdaIndex); int mcCollisionIndex = -1; @@ -3408,7 +3408,7 @@ struct sigma0builder { pi0Candidates.push_back(pi0Info); } } - + // second build Xi0 candidates for (size_t i = 0; i < bestLambdasArray.size(); ++i) { // loop over Lambda auto lambda = fullV0s.rawIteratorAt(bestLambdasArray[i]); From 178e4e041521eb3eb4e868b82972f81d69a64907 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 00:14:02 +0200 Subject: [PATCH 07/13] Fix sigma0builder --- .../Strangeness/sigma0builder.cxx | 31 ++++++++++++------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index bdac95866ef..d7710bb15e3 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -42,6 +42,9 @@ #include #include #include +#include +#include +#include #include // IWYU pragma: keep (do not replace with Math/Vector3Dfwd.h) #include @@ -1094,8 +1097,9 @@ struct sigma0builder { auto v02MC = v02.template v0MCCore_as>(); // Sanity check: Is V0Pair <-> Mother assignment correct? + int expectedPairPDG = doLambdaStar ? 3124 : PDG_t::kSigma0; bool fIsSigma0 = false; - if ((v01MC.pdgCode() == PDG_t::kGamma) && (v01MC.pdgCodeMother() == PDG_t::kSigma0) && (v02MC.pdgCode() == PDG_t::kLambda0) && (v02MC.pdgCodeMother() == PDG_t::kSigma0) && (v01.motherMCPartId() == v02.motherMCPartId())) + if ((v01MC.pdgCode() == PDG_t::kGamma) && (v01MC.pdgCodeMother() == expectedPairPDG) && (v02MC.pdgCode() == PDG_t::kLambda0) && (v02MC.pdgCodeMother() == expectedPairPDG) && (v01.motherMCPartId() == v02.motherMCPartId())) fIsSigma0 = true; bool fIsKStar = false; @@ -1246,11 +1250,11 @@ struct sigma0builder { } } // Check association correctness - if (fIsSigma0 && (MCinfo.V0PairPDGCode == PDG_t::kSigma0)) + if (fIsSigma0 && (MCinfo.V0PairPDGCode == expectedPairPDG)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 1); // match - if (fIsSigma0 && !(MCinfo.V0PairPDGCode == PDG_t::kSigma0)) + if (fIsSigma0 && !(MCinfo.V0PairPDGCode == expectedPairPDG)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 2); // mismatch - if (!fIsSigma0 && (MCinfo.V0PairPDGCode == PDG_t::kSigma0)) + if (!fIsSigma0 && (MCinfo.V0PairPDGCode == expectedPairPDG)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 3); // mismatch // Check association correctness @@ -1905,6 +1909,8 @@ struct sigma0builder { GenInfo.IsPi0 = mcParticle.pdgCode() == PDG_t::kPi0; // 111; GenInfo.IsSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0; // PDG_t::kSigma0 GenInfo.IsAntiSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0Bar; //-3212 + GenInfo.IsLambdaStar = mcParticle.pdgCode() == 3124; // 102134 (PYTHIA8) + GenInfo.IsAntiLambdaStar = mcParticle.pdgCode() == -3124; // -102134 GenInfo.IsKStar = std::abs(mcParticle.pdgCode()) == o2::constants::physics::Pdg::kK0Star892; // 313; GenInfo.IsProducedByGenerator = mcParticle.producedByGenerator(); GenInfo.MCProcess = mcParticle.getProcess(); @@ -1916,7 +1922,7 @@ struct sigma0builder { GenInfo.MCCollId = mcParticle.mcCollisionId(); // save this reference, please // Checking decay mode if sigma0 or pi0 (it is easier here) - if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0 || GenInfo.IsPi0 || GenInfo.IsKStar) { + if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0 || GenInfo.IsPi0 || GenInfo.IsKStar || GenInfo.IsLambdaStar || GenInfo.IsAntiLambdaStar) { // This is a costly operation, so we do it only for pi0s and sigma0s auto const& daughters = mcParticle.template daughters_as(); @@ -1926,7 +1932,7 @@ struct sigma0builder { auto const& GenMothersList = mcParticle.template mothers_as(); GenInfo.PDGCodeMother = (!GenMothersList.empty()) ? GenMothersList.front().pdgCode() : 0; - if ((GenInfo.IsSigma0 || GenInfo.IsAntiSigma0) && genSelections.doQA) { + if ((doLambdaStar ? (GenInfo.IsLambdaStar || GenInfo.IsAntiLambdaStar) : (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0)) && genSelections.doQA) { histos.fill(HIST("GenQA/h2dSigma0MCSourceVsPDGMother"), GenInfo.IsProducedByGenerator, GenInfo.PDGCodeMother); // Checking decay modes and getting daughter pTs @@ -1989,7 +1995,7 @@ struct sigma0builder { histos.fill(HIST("GenQA/hGenSpeciesKStar"), 0); // Checking decay mode - if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0) { + if (doLambdaStar ? (GenInfo.IsLambdaStar || GenInfo.IsAntiLambdaStar) : (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0)) { histos.fill(HIST("GenQA/hSigma0NDau"), GenInfo.NDaughters); histos.fill(HIST("GenQA/h2dSigma0NDauVsProcess"), GenInfo.NDaughters, GenInfo.MCProcess); @@ -2009,7 +2015,7 @@ struct sigma0builder { histos.fill(HIST("GenQA/h2DGenSigma0TypeVsProducedByGen"), typeIndex, genIndex); // Fill histograms - if (GenInfo.IsSigma0) { + if (doLambdaStar ? GenInfo.IsLambdaStar : GenInfo.IsSigma0) { histos.fill(HIST("GenQA/hGenSpecies"), 2); histos.fill(HIST("GenQA/hGenSigma0"), GenInfo.MCPt); histos.fill(HIST("GenQA/h3dGenSigma0_pTMap"), GenInfo.MCPt, GenInfo.MCDau1Pt, GenInfo.MCDau2Pt); @@ -2018,7 +2024,7 @@ struct sigma0builder { if (GenInfo.IsPrimary) histos.fill(HIST("GenQA/hPrimarySigma0s"), 1); } - if (GenInfo.IsAntiSigma0) { + if (doLambdaStar ? GenInfo.IsAntiLambdaStar : GenInfo.IsAntiSigma0) { histos.fill(HIST("GenQA/hGenSpecies"), 3); histos.fill(HIST("GenQA/hGenAntiSigma0"), GenInfo.MCPt); histos.fill(HIST("GenQA/h3dGenASigma0_pTMap"), GenInfo.MCPt, GenInfo.MCDau1Pt, GenInfo.MCDau2Pt); @@ -2095,9 +2101,10 @@ struct sigma0builder { pi0GenCollRefs(MCGenInfo.MCCollId); // link to stramccollision table } - // Sigma0/ASigma0 - if (fillSigma0Tables && (MCGenInfo.IsSigma0 || MCGenInfo.IsAntiSigma0)) { - sigma0Gens(MCGenInfo.IsSigma0, MCGenInfo.IsProducedByGenerator, MCGenInfo.MCPt, mcParticle.y()); + // Sigma0/ASigma0 (Lambda(1520)/ALambda(1520)) + bool fIsGenSigma0Like = doLambdaStar ? (MCGenInfo.IsLambdaStar || MCGenInfo.IsAntiLambdaStar) : (MCGenInfo.IsSigma0 || MCGenInfo.IsAntiSigma0); + if (fillSigma0Tables && fIsGenSigma0Like) { + sigma0Gens(doLambdaStar ? MCGenInfo.IsLambdaStar : MCGenInfo.IsSigma0, MCGenInfo.IsProducedByGenerator, MCGenInfo.MCPt, mcParticle.y()); sigma0GenCollRefs(MCGenInfo.MCCollId); // link to stramccollision table } From 111a5a5b680f81963faa37b4c56241232cfc5251 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 28 Sep 2026 22:15:06 +0000 Subject: [PATCH 08/13] Please consider the following formatting changes --- PWGLF/TableProducer/Strangeness/sigma0builder.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index 8668fe1abcd..9feca132488 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -31,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -42,9 +43,8 @@ #include #include #include -#include #include -#include +#include #include // IWYU pragma: keep (do not replace with Math/Vector3Dfwd.h) #include From b6834b717cb197a7e63978ab9a66f0485921ff78 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 00:16:24 +0200 Subject: [PATCH 09/13] Remove and re-add LStarAlpha and LStarQtArm properties --- PWGLF/DataModel/LFSigmaTables.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index f36c89c38c2..4428b3b2dff 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -179,8 +179,6 @@ DECLARE_SOA_TABLE(Sigma0Cores, "AOD", "SIGMA0CORES", // Basic properties sigma0Core::X, sigma0Core::Y, sigma0Core::Z, sigma0Core::DCADaughters, - sigma0Core::LStarAlpha, - sigma0Core::LStarQtArm, sigma0Core::PhotonPx, sigma0Core::PhotonPy, sigma0Core::PhotonPz, sigma0Core::PhotonMass, sigma0Core::LambdaPx, sigma0Core::LambdaPy, sigma0Core::LambdaPz, sigma0Core::LambdaMass, sigma0Core::AntiLambdaMass, @@ -197,6 +195,8 @@ DECLARE_SOA_TABLE(Sigma0Cores, "AOD", "SIGMA0CORES", sigma0Core::Eta, sigma0Core::Radius, sigma0Core::OPAngle, + sigma0Core::LStarAlpha, + sigma0Core::LStarQtArm, sigma0Core::PhotonPt, sigma0Core::PhotonP, From e65ac6958a5578ba422a7f6ec5403a2ed1b82434 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 00:17:47 +0200 Subject: [PATCH 10/13] Add flags for LambdaStar and AntiLambdaStar --- PWGLF/TableProducer/Strangeness/sigma0builder.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index 9feca132488..96e80da2bb7 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -951,6 +951,8 @@ struct sigma0builder { bool IsSigma0 = false; bool IsAntiSigma0 = false; bool IsKStar = false; + bool IsLambdaStar = false; + bool IsAntiLambdaStar = false; bool IsProducedByGenerator = false; bool IsSterile = false; int MCProcess = -1; From a616bea58705a701b60d85a55c1bdf30300a5332 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 08:14:31 +0200 Subject: [PATCH 11/13] Remove unused variable --- PWGLF/TableProducer/Strangeness/sigma0builder.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index 96e80da2bb7..c9df9aa31d5 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -1432,8 +1432,8 @@ struct sigma0builder { return MCinfo; } - template - Xi0MCInfo getXi0MCInfo(TV0 const& gamma1, TV0 const& gamma2, TV0 const& lambda, TCollision const& collision, TMCParticles const& mcparticles) + template + Xi0MCInfo getXi0MCInfo(TV0 const& gamma1, TV0 const& gamma2, TV0 const& lambda, TMCParticles const& mcparticles) { Xi0MCInfo MCinfo; @@ -2843,7 +2843,7 @@ struct sigma0builder { // Check if MC data and populate corresponding table if constexpr (requires { collision.straMCCollisionId(); gamma1.motherMCPartId(); gamma2.motherMCPartId(); lambda.motherMCPartId(); }) { - auto xi0MCInfo = getXi0MCInfo(gamma1, gamma2, lambda, collision, mcparticles); + auto xi0MCInfo = getXi0MCInfo(gamma1, gamma2, lambda, mcparticles); xi0mccores(xi0MCInfo.gamma1MCpx, xi0MCInfo.gamma1MCpy, xi0MCInfo.gamma1MCpz, xi0MCInfo.gamma1IsPhysicalPrimary, xi0MCInfo.gamma1PDGCodePos, xi0MCInfo.gamma1PDGCodeNeg, xi0MCInfo.gamma1PDGCode, xi0MCInfo.gamma1PDGCodeMother, From d3430a090a2f56c1ed82956276fc44948bb0595e Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 08:15:34 +0200 Subject: [PATCH 12/13] Fix typo in CMakeLists.txt --- PWGLF/Tasks/Strangeness/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index 9ddda410938..7a757d05e88 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -70,7 +70,7 @@ o2physics_add_dpl_workflow(derivedcascadeanalysis COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(derivedxi0analysis - SOURCES derivedXi0Analysis.cxx + SOURCES derivedXi0analysis.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB COMPONENT_NAME Analysis) From 10d1e18e1011cd7d1dd4750821af11675dd204d9 Mon Sep 17 00:00:00 2001 From: SCHOTTER Romain <47983209+romainschotter@users.noreply.github.com> Date: Tue, 29 Sep 2026 08:27:37 +0200 Subject: [PATCH 13/13] Fix Xi0 mass --- PWGLF/TableProducer/Strangeness/sigma0builder.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index c9df9aa31d5..cb51f93b6d3 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -2812,7 +2812,7 @@ struct sigma0builder { std::array momPi0{gamma1.px() + gamma2.px(), gamma1.py() + gamma2.py(), gamma1.pz() + gamma2.pz()}; auto arrMom = std::array{momPi0, momLambda}; - cascade.xi0Mass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassXi0}); + cascade.xi0Mass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassLambda0}); // Calculate position covariance matrix auto covVtxV = fitter3Prongs.calcPCACovMatrix(0);