diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index e11567b2d7a..de5d0f46acd 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -1395,6 +1395,530 @@ 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 +{ +DECLARE_SOA_COLUMN(Photon1Index, photon1Index, int); +DECLARE_SOA_COLUMN(Photon2Index, photon2Index, int); +DECLARE_SOA_COLUMN(LambdaIndex, lambdaIndex, int); +DECLARE_SOA_COLUMN(X, x, float); +DECLARE_SOA_COLUMN(Y, y, float); +DECLARE_SOA_COLUMN(Z, z, float); +DECLARE_SOA_COLUMN(XPi0, xPi0, float); +DECLARE_SOA_COLUMN(YPi0, yPi0, float); +DECLARE_SOA_COLUMN(ZPi0, zPi0, float); +DECLARE_SOA_COLUMN(XLambda, xLambda, float); +DECLARE_SOA_COLUMN(YLambda, yLambda, float); +DECLARE_SOA_COLUMN(ZLambda, zLambda, float); +DECLARE_SOA_COLUMN(DCADaughters, dcadaughters, float); +DECLARE_SOA_COLUMN(DCADaughtersPi0, dcadaughtersPi0, float); +DECLARE_SOA_COLUMN(DCAXYCascToPV, dcaXYCascToPV, float); +DECLARE_SOA_COLUMN(DCAZCascToPV, dcaZCascToPV, float); + +DECLARE_SOA_COLUMN(Photon1Px, photon1Px, float); +DECLARE_SOA_COLUMN(Photon1Py, photon1Py, float); +DECLARE_SOA_COLUMN(Photon1Pz, photon1Pz, float); + +DECLARE_SOA_COLUMN(Photon2Px, photon2Px, float); +DECLARE_SOA_COLUMN(Photon2Py, photon2Py, float); +DECLARE_SOA_COLUMN(Photon2Pz, photon2Pz, float); + +DECLARE_SOA_COLUMN(LambdaPx, lambdaPx, float); +DECLARE_SOA_COLUMN(LambdaPy, lambdaPy, float); +DECLARE_SOA_COLUMN(LambdaPz, lambdaPz, float); + +//______________________________________________________ +// DYNAMIC COLUMNS +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, //! Pi0 px + [](float photon1Px, float photon2Px, float lambdaPx) -> float { return photon1Px + photon2Px + lambdaPx; }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, //! Pi0 py + [](float photon1Py, float photon2Py, float lambdaPy) -> float { return photon1Py + photon2Py + lambdaPy; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! Pi0 pz + [](float photon1Pz, float photon2Pz, float lambdaPz) -> float { return photon1Pz + photon2Pz + lambdaPz; }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, + [](float photon1Px, float photon1Py, float photon2Px, float photon2Py, float lambdaPx, float lambdaPy) -> float { + return RecoDecay::pt(std::array{photon1Px + photon2Px + lambdaPx, photon1Py + photon2Py + lambdaPy}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Total momentum in GeV/c + [](float photon1Px, float photon1Py, float photon1Pz, float photon2Px, float photon2Py, float photon2Pz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { + return RecoDecay::sqrtSumOfSquares(photon1Px + photon2Px + lambdaPx, photon1Py + photon2Py + lambdaPy, photon1Pz + photon2Pz + lambdaPz); + }); + +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 pVecLambda{lambdaPx, lambdaPy, lambdaPz}; + auto arrMom = std::array{pVecPi0, pVecLambda}; + return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPi0, o2::constants::physics::MassXi0}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PxPi0, pxPi0, //! Pi0 px + [](float photon1Px, float photon2Px) -> float { return photon1Px + photon2Px; }); +DECLARE_SOA_DYNAMIC_COLUMN(PyPi0, pyPi0, //! Pi0 py + [](float photon1Py, float photon2Py) -> float { return photon1Py + photon2Py; }); +DECLARE_SOA_DYNAMIC_COLUMN(PzPi0, pzPi0, //! Pi0 pz + [](float photon1Pz, float photon2Pz) -> float { return photon1Pz + photon2Pz; }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0Mass, pi0Mass, + [](float photon1Px, float photon1Py, float photon1Pz, float photon2Px, float photon2Py, float photon2Pz) -> float { + std::array pVecPhoton1{photon1Px, photon1Py, photon1Pz}; + std::array pVecPhoton2{photon2Px, photon2Py, photon2Pz}; + auto arrMom = std::array{pVecPhoton1, pVecPhoton2}; + return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassPhoton}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Rapidity, rapidity, + [](float photon1Px, float photon1Py, float photon1Pz, float photon2Px, float photon2Py, float photon2Pz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { + return RecoDecay::y(std::array{photon1Px + photon2Px + lambdaPx, photon1Py + photon2Py + lambdaPy, photon1Pz + photon2Pz + lambdaPz}, o2::constants::physics::MassXi0); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, //! Phi in the range [0, 2pi) + [](float photon1Px, float photon1Py, float photon2Px, float photon2Py, float lambdaPx, float lambdaPy) -> float { return RecoDecay::phi(photon1Px + photon2Px + lambdaPx, photon1Py + photon2Py + lambdaPy); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, //! Pseudorapidity + [](float photon1Px, float photon1Py, float photon1Pz, float photon2Px, float photon2Py, float photon2Pz, float lambdaPx, float lambdaPy, float lambdaPz) -> float { + return RecoDecay::eta(std::array{photon1Px + photon2Px + lambdaPx, photon1Py + photon2Py + lambdaPy, photon1Pz + photon2Pz + lambdaPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Radius, radius, //! Xi0 decay radius (2D, centered at zero) + [](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); }); + +DECLARE_SOA_DYNAMIC_COLUMN(RadiusLambda, radiusLambda, //! Pi0 decay radius (2D, centered at zero) + [](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 { + return RecoDecay::cpa(std::array{pvX, pvY, pvZ}, std::array{xPi0, yPi0, zPi0}, std::array{pxGamma1 + pxGamma2, pyGamma1 + pyGamma2, pzGamma1 + pzGamma2}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaCosPA, lambdaCosPA, //! + [](float xLambda, float yLambda, float zLambda, float pxLambda, float pyLambda, float pzLambda, float pvX, float pvY, float pvZ) -> float { + return RecoDecay::cpa(std::array{pvX, pvY, pvZ}, std::array{xLambda, yLambda, zLambda}, std::array{pxLambda, pyLambda, pzLambda}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(CascCosPA, cascCosPA, //! + [](float X, float Y, float Z, float PxGamma1, float PxGamma2, float PxLambda, float PyGamma1, float PyGamma2, float PyLambda, float PzGamma1, float PzGamma2, float PzLambda, float pvX, float pvY, float pvZ) -> float { return RecoDecay::cpa(std::array{pvX, pvY, pvZ}, std::array{X, Y, Z}, std::array{PxGamma1 + PxGamma2 + PxLambda, PyGamma1 + PyGamma2 + PyLambda, PzGamma1 + PzGamma2 + PzLambda}); }); + +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 { + return RecoDecay::sqrtSumOfSquares(photon1Px, photon1Py); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1P, photon1p, //! Total momentum in GeV/c + [](float photon1Px, float photon1Py, float photon1Pz) -> float { + return RecoDecay::sqrtSumOfSquares(photon1Px, photon1Py, photon1Pz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1Eta, photon1Eta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float photon1Px, float photon1Py, float photon1Pz) -> float { + return RecoDecay::eta(std::array{photon1Px, photon1Py, photon1Pz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1Y, photon1Y, //! Rapidity + [](float photon1Px, float photon1Py, float photon1Pz) -> float { + return RecoDecay::y(std::array{photon1Px, photon1Py, photon1Pz}, o2::constants::physics::MassGamma); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1Phi, photon1Phi, //! Phi in the range [0, 2pi) + [](float photon1Px, float photon1Py) -> float { return RecoDecay::phi(photon1Px, photon1Py); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2Pt, photon2Pt, //! Transverse momentum in GeV/c + [](float photon2Px, float photon2Py) -> float { + return RecoDecay::sqrtSumOfSquares(photon2Px, photon2Py); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2P, photon2p, //! Total momentum in GeV/c + [](float photon2Px, float photon2Py, float photon2Pz) -> float { + return RecoDecay::sqrtSumOfSquares(photon2Px, photon2Py, photon2Pz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2Eta, photon2Eta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float photon2Px, float photon2Py, float photon2Pz) -> float { + return RecoDecay::eta(std::array{photon2Px, photon2Py, photon2Pz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2Y, photon2Y, //! Rapidity + [](float photon2Px, float photon2Py, float photon2Pz) -> float { + return RecoDecay::y(std::array{photon2Px, photon2Py, photon2Pz}, o2::constants::physics::MassGamma); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2Phi, photon2Phi, //! Phi in the range [0, 2pi) + [](float photon2Px, float photon2Py) -> float { return RecoDecay::phi(photon2Px, photon2Py); }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaPt, lambdaPt, //! Transverse momentum in GeV/c + [](float lambdaPx, float lambdaPy) -> float { + return RecoDecay::sqrtSumOfSquares(lambdaPx, lambdaPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaP, lambdaP, //! Total momentum in GeV/c + [](float lambdaPx, float lambdaPy, float lambdaPz) -> float { + return RecoDecay::sqrtSumOfSquares(lambdaPx, lambdaPy, lambdaPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaEta, lambdaEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float lambdaPx, float lambdaPy, float lambdaPz) -> float { + return RecoDecay::eta(std::array{lambdaPx, lambdaPy, lambdaPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaY, lambdaY, //! Rapidity + [](float lambdaPx, float lambdaPy, float lambdaPz) -> float { + return RecoDecay::y(std::array{lambdaPx, lambdaPy, lambdaPz}, o2::constants::physics::MassLambda); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaPhi, lambdaPhi, //! Phi in the range [0, 2pi) + [](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::XPi0, Xi0Core::YPi0, Xi0Core::ZPi0, + Xi0Core::XLambda, Xi0Core::YLambda, Xi0Core::ZLambda, + Xi0Core::DCADaughters, Xi0Core::DCADaughtersPi0, Xi0Core::DCAXYCascToPV, Xi0Core::DCAZCascToPV, + + // Xi0 daughter momenta + Xi0Core::Photon1Px, Xi0Core::Photon1Py, Xi0Core::Photon1Pz, + Xi0Core::Photon2Px, Xi0Core::Photon2Py, Xi0Core::Photon2Pz, + Xi0Core::LambdaPx, Xi0Core::LambdaPy, Xi0Core::LambdaPz, + + // Dynamic columns + Xi0Core::Px, + Xi0Core::Py, + Xi0Core::Pz, + Xi0Core::Pt, + Xi0Core::P, + Xi0Core::Rapidity, + Xi0Core::Phi, + Xi0Core::Eta, + Xi0Core::Radius, + Xi0Core::CascCosPA, + + // Dynamic columns for pi0 + Xi0Core::PxPi0, + Xi0Core::PyPi0, + Xi0Core::PzPi0, + Xi0Core::Pi0Mass, + Xi0Core::RadiusPi0, + Xi0Core::Pi0CosPA, + Xi0Core::DCAPi0ToPV, + + // Dynamic columns for Lambda + Xi0Core::RadiusLambda, + Xi0Core::LambdaCosPA, + Xi0Core::DCALambdaToPV, + Xi0Core::LambdaPt, + Xi0Core::LambdaP, + Xi0Core::LambdaEta, + Xi0Core::LambdaY, + Xi0Core::LambdaPhi, + + // Dynamic columns for Photon1 + Xi0Core::Photon1Pt, + Xi0Core::Photon1P, + Xi0Core::Photon1Eta, + Xi0Core::Photon1Y, + Xi0Core::Photon1Phi, + + // Dynamic columns for Photon2 + Xi0Core::Photon2Pt, + Xi0Core::Photon2P, + Xi0Core::Photon2Eta, + Xi0Core::Photon2Y, + Xi0Core::Photon2Phi); + +DECLARE_SOA_TABLE(Xi0CollRefs, "AOD", "XI0COLLREFS", //! optional table to refer back to a collision + 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); + +// 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(Photon1MCPx, photon1mcpx, float); +DECLARE_SOA_COLUMN(Photon1MCPy, photon1mcpy, float); +DECLARE_SOA_COLUMN(Photon1MCPz, photon1mcpz, float); +DECLARE_SOA_COLUMN(IsPhoton1Primary, isPhoton1Primary, bool); +DECLARE_SOA_COLUMN(Photon1PDGCodePositive, photon1PDGCodePositive, int); +DECLARE_SOA_COLUMN(Photon1PDGCodeNegative, photon1PDGCodeNegative, int); +DECLARE_SOA_COLUMN(Photon1PDGCode, photon1PDGCode, int); +DECLARE_SOA_COLUMN(Photon1PDGCodeMother, photon1PDGCodeMother, int); + +DECLARE_SOA_COLUMN(Photon2MCPx, photon2mcpx, float); +DECLARE_SOA_COLUMN(Photon2MCPy, photon2mcpy, float); +DECLARE_SOA_COLUMN(Photon2MCPz, photon2mcpz, float); +DECLARE_SOA_COLUMN(IsPhoton2Primary, isPhoton2Primary, bool); +DECLARE_SOA_COLUMN(Photon2PDGCodePositive, photon2PDGCodePositive, int); +DECLARE_SOA_COLUMN(Photon2PDGCodeNegative, photon2PDGCodeNegative, int); +DECLARE_SOA_COLUMN(Photon2PDGCode, photon2PDGCode, int); +DECLARE_SOA_COLUMN(Photon2PDGCodeMother, photon2PDGCodeMother, int); + +DECLARE_SOA_COLUMN(LambdaMCPx, lambdamcpx, float); +DECLARE_SOA_COLUMN(LambdaMCPy, lambdamcpy, float); +DECLARE_SOA_COLUMN(LambdaMCPz, lambdamcpz, float); +DECLARE_SOA_COLUMN(IsLambdaPrimary, isLambdaPrimary, bool); +DECLARE_SOA_COLUMN(LambdaPDGCodePositive, lambdaPDGCodePositive, int); +DECLARE_SOA_COLUMN(LambdaPDGCodeNegative, lambdaPDGCodeNegative, int); +DECLARE_SOA_COLUMN(LambdaPDGCode, lambdaPDGCode, int); +DECLARE_SOA_COLUMN(LambdaPDGCodeMother, lambdaPDGCodeMother, int); + +DECLARE_SOA_COLUMN(Pi0MCx, pi0mcx, float); +DECLARE_SOA_COLUMN(Pi0MCy, pi0mcy, float); +DECLARE_SOA_COLUMN(Pi0MCz, pi0mcz, float); +DECLARE_SOA_COLUMN(IsPi0Primary, isPi0Primary, bool); +DECLARE_SOA_COLUMN(Pi0MCprocess, pi0mcprocess, int); +DECLARE_SOA_COLUMN(Pi0PDGCode, pi0PDGCode, int); + +DECLARE_SOA_COLUMN(MCx, mcx, float); +DECLARE_SOA_COLUMN(MCy, mcy, float); +DECLARE_SOA_COLUMN(MCz, mcz, float); +DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, bool); +DECLARE_SOA_COLUMN(MCprocess, mcprocess, int); +DECLARE_SOA_COLUMN(PDGCode, pdgCode, int); + +DECLARE_SOA_COLUMN(PDGCodeMother, pdgCodeMother, int); + +DECLARE_SOA_DYNAMIC_COLUMN(MCPx, mcpx, //! Xi0 MC px + [](float photon1MCPx, float photon2MCPx, float lambdaMCPx) -> float { return photon1MCPx + photon2MCPx + lambdaMCPx; }); +DECLARE_SOA_DYNAMIC_COLUMN(MCPy, mcpy, //! Xi0 MC py + [](float photon1MCPy, float photon2MCPy, float lambdaMCPy) -> float { return photon1MCPy + photon2MCPy + lambdaMCPy; }); +DECLARE_SOA_DYNAMIC_COLUMN(MCPz, mcpz, //! Xi0 MC pz + [](float photon1MCPz, float photon2MCPz, float lambdaMCPz) -> float { return photon1MCPz + photon2MCPz + lambdaMCPz; }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCPt, mcpt, + [](float photon1MCPx, float photon1MCPy, float photon2MCPx, float photon2MCPy, float lambdaMCPx, float lambdaMCPy) -> float { + return RecoDecay::pt(std::array{photon1MCPx + photon2MCPx + lambdaMCPx, photon1MCPy + photon2MCPy + lambdaMCPy}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCP, mcp, //! Total momentum in GeV/c + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz, float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { + return RecoDecay::sqrtSumOfSquares(photon1MCPx + photon2MCPx + lambdaMCPx, photon1MCPy + photon2MCPy + lambdaMCPy, photon1MCPz + photon2MCPz + lambdaMCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(RapidityMC, rapidityMC, + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz, float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { + return RecoDecay::y(std::array{photon1MCPx + photon2MCPx + lambdaMCPx, photon1MCPy + photon2MCPy + lambdaMCPy, photon1MCPz + photon2MCPz + lambdaMCPz}, o2::constants::physics::MassXi0); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcPhi, //! Phi in the range [0, 2pi) + [](float photon1MCPx, float photon1MCPy, float photon2MCPx, float photon2MCPy, float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(photon1MCPx + photon2MCPx + lambdaMCPx, photon1MCPy + photon2MCPy + lambdaMCPy); }); + +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); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPx, pi0mcpx, //! Pi0 MC px + [](float photon1MCPx, float photon2MCPx) -> float { return photon1MCPx + photon2MCPx; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPy, pi0mcpy, //! Pi0 MC py + [](float photon1MCPy, float photon2MCPy) -> float { return photon1MCPy + photon2MCPy; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPz, pi0mcpz, //! Pi0 MC pz + [](float photon1MCPz, float photon2MCPz) -> float { return photon1MCPz + photon2MCPz; }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPt, pi0mcpt, + [](float photon1MCPx, float photon1MCPy, float photon2MCPx, float photon2MCPy) -> float { + return RecoDecay::pt(std::array{photon1MCPx + photon2MCPx, photon1MCPy + photon2MCPy}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCP, pi0mcp, //! Total momentum in GeV/c + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + return RecoDecay::sqrtSumOfSquares(photon1MCPx + photon2MCPx, photon1MCPy + photon2MCPy, photon1MCPz + photon2MCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCMass, pi0MCMass, + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + std::array pVecPhoton1{photon1MCPx, photon1MCPy, photon1MCPz}; + std::array pVecPhoton2{photon2MCPx, photon2MCPy, photon2MCPz}; + auto arrMom = std::array{pVecPhoton1, pVecPhoton2}; + return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassPhoton}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0RapidityMC, pi0RapidityMC, + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + return RecoDecay::y(std::array{photon1MCPx + photon2MCPx, photon1MCPy + photon2MCPy, photon1MCPz + photon2MCPz}, o2::constants::physics::MassPi0); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCPhi, pi0mcphi, //! Phi in the range [0, 2pi) + [](float photon1MCPx, float photon1MCPy, float photon2MCPx, float photon2MCPy) -> float { return RecoDecay::phi(photon1MCPx + photon2MCPx, photon1MCPy + photon2MCPy); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCEta, pi0mceta, //! Pseudorapidity + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + return RecoDecay::eta(std::array{photon1MCPx + photon2MCPx, photon1MCPy + photon2MCPy, photon1MCPz + photon2MCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pi0MCOPAngle, pi0mcopAngle, + [](float photon1MCPx, float photon1MCPy, float photon1MCPz, float photon2MCPx, float photon2MCPy, float photon2MCPz) { + TVector3 v1(photon1MCPx, photon1MCPy, photon1MCPz); + TVector3 v2(photon2MCPx, photon2MCPy, photon2MCPz); + return v1.Angle(v2); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1MCPt, photon1MCPt, //! Transverse momentum in GeV/c + [](float photon1MCPx, float photon1MCPy) -> float { + return RecoDecay::sqrtSumOfSquares(photon1MCPx, photon1MCPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1MCP, photon1MCp, //! Total momentum in GeV/c + [](float photon1MCPx, float photon1MCPy, float photon1MCPz) -> float { + return RecoDecay::sqrtSumOfSquares(photon1MCPx, photon1MCPy, photon1MCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1MCEta, photon1MCEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float photon1MCPx, float photon1MCPy, float photon1MCPz) -> float { + return RecoDecay::eta(std::array{photon1MCPx, photon1MCPy, photon1MCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1MCY, photon1MCY, //! Rapidity + [](float photon1MCPx, float photon1MCPy, float photon1MCPz) -> float { + return RecoDecay::y(std::array{photon1MCPx, photon1MCPy, photon1MCPz}, o2::constants::physics::MassGamma); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon1MCPhi, photon1MCPhi, //! Phi in the range [0, 2pi) + [](float photon1MCPx, float photon1MCPy) -> float { return RecoDecay::phi(photon1MCPx, photon1MCPy); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2MCPt, photon2MCPt, //! Transverse momentum in GeV/c + [](float photon2MCPx, float photon2MCPy) -> float { + return RecoDecay::sqrtSumOfSquares(photon2MCPx, photon2MCPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2MCP, photon2MCp, //! Total momentum in GeV/c + [](float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + return RecoDecay::sqrtSumOfSquares(photon2MCPx, photon2MCPy, photon2MCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2MCEta, photon2MCEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + return RecoDecay::eta(std::array{photon2MCPx, photon2MCPy, photon2MCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2MCY, photon2MCY, //! Rapidity + [](float photon2MCPx, float photon2MCPy, float photon2MCPz) -> float { + return RecoDecay::y(std::array{photon2MCPx, photon2MCPy, photon2MCPz}, o2::constants::physics::MassGamma); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Photon2MCPhi, photon2MCPhi, //! Phi in the range [0, 2pi) + [](float photon2MCPx, float photon2MCPy) -> float { return RecoDecay::phi(photon2MCPx, photon2MCPy); }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCPt, lambdaMCPt, //! Transverse momentum in GeV/c + [](float lambdaMCPx, float lambdaMCPy) -> float { + return RecoDecay::sqrtSumOfSquares(lambdaMCPx, lambdaMCPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCP, lambdaMCp, //! Total momentum in GeV/c + [](float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { + return RecoDecay::sqrtSumOfSquares(lambdaMCPx, lambdaMCPy, lambdaMCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCEta, lambdaMCEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { + return RecoDecay::eta(std::array{lambdaMCPx, lambdaMCPy, lambdaMCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCY, lambdaMCY, //! Rapidity + [](float lambdaMCPx, float lambdaMCPy, float lambdaMCPz) -> float { + return RecoDecay::y(std::array{lambdaMCPx, lambdaMCPy, lambdaMCPz}, o2::constants::physics::MassLambda); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCPhi, lambdaMCPhi, //! Phi in the range [0, 2pi) + [](float lambdaMCPx, float lambdaMCPy) -> float { return RecoDecay::phi(lambdaMCPx, lambdaMCPy); }); + +} // namespace Xi0MCCore + +DECLARE_SOA_TABLE(Xi0MCCores, "AOD", "XI0MCCORES", + // Basic properties + Xi0MCCore::Photon1MCPx, Xi0MCCore::Photon1MCPy, Xi0MCCore::Photon1MCPz, + Xi0MCCore::IsPhoton1Primary, Xi0MCCore::Photon1PDGCodePositive, Xi0MCCore::Photon1PDGCodeNegative, Xi0MCCore::Photon1PDGCode, Xi0MCCore::Photon1PDGCodeMother, + + Xi0MCCore::Photon2MCPx, Xi0MCCore::Photon2MCPy, Xi0MCCore::Photon2MCPz, + Xi0MCCore::IsPhoton2Primary, Xi0MCCore::Photon2PDGCodePositive, Xi0MCCore::Photon2PDGCodeNegative, Xi0MCCore::Photon2PDGCode, Xi0MCCore::Photon2PDGCodeMother, + + Xi0MCCore::LambdaMCPx, Xi0MCCore::LambdaMCPy, Xi0MCCore::LambdaMCPz, + Xi0MCCore::IsLambdaPrimary, Xi0MCCore::LambdaPDGCodePositive, Xi0MCCore::LambdaPDGCodeNegative, Xi0MCCore::LambdaPDGCode, Xi0MCCore::LambdaPDGCodeMother, + + Xi0MCCore::Pi0MCx, Xi0MCCore::Pi0MCy, Xi0MCCore::Pi0MCz, + Xi0MCCore::IsPi0Primary, Xi0MCCore::Pi0MCprocess, Xi0MCCore::Pi0PDGCode, + + Xi0MCCore::MCx, Xi0MCCore::MCy, Xi0MCCore::MCz, + Xi0MCCore::IsPhysicalPrimary, Xi0MCCore::MCprocess, Xi0MCCore::PDGCode, Xi0MCCore::PDGCodeMother, + + // Dynamic columns + Xi0MCCore::MCPx, + Xi0MCCore::MCPy, + Xi0MCCore::MCPz, + Xi0MCCore::MCPt, + Xi0MCCore::MCP, + Xi0MCCore::RapidityMC, + Xi0MCCore::MCPhi, + Xi0MCCore::MCEta, + + Xi0MCCore::MCRadius, + Xi0MCCore::MCRadiusPi0, + + Xi0MCCore::MCPx, + Xi0MCCore::MCPy, + Xi0MCCore::MCPz, + Xi0MCCore::MCPt, + Xi0MCCore::MCP, + Xi0MCCore::Pi0MCMass, + Xi0MCCore::Pi0RapidityMC, + Xi0MCCore::Pi0MCPhi, + Xi0MCCore::Pi0MCEta, + Xi0MCCore::Pi0MCOPAngle, + + Xi0MCCore::Photon1MCPt, + Xi0MCCore::Photon1MCP, + Xi0MCCore::Photon1MCEta, + Xi0MCCore::Photon1MCY, + Xi0MCCore::Photon1MCPhi, + + Xi0MCCore::Photon2MCPt, + Xi0MCCore::Photon2MCP, + Xi0MCCore::Photon2MCEta, + Xi0MCCore::Photon2MCY, + Xi0MCCore::Photon2MCPhi, + + Xi0MCCore::LambdaMCPt, + Xi0MCCore::LambdaMCP, + Xi0MCCore::LambdaMCEta, + Xi0MCCore::LambdaMCY, + Xi0MCCore::LambdaMCPhi); + +DECLARE_SOA_TABLE(Xi0MCCollRefs, "AOD", "XI0MCCOLLREFS", //! optional table to refer back to a collision + 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); + } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFSIGMATABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 643c63aa9c0..a61a01e6a56 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -149,7 +149,7 @@ o2physics_add_dpl_workflow(cascademlselection o2physics_add_dpl_workflow(sigma0builder SOURCES sigma0builder.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2Physics::AnalysisCore O2Physics::MLCore O2Physics::AnalysisCCDB + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2Physics::AnalysisCore O2Physics::MLCore O2Physics::AnalysisCCDB O2::DCAFitter COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(lambdajetpolarizationbuilder # o2-linter: disable=name/o2-workflow (Keep historical workflow name for backward compatibility.) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index 214259f75c4..090103f41db 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -25,10 +25,13 @@ #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 @@ -40,6 +43,7 @@ #include #include #include +#include #include // IWYU pragma: keep (do not replace with Math/Vector3Dfwd.h) #include @@ -73,7 +77,10 @@ static const std::vector DirList = {"V0BeforeSel", "PhotonSel", "La static const std::vector DirList2 = {"EMCalPhotonBeforeSel", "EMCalPhotonSel"}; struct sigma0builder { + o2::common::StandardCCDBLoaderConfigurables standardCCDBLoaderConfigurables; + o2::common::StandardCCDBLoader ccdbLoader; Service ccdb; + int mRunNumber = 0; ctpRateFetcher rateFetcher; //___________________________________________________ @@ -105,6 +112,15 @@ struct sigma0builder { Produces pi0Gens; // Generated pi0s Produces pi0GenCollRefs; // references collisions from pi0Gens + //__________________________________________________ + // Xi0 specific + Produces xi0cores; // xi0 candidates info for analysis + Produces xi0collRefs; // references to straCollision + Produces xi0Indices; // references to V0Cores + Produces xi0mccores; // Reco xi0 MC properties + Produces xi0mccollRefs; // references to straMCCollision + Produces xi0mcIndices; // references to V0MCCores + //__________________________________________________ // pack track quality but separte also afterburner // dynamic range: 0-31 @@ -176,6 +192,7 @@ struct sigma0builder { Configurable fillPi0Tables{"fillPi0Tables", false, "fill pi0 tables for QA"}; Configurable fillSigma0Tables{"fillSigma0Tables", true, "fill sigma0 tables for analysis"}; Configurable fillKStarTables{"fillKStarTables", true, "fill kstar tables for analysis"}; + Configurable fillXi0Tables{"fillXi0Tables", false, "fill xi0 tables for analysis"}; // For ML Selection Configurable useMLScores{"useMLScores", false, "use ML scores to select candidates"}; @@ -284,6 +301,19 @@ struct sigma0builder { Configurable KShortMaxTPCNSigmas{"KShortMaxTPCNSigmas", 1e+9, "Max |TPC NSigma| (pion hypothesis) for K0S daughters"}; } kshortSelections; + // Xi0 criteria: + struct : ConfigurableGroup { + std::string prefix = "cascadeSelections"; // JSON group name + Configurable maxPi0Rapidity{"maxPi0Rapidity", 1.0, "Maximum rapidity of pi0"}; + Configurable massPi0Window{"massPi0Window", 0.115, "Maximum pi0 mass (GeV/c)"}; + Configurable dcaPi0daughters{"dcaPi0daughters", 1., "Maximum DCA between pi0 daughters"}; + Configurable radiusPi0{"radiusPi0", 2.5, "Max DCA V0 Daughters (cm)"}; + Configurable dcaPi0ToPV{"dcaPi0ToPV", 0.0, "Minimum Pi0 DCA to PV(cm)"}; + Configurable dcaCascadeDaughters{"dcaCascadeDaughters", 1., "Maximum DCA between cascade daughters"}; + Configurable radiusCascade{"radiusCascade", 2.5, "Minimum cascade decay radius (cm)"}; + Configurable cosPACascade{"cosPACascade", 0.95, "Min cascade CosPA"}; + } cascadeSelections; + // KStar criteria: Configurable KStarWindow{"KStarWindow", 0.1, "Mass window around expected (in GeV/c2)"}; Configurable KStarMaxRap{"KStarMaxRap", 0.8, "Max kstar rapidity"}; @@ -315,11 +345,13 @@ struct sigma0builder { ConfigurableAxis axisNch{"axisNch", {300, 0.0f, 3000.0f}, "N_{ch}"}; // Invariant Mass - ConfigurableAxis axisSigmaMass{"axisSigmaMass", {500, 1.10f, 1.30f}, "M_{#Sigma^{0}} (GeV/c^{2})"}; - ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, "M_{#Lambda} (GeV/c^{2})"}; - ConfigurableAxis axisPhotonMass{"axisPhotonMass", {200, 0.0f, 0.3f}, "M_{#Gamma}"}; - ConfigurableAxis axisK0SMass{"axisK0SMass", {200, 0.4f, 0.6f}, "M_{K^{0}}"}; - ConfigurableAxis axisKStarMass{"axisKStarMass", {500, 0.6f, 1.6f}, "M_{K^{*}} (GeV/c^{2})"}; + ConfigurableAxis axisSigmaMass{"axisSigmaMass", {500, 1.10f, 1.30f}, "#it{M}_{#Sigma^{0}} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, "#it{M}_{#Lambda} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisPhotonMass{"axisPhotonMass", {200, 0.0f, 0.3f}, "#it{M}_{#Gamma} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisK0SMass{"axisK0SMass", {200, 0.4f, 0.6f}, "#it{M}_{K^{0}_{S}} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisKStarMass{"axisKStarMass", {500, 0.6f, 1.6f}, "#it{M}_{K^{*}} (GeV/#it{c}}^{2})"}; + ConfigurableAxis axisPi0Mass{"axisPi0Mass", {100, 0.08f, 0.18f}, "#it{M}_{#pi^{0}} (GeV/#it{c}^{2})"}; + ConfigurableAxis axisXi0Mass{"axisXi0Mass", {100, 1.27f, 1.37f}, "#it{M}_{#Xi^{0}} (GeV/#it{c}^{2})"}; // AP plot axes ConfigurableAxis axisAPAlpha{"axisAPAlpha", {220, -1.1f, 1.1f}, "V0 AP alpha"}; @@ -350,6 +382,30 @@ struct sigma0builder { ConfigurableAxis axisClrShape{"axisClrShape", {100, 0.0, 1.0}, "cluster shape"}; } axisConfig; + struct : ConfigurableGroup { + std::string prefix = "fitterConfiguration"; // JSON group name + Configurable propagateToPCA{"propagateToPCA", true, "Propagate to PCA?"}; + Configurable minParamChange{"minParamChange", 4., "Stop minimization iterations if largest change of any X is smaller than this."}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "Stop iterations is chi2/chi2old > this"}; + Configurable maxR{"maxR", 200., "Don't consider as a seed (circles intersection) if its R exceeds this"}; + Configurable maxDXYIni{"maxDXYIni", 4.0f, "Don't consider as a seed (circles intersection) if XY distance exceeds this"}; + Configurable maxDZIni{"maxDZIni", 1e9, "Don't consider as a seed (circles intersection) if Z distance exceeds this"}; + Configurable maxChi2{"maxChi2", 1e9, "Max dca from prongs to vertex"}; + Configurable useAbsDCA{"useAbsDCA", true, "Use abs dca minimization"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Max generated particle rapidity"}; + } fitterConfiguration; + + // CCDB options + struct : ConfigurableGroup { + // manual + Configurable useCustomRunNumber{"ccdbConfigurations.useCustomRunNumber", false, "Use custom run number"}; + Configurable customRunNumber{"ccdbConfigurations.customRunNumber", 544122, "Manually set the run number "}; + Configurable ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + } ccdbConfigurations; + + o2::vertexing::DCAFitterN<2> fitter2Prongs; // 2-prong o2 dca fitter + o2::vertexing::DCAFitterN<3> fitter3Prongs; // 3-prong o2 dca fitter + void init(InitContext const&) { LOGF(info, "Initializing now: cross-checking correctness..."); @@ -366,9 +422,9 @@ struct sigma0builder { } // setting CCDB service - ccdb->setURL("http://alice-ccdb.cern.ch"); ccdb->setCaching(true); - ccdb->setFatalWhenNull(false); + ccdb->setLocalObjectValidityChecking(); + ccdb->setURL(ccdbConfigurations.ccdburl.value); histos.add("hEventCentrality", "hEventCentrality", kTH1D, {axisConfig.axisCentrality}); @@ -562,6 +618,11 @@ struct sigma0builder { histos.add("KStarSel/hKStarMassSelected", "hKStarMassSelected", kTH1F, {axisConfig.axisKStarMass}); } + if (fillXi0Tables) { + histos.add("Xi0Sel/hPi0Mass", "hPi0Mass", kTH1F, {axisConfig.axisPi0Mass}); + histos.add("Xi0Sel/hXi0Mass", "hXi0Mass", kTH1F, {axisConfig.axisXi0Mass}); + } + if (doAssocStudy && (doprocessMonteCarlo || doprocessMonteCarloWithTOF)) { histos.add("V0AssoQA/h2dIRVsPt_TrueGamma", "h2dIRVsPt_TrueGamma", kTH2F, {axisConfig.axisIRBinning, axisConfig.axisPt}); histos.add("V0AssoQA/h3dPAVsIRVsPt_TrueGamma", "h3dPAVsIRVsPt_TrueGamma", kTH3F, {axisConfig.axisPA, axisConfig.axisIRBinning, axisConfig.axisPt}); @@ -782,18 +843,66 @@ struct sigma0builder { histos.add("V0QA/h2dGenKShortVsMultMC", "h2dGenKShortVsMultMC", kTH2D, {axisConfig.axisNch, axisConfig.axisPt}); } + fitter2Prongs.setPropagateToPCA(fitterConfiguration.propagateToPCA); + fitter2Prongs.setMaxR(fitterConfiguration.maxR); + fitter2Prongs.setMinParamChange(fitterConfiguration.minParamChange); + fitter2Prongs.setMinRelChi2Change(fitterConfiguration.minRelChi2Change); + fitter2Prongs.setMaxDZIni(fitterConfiguration.maxDZIni); + fitter2Prongs.setMaxDXYIni(fitterConfiguration.maxDXYIni); + fitter2Prongs.setMaxChi2(fitterConfiguration.maxChi2); + fitter2Prongs.setUseAbsDCA(fitterConfiguration.useAbsDCA); + fitter2Prongs.setWeightedFinalPCA(fitterConfiguration.useWeightedFinalPCA); + + fitter3Prongs.setPropagateToPCA(fitterConfiguration.propagateToPCA); + fitter3Prongs.setMaxR(fitterConfiguration.maxR); + fitter3Prongs.setMinParamChange(fitterConfiguration.minParamChange); + fitter3Prongs.setMinRelChi2Change(fitterConfiguration.minRelChi2Change); + fitter3Prongs.setMaxDZIni(fitterConfiguration.maxDZIni); + fitter3Prongs.setMaxDXYIni(fitterConfiguration.maxDXYIni); + fitter3Prongs.setMaxChi2(fitterConfiguration.maxChi2); + fitter3Prongs.setUseAbsDCA(fitterConfiguration.useAbsDCA); + fitter3Prongs.setWeightedFinalPCA(fitterConfiguration.useWeightedFinalPCA); + + fitter2Prongs.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrLUT); + fitter3Prongs.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrLUT); + + // mag field has to be set later + fitter2Prongs.setBz(-999.9f); // will NOT make sense if not changed + fitter3Prongs.setBz(-999.9f); // will NOT make sense if not changed + // inspect histogram sizes, please histos.print(); } + template // TCollision should be of the type: soa::Join::iterator or so + void initCCDB(TCollision const& collision) + { + if (mRunNumber == collision.runNumber() || (ccdbConfigurations.useCustomRunNumber && mRunNumber == ccdbConfigurations.customRunNumber)) { + return; + } + + mRunNumber = ccdbConfigurations.useCustomRunNumber ? ccdbConfigurations.customRunNumber : collision.runNumber(); + + ccdbLoader.initCCDB(standardCCDBLoaderConfigurables, ccdb, mRunNumber); + + auto magneticField = o2::base::Propagator::Instance()->getNominalBz(); + // Set magnetic field value once known + fitter2Prongs.setBz(magneticField); + fitter3Prongs.setBz(magneticField); + } + // ______________________________________________________ // Struct to store V0Pair properties struct V0PairTopoInfo { - float X = -999.f; - float Y = -999.f; - float Z = -999.f; - float DCADau = -999.f; - float CosPA = -1.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; + + int v01Index = 0; // index to de-reference V0Cores table + int v02Index = 0; // index to de-reference V0Cores table }; // ______________________________________________________ @@ -856,6 +965,79 @@ struct sigma0builder { float MCvy = 999.f; }; + // ______________________________________________________ + // 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.}; + 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.}; + }; + + 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; + float gamma1MCpz = -999.f; + float gamma2MCpx = -999.f; + float gamma2MCpy = -999.f; + float gamma2MCpz = -999.f; + float lambdaMCpx = -999.f; + float lambdaMCpy = -999.f; + float lambdaMCpz = -999.f; + + // MC association info + bool gamma1IsPhysicalPrimary = false; + bool gamma2IsPhysicalPrimary = false; + bool lambdaIsPhysicalPrimary = false; + int gamma1PDGCodePos = 0; + int gamma1PDGCodeNeg = 0; + int gamma2PDGCodePos = 0; + int gamma2PDGCodeNeg = 0; + int lambdaPDGCodePos = 0; + int lambdaPDGCodeNeg = 0; + int gamma1PDGCode = 0; + int gamma2PDGCode = 0; + int lambdaPDGCode = 0; + int gamma1PDGCodeMother = 0; + int gamma2PDGCodeMother = 0; + int lambdaPDGCodeMother = 0; + + int pi0PDGCode = 0; + int pi0MCProcess = 0; + bool pi0IsPhysicalPrimary = false; + float pi0X = -999.f; + float pi0Y = -999.f; + float pi0Z = -999.f; + + int xi0PDGCode = 0; + int xi0MCProcess = 0; + bool xi0IsPhysicalPrimary = false; + float xi0X = -999.f; + float xi0Y = -999.f; + float xi0Z = -999.f; + + int xi0PDGCodeMother = 0; + }; + template V0PairTopoInfo propagateV0PairToDCA(TV01 const& v01, TV02 const& v02) { @@ -886,18 +1068,18 @@ struct sigma0builder { ROOT::Math::XYZVector PCA = 0.5 * (pointOn1 + pointOn2); // Calculate properties and fill struct - info.DCADau = (cross.Mag2() > 0) ? std::abs(posdiff.Dot(cross)) / cross.R() : 999.f; - info.CosPA = v01momentumNorm.Dot(v02momentumNorm); + info.daughterDCA = (cross.Mag2() > 0) ? std::abs(posdiff.Dot(cross)) / cross.R() : 999.f; + info.cosPA = v01momentumNorm.Dot(v02momentumNorm); float Min_threshold = 1e-5f; // Threshold to consider lines as parallel, can be tuned if (d < Min_threshold) { // Parallel or nearly parallel lines - info.X = info.Y = info.Z = 0.f; // should we use another dummy value? Perhaps 999.f? + info.position[0] = info.position[1] = info.position[2] = 0.f; // should we use another dummy value? Perhaps 999.f? return info; } - info.X = PCA.X(); - info.Y = PCA.Y(); - info.Z = PCA.Z(); + info.position[0] = PCA.X(); + info.position[1] = PCA.Y(); + info.position[2] = PCA.Z(); return info; } @@ -1249,6 +1431,105 @@ struct sigma0builder { return MCinfo; } + template + Xi0MCInfo getXi0MCInfo(TV0 const& gamma1, TV0 const& gamma2, TV0 const& lambda, TMCParticles const& mcparticles) + { + Xi0MCInfo MCinfo; + + if (!gamma1.has_v0MCCore() || !gamma2.has_v0MCCore() || !lambda.has_v0MCCore()) { + return MCinfo; + } + + auto gamma1MC = gamma1.template v0MCCore_as>(); + auto gamma2MC = gamma2.template v0MCCore_as>(); + auto lambdaMC = lambda.template v0MCCore_as>(); + + // Basic kinematic info + MCinfo.gamma1MCpx = gamma1MC.pxMC(); + MCinfo.gamma1MCpy = gamma1MC.pyMC(); + MCinfo.gamma1MCpz = gamma1MC.pzMC(); + MCinfo.gamma2MCpx = gamma2MC.pxMC(); + MCinfo.gamma2MCpy = gamma2MC.pyMC(); + MCinfo.gamma2MCpz = gamma2MC.pzMC(); + MCinfo.lambdaMCpx = lambdaMC.pxMC(); + MCinfo.lambdaMCpy = lambdaMC.pyMC(); + MCinfo.lambdaMCpz = lambdaMC.pzMC(); + + // MC association info + MCinfo.gamma1Index = gamma1MC.globalIndex(); + MCinfo.gamma2Index = gamma2MC.globalIndex(); + MCinfo.lambdaIndex = lambdaMC.globalIndex(); + MCinfo.gamma1IsPhysicalPrimary = gamma1MC.isPhysicalPrimary(); + MCinfo.gamma2IsPhysicalPrimary = gamma2MC.isPhysicalPrimary(); + MCinfo.lambdaIsPhysicalPrimary = lambdaMC.isPhysicalPrimary(); + MCinfo.gamma1PDGCodePos = gamma1MC.pdgCodePositive(); + MCinfo.gamma1PDGCodeNeg = gamma1MC.pdgCodeNegative(); + MCinfo.gamma2PDGCodePos = gamma2MC.pdgCodePositive(); + MCinfo.gamma2PDGCodeNeg = gamma2MC.pdgCodeNegative(); + MCinfo.lambdaPDGCodePos = lambdaMC.pdgCodePositive(); + MCinfo.lambdaPDGCodeNeg = lambdaMC.pdgCodeNegative(); + MCinfo.gamma1PDGCode = gamma1MC.pdgCode(); + MCinfo.gamma2PDGCode = gamma2MC.pdgCode(); + MCinfo.lambdaPDGCode = lambdaMC.pdgCode(); + + // Get corresponding entries in MCParticles table + auto MCParticle_gamma1 = mcparticles.rawIteratorAt(gamma1MC.particleIdMC()); + 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) && + (std::abs(MCParticle_lambda.pdgCode()) == PDG_t::kLambda0)) { + + // Get MC Mothers + auto const& MCMothersList_gamma1 = MCParticle_gamma1.template mothers_as(); + auto const& MCMothersList_gamma2 = MCParticle_gamma2.template mothers_as(); + 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 + + MCinfo.gamma1PDGCodeMother = MCMother_gamma1.pdgCode(); + MCinfo.gamma2PDGCodeMother = MCMother_gamma2.pdgCode(); + MCinfo.lambdaPDGCodeMother = MCMother_lambda.pdgCode(); + + if (MCMother_gamma1.globalIndex() == MCMother_gamma2.globalIndex() && std::abs(MCinfo.gamma2PDGCodeMother) == PDG_t::kPi0) { // check that gamma1 and gamma2 have the same pi0 mother + + MCinfo.pi0PDGCode = MCMother_gamma1.pdgCode(); + MCinfo.pi0MCProcess = MCMother_gamma1.getProcess(); + MCinfo.pi0IsPhysicalPrimary = MCMother_gamma1.isPhysicalPrimary(); + MCinfo.pi0X = MCMother_gamma1.vx(); + MCinfo.pi0Y = MCMother_gamma1.vy(); + 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 (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.xi0X = MCMother_pi0.vx(); + MCinfo.xi0Y = MCMother_pi0.vy(); + MCinfo.xi0Z = MCMother_pi0.vz(); + + auto const& v0pairmothers = MCMother_pi0.template mothers_as(); // Get mothers + if (!v0pairmothers.empty()) { + auto& v0PairMother = v0pairmothers.front(); // V0Pair mother, V0s grandmother + MCinfo.xi0PDGCodeMother = v0PairMother.pdgCode(); + } + } + } + } + } + } + + return MCinfo; + } + // ______________________________________________________ // Check whether the collision passes our collision selections // Should work with collisions, mccollisions, stracollisions and stramccollisions tables! @@ -1629,9 +1910,9 @@ 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.IsKStar = std::abs(mcParticle.pdgCode()) == o2::constants::physics::Pdg::kK0Star892; // 313; 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(); GenInfo.MCPt = mcParticle.pt(); @@ -2310,7 +2591,7 @@ struct sigma0builder { pi0MCInfo.fIsV02Primary, pi0MCInfo.V02PDGCode, pi0MCInfo.V02PDGCodeMother, pi0MCInfo.fIsV02CorrectlyAssign); } - pi0cores(pi0TopoInfo.X, pi0TopoInfo.Y, pi0TopoInfo.Z, pi0TopoInfo.DCADau, pi0TopoInfo.CosPA, + pi0cores(pi0TopoInfo.position[0], pi0TopoInfo.position[1], pi0TopoInfo.position[2], pi0TopoInfo.daughterDCA, pi0TopoInfo.cosPA, gamma1.px(), gamma1.py(), gamma1.pz(), gamma1.mGamma(), gamma1.qtarm(), gamma1.alpha(), gamma1.dcapostopv(), gamma1.dcanegtopv(), gamma1.dcaV0daughters(), gamma1.negativeeta(), gamma1.positiveeta(), gamma1.v0cosPA(), gamma1.v0radius(), gamma1.z(), @@ -2325,6 +2606,274 @@ struct sigma0builder { return true; } + //_______________________________________________ + // Build pi0 candidate for QA + template + bool buildPi0ForXi0(TV0Object const& gamma1, TV0Object const& gamma2, TCollision const& collision, V0PairTopoInfo& info) + { + //_______________________________________________ + // Check if both V0s are made of the same tracks + if (gamma1.posTrackExtraId() == gamma2.posTrackExtraId() || + gamma1.negTrackExtraId() == gamma2.negTrackExtraId()) { + return false; + } + + //_______________________________________________ + // Calculate pi0 properties + const std::array vtxGamma1 = {gamma1.x(), gamma1.y(), gamma1.z()}; + const std::array vtxGamma2 = {gamma2.x(), gamma2.y(), gamma2.z()}; + + const std::array momGamma1 = {gamma1.px(), gamma1.py(), gamma1.pz()}; + const std::array momGamma2 = {gamma2.px(), gamma2.py(), gamma2.pz()}; + + const std::array covGamma1 = {999.}; + const std::array covGamma2 = {999.}; + + std::array pVecGamma1{gamma1.px(), gamma1.py(), gamma1.pz()}; + std::array pVecGamma2{gamma2.px(), gamma2.py(), gamma2.pz()}; + std::array arrpi0{pVecGamma1, pVecGamma2}; + info.pi0Mass = RecoDecay::m(arrpi0, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassPhoton}); + info.pi0Y = RecoDecay::y(std::array{gamma1.px() + gamma2.px(), gamma1.py() + gamma2.py(), gamma1.pz() + gamma2.pz()}, o2::constants::physics::MassPi0); + + info.v01Index = gamma1.globalIndex(); + info.v02Index = gamma2.globalIndex(); + + //_______________________________________________ + // Pi0-specific selections: + if (std::abs(info.pi0Y) > cascadeSelections.maxPi0Rapidity) + return false; + + if (std::abs(info.pi0Mass - o2::constants::physics::MassPi0) > cascadeSelections.massPi0Window) + return false; + + o2::track::TrackParCov gammaTrack1(vtxGamma1, momGamma1, covGamma1, 0, true); + gammaTrack1.setAbsCharge(0); + gammaTrack1.setPID(o2::track::PID::Photon); + o2::track::TrackParCov gammaTrack2(vtxGamma2, momGamma2, covGamma2, 0, true); + gammaTrack2.setAbsCharge(0); + gammaTrack2.setPID(o2::track::PID::Photon); + + // First build Pi0 to select secondary pi0s + int nCandPi0 = 0; + try { + nCandPi0 = fitter2Prongs.process(gammaTrack1, gammaTrack2); + } catch (...) { + return false; + } + if (nCandPi0 == 0) { + return false; + } + + fitter2Prongs.propagateTracksToVertex(); // propagate e and K to D vertex + if (!fitter2Prongs.isPropagateTracksToVertexDone()) { + return false; + } + + info.daughterDCA = TMath::Sqrt(fitter2Prongs.getChi2AtPCACandidate()); + + if (info.daughterDCA > cascadeSelections.dcaPi0daughters) { + return false; + } + + // get decay vertex coordinates + const auto& vtxPi0 = fitter2Prongs.getPCACandidate(); + for (int i = 0; i < 3; i++) { + info.position[i] = vtxPi0[i]; + } + 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]; + float px = pVecGamma1[0] + pVecGamma2[0]; + float py = pVecGamma1[1] + pVecGamma2[1]; + float pz = pVecGamma1[2] + pVecGamma2[2]; + info.dcaPi0ToPV = std::sqrt((std::pow((collision.posY() - y) * pz - (collision.posZ() - z) * py, 2) + std::pow((collision.posX() - x) * pz - (collision.posZ() - z) * px, 2) + std::pow((collision.posX() - x) * py - (collision.posY() - y) * px, 2)) / (px * px + py * py + pz * pz)); + if (info.dcaPi0ToPV < cascadeSelections.dcaPi0ToPV) { + return false; + } + + histos.fill(HIST("Xi0Sel/hPi0Mass"), info.pi0Mass); + + return true; + } + + //_______________________________________________ + // Build Xi0 candidate for analysis + template + bool buildXi0(TV0Object const& gamma1, TV0Object const& gamma2, TV0Object const& lambda, TCollision const& collision, TMCParticles const& mcparticles) + { + Xi0Info cascade; + //_______________________________________________ + // Check if the 3 V0s are made of the same tracks + if (gamma1.posTrackExtraId() == gamma2.posTrackExtraId() || + gamma1.negTrackExtraId() == gamma2.negTrackExtraId() || + gamma1.posTrackExtraId() == lambda.posTrackExtraId() || + gamma1.negTrackExtraId() == lambda.negTrackExtraId() || + gamma2.posTrackExtraId() == lambda.posTrackExtraId() || + gamma2.negTrackExtraId() == lambda.negTrackExtraId()) { + return false; + } + + const std::array vtxGamma1 = {gamma1.x(), gamma1.y(), gamma1.z()}; + const std::array vtxGamma2 = {gamma2.x(), gamma2.y(), gamma2.z()}; + const std::array vtxLambda = {lambda.x(), lambda.y(), lambda.z()}; + + const std::array momGamma1 = {gamma1.px(), gamma1.py(), gamma1.pz()}; + const std::array momGamma2 = {gamma2.px(), gamma2.py(), gamma2.pz()}; + const std::array momLambda = {lambda.px(), lambda.py(), lambda.pz()}; + + const std::array covGamma1 = {999.}; + const std::array covGamma2 = {999.}; + const std::array covLambda = {999.}; + + o2::track::TrackParCov gammaTrack1(vtxGamma1, momGamma1, covGamma1, 0, true); + gammaTrack1.setAbsCharge(0); + gammaTrack1.setPID(o2::track::PID::Photon); + o2::track::TrackParCov gammaTrack2(vtxGamma2, momGamma2, covGamma2, 0, true); + gammaTrack2.setAbsCharge(0); + gammaTrack2.setPID(o2::track::PID::Photon); + o2::track::TrackParCov lambdaTrack(vtxLambda, momLambda, covLambda, 0, true); + lambdaTrack.setAbsCharge(0); + lambdaTrack.setPID(o2::track::PID::Lambda); + + // build Xi0 --> pi0 Lambda --> gamma gamma Lambda + int nCandXi0 = 0; + try { + nCandXi0 = fitter3Prongs.process(gammaTrack1, gammaTrack2, lambdaTrack); + } catch (...) { + return false; + } + if (nCandXi0 == 0) { + return false; + } + + fitter3Prongs.propagateTracksToVertex(); // propagate e and K to D vertex + if (!fitter3Prongs.isPropagateTracksToVertexDone()) { + return false; + } + + auto lGamma1Track = fitter3Prongs.getTrack(0); + auto lGamma2Track = fitter3Prongs.getTrack(1); + auto lLambdaTrack = fitter3Prongs.getTrack(2); + + cascade.cascadeDaughterDCA = TMath::Sqrt(fitter3Prongs.getChi2AtPCACandidate()); + if (cascade.cascadeDaughterDCA > cascadeSelections.dcaCascadeDaughters) { + return false; + } + + lGamma1Track.getPxPyPzGlo(cascade.gamma1Momentum); + lGamma2Track.getPxPyPzGlo(cascade.gamma2Momentum); + lLambdaTrack.getPxPyPzGlo(cascade.lambdaMomentum); + // get decay vertex coordinates + const auto& vtx = fitter3Prongs.getPCACandidate(); + for (int i = 0; i < 3; i++) { + cascade.cascadePosition[i] = vtx[i]; + } + if (std::hypot(cascade.cascadePosition[0], cascade.cascadePosition[1]) < cascadeSelections.radiusCascade) { + return false; + } + + double cosPA = RecoDecay::cpa( + std::array{collision.posX(), collision.posY(), collision.posZ()}, + std::array{cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2]}, + std::array{cascade.gamma1Momentum[0] + cascade.gamma2Momentum[0] + cascade.lambdaMomentum[0], + cascade.gamma1Momentum[1] + cascade.gamma2Momentum[1] + cascade.lambdaMomentum[1], + cascade.gamma1Momentum[2] + cascade.gamma2Momentum[2] + cascade.lambdaMomentum[2]}); + if (cosPA < cascadeSelections.cosPACascade) { + return false; + } + + // Calculate DCAxy of the cascade (with bending) + auto lCascadeTrack = fitter3Prongs.createParentTrackParCov(); + lCascadeTrack.setAbsCharge(0); // to be sure + lCascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas + 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]; + cascade.cascadeDCAz = dcaInfo[1]; + + // Populate information + // cascadecandidate.v0Id = v0index.globalIndex(); + cascade.collisionId = collision.globalIndex(); + cascade.gamma1Index = gamma1.globalIndex(); + cascade.gamma2Index = gamma2.globalIndex(); + cascade.lambdaIndex = lambda.globalIndex(); + + 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; + float positionCovariance[6]; + positionCovariance[0] = covVtxV(0, 0); + positionCovariance[1] = covVtxV(1, 0); + positionCovariance[2] = covVtxV(1, 1); + positionCovariance[3] = covVtxV(2, 0); + positionCovariance[4] = covVtxV(2, 1); + positionCovariance[5] = covVtxV(2, 2); + // store momentum covariance matrix + std::array covTgamma1 = {0.}; + std::array covTgamma2 = {0.}; + std::array covTlambda = {0.}; + // std::array momentumCovariance; + lGamma1Track.getCovXYZPxPyPzGlo(covTgamma1); + lGamma2Track.getCovXYZPxPyPzGlo(covTgamma2); + lLambdaTrack.getCovXYZPxPyPzGlo(covTlambda); + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 21; i++) { + cascade.covariance[i] = 0.0f; + } + for (int i = 0; i < 6; i++) { + cascade.covariance[i] = positionCovariance[i]; + cascade.covariance[MomInd[i]] = covTgamma1[MomInd[i]] + covTgamma2[MomInd[i]] + covTlambda[MomInd[i]]; + } + + // 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); + + xi0mccores(xi0MCInfo.gamma1MCpx, xi0MCInfo.gamma1MCpy, xi0MCInfo.gamma1MCpz, + xi0MCInfo.gamma1IsPhysicalPrimary, xi0MCInfo.gamma1PDGCodePos, xi0MCInfo.gamma1PDGCodeNeg, xi0MCInfo.gamma1PDGCode, xi0MCInfo.gamma1PDGCodeMother, + xi0MCInfo.gamma2MCpx, xi0MCInfo.gamma2MCpy, xi0MCInfo.gamma2MCpz, + xi0MCInfo.gamma2IsPhysicalPrimary, xi0MCInfo.gamma2PDGCodePos, xi0MCInfo.gamma2PDGCodeNeg, xi0MCInfo.gamma2PDGCode, xi0MCInfo.gamma2PDGCodeMother, + xi0MCInfo.lambdaMCpx, xi0MCInfo.lambdaMCpy, xi0MCInfo.lambdaMCpz, + 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; + if (collision.has_straMCCollision()) { + auto mcCollision = collision.template straMCCollision_as>(); + mcCollisionIndex = mcCollision.globalIndex(); + } + xi0mccollRefs(mcCollisionIndex); + } + + xi0cores(cascade.cascadePosition[0], cascade.cascadePosition[1], cascade.cascadePosition[2], + cascade.pi0Position[0], cascade.pi0Position[1], cascade.pi0Position[2], + vtxLambda[0], vtxLambda[1], vtxLambda[2], + cascade.cascadeDaughterDCA, cascade.pi0DaughterDCA, cascade.cascadeDCAxy, cascade.cascadeDCAz, + cascade.gamma1Momentum[0], cascade.gamma1Momentum[1], cascade.gamma1Momentum[2], + cascade.gamma2Momentum[0], cascade.gamma2Momentum[1], cascade.gamma2Momentum[2], + cascade.lambdaMomentum[0], cascade.lambdaMomentum[1], cascade.lambdaMomentum[2]); + + xi0Indices(cascade.gamma1Index, cascade.gamma2Index, cascade.lambdaIndex); + + xi0collRefs(collision.globalIndex()); + + histos.fill(HIST("Xi0Sel/hXi0Mass"), cascade.xi0Mass); + + return true; + } + //_______________________________________________ // Build sigma0 candidate with PCM photons template @@ -2374,7 +2923,7 @@ struct sigma0builder { // Sigma0 topological info auto sigma0TopoInfo = propagateV0PairToDCA(gamma, lambda); - sigma0cores(gamma.globalIndex(), lambda.globalIndex(), sigma0TopoInfo.X, sigma0TopoInfo.Y, sigma0TopoInfo.Z, sigma0TopoInfo.DCADau, + sigma0cores(gamma.globalIndex(), lambda.globalIndex(), sigma0TopoInfo.position[0], sigma0TopoInfo.position[1], sigma0TopoInfo.position[2], sigma0TopoInfo.daughterDCA, gamma.px(), gamma.py(), gamma.pz(), gamma.mGamma(), lambda.px(), lambda.py(), lambda.pz(), lambda.mLambda(), lambda.mAntiLambda()); // MC properties @@ -2601,7 +3150,7 @@ struct sigma0builder { auto kstarTopoInfo = propagateV0PairToDCA(gamma, kshort); - kstarcores(gamma.globalIndex(), kshort.globalIndex(), kstarTopoInfo.X, kstarTopoInfo.Y, kstarTopoInfo.Z, kstarTopoInfo.DCADau, + kstarcores(gamma.globalIndex(), kshort.globalIndex(), kstarTopoInfo.position[0], kstarTopoInfo.position[1], kstarTopoInfo.position[2], kstarTopoInfo.daughterDCA, gamma.px(), gamma.py(), gamma.pz(), gamma.mGamma(), kshort.px(), kshort.py(), kshort.pz(), kshort.mK0Short()); // MC properties @@ -2720,6 +3269,8 @@ struct sigma0builder { //_______________________________________________ // Collisions loop for (const auto& coll : collisions) { + initCCDB(coll); + // Event selection if (eventSelections.fUseEventSelection) { if (!IsEventAccepted(coll, true)) @@ -2842,6 +3393,36 @@ struct sigma0builder { } } } + + //_______________________________________________ + // Xi0 loop + if (fillXi0Tables) { + // First find pi0 candidates + std::vector pi0Candidates; + for (size_t i = 0; i < bestGammasArray.size(); ++i) { // loop over photons + auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); + for (size_t j = i + 1; j < bestGammasArray.size(); ++j) { // loop over photons + auto gamma2 = fullV0s.rawIteratorAt(bestGammasArray[j]); + V0PairTopoInfo pi0Info; + if (!buildPi0ForXi0(gamma1, gamma2, coll, pi0Info)) { + continue; + } + 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]); + for (size_t j = 0; j < pi0Candidates.size(); ++j) { // loop over pi0 + 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)) + continue; + } + } + } } }