diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index 521759750a7..9d728bab596 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -193,6 +193,261 @@ DECLARE_SOA_TABLE(Sigma0Cores, "AOD", "SIGMA0CORES", sigma0Core::LambdaY, sigma0Core::LambdaPhi); +// KSTAR core info +namespace kstarCore +{ +DECLARE_SOA_COLUMN(X, x, float); +DECLARE_SOA_COLUMN(Y, y, float); +DECLARE_SOA_COLUMN(Z, z, float); +DECLARE_SOA_COLUMN(DCADaughters, dcadaughters, float); + +DECLARE_SOA_COLUMN(PhotonPx, photonPx, float); +DECLARE_SOA_COLUMN(PhotonPy, photonPy, float); +DECLARE_SOA_COLUMN(PhotonPz, photonPz, float); +DECLARE_SOA_COLUMN(PhotonMass, photonMass, float); + +DECLARE_SOA_COLUMN(KShortPx, kshortPx, float); +DECLARE_SOA_COLUMN(KShortPy, kshortPy, float); +DECLARE_SOA_COLUMN(KShortPz, kshortPz, float); +DECLARE_SOA_COLUMN(KShortMass, kshortMass, float); + +// Dynamic Columns + +// KStar +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, //! KStar px + [](float photonPx, float kshortPx) -> float { return photonPx + kshortPx; }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, //! KStar py + [](float photonPy, float kshortPy) -> float { return photonPy + kshortPy; }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! KStar pz + [](float photonPz, float kshortPz) -> float { return photonPz + kshortPz; }); + +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, + [](float photonPx, float photonPy, float kshortPx, float kshortPy) -> float { + return RecoDecay::pt(array{photonPx + kshortPx, photonPy + kshortPy}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Total momentum in GeV/c + [](float photonPx, float photonPy, float photonPz, float kshortPx, float kshortPy, float kshortPz) -> float { + return RecoDecay::sqrtSumOfSquares(photonPx + kshortPx, photonPy + kshortPy, photonPz + kshortPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KStarMass, kstarMass, + [](float photonPx, float photonPy, float photonPz, float kshortPx, float kshortPy, float kshortPz) -> float { + std::array pVecPhotons{photonPx, photonPy, photonPz}; + std::array pVecKShort{kshortPx, kshortPy, kshortPz}; + auto arrMom = std::array{pVecPhotons, pVecKShort}; + return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassKaonNeutral}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KStarY, kstarY, + [](float photonPx, float photonPy, float photonPz, float kshortPx, float kshortPy, float kshortPz) -> float { + return RecoDecay::y(std::array{photonPx + kshortPx, photonPy + kshortPy, photonPz + kshortPz}, o2::constants::physics::MassK0Star892); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, //! Phi in the range [0, 2pi) + [](float photonPx, float photonPy, float kshortPx, float kshortPy) -> float { return RecoDecay::phi(photonPx + kshortPx, photonPy + kshortPy); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, //! Pseudorapidity + [](float photonPx, float photonPy, float photonPz, float kshortPx, float kshortPy, float kshortPz) -> float { + return RecoDecay::eta(std::array{photonPx + kshortPx, photonPy + kshortPy, photonPz + kshortPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(Radius, radius, //! Sigma0 decay radius (2D, centered at zero) + [](float x, float y) -> float { return RecoDecay::sqrtSumOfSquares(x, y); }); + +DECLARE_SOA_DYNAMIC_COLUMN(OPAngle, opAngle, + [](float photonPx, float photonPy, float photonPz, float kshortPx, float kshortPy, float kshortPz) { + TVector3 v1(photonPx, photonPy, photonPz); + TVector3 v2(kshortPx, kshortPy, kshortPz); + return v1.Angle(v2); + }); + +// Photon +DECLARE_SOA_DYNAMIC_COLUMN(PhotonPt, photonPt, //! Transverse momentum in GeV/c + [](float photonPx, float photonPy) -> float { + return RecoDecay::sqrtSumOfSquares(photonPx, photonPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonP, photonp, //! Total momentum in GeV/c + [](float photonPx, float photonPy, float photonPz) -> float { + return RecoDecay::sqrtSumOfSquares(photonPx, photonPy, photonPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonEta, photonEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float photonPx, float photonPy, float photonPz) -> float { + return RecoDecay::eta(std::array{photonPx, photonPy, photonPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonY, photonY, //! Rapidity + [](float photonPx, float photonPy, float photonPz) -> float { + return RecoDecay::y(std::array{photonPx, photonPy, photonPz}, o2::constants::physics::MassGamma); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonPhi, photonPhi, //! Phi in the range [0, 2pi) + [](float photonPx, float photonPy) -> float { return RecoDecay::phi(photonPx, photonPy); }); + +// KShort +DECLARE_SOA_DYNAMIC_COLUMN(KShortPt, kshortPt, //! Transverse momentum in GeV/c + [](float kshortPx, float kshortPy) -> float { + return RecoDecay::sqrtSumOfSquares(kshortPx, kshortPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortP, kshortp, //! Total momentum in GeV/c + [](float kshortPx, float kshortPy, float kshortPz) -> float { + return RecoDecay::sqrtSumOfSquares(kshortPx, kshortPy, kshortPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortEta, kshortEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float kshortPx, float kshortPy, float kshortPz) -> float { + return RecoDecay::eta(std::array{kshortPx, kshortPy, kshortPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortY, kshortY, //! Rapidity + [](float kshortPx, float kshortPy, float kshortPz) -> float { + return RecoDecay::y(std::array{kshortPx, kshortPy, kshortPz}, o2::constants::physics::MassKaonNeutral); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortPhi, kshortPhi, //! Phi in the range [0, 2pi) + [](float kshortPx, float kshortPy) -> float { return RecoDecay::phi(kshortPx, kshortPy); }); + +} // namespace kstarCore + +DECLARE_SOA_TABLE(KStarCores, "AOD", "KSTARCORES", + kstarCore::X, kstarCore::Y, kstarCore::Z, kstarCore::DCADaughters, + kstarCore::PhotonPx, kstarCore::PhotonPy, kstarCore::PhotonPz, kstarCore::PhotonMass, + kstarCore::KShortPx, kstarCore::KShortPy, kstarCore::KShortPz, kstarCore::KShortMass, + + // Dynamic columns + kstarCore::Px, + kstarCore::Py, + kstarCore::Pz, + kstarCore::Pt, + kstarCore::P, + kstarCore::KStarMass, + kstarCore::KStarY, + kstarCore::Phi, + kstarCore::Eta, + kstarCore::Radius, + kstarCore::OPAngle, + + kstarCore::PhotonPt, + kstarCore::PhotonP, + kstarCore::PhotonEta, + kstarCore::PhotonY, + kstarCore::PhotonPhi, + + kstarCore::KShortPt, + kstarCore::KShortP, + kstarCore::KShortEta, + kstarCore::KShortY, + kstarCore::KShortPhi); + +// For KStar extra info +namespace kshortExtra +{ +DECLARE_SOA_COLUMN(KShortQt, kshortQt, float); +DECLARE_SOA_COLUMN(KShortAlpha, kshortAlpha, float); +DECLARE_SOA_COLUMN(KShortLifeTime, kshortLifeTime, float); +DECLARE_SOA_COLUMN(KShortRadius, kshortRadius, float); +DECLARE_SOA_COLUMN(KShortCosPA, kshortCosPA, float); +DECLARE_SOA_COLUMN(KShortDCADau, kshortDCADau, float); +DECLARE_SOA_COLUMN(KShortDCANegPV, kshortDCANegPV, float); +DECLARE_SOA_COLUMN(KShortDCAPosPV, kshortDCAPosPV, float); +DECLARE_SOA_COLUMN(KShortPosPiTPCNSigma, kshortPosPiTPCNSigma, float); +DECLARE_SOA_COLUMN(KShortNegPiTPCNSigma, kshortNegPiTPCNSigma, float); +// // DECLARE_SOA_COLUMN(KShortPosPiTOFNSigma, kshortPosPiTOFNSigma, float); +// DECLARE_SOA_COLUMN(KShortNegPiTOFNSigma, kshortNegPiTOFNSigma, float); +DECLARE_SOA_COLUMN(KShortPiTOFNSigma, kshortPiTOFNSigma, float); +DECLARE_SOA_COLUMN(KShortPosTPCCrossedRows, kshortPosTPCCrossedRows, uint8_t); +DECLARE_SOA_COLUMN(KShortNegTPCCrossedRows, kshortNegTPCCrossedRows, uint8_t); +DECLARE_SOA_COLUMN(KShortPosEta, kshortPosEta, float); +DECLARE_SOA_COLUMN(KShortNegEta, kshortNegEta, float); +DECLARE_SOA_COLUMN(KShortPosITSCls, kshortPosITSCls, int); +DECLARE_SOA_COLUMN(KShortNegITSCls, kshortNegITSCls, int); +DECLARE_SOA_COLUMN(KShortPosITSChi2PerNcl, kshortPosChi2PerNcl, float); +DECLARE_SOA_COLUMN(KShortNegITSChi2PerNcl, kshortNegChi2PerNcl, float); +DECLARE_SOA_COLUMN(KShortPosTrackCode, kshortPosTrackCode, uint8_t); +DECLARE_SOA_COLUMN(KShortNegTrackCode, kshortNegTrackCode, uint8_t); +DECLARE_SOA_COLUMN(KShortV0Type, kshortV0Type, uint8_t); +} // namespace kshortExtra + +DECLARE_SOA_TABLE(KShortExtras, "AOD", "KSHORT", + kshortExtra::KShortQt, + kshortExtra::KShortAlpha, + kshortExtra::KShortLifeTime, + kshortExtra::KShortRadius, + kshortExtra::KShortCosPA, + kshortExtra::KShortDCADau, + kshortExtra::KShortDCANegPV, + kshortExtra::KShortDCAPosPV, + kshortExtra::KShortPosPiTPCNSigma, + kshortExtra::KShortNegPiTPCNSigma, + // kshortExtra::KShortPosPiTOFNSigma, + // kshortExtra::KShortNegPiTOFNSigma, + kshortExtra::KShortPiTOFNSigma, + kshortExtra::KShortPosTPCCrossedRows, + kshortExtra::KShortNegTPCCrossedRows, + kshortExtra::KShortPosEta, + kshortExtra::KShortNegEta, + kshortExtra::KShortPosITSCls, + kshortExtra::KShortNegITSCls, + kshortExtra::KShortPosITSChi2PerNcl, + kshortExtra::KShortNegITSChi2PerNcl, + kshortExtra::KShortPosTrackCode, + kshortExtra::KShortNegTrackCode, + kshortExtra::KShortV0Type); + +// For KStar extra info +namespace kstarPhotonExtra +{ +DECLARE_SOA_COLUMN(PhotonQt, photonQt, float); +DECLARE_SOA_COLUMN(PhotonAlpha, photonAlpha, float); +DECLARE_SOA_COLUMN(PhotonCosPA, photonCosPA, float); +DECLARE_SOA_COLUMN(PhotonDCADau, photonDCADau, float); +DECLARE_SOA_COLUMN(PhotonDCANegPV, photonDCANegPV, float); +DECLARE_SOA_COLUMN(PhotonDCAPosPV, photonDCAPosPV, float); +DECLARE_SOA_COLUMN(PhotonRadius, photonRadius, float); +DECLARE_SOA_COLUMN(PhotonZconv, photonZconv, float); +DECLARE_SOA_COLUMN(PhotonPosTPCNSigmaEl, photonPosTPCNSigmaEl, float); +DECLARE_SOA_COLUMN(PhotonNegTPCNSigmaEl, photonNegTPCNSigmaEl, float); +DECLARE_SOA_COLUMN(PhotonPosTPCCrossedRows, photonPosTPCCrossedRows, uint8_t); +DECLARE_SOA_COLUMN(PhotonNegTPCCrossedRows, photonNegTPCCrossedRows, uint8_t); +DECLARE_SOA_COLUMN(PhotonPosEta, photonPosEta, float); +DECLARE_SOA_COLUMN(PhotonNegEta, photonNegEta, float); +DECLARE_SOA_COLUMN(PhotonPsiPair, photonPsiPair, float); +DECLARE_SOA_COLUMN(PhotonPosITSCls, photonPosITSCls, int); +DECLARE_SOA_COLUMN(PhotonNegITSCls, photonNegITSCls, int); +DECLARE_SOA_COLUMN(PhotonPosITSChi2PerNcl, photonPosITSChi2PerNcl, float); +DECLARE_SOA_COLUMN(PhotonNegITSChi2PerNcl, photonNegITSChi2PerNcl, float); +DECLARE_SOA_COLUMN(PhotonPosTrackCode, photonPosTrackCode, uint8_t); +DECLARE_SOA_COLUMN(PhotonNegTrackCode, photonNegTrackCode, uint8_t); +DECLARE_SOA_COLUMN(PhotonV0Type, photonV0Type, uint8_t); +} // namespace kstarPhotonExtra + +DECLARE_SOA_TABLE(KStarPhotonExtras, "AOD", "KSTARPHOTON", + kstarPhotonExtra::PhotonQt, + kstarPhotonExtra::PhotonAlpha, + kstarPhotonExtra::PhotonCosPA, + kstarPhotonExtra::PhotonDCADau, + kstarPhotonExtra::PhotonDCANegPV, + kstarPhotonExtra::PhotonDCAPosPV, + kstarPhotonExtra::PhotonRadius, + kstarPhotonExtra::PhotonZconv, + kstarPhotonExtra::PhotonPosTPCNSigmaEl, + kstarPhotonExtra::PhotonNegTPCNSigmaEl, + kstarPhotonExtra::PhotonPosTPCCrossedRows, + kstarPhotonExtra::PhotonNegTPCCrossedRows, + kstarPhotonExtra::PhotonPosEta, + kstarPhotonExtra::PhotonNegEta, + kstarPhotonExtra::PhotonPsiPair, + kstarPhotonExtra::PhotonPosITSCls, + kstarPhotonExtra::PhotonNegITSCls, + kstarPhotonExtra::PhotonPosITSChi2PerNcl, + kstarPhotonExtra::PhotonNegITSChi2PerNcl, + kstarPhotonExtra::PhotonPosTrackCode, + kstarPhotonExtra::PhotonNegTrackCode, + kstarPhotonExtra::PhotonV0Type); + // For Photon extra info namespace sigma0PhotonExtra { @@ -334,11 +589,11 @@ DECLARE_SOA_COLUMN(LambdaPDGCode, lambdaPDGCode, int); DECLARE_SOA_COLUMN(LambdaPDGCodeMother, lambdaPDGCodeMother, int); DECLARE_SOA_COLUMN(LambdaIsCorrectlyAssoc, lambdaIsCorrectlyAssoc, bool); -DECLARE_SOA_DYNAMIC_COLUMN(IsSigma0, isSigma0, //! IsSigma0 - [](int pdgCode) -> bool { return pdgCode == 3212; }); +DECLARE_SOA_DYNAMIC_COLUMN(IsSigma0, isSigma0, //! IsSigma0 + [](int pdgCode) -> bool { return pdgCode == PDG_t::kSigma0; }); // 3212 -DECLARE_SOA_DYNAMIC_COLUMN(IsAntiSigma0, isAntiSigma0, //! IsASigma0 - [](int pdgCode) -> bool { return pdgCode == -3212; }); +DECLARE_SOA_DYNAMIC_COLUMN(IsAntiSigma0, isAntiSigma0, //! IsASigma0 + [](int pdgCode) -> bool { return pdgCode == PDG_t::kSigma0Bar; }); //-3212 DECLARE_SOA_DYNAMIC_COLUMN(MCPx, mcpx, //! Sigma0 px [](float photonMCPx, float lambdaMCPx) -> float { return photonMCPx + lambdaMCPx; }); @@ -477,6 +732,166 @@ namespace sigma0MCCore { DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); //! MC particle for Sigma0 } + +// for MC +namespace kstarMCCore +{ +DECLARE_SOA_COLUMN(MCradius, mcradius, float); +DECLARE_SOA_COLUMN(PDGCode, pdgCode, int); +DECLARE_SOA_COLUMN(PDGCodeMother, pdgCodeMother, int); +DECLARE_SOA_COLUMN(MCprocess, mcprocess, int); +DECLARE_SOA_COLUMN(IsProducedByGenerator, isProducedByGenerator, bool); + +DECLARE_SOA_COLUMN(PhotonMCPx, photonmcpx, float); +DECLARE_SOA_COLUMN(PhotonMCPy, photonmcpy, float); +DECLARE_SOA_COLUMN(PhotonMCPz, photonmcpz, float); +DECLARE_SOA_COLUMN(IsPhotonPrimary, isPhotonPrimary, bool); +DECLARE_SOA_COLUMN(PhotonPDGCode, photonPDGCode, int); +DECLARE_SOA_COLUMN(PhotonPDGCodeMother, photonPDGCodeMother, int); +DECLARE_SOA_COLUMN(PhotonIsCorrectlyAssoc, photonIsCorrectlyAssoc, bool); + +DECLARE_SOA_COLUMN(KShortMCPx, kshortmcpx, float); +DECLARE_SOA_COLUMN(KShortMCPy, kshortmcpy, float); +DECLARE_SOA_COLUMN(KShortMCPz, kshortmcpz, float); +DECLARE_SOA_COLUMN(IsKShortPrimary, isKShortPrimary, bool); +DECLARE_SOA_COLUMN(KShortPDGCode, kshortPDGCode, int); +DECLARE_SOA_COLUMN(KShortPDGCodeMother, kshortPDGCodeMother, int); +DECLARE_SOA_COLUMN(KShortIsCorrectlyAssoc, kshortIsCorrectlyAssoc, bool); + +DECLARE_SOA_DYNAMIC_COLUMN(IsKStar, isKStar, //! IsSigma0 + [](int pdgCode) -> bool { return pdgCode == o2::constants::physics::Pdg::kK0Star892; }); // 313 + +DECLARE_SOA_DYNAMIC_COLUMN(MCPx, mcpx, //! Sigma0 px + [](float photonMCPx, float kshortMCPx) -> float { return photonMCPx + kshortMCPx; }); +DECLARE_SOA_DYNAMIC_COLUMN(MCPy, mcpy, //! Sigma0 py + [](float photonMCPy, float kshortMCPy) -> float { return photonMCPy + kshortMCPy; }); +DECLARE_SOA_DYNAMIC_COLUMN(MCPz, mcpz, //! Sigma0 pz + [](float photonMCPz, float kshortMCPz) -> float { return photonMCPz + kshortMCPz; }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCPt, mcpt, + [](float photonMCPx, float photonMCPy, float kshortMCPx, float kshortMCPy) -> float { + return RecoDecay::pt(array{photonMCPx + kshortMCPx, photonMCPy + kshortMCPy}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCP, mcp, //! Total momentum in GeV/c + [](float photonMCPx, float photonMCPy, float photonMCPz, float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + return RecoDecay::sqrtSumOfSquares(photonMCPx + kshortMCPx, photonMCPy + kshortMCPy, photonMCPz + kshortMCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KStarMCMass, kstarMCMass, + [](float photonMCPx, float photonMCPy, float photonMCPz, float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + std::array pVecPhotons{photonMCPx, photonMCPy, photonMCPz}; + std::array pVecKShort{kshortMCPx, kshortMCPy, kshortMCPz}; + auto arrMom = std::array{pVecPhotons, pVecKShort}; + return RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassKaonNeutral}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KStarMCY, kstarMCY, + [](float photonMCPx, float photonMCPy, float photonMCPz, float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + return RecoDecay::y(std::array{photonMCPx + kshortMCPx, photonMCPy + kshortMCPy, photonMCPz + kshortMCPz}, o2::constants::physics::MassKaonNeutral); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcphi, //! Phi in the range [0, 2pi) + [](float photonMCPx, float photonMCPy, float kshortMCPx, float kshortMCPy) -> float { return RecoDecay::phi(photonMCPx + kshortMCPx, photonMCPy + kshortMCPy); }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCEta, mceta, //! Pseudorapidity + [](float photonMCPx, float photonMCPy, float photonMCPz, float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + return RecoDecay::eta(std::array{photonMCPx + kshortMCPx, photonMCPy + kshortMCPy, photonMCPz + kshortMCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(MCOPAngle, mcopAngle, + [](float photonMCPx, float photonMCPy, float photonMCPz, float kshortMCPx, float kshortMCPy, float kshortMCPz) { + TVector3 v1(photonMCPx, photonMCPy, photonMCPz); + TVector3 v2(kshortMCPx, kshortMCPy, kshortMCPz); + return v1.Angle(v2); + }); + +// Photon +DECLARE_SOA_DYNAMIC_COLUMN(PhotonMCPt, photonmcpt, //! Transverse momentum in GeV/c + [](float photonMCPx, float photonMCPy) -> float { + return RecoDecay::sqrtSumOfSquares(photonMCPx, photonMCPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonMCP, photonmcp, //! Total momentum in GeV/c + [](float photonMCPx, float photonMCPy, float photonMCPz) -> float { + return RecoDecay::sqrtSumOfSquares(photonMCPx, photonMCPy, photonMCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonMCEta, photonMCEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float photonMCPx, float photonMCPy, float photonMCPz) -> float { + return RecoDecay::eta(std::array{photonMCPx, photonMCPy, photonMCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonMCY, photonMCY, //! Rapidity + [](float photonMCPx, float photonMCPy, float photonMCPz) -> float { + return RecoDecay::y(std::array{photonMCPx, photonMCPy, photonMCPz}, o2::constants::physics::MassGamma); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(PhotonMCPhi, photonMCPhi, //! Phi in the range [0, 2pi) + [](float photonMCPx, float photonMCPy) -> float { return RecoDecay::phi(photonMCPx, photonMCPy); }); + +// KShort +DECLARE_SOA_DYNAMIC_COLUMN(KShortMCPt, kshortmcpt, //! Transverse momentum in GeV/c + [](float kshortMCPx, float kshortMCPy) -> float { + return RecoDecay::sqrtSumOfSquares(kshortMCPx, kshortMCPy); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortMCP, kshortmcp, //! Total momentum in GeV/c + [](float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + return RecoDecay::sqrtSumOfSquares(kshortMCPx, kshortMCPy, kshortMCPz); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortMCEta, kshortMCEta, //! Pseudorapidity, conditionally defined to avoid FPEs + [](float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + return RecoDecay::eta(std::array{kshortMCPx, kshortMCPy, kshortMCPz}); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortMCY, kshortMCY, //! Rapidity + [](float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { + return RecoDecay::y(std::array{kshortMCPx, kshortMCPy, kshortMCPz}, o2::constants::physics::MassKaonNeutral); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(KShortMCPhi, kshortMCPhi, //! Phi in the range [0, 2pi) + [](float kshortMCPx, float kshortMCPy) -> float { return RecoDecay::phi(kshortMCPx, kshortMCPy); }); + +} // namespace kstarMCCore + +DECLARE_SOA_TABLE(KStarMCCores, "AOD", "KSTARMCCORES", + // Basic properties + kstarMCCore::MCradius, kstarMCCore::PDGCode, kstarMCCore::PDGCodeMother, kstarMCCore::MCprocess, kstarMCCore::IsProducedByGenerator, + + kstarMCCore::PhotonMCPx, kstarMCCore::PhotonMCPy, kstarMCCore::PhotonMCPz, + kstarMCCore::IsPhotonPrimary, kstarMCCore::PhotonPDGCode, kstarMCCore::PhotonPDGCodeMother, kstarMCCore::PhotonIsCorrectlyAssoc, + + kstarMCCore::KShortMCPx, kstarMCCore::KShortMCPy, kstarMCCore::KShortMCPz, + kstarMCCore::IsKShortPrimary, kstarMCCore::KShortPDGCode, kstarMCCore::KShortPDGCodeMother, kstarMCCore::KShortIsCorrectlyAssoc, + + // Dynamic columns + kstarMCCore::IsKStar, + + kstarMCCore::MCPx, + kstarMCCore::MCPy, + kstarMCCore::MCPz, + kstarMCCore::MCPt, + kstarMCCore::MCP, + kstarMCCore::KStarMCMass, + kstarMCCore::KStarMCY, + kstarMCCore::MCPhi, + kstarMCCore::MCEta, + kstarMCCore::MCOPAngle, + + kstarMCCore::PhotonMCPt, + kstarMCCore::PhotonMCP, + kstarMCCore::PhotonMCEta, + kstarMCCore::PhotonMCY, + kstarMCCore::PhotonMCPhi, + + kstarMCCore::KShortMCPt, + kstarMCCore::KShortMCP, + kstarMCCore::KShortMCEta, + kstarMCCore::KShortMCY, + kstarMCCore::KShortMCPhi); + namespace sigma0Gen { DECLARE_SOA_COLUMN(IsSigma0, isSigma0, bool); // true: sigma0, false: antisigma0 @@ -492,6 +907,19 @@ DECLARE_SOA_TABLE(Sigma0Gens, "AOD", "SIGMA0GENS", sigma0Gen::MCPt, sigma0Gen::MCY); +namespace kstarGen +{ +DECLARE_SOA_COLUMN(IsKStar, isKStar, bool); // true: sigma0, false: antisigma0 +DECLARE_SOA_COLUMN(ProducedByGenerator, producedByGenerator, bool); +DECLARE_SOA_COLUMN(KStarMCPt, kstarMCPt, float); // MC pT + +} // namespace kstarGen + +DECLARE_SOA_TABLE(KStarGens, "AOD", "KSTARGENS", + kstarGen::IsKStar, + kstarGen::ProducedByGenerator, + kstarGen::KStarMCPt); + DECLARE_SOA_TABLE(SigmaCollRef, "AOD", "SIGMACOLLREF", //! optional table to refer back to a collision o2::soa::Index<>, v0data::StraCollisionId); @@ -504,6 +932,12 @@ DECLARE_SOA_TABLE(SigmaMCLabels, "AOD", "SIGMAMCLABEL", //! optional table to re DECLARE_SOA_TABLE(SigmaGenCollRef, "AOD", "SIGMAGENCOLLREF", //! optional table to refer back to a collision o2::soa::Index<>, v0data::StraMCCollisionId); +DECLARE_SOA_TABLE(KStarCollRef, "AOD", "KSTARCOLLREF", + o2::soa::Index<>, v0data::StraCollisionId); + +DECLARE_SOA_TABLE(KStarGenCollRef, "AOD", "KSTARGENCOLLREF", + o2::soa::Index<>, v0data::StraMCCollisionId); + // ___________________________________________________________________________ // pi0 QA namespace Pi0Core @@ -719,11 +1153,11 @@ DECLARE_SOA_COLUMN(Photon2PDGCode, photon2PDGCode, int); DECLARE_SOA_COLUMN(Photon2PDGCodeMother, photon2PDGCodeMother, int); DECLARE_SOA_COLUMN(Photon2IsCorrectlyAssoc, photon2IsCorrectlyAssoc, bool); -DECLARE_SOA_DYNAMIC_COLUMN(IsPi0, isPi0, //! IsPi0 - [](int pdgCode) -> bool { return pdgCode == 111; }); +DECLARE_SOA_DYNAMIC_COLUMN(IsPi0, isPi0, //! IsPi0 + [](int pdgCode) -> bool { return pdgCode == PDG_t::kPi0; }); // 111 -DECLARE_SOA_DYNAMIC_COLUMN(IsFromXi0, isFromXi0, //! Pi0 from Xi0 - [](int pdgCodeMother) -> bool { return pdgCodeMother == 3322; }); +DECLARE_SOA_DYNAMIC_COLUMN(IsFromXi0, isFromXi0, //! Pi0 from Xi0 + [](int pdgCodeMother) -> bool { return pdgCodeMother == o2::constants::physics::Pdg::kXi0; }); // 3322 DECLARE_SOA_DYNAMIC_COLUMN(MCPx, mcpx, //! Pi0 MC px [](float photon1MCPx, float photon2MCPx) -> float { return photon1MCPx + photon2MCPx; }); diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index dcbf8f8f31c..973fef69a8b 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -74,6 +74,17 @@ struct sigma0builder { Service ccdb; ctpRateFetcher rateFetcher; + //___________________________________________________ + // KStar Specific + + Produces kstarcores; // kstar candidates info for analysis + Produces kshortExtras; // lambdas from sigma0 candidates info + Produces kstarPhotonExtras; // photons from kstar candidates info + Produces kstarCollRefs; // references collisions from kstarcores + Produces kstarmccores; // Reco sigma0 MC properties + Produces kstarGens; // Generated sigma0s + Produces kstarGenCollRefs; // references collisions from sigma0Gens + //__________________________________________________ // Sigma0 specific Produces sigmaIndices; // references V0Cores from sigma0Gens @@ -83,7 +94,6 @@ struct sigma0builder { Produces sigma0CollRefs; // references collisions from Sigma0Cores Produces sigma0mccores; // Reco sigma0 MC properties Produces sigma0mclabel; // Link of reco sigma0 to mcparticles - Produces sigma0Gens; // Generated sigma0s Produces sigma0GenCollRefs; // references collisions from sigma0Gens[ @@ -156,6 +166,7 @@ struct sigma0builder { // Tables to fill 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"}; // For ML Selection Configurable useMLScores{"useMLScores", false, "use ML scores to select candidates"}; @@ -221,6 +232,20 @@ struct sigma0builder { Configurable PhotonPhiMax2{"PhotonPhiMax2", -1, "Phi min value to reject photons, region 2 (leave negative if no selection desired)"}; } photonSelections; + /// KShort criteria: + Configurable V0Rapidity{"V0Rapidity", 0.5, "v0 rapidity"}; + Configurable K0ShortDauPseudoRap{"K0SDauPseudoRap", 1.5, "Max pseudorapidity of daughter tracks"}; + Configurable K0ShortMinDCANegToPv{"K0SMinDCANegToPv", 0.0, "min DCA Neg To PV (cm)"}; + Configurable K0ShortMinDCAPosToPv{"K0SMinDCAPosToPv", 0.0, "min DCA Pos To PV (cm)"}; + Configurable K0ShortMaxDCAV0Dau{"K0SMaxDCAV0Dau", 3.5, "Max DCA V0 Daughters (cm)"}; + Configurable K0ShortMinv0radius{"K0SMinv0radius", 0.0, "Min V0 radius (cm)"}; + Configurable K0ShortMaxv0radius{"K0SMaxv0radius", 60, "Max V0 radius (cm)"}; + Configurable K0ShortWindow{"K0SWindow", 0.1, "Mass window around expected (in GeV/c2)"}; + + // KStar criteria: + Configurable KStarWindow{"KStarWindow", 0.1, "Mass window around expected (in GeV/c2)"}; + Configurable KStarMaxRap{"KStarMaxRap", 0.8, "Max kstar rapidity"}; + //// Sigma0 criteria: Configurable Sigma0Window{"Sigma0Window", 0.1, "Mass window around expected (in GeV/c2)"}; Configurable SigmaMaxRap{"SigmaMaxRap", 0.8, "Max sigma0 rapidity"}; @@ -251,7 +276,7 @@ struct sigma0builder { 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})"}; // AP plot axes ConfigurableAxis axisAPAlpha{"axisAPAlpha", {220, -1.1f, 1.1f}, "V0 AP alpha"}; ConfigurableAxis axisAPQt{"axisAPQt", {220, 0.0f, 0.5f}, "V0 AP alpha"}; @@ -330,6 +355,31 @@ struct sigma0builder { } } + histos.add("K0ShortSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); + histos.get(HIST("K0ShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); + histos.get(HIST("K0ShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "K0Short Mass Cut"); + histos.get(HIST("K0ShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "K0Short Eta/Y Cut"); + histos.get(HIST("K0ShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(4, "K0Short DCAToPV Cut"); + histos.get(HIST("K0ShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(5, "K0Short Radius Cut"); + histos.get(HIST("K0ShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(6, "K0Short DCADau Cut"); + + histos.add("K0ShortSel/hK0ShortMass", "hK0ShortMass", kTH1F, {axisK0SMass}); + histos.add("K0ShortSel/hK0ShortNegEta", "hK0ShortNegEta", kTH1F, {axisRapidity}); + histos.add("K0ShortSel/hK0ShortPosEta", "hK0ShortPosEta", kTH1F, {axisRapidity}); + histos.add("K0ShortSel/hK0ShortY", "hK0ShortY", kTH1F, {axisRapidity}); + histos.add("K0ShortSel/hK0ShortDCANegToPV", "hK0ShortDCANegToPV", kTH1F, {axisDCAtoPV}); + histos.add("K0ShortSel/hK0ShortDCAPosToPV", "hK0ShortDCAPosToPV", kTH1F, {axisDCAtoPV}); + histos.add("K0ShortSel/hK0ShortDCADau", "hK0ShortDCADau", kTH1F, {axisDCAdau}); + histos.add("K0ShortSel/hK0ShortRadius", "hK0ShortRadius", kTH1F, {axisRadius}); + histos.add("K0ShortSel/h3dK0ShortMass", "h3dK0ShortMass", kTH3D, {axisCentrality, axisPt, axisK0SMass}); + + histos.add("KStarSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); + histos.get(HIST("KStarSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); + histos.get(HIST("KStarSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "KStar Mass Window"); + histos.get(HIST("KStarSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "KStar Y Window"); + + histos.add("KStarSel/hKStarMassSelected", "hKStarMassSelected", kTH1F, {axisKStarMass}); + for (const auto& histodir : DirList) { if ((histodir == "V0BeforeSel" && !fFillNoSelV0Histos) || (histodir == "PhotonSel" && !fFillSelPhotonHistos) || @@ -494,6 +544,35 @@ struct sigma0builder { h2DGenSigma0TypeVsProducedByGen->GetXaxis()->SetBinLabel(2, "Non-Sterile"); h2DGenSigma0TypeVsProducedByGen->GetYaxis()->SetBinLabel(1, "Generator"); h2DGenSigma0TypeVsProducedByGen->GetYaxis()->SetBinLabel(2, "Transport"); + + // ______________________________________________________ + // KStar + histos.add("GenQA/hGenKStar", "hGenKStar", kTH1D, {axisPt}); + + histos.add("GenQA/h2dGenKStarxy_Generator", "hGenKStarxy_Generator", kTH2D, {axisXY, axisXY}); + histos.add("GenQA/h2dGenKStarxy_Transport", "hGenKStarxy_Transport", kTH2D, {axisXY, axisXY}); + histos.add("GenQA/hGenKStarRadius_Generator", "hGenKStarRadius_Generator", kTH1D, {axisRadius}); + histos.add("GenQA/hGenKStarRadius_Transport", "hGenKStarRadius_Transport", kTH1D, {axisRadius}); + + histos.add("GenQA/h2dKStarMCSourceVsPDGMother", "h2dKStarMCSourceVsPDGMother", kTHnSparseD, {{2, -0.5f, 1.5f}, {10001, -5000.5f, +5000.5f}}); + histos.add("GenQA/h2dKStarNDaughtersVsPDG", "h2dKStarNDaughtersVsPDG", kTHnSparseD, {{10, -0.5f, +9.5f}, {10001, -5000.5f, +5000.5f}}); + + auto hPrimaryKStars = histos.add("GenQA/hPrimaryKStars", "hPrimaryKStars", kTH1D, {{2, -0.5f, 1.5f}}); + hPrimaryKStars->GetXaxis()->SetBinLabel(1, "All KStars"); + hPrimaryKStars->GetXaxis()->SetBinLabel(2, "Primary KStars"); + + auto hGenSpeciesKStar = histos.add("GenQA/hGenSpeciesKStar", "hGenSpeciesKStar", kTH1D, {{4, -0.5f, 3.5f}}); + hGenSpeciesKStar->GetXaxis()->SetBinLabel(1, "All Prim. KShort"); + hGenSpeciesKStar->GetXaxis()->SetBinLabel(5, "All KStars"); + + histos.add("GenQA/hKStarNDau", "hKStarNDau", kTH1D, {{10, -0.5f, +9.5f}}); + histos.add("GenQA/h2dKStarNDauVsProcess", "h2dKStarNDauVsProcess", kTH2D, {{10, -0.5f, +9.5f}, {50, -0.5f, 49.5f}}); + + auto h2DGenKStarTypeVsProducedByGen = histos.add("GenQA/h2DGenKStarTypeVsProducedByGen", "h2DGenKStarTypeVsProducedByGen", kTH2D, {{2, -0.5f, 1.5f}, {2, -0.5f, 1.5f}}); + h2DGenKStarTypeVsProducedByGen->GetXaxis()->SetBinLabel(1, "Sterile"); + h2DGenKStarTypeVsProducedByGen->GetXaxis()->SetBinLabel(2, "Non-Sterile"); + h2DGenKStarTypeVsProducedByGen->GetYaxis()->SetBinLabel(1, "Generator"); + h2DGenKStarTypeVsProducedByGen->GetYaxis()->SetBinLabel(2, "Transport"); } if (doprocessPhotonLambdaQA || doprocessPhotonLambdaMCQA) { @@ -598,9 +677,11 @@ struct sigma0builder { bool IsPrimary = false; bool IsV0Lambda = false; bool IsV0AntiLambda = false; + bool IsV0KShort = false; bool IsPi0 = false; bool IsSigma0 = false; bool IsAntiSigma0 = false; + bool IsKStar = false; bool IsProducedByGenerator = false; bool IsSterile = false; int MCProcess = -1; @@ -674,7 +755,7 @@ struct sigma0builder { // Sanity check: Is V0Pair <-> Mother assignment correct? bool fIsSigma0 = false; - if ((v01MC.pdgCode() == 22) && (v01MC.pdgCodeMother() == 3212) && (v02MC.pdgCode() == 3122) && (v02MC.pdgCodeMother() == 3212) && (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; // Check collision assignment @@ -728,7 +809,7 @@ struct sigma0builder { // MC QA histograms // Parenthood check for sigma0-like candidate - if (MCParticle_v01.pdgCode() == 22 && TMath::Abs(MCParticle_v02.pdgCode()) == 3122) { + if (MCParticle_v01.pdgCode() == PDG_t::kGamma && TMath::Abs(MCParticle_v02.pdgCode()) == PDG_t::kLambda0) { for (const auto& mother1 : MCMothersList_v01) { // Photon mothers histos.fill(HIST("MCQA/h2dPhotonNMothersVsPDG"), MCMothersList_v01.size(), mother1.pdgCode()); histos.fill(HIST("MCQA/h2dPhotonNMothersVsMCProcess"), MCMothersList_v01.size(), mother1.getProcess()); @@ -753,11 +834,11 @@ struct sigma0builder { } // Check association correctness - if (fIsSigma0 && (MCinfo.V0PairPDGCode == 3212)) + if (fIsSigma0 && (MCinfo.V0PairPDGCode == PDG_t::kSigma0)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 1); // match - if (fIsSigma0 && !(MCinfo.V0PairPDGCode == 3212)) + if (fIsSigma0 && !(MCinfo.V0PairPDGCode == PDG_t::kSigma0)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 2); // mismatch - if (!fIsSigma0 && (MCinfo.V0PairPDGCode == 3212)) + if (!fIsSigma0 && (MCinfo.V0PairPDGCode == PDG_t::kSigma0)) histos.fill(HIST("MCQA/hSigma0MCCheck"), 3); // mismatch } } @@ -1080,7 +1161,7 @@ struct sigma0builder { bool fIsV0CorrectlyAssigned = (v0MC.straMCCollisionId() == v0MCCollision.globalIndex()); bool isPrimary = v0MC.isPhysicalPrimary(); - if ((v0MC.pdgCode() == 22) && isPhotonAnalysis && isPrimary) { // True Gamma + if ((v0MC.pdgCode() == PDG_t::kGamma) && isPhotonAnalysis && isPrimary) { // True Gamma histos.fill(HIST("V0AssoQA/h2dIRVsPt_TrueGamma"), IR, V0MCpT); histos.fill(HIST("V0AssoQA/h3dPAVsIRVsPt_TrueGamma"), V0PA, IR, V0MCpT); @@ -1089,7 +1170,7 @@ struct sigma0builder { histos.fill(HIST("V0AssoQA/h3dPAVsIRVsPt_TrueGamma_BadCollAssig"), V0PA, IR, V0MCpT); } } - if ((v0MC.pdgCode() == 3122) && !isPhotonAnalysis && isPrimary) { // True Lambda + if ((v0MC.pdgCode() == PDG_t::kLambda0) && !isPhotonAnalysis && isPrimary) { // True Lambda histos.fill(HIST("V0AssoQA/h2dIRVsPt_TrueLambda"), IR, V0MCpT); histos.fill(HIST("V0AssoQA/h3dPAVsIRVsPt_TrueLambda"), V0PA, IR, V0MCpT); @@ -1108,11 +1189,13 @@ struct sigma0builder { // Fill with properties GenInfo.IsPrimary = mcParticle.isPhysicalPrimary(); - GenInfo.IsV0Lambda = mcParticle.pdgCode() == 3122; - GenInfo.IsV0AntiLambda = mcParticle.pdgCode() == -3122; - GenInfo.IsPi0 = mcParticle.pdgCode() == 111; - GenInfo.IsSigma0 = mcParticle.pdgCode() == 3212; - GenInfo.IsAntiSigma0 = mcParticle.pdgCode() == -3212; + GenInfo.IsV0Lambda = mcParticle.pdgCode() == PDG_t::kLambda0; // 3122 + GenInfo.IsV0AntiLambda = mcParticle.pdgCode() == PDG_t::kLambda0Bar; //-3122 + GenInfo.IsV0KShort = mcParticle.pdgCode() == PDG_t::kK0Short; // 310 + 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 = mcParticle.pdgCode() == o2::constants::physics::Pdg::kK0Star892; // 313; GenInfo.IsProducedByGenerator = mcParticle.producedByGenerator(); GenInfo.MCProcess = mcParticle.getProcess(); GenInfo.MCPt = mcParticle.pt(); @@ -1123,7 +1206,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) { + 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(); @@ -1141,15 +1224,21 @@ struct sigma0builder { histos.fill(HIST("GenQA/h2dSigma0NDaughtersVsPDG"), daughters.size(), daughter.pdgCode()); if (GenInfo.NDaughters == 2) { - if (daughter.pdgCode() == 22) + if (daughter.pdgCode() == PDG_t::kGamma) GenInfo.MCDau1Pt = daughter.pt(); - if (TMath::Abs(daughter.pdgCode()) == 3122) + if (TMath::Abs(daughter.pdgCode()) == PDG_t::kLambda0) GenInfo.MCDau2Pt = daughter.pt(); } } } + if ((GenInfo.IsKStar) && genSelections.doQA) { + histos.fill(HIST("GenQA/h2dKStarMCSourceVsPDGMother"), GenInfo.IsProducedByGenerator, GenInfo.PDGCodeMother); + for (auto& daughter : daughters) // checking decay modes + histos.fill(HIST("GenQA/h2dKStarNDaughtersVsPDG"), daughters.size(), daughter.pdgCode()); + } + if (GenInfo.IsPi0 && genSelections.doQA) { histos.fill(HIST("GenQA/h2dPi0MCSourceVsPDGMother"), GenInfo.IsProducedByGenerator, GenInfo.PDGCodeMother); for (auto& daughter : daughters) // checking decay modes @@ -1186,6 +1275,8 @@ struct sigma0builder { histos.fill(HIST("GenQA/hGenSpecies"), 0); if (GenInfo.IsV0AntiLambda && GenInfo.IsPrimary) histos.fill(HIST("GenQA/hGenSpecies"), 1); + if (GenInfo.IsV0KShort && GenInfo.IsPrimary) + histos.fill(HIST("GenQA/hGenSpeciesKStar"), 0); // Checking decay mode if (GenInfo.IsSigma0 || GenInfo.IsAntiSigma0) { @@ -1223,6 +1314,36 @@ struct sigma0builder { histos.fill(HIST("GenQA/h3dGenASigma0_pTMap"), GenInfo.MCPt, GenInfo.MCDau1Pt, GenInfo.MCDau2Pt); } } + + // Checking decay mode + if (GenInfo.IsKStar) { + histos.fill(HIST("GenQA/hKStarNDau"), GenInfo.NDaughters); + histos.fill(HIST("GenQA/h2dKStarNDauVsProcess"), GenInfo.NDaughters, GenInfo.MCProcess); + + const auto radius = std::hypot(GenInfo.MCvx, GenInfo.MCvy); + // KStar XY and radius (separate histos for Gen/Transport) + if (GenInfo.IsProducedByGenerator) { + histos.fill(HIST("GenQA/h2dGenKStarxy_Generator"), GenInfo.MCvx, GenInfo.MCvy); + histos.fill(HIST("GenQA/hGenKStarRadius_Generator"), radius); + } else { + histos.fill(HIST("GenQA/h2dGenKStarxy_Transport"), GenInfo.MCvx, GenInfo.MCvy); + histos.fill(HIST("GenQA/hGenKStarRadius_Transport"), radius); + } + + // kstar type vs origin (single 2D histo) + const int genIndex = GenInfo.IsProducedByGenerator ? 0 : 1; // 0 = Generator, 1 = Transport + const int typeIndex = GenInfo.IsSterile ? 0 : 1; // 0 = Sterile, 1 = Normal + histos.fill(HIST("GenQA/h2DGenKStarTypeVsProducedByGen"), typeIndex, genIndex); + + // Fill histograms + if (GenInfo.IsKStar) { + histos.fill(HIST("GenQA/hGenSpeciesKStar"), 2); + histos.fill(HIST("GenQA/hGenKStar"), GenInfo.MCPt); + histos.fill(HIST("GenQA/hPrimaryKStars"), 0); + if (GenInfo.IsPrimary) + histos.fill(HIST("GenQA/hPrimaryKStars"), 1); + } + } } // ______________________________________________________ @@ -1269,6 +1390,12 @@ struct sigma0builder { sigma0Gens(MCGenInfo.IsSigma0, MCGenInfo.IsProducedByGenerator, MCGenInfo.MCPt, mcParticle.y()); sigma0GenCollRefs(MCGenInfo.MCCollId); // link to stramccollision table } + + // KStar + if (fillKStarTables && MCGenInfo.IsKStar) { + kstarGens(MCGenInfo.IsKStar, MCGenInfo.IsProducedByGenerator, MCGenInfo.MCPt); + kstarGenCollRefs(MCGenInfo.MCCollId); // link to stramccollision table + } } } @@ -1339,7 +1466,7 @@ struct sigma0builder { if constexpr (requires { gamma.motherMCPartId(); }) { if (gamma.has_v0MCCore()) { auto gammaMC = gamma.template v0MCCore_as>(); - if (gammaMC.pdgCode() != 22) + if (gammaMC.pdgCode() != PDG_t::kGamma) return false; } } @@ -1436,7 +1563,7 @@ struct sigma0builder { if constexpr (requires { lambda.motherMCPartId(); }) { if (lambda.has_v0MCCore()) { auto lambdaMC = lambda.template v0MCCore_as>(); - if (TMath::Abs(lambdaMC.pdgCode()) != 3122) + if (TMath::Abs(lambdaMC.pdgCode()) != PDG_t::kLambda0) return false; } } @@ -1522,6 +1649,58 @@ struct sigma0builder { return true; } + template + bool processKShortCandidate(TV0Object const& kshort, TCollision const& collision) + { + + if (kshort.v0Type() != 1) + return false; + + float centrality = doPPAnalysis ? collision.centFT0M() : collision.centFT0C(); + + histos.fill(HIST("K0ShortSel/h3dK0ShortMass"), centrality, kshort.pt(), kshort.mK0Short()); + + histos.fill(HIST("K0ShortSel/hSelectionStatistics"), 1.); + histos.fill(HIST("K0ShortSel/hK0ShortMass"), kshort.mK0Short()); + + if ((TMath::Abs(kshort.mK0Short() - o2::constants::physics::MassK0Short) > K0ShortWindow)) + return false; + histos.fill(HIST("K0ShortSel/hK0ShortNegEta"), kshort.negativeeta()); + histos.fill(HIST("K0ShortSel/hK0ShortPosEta"), kshort.positiveeta()); + histos.fill(HIST("K0ShortSel/hK0ShortY"), kshort.yK0Short()); + histos.fill(HIST("K0ShortSel/hSelectionStatistics"), 2.); + + if ((TMath::Abs(kshort.yK0Short()) > V0Rapidity) || (TMath::Abs(kshort.negativeeta()) > K0ShortDauPseudoRap) || (TMath::Abs(kshort.positiveeta()) > K0ShortDauPseudoRap)) + return false; + histos.fill(HIST("K0ShortSel/hK0ShortDCANegToPV"), kshort.dcanegtopv()); + histos.fill(HIST("K0ShortSel/hK0ShortDCAPosToPV"), kshort.dcapostopv()); + histos.fill(HIST("K0ShortSel/hSelectionStatistics"), 3.); + + if ((TMath::Abs(kshort.dcapostopv()) < K0ShortMinDCAPosToPv) || (TMath::Abs(kshort.dcanegtopv()) < K0ShortMinDCANegToPv)) + return false; + histos.fill(HIST("K0ShortSel/hK0ShortRadius"), kshort.v0radius()); + histos.fill(HIST("K0ShortSel/hSelectionStatistics"), 4.); + + if ((kshort.v0radius() < K0ShortMinv0radius) || (kshort.v0radius() > K0ShortMaxv0radius)) + return false; + histos.fill(HIST("K0ShortSel/hK0ShortDCADau"), kshort.dcaV0daughters()); + histos.fill(HIST("K0ShortSel/hSelectionStatistics"), 5.); + + if (TMath::Abs(kshort.dcaV0daughters()) > K0ShortMaxDCAV0Dau) + return false; + histos.fill(HIST("K0ShortSel/hSelectionStatistics"), 6.); + // MC Processing (if available) + // if constexpr (requires { kshort.motherMCPartId(); }) { + // if (kshort.has_v0MCCore()) { + // auto kshortMC = kshort.template v0MCCore_as>(); + // // if (kshortMC.pdgCode() == 310) // Is K0Short + // // histos.fill(HIST("MC/h2dPtVsCentrality_MCAssocKShort"), centrality, kshort.pt()); + // } + // } + + return true; + } + //_______________________________________________ // Build pi0 candidate for QA template @@ -1712,6 +1891,127 @@ struct sigma0builder { return true; } + //_______________________________________________ + // Build kstar candidate + template + bool buildKStar(TV0Object const& kshort, TV0Object const& gamma, TCollision const& collision, TMCParticles const& mcparticles) + { + + //_______________________________________________ + // Checking if both V0s are made of the very same tracks + if (gamma.posTrackExtraId() == kshort.posTrackExtraId() || + gamma.negTrackExtraId() == kshort.negTrackExtraId()) { + return false; + } + + //_______________________________________________ + // Kstar pre-selections + std::array pVecPhotons{gamma.px(), gamma.py(), gamma.pz()}; + std::array pVecKShort{kshort.px(), kshort.py(), kshort.pz()}; + + auto arrMom = std::array{pVecPhotons, pVecKShort}; + float kstarMass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassK0Short}); + float kstarY = -999.f; + + if constexpr (requires { gamma.pxMC(); kshort.pxMC(); }) // If MC + kstarY = RecoDecay::y(std::array{gamma.pxMC() + kshort.pxMC(), gamma.pyMC() + kshort.pyMC(), gamma.pzMC() + kshort.pzMC()}, o2::constants::physics::MassK0Star892); + else // If DATA + kstarY = RecoDecay::y(std::array{gamma.px() + kshort.px(), gamma.py() + kshort.py(), gamma.pz() + kshort.pz()}, o2::constants::physics::MassK0Star892); + + histos.fill(HIST("KStarSel/hSelectionStatistics"), 1.); + if (TMath::Abs(kstarMass - o2::constants::physics::MassK0Star892) > KStarWindow) + return false; + + histos.fill(HIST("KStarSel/hSelectionStatistics"), 2.); + if (TMath::Abs(kstarY) > KStarMaxRap) + return false; + + histos.fill(HIST("KStarSel/hKStarMassSelected"), kstarMass); + histos.fill(HIST("KStarSel/hSelectionStatistics"), 3.); + + auto kstarTopoInfo = propagateV0PairToDCA(gamma, kshort); + + kstarcores(kstarTopoInfo.X, kstarTopoInfo.Y, kstarTopoInfo.Z, kstarTopoInfo.DCADau, + gamma.px(), gamma.py(), gamma.pz(), gamma.mGamma(), kshort.px(), kshort.py(), kshort.pz(), kshort.mK0Short()); + + // MC properties + if constexpr (requires { gamma.motherMCPartId(); kshort.motherMCPartId(); }) { + auto kstarMCInfo = getV0PairMCInfo(gamma, kshort, collision, mcparticles); + + kstarmccores(kstarMCInfo.V0PairMCRadius, kstarMCInfo.V0PairPDGCode, kstarMCInfo.V0PairPDGCodeMother, kstarMCInfo.V0PairMCProcess, kstarMCInfo.fV0PairProducedByGenerator, + kstarMCInfo.V01MCpx, kstarMCInfo.V01MCpy, kstarMCInfo.V01MCpz, + kstarMCInfo.fIsV01Primary, kstarMCInfo.V01PDGCode, kstarMCInfo.V01PDGCodeMother, kstarMCInfo.fIsV01CorrectlyAssign, + kstarMCInfo.V02MCpx, kstarMCInfo.V02MCpy, kstarMCInfo.V02MCpz, + kstarMCInfo.fIsV02Primary, kstarMCInfo.V02PDGCode, kstarMCInfo.V02PDGCodeMother, kstarMCInfo.fIsV02CorrectlyAssign); + } + + // KStar -> stracollisions link + kstarCollRefs(collision.globalIndex()); + + //_______________________________________________ + // Photon extra properties + + auto posTrackGamma = gamma.template posTrackExtra_as(); + auto negTrackGamma = gamma.template negTrackExtra_as(); + + uint8_t fPhotonPosTrackCode = ((uint8_t(posTrackGamma.hasTPC()) << hasTPC) | + (uint8_t(posTrackGamma.hasITSTracker()) << hasITSTracker) | + (uint8_t(posTrackGamma.hasITSAfterburner()) << hasITSAfterburner) | + (uint8_t(posTrackGamma.hasTRD()) << hasTRD) | + (uint8_t(posTrackGamma.hasTOF()) << hasTOF)); + + uint8_t fPhotonNegTrackCode = ((uint8_t(negTrackGamma.hasTPC()) << hasTPC) | + (uint8_t(negTrackGamma.hasITSTracker()) << hasITSTracker) | + (uint8_t(negTrackGamma.hasITSAfterburner()) << hasITSAfterburner) | + (uint8_t(negTrackGamma.hasTRD()) << hasTRD) | + (uint8_t(negTrackGamma.hasTOF()) << hasTOF)); + + kstarPhotonExtras(gamma.qtarm(), gamma.alpha(), gamma.v0cosPA(), gamma.dcaV0daughters(), gamma.dcanegtopv(), gamma.dcapostopv(), gamma.v0radius(), gamma.z(), + posTrackGamma.tpcNSigmaEl(), negTrackGamma.tpcNSigmaEl(), posTrackGamma.tpcCrossedRows(), negTrackGamma.tpcCrossedRows(), + gamma.positiveeta(), gamma.negativeeta(), gamma.psipair(), posTrackGamma.itsNCls(), negTrackGamma.itsNCls(), posTrackGamma.itsChi2PerNcl(), negTrackGamma.itsChi2PerNcl(), + fPhotonPosTrackCode, fPhotonNegTrackCode, gamma.v0Type()); + + //_______________________________________________ + // KShort extra properties + + auto posTrackKShort = kshort.template posTrackExtra_as(); + auto negTrackKShort = kshort.template negTrackExtra_as(); + + uint8_t fKShortPosTrackCode = ((uint8_t(posTrackKShort.hasTPC()) << hasTPC) | + (uint8_t(posTrackKShort.hasITSTracker()) << hasITSTracker) | + (uint8_t(posTrackKShort.hasITSAfterburner()) << hasITSAfterburner) | + (uint8_t(posTrackKShort.hasTRD()) << hasTRD) | + (uint8_t(posTrackKShort.hasTOF()) << hasTOF)); + + uint8_t fKShortNegTrackCode = ((uint8_t(negTrackKShort.hasTPC()) << hasTPC) | + (uint8_t(negTrackKShort.hasITSTracker()) << hasITSTracker) | + (uint8_t(negTrackKShort.hasITSAfterburner()) << hasITSAfterburner) | + (uint8_t(negTrackKShort.hasTRD()) << hasTRD) | + (uint8_t(negTrackKShort.hasTOF()) << hasTOF)); + + float fKShortLifeTime = kshort.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0Short; + float fKShortPiTOFNSigma = -999.f; + + // float postrackTofNSigmaPi = -999.f; + // float negtrackTofNSigmaPi = -999.f; + + if constexpr (requires { kshort.tofNSigmaK0Pi(); }) { // If TOF info avaiable + + // postrackTofNSigmaPi = posTrackKShort.tofNSigmaPi(); + // negtrackTofNSigmaPi = negTrackKShort.tofNSigmaPi(); + fKShortPiTOFNSigma = kshort.tofNSigmaLaPi(); + } + + kshortExtras(kshort.qtarm(), kshort.alpha(), fKShortLifeTime, kshort.v0radius(), kshort.v0cosPA(), kshort.dcaV0daughters(), kshort.dcanegtopv(), kshort.dcapostopv(), + posTrackKShort.tpcNSigmaPi(), negTrackKShort.tpcNSigmaPi(), fKShortPiTOFNSigma, + posTrackKShort.tpcCrossedRows(), negTrackKShort.tpcCrossedRows(), + kshort.positiveeta(), kshort.negativeeta(), + posTrackKShort.itsNCls(), negTrackKShort.itsNCls(), posTrackKShort.itsChi2PerNcl(), negTrackKShort.itsChi2PerNcl(), + fKShortPosTrackCode, fKShortNegTrackCode, kshort.v0Type()); + + return true; + } + // Process photon and lambda candidates to build sigma0 candidates template void dataProcess(TCollision const& collisions, TV0s const& fullV0s, TMCParticles const& mcparticles) @@ -1722,6 +2022,9 @@ struct sigma0builder { std::vector bestGammasArray; std::vector bestLambdasArray; + std::vector bestKStarGammasArray; + std::vector bestKShortsArray; + // Custom grouping std::vector> v0grouped(collisions.size()); @@ -1742,6 +2045,9 @@ struct sigma0builder { bestGammasArray.clear(); bestLambdasArray.clear(); + bestKStarGammasArray.clear(); + bestKShortsArray.clear(); + float centrality = doPPAnalysis ? coll.centFT0M() : coll.centFT0C(); histos.fill(HIST("hEventCentrality"), centrality); @@ -1764,6 +2070,9 @@ struct sigma0builder { fillHistos<2>(v0, coll); // QA histos bestLambdasArray.push_back(v0.globalIndex()); // Save indices of best lambda candidates } + + if (processKShortCandidate(v0, coll)) // selecting kshorts + bestKShortsArray.push_back(v0.globalIndex()); // Save indices of best kshort candidates } //_______________________________________________ @@ -1792,6 +2101,18 @@ struct sigma0builder { } } + //_______________________________________________ + // KStar loop + if (fillKStarTables) { + for (size_t j = 0; j < bestKShortsArray.size(); ++j) { + auto kshort = fullV0s.rawIteratorAt(bestKShortsArray[j]); + + // Building kstar candidate & filling tables + if (!buildKStar(kshort, gamma1, coll, mcparticles)) + continue; + } + } + //_______________________________________________ // pi0 loop if (fillPi0Tables) { @@ -1878,25 +2199,25 @@ struct sigma0builder { auto v0MC = v0.template v0MCCore_as>(); // True Photon - if (v0MC.pdgCode() == 22 && fPassPhotonSel) { + if (v0MC.pdgCode() == PDG_t::kGamma && fPassPhotonSel) { histos.fill(HIST("PhotonLambdaQA/h3dTruePhotonMass"), centrality, v0MC.ptMC(), v0.mGamma()); - if (TMath::Abs(v0MC.pdgCodeMother()) == 3212) { // If from sigma0 or ASigma0 + if (TMath::Abs(v0MC.pdgCodeMother()) == PDG_t::kSigma0) { // If from sigma0 or ASigma0 histos.fill(HIST("PhotonLambdaQA/h2dTrueSigma0PhotonMass"), v0MC.ptMC(), v0.mGamma()); } } // True Lambda - if (v0MC.pdgCode() == 3122 && fPassLambdaSel) { + if (v0MC.pdgCode() == PDG_t::kLambda0 && fPassLambdaSel) { histos.fill(HIST("PhotonLambdaQA/h3dTrueLambdaMass"), centrality, v0MC.ptMC(), v0.mLambda()); - if (v0MC.pdgCodeMother() == 3212) { // If from sigma0 + if (v0MC.pdgCodeMother() == PDG_t::kSigma0) { // If from sigma0 histos.fill(HIST("PhotonLambdaQA/h2dTrueSigma0LambdaMass"), v0MC.ptMC(), v0.mLambda()); } } // True ALambda - if (v0MC.pdgCode() == -3122 && fPassLambdaSel) { + if (v0MC.pdgCode() == PDG_t::kLambda0Bar && fPassLambdaSel) { histos.fill(HIST("PhotonLambdaQA/h3dTrueALambdaMass"), centrality, v0MC.ptMC(), v0.mAntiLambda()); - if (v0MC.pdgCodeMother() == -3212) { // If from asigma0 + if (v0MC.pdgCodeMother() == PDG_t::kSigma0Bar) { // If from asigma0 histos.fill(HIST("PhotonLambdaQA/h2dTrueASigma0ALambdaMass"), v0MC.ptMC(), v0.mAntiLambda()); } } diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index 660eb015a1e..e4bc31fc311 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -292,3 +292,7 @@ o2physics_add_dpl_workflow(phitutorial-step3 SOURCES phitutorial_step3.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(k892hadronphoton + SOURCES k892hadronphoton.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB + COMPONENT_NAME Analysis) diff --git a/PWGLF/Tasks/Resonances/k892hadronphoton.cxx b/PWGLF/Tasks/Resonances/k892hadronphoton.cxx new file mode 100644 index 00000000000..e83ac55abd8 --- /dev/null +++ b/PWGLF/Tasks/Resonances/k892hadronphoton.cxx @@ -0,0 +1,1272 @@ +// 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. +// +// This is a task that reads kstar tables (from sigma0builder) to perform analysis. +// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// k892hadronphoton analysis task +// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// +// Comments, questions, complaints, suggestions? +// Please write to: +// oussama.benchikhi@cern.ch +// gianni.shigeru.setoue.liveraro@cern.ch +// + +#include "PWGLF/DataModel/LFSigmaTables.h" +#include "PWGLF/DataModel/LFStrangenessMLTables.h" +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/CCDB/ctpRateFetcher.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "Framework/ASoA.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +using std::array; +using KStars = soa::Join; +using MCKStars = soa::Join; + +static const std::vector PhotonSels = {"NoSel", "V0Type", "DCADauToPV", + "DCADau", "DauTPCCR", "TPCNSigmaEl", "V0pT", + "Y", "V0Radius", "RZCut", "Armenteros", "CosPA", + "PsiPair", "Phi", "Mass"}; + +static const std::vector KShortSels = {"NoSel", "V0Radius", "DCADau", "Armenteros", + "CosPA", "Y", "TPCCR", "DauITSCls", "Lifetime", + "TPCTOFPID", "DCADauToPV", "Mass"}; + +static const std::vector DirList = {"BeforeSel", "AfterSel"}; + +struct k892hadronphoton { + Service ccdb; + ctpRateFetcher rateFetcher; + + //__________________________________________________ + // For manual sliceBy + // SliceCache cache; + PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; + + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // Event level + Configurable doPPAnalysis{"doPPAnalysis", true, "if in pp, set to true"}; + 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 { + 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", true, "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", false, "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 requireNoCollInTimeRangeVzDep{"requireNoCollInTimeRangeVzDep", false, "reject collisions corrupted by the cannibalism, with other collisions with pvZ of drifting TPC tracks from past/future collisions within 2.5 cm the current pvZ"}; + 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", false, "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)"}; + // fast check on occupancy + Configurable minOccupancy{"minOccupancy", -1, "minimum occupancy from neighbouring collisions"}; + Configurable maxOccupancy{"maxOccupancy", -1, "maximum occupancy from neighbouring collisions"}; + + // fast check on interaction rate + Configurable minIR{"minIR", -1, "minimum IR collisions"}; + Configurable maxIR{"maxIR", -1, "maximum IR collisions"}; + + } eventSelections; + + // generated + Configurable mc_keepOnlyFromGenerator{"mc_keepOnlyFromGenerator", true, "if true, consider only particles from generator to calculate efficiency."}; + + // QA + Configurable fillBkgQAhistos{"fillBkgQAhistos", false, "if true, fill MC QA histograms for Bkg study. Only works with MC."}; + Configurable fillResoQAhistos{"fillResoQAhistos", false, "if true, fill MC QA histograms for pT resolution study. Only works with MC."}; + + // Analysis strategy: + Configurable doMCAssociation{"doMCAssociation", false, "Flag to process only signal candidates. Use only with processMonteCarlo!"}; + Configurable selRecoFromGenerator{"selRecoFromGenerator", false, "Flag to process only signal candidates from generator"}; + Configurable doPhotonKShortSelQA{"doPhotonKShortSelQA", false, "Flag to fill photon and kshort QA histos!"}; + + //// K0Short criteria: + struct : ConfigurableGroup { + Configurable KShort_MLThreshold{"KShort_MLThreshold", 0.1, "Decision Threshold value to select kshorts"}; + Configurable KShortMinDCANegToPv{"KShortMinDCANegToPv", .05, "min DCA Neg To PV (cm)"}; + Configurable KShortMinDCAPosToPv{"KShortMinDCAPosToPv", .05, "min DCA Pos To PV (cm)"}; + Configurable KShortMaxDCAV0Dau{"KShortMaxDCAV0Dau", 2.5, "Max DCA V0 Daughters (cm)"}; + Configurable KShortMinv0radius{"KShortMinv0radius", 0.0, "Min V0 radius (cm)"}; + Configurable KShortMaxv0radius{"KShortMaxv0radius", 40, "Max V0 radius (cm)"}; + Configurable KShortMinQt{"KShortMinQt", 0.1, "Min kshort qt value (AP plot) (GeV/c)"}; + Configurable KShortMaxQt{"KShortMaxQt", 2.5, "Max kshort qt value (AP plot) (GeV/c)"}; + Configurable KShortMinAlpha{"KShortMinAlpha", 0.0, "Min kshort alpha absolute value (AP plot)"}; + Configurable KShortMaxAlpha{"KShortMaxAlpha", 1.0, "Max kshort alpha absolute value (AP plot)"}; + Configurable KShortMinv0cospa{"KShortMinv0cospa", 0.95, "Min V0 CosPA"}; + Configurable KShortMaxLifeTime{"KShortMaxLifeTime", 30, "Max lifetime"}; + Configurable KShortWindow{"KShortWindow", 0.015, "Mass window around expected (in GeV/c2)"}; + Configurable KShortMaxRap{"KShortMaxRap", 0.8, "Max kshort rapidity"}; + Configurable KShortMaxDauEta{"KShortMaxDauEta", 0.8, "Max pseudorapidity of daughter tracks"}; + Configurable fselKShortTPCPID{"fselKShortTPCPID", true, "Flag to select kshort-like candidates using TPC NSigma."}; + Configurable fselKShortTOFPID{"fselKShortTOFPID", false, "Flag to select kshort-like candidates using TOF NSigma."}; + Configurable KShortMaxTPCNSigmas{"KShortMaxTPCNSigmas", 1e+9, "Max TPC NSigmas for daughters"}; + // Configurable KShortPrMaxTOFNSigmas{"KShortPrMaxTOFNSigmas", 1e+9, "Max TOF NSigmas for daughters"}; + Configurable KShortPiMaxTOFNSigmas{"KShortPiMaxTOFNSigmas", 1e+9, "Max TOF NSigmas for daughters"}; + Configurable KShortMinTPCCrossedRows{"KShortMinTPCCrossedRows", 50, "Min daughter TPC Crossed Rows"}; + Configurable KShortMinITSclusters{"KShortMinITSclusters", 1, "minimum ITS clusters"}; + Configurable KShortRejectPosITSafterburner{"KShortRejectPosITSafterburner", false, "reject positive track formed out of afterburner ITS tracks"}; + Configurable KShortRejectNegITSafterburner{"KShortRejectNegITSafterburner", false, "reject negative track formed out of afterburner ITS tracks"}; + } kshortSelections; + + //// Photon criteria: + struct : ConfigurableGroup { + Configurable Gamma_MLThreshold{"Gamma_MLThreshold", 0.1, "Decision Threshold value to select gammas"}; + Configurable Photonv0TypeSel{"Photonv0TypeSel", 7, "select on a certain V0 type (leave negative if no selection desired)"}; + Configurable PhotonMinDCADauToPv{"PhotonMinDCADauToPv", 0.0, "Min DCA daughter To PV (cm)"}; + Configurable PhotonMaxDCAV0Dau{"PhotonMaxDCAV0Dau", 3.5, "Max DCA V0 Daughters (cm)"}; + Configurable PhotonMinTPCCrossedRows{"PhotonMinTPCCrossedRows", 30, "Min daughter TPC Crossed Rows"}; + Configurable PhotonMinTPCNSigmas{"PhotonMinTPCNSigmas", -7, "Min TPC NSigmas for daughters"}; + Configurable PhotonMaxTPCNSigmas{"PhotonMaxTPCNSigmas", 7, "Max TPC NSigmas for daughters"}; + Configurable PhotonMinPt{"PhotonMinPt", 0.0, "Min photon pT (GeV/c)"}; + Configurable PhotonMaxPt{"PhotonMaxPt", 50.0, "Max photon pT (GeV/c)"}; + Configurable PhotonMaxRap{"PhotonMaxRap", 0.5, "Max photon rapidity"}; + Configurable PhotonMinRadius{"PhotonMinRadius", 3.0, "Min photon conversion radius (cm)"}; + Configurable PhotonMaxRadius{"PhotonMaxRadius", 115, "Max photon conversion radius (cm)"}; + Configurable PhotonMaxZ{"PhotonMaxZ", 240, "Max photon conversion point z value (cm)"}; + Configurable PhotonMaxQt{"PhotonMaxQt", 0.05, "Max photon qt value (AP plot) (GeV/c)"}; + Configurable PhotonMaxAlpha{"PhotonMaxAlpha", 0.95, "Max photon alpha absolute value (AP plot)"}; + Configurable PhotonMinV0cospa{"PhotonMinV0cospa", 0.80, "Min V0 CosPA"}; + Configurable PhotonMaxMass{"PhotonMaxMass", 0.10, "Max photon mass (GeV/c^{2})"}; + Configurable PhotonPsiPairMax{"PhotonPsiPairMax", 1e+9, "maximum psi angle of the track pair"}; + Configurable PhotonMaxDauEta{"PhotonMaxDauEta", 0.8, "Max pseudorapidity of daughter tracks"}; + Configurable PhotonLineCutZ0{"PhotonLineCutZ0", 7.0, "The offset for the linecute used in the Z vs R plot"}; + Configurable PhotonPhiMin1{"PhotonPhiMin1", -1, "Phi min value to reject photons, region 1 (leave negative if no selection desired)"}; + Configurable PhotonPhiMax1{"PhotonPhiMax1", -1, "Phi max value to reject photons, region 1 (leave negative if no selection desired)"}; + Configurable PhotonPhiMin2{"PhotonPhiMin2", -1, "Phi max value to reject photons, region 2 (leave negative if no selection desired)"}; + Configurable PhotonPhiMax2{"PhotonPhiMax2", -1, "Phi min value to reject photons, region 2 (leave negative if no selection desired)"}; + } photonSelections; + + struct : ConfigurableGroup { + Configurable KStarMaxRap{"KStarMaxRap", 0.5, "Max kstar rapidity"}; + Configurable KStarMaxRadius{"KStarMaxRadius", 200, "Max kstar decay radius"}; + Configurable KStarMaxDCADau{"KStarMaxDCADau", 50, "Max kstar DCA between daughters"}; + Configurable KStarMaxOPAngle{"KStarMaxOPAngle", 7, "Max kstar OP Angle between daughters"}; + } kstarSelections; + + // Axis + // base properties + ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "Centrality"}; + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "p_{T} (GeV/c)"}; + ConfigurableAxis axisInvPt{"axisInvPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 2.0, 5.0, 10.0, 20.0, 50.0}, ""}; + ConfigurableAxis axisDeltaPt{"axisDeltaPt", {400, -50.0, 50.0}, ""}; + ConfigurableAxis axisRapidity{"axisRapidity", {100, -2.0f, 2.0f}, "Rapidity"}; + ConfigurableAxis axisIRBinning{"axisIRBinning", {150, 0, 1500}, "Binning for the interaction rate (kHz)"}; + ConfigurableAxis axisNch{"axisNch", {300, 0.0f, 3000.0f}, "N_{ch}"}; + ConfigurableAxis axisGeneratorIds{"axisGeneratorIds", {256, -0.5f, 255.5f}, "axis for generatorIds"}; + + // invariant mass + ConfigurableAxis axisPhotonMass{"axisPhotonMass", {200, -0.1f, 0.5f}, "M_{#Gamma}"}; + ConfigurableAxis axisKShortMass{"axisKShortMass", {200, 0.3f, 0.6f}, "M_{#K_{s}^{0}}"}; + ConfigurableAxis axisKStarMass{"axisKStarMass", {200, 0.08f, 1.5f}, "M_{#K^{*}}"}; + + // AP plot axes + ConfigurableAxis axisAPAlpha{"axisAPAlpha", {220, -1.1f, 1.1f}, "V0 AP alpha"}; + ConfigurableAxis axisAPQt{"axisAPQt", {220, 0.0f, 0.5f}, "V0 AP alpha"}; + + // Track quality, PID and other axes + ConfigurableAxis axisTPCrows{"axisTPCrows", {160, 0.0f, 160.0f}, "N TPC rows"}; + ConfigurableAxis axisNCls{"axisNCls", {8, -0.5, 7.5}, "NCls"}; + ConfigurableAxis axisChi2PerNcl{"axisChi2PerNcl", {80, -40, 40}, "Chi2 Per Ncl"}; + ConfigurableAxis axisTPCNSigma{"axisTPCNSigma", {120, -30, 30}, "TPC NSigma"}; + ConfigurableAxis axisTOFNSigma{"axisTOFNSigma", {120, -30, 30}, "TOF NSigma"}; + ConfigurableAxis axisLifetime{"axisLifetime", {100, 0, 100}, "Chi2 Per Ncl"}; + + // topological variable QA axes + ConfigurableAxis axisV0Radius{"axisV0Radius", {240, 0.0f, 120.0f}, "V0 radius (cm)"}; + ConfigurableAxis axisV0PairRadius{"axisV0PairRadius", {200, 0.0f, 20.0f}, "V0Pair radius (cm)"}; + ConfigurableAxis axisDCAtoPV{"axisDCAtoPV", {500, 0.0f, 50.0f}, "DCA (cm)"}; + ConfigurableAxis axisDCAdau{"axisDCAdau", {50, 0.0f, 5.0f}, "DCA (cm)"}; + ConfigurableAxis axisCosPA{"axisCosPA", {200, 0.5f, 1.0f}, "Cosine of pointing angle"}; + ConfigurableAxis axisPA{"axisPA", {100, 0.0f, 1}, "Pointing angle"}; + ConfigurableAxis axisPsiPair{"axisPsiPair", {250, -5.0f, 5.0f}, "Psipair for photons"}; + ConfigurableAxis axisPhi{"axisPhi", {200, 0, 2 * o2::constants::math::PI}, "Phi for photons"}; + ConfigurableAxis axisZ{"axisZ", {120, -120.0f, 120.0f}, "V0 Z position (cm)"}; + + ConfigurableAxis axisCandSel{"axisCandSel", {20, 0.5f, +20.5f}, "Candidate Selection"}; + + // ML + ConfigurableAxis MLProb{"MLOutput", {100, 0.0f, 1.0f}, ""}; + + void init(InitContext const&) + { + LOGF(info, "Initializing now: cross-checking correctness..."); + if (doprocessRealData + doprocessMonteCarlo > 1) { + LOGF(fatal, "You have enabled more than one process function. Please check your configuration! Aborting now."); + } + + // setting CCDB service + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setFatalWhenNull(false); + + // Event Counters + histos.add("hEventCentrality", "hEventCentrality", kTH1D, {axisCentrality}); + + histos.add("hEventSelection", "hEventSelection", kTH1D, {{21, -0.5f, +20.5f}}); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(2, "sel8 cut"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(3, "kIsTriggerTVX"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(4, "kNoITSROFrameBorder"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(5, "kNoTimeFrameBorder"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(6, "posZ cut"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(7, "kIsVertexITSTPC"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(8, "kIsGoodZvtxFT0vsPV"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(9, "kIsVertexTOFmatched"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(10, "kIsVertexTRDmatched"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(11, "kNoSameBunchPileup"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(12, "kNoCollInTimeRangeStd"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(13, "kNoCollInTimeRangeStrict"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(14, "kNoCollInTimeRangeNarrow"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(15, "kNoCollInRofStd"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(16, "kNoCollInRofStrict"); + if (doPPAnalysis) { + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(17, "INEL>0"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(18, "INEL>1"); + } else { + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(17, "Below min occup."); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(18, "Above max occup."); + } + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(19, "Below min IR"); + histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(20, "Above max IR"); + + if (fGetIR) { + histos.add("GeneralQA/hRunNumberNegativeIR", "", kTH1D, {{1, 0., 1.}}); + histos.add("GeneralQA/hInteractionRate", "hInteractionRate", kTH1D, {axisIRBinning}); + histos.add("GeneralQA/hCentralityVsInteractionRate", "hCentralityVsInteractionRate", kTH2D, {axisCentrality, axisIRBinning}); + } + + if (doprocessRealData || doprocessMonteCarlo) { + for (const auto& histodir : DirList) { + + histos.add(histodir + "/Photon/hTrackCode", "hTrackCode", kTH1D, {{11, 0.5f, 11.5f}}); + histos.add(histodir + "/Photon/hV0Type", "hV0Type", kTH1D, {{8, 0.5f, 8.5f}}); + histos.add(histodir + "/Photon/hDCANegToPV", "hDCANegToPV", kTH1D, {axisDCAtoPV}); + histos.add(histodir + "/Photon/hDCAPosToPV", "hDCAPosToPV", kTH1D, {axisDCAtoPV}); + histos.add(histodir + "/Photon/hDCADau", "hDCADau", kTH1D, {axisDCAdau}); + histos.add(histodir + "/Photon/hPosTPCCR", "hPosTPCCR", kTH1D, {axisTPCrows}); + histos.add(histodir + "/Photon/hNegTPCCR", "hNegTPCCR", kTH1D, {axisTPCrows}); + histos.add(histodir + "/Photon/hPosTPCNSigmaEl", "hPosTPCNSigmaEl", kTH1D, {axisTPCNSigma}); + histos.add(histodir + "/Photon/hNegTPCNSigmaEl", "hNegTPCNSigmaEl", kTH1D, {axisTPCNSigma}); + + histos.add(histodir + "/Photon/hpT", "hpT", kTH1D, {axisPt}); + histos.add(histodir + "/Photon/hY", "hY", kTH1D, {axisRapidity}); + histos.add(histodir + "/Photon/hPosEta", "hPosEta", kTH1D, {axisRapidity}); + histos.add(histodir + "/Photon/hNegEta", "hNegEta", kTH1D, {axisRapidity}); + histos.add(histodir + "/Photon/hRadius", "hRadius", kTH1D, {axisV0Radius}); + histos.add(histodir + "/Photon/hZ", "hZ", kTH1D, {axisZ}); + histos.add(histodir + "/Photon/h2dRZCut", "h2dRZCut", kTH2D, {axisZ, axisV0Radius}); + histos.add(histodir + "/Photon/h2dRZPlane", "h2dRZPlane", kTH2D, {axisZ, axisV0Radius}); + histos.add(histodir + "/Photon/hCosPA", "hCosPA", kTH1D, {axisCosPA}); + histos.add(histodir + "/Photon/hPsiPair", "hPsiPair", kTH1D, {axisPsiPair}); + histos.add(histodir + "/Photon/hPhi", "hPhi", kTH1D, {axisPhi}); + histos.add(histodir + "/Photon/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisPhotonMass}); + histos.add(histodir + "/Photon/hMass", "hMass", kTH1D, {axisPhotonMass}); + + histos.add(histodir + "/KShort/hTrackCode", "hTrackCode", kTH1D, {{11, 0.5f, 11.5f}}); + histos.add(histodir + "/KShort/hRadius", "hRadius", kTH1D, {axisV0Radius}); + histos.add(histodir + "/KShort/hDCADau", "hDCADau", kTH1D, {axisDCAdau}); + histos.add(histodir + "/KShort/hCosPA", "hCosPA", kTH1D, {axisCosPA}); + histos.add(histodir + "/KShort/hY", "hY", kTH1D, {axisRapidity}); + histos.add(histodir + "/KShort/hPosEta", "hPosEta", kTH1D, {axisRapidity}); + histos.add(histodir + "/KShort/hNegEta", "hNegEta", kTH1D, {axisRapidity}); + histos.add(histodir + "/KShort/hPosTPCCR", "hPosTPCCR", kTH1D, {axisTPCrows}); + histos.add(histodir + "/KShort/hNegTPCCR", "hNegTPCCR", kTH1D, {axisTPCrows}); + histos.add(histodir + "/KShort/hPosITSCls", "hPosITSCls", kTH1D, {axisNCls}); + histos.add(histodir + "/KShort/hNegITSCls", "hNegITSCls", kTH1D, {axisNCls}); + histos.add(histodir + "/KShort/hPosChi2PerNc", "hPosChi2PerNc", kTH1D, {axisChi2PerNcl}); + histos.add(histodir + "/KShort/hNegChi2PerNc", "hNegChi2PerNc", kTH1D, {axisChi2PerNcl}); + histos.add(histodir + "/KShort/hLifeTime", "hLifeTime", kTH1D, {axisLifetime}); + histos.add(histodir + "/KShort/h2dTPCvsTOFNSigma_KShortPi", "h2dTPCvsTOFNSigma_KShortPi", kTH2D, {axisTPCNSigma, axisTOFNSigma}); + histos.add(histodir + "/KShort/hKShortDCANegToPV", "hKShortDCANegToPV", kTH1D, {axisDCAtoPV}); + histos.add(histodir + "/KShort/hKShortDCAPosToPV", "hKShortDCAPosToPV", kTH1D, {axisDCAtoPV}); + histos.add(histodir + "/KShort/hKShortpT", "hKShortpT", kTH1D, {axisPt}); + histos.add(histodir + "/KShort/hKShortMass", "hKShortMass", kTH1D, {axisKShortMass}); + histos.add(histodir + "/KShort/h3dKShortMass", "h3dKShortMass", kTH3D, {axisCentrality, axisPt, axisKShortMass}); + + histos.add(histodir + "/h2dArmenteros", "h2dArmenteros", kTH2D, {axisAPAlpha, axisAPQt}); + + histos.add(histodir + "/KStar/hMass", "hMass", kTH1D, {axisKStarMass}); + histos.add(histodir + "/KStar/hPt", "hPt", kTH1D, {axisPt}); + histos.add(histodir + "/KStar/hY", "hY", kTH1D, {axisRapidity}); + histos.add(histodir + "/KStar/hRadius", "hRadius", kTH1D, {axisV0PairRadius}); + histos.add(histodir + "/KStar/h2dRadiusVspT", "h2dRadiusVspT", kTH2D, {axisV0PairRadius, axisPt}); + histos.add(histodir + "/KStar/hDCAPairDau", "hDCAPairDau", kTH1D, {axisDCAdau}); + histos.add(histodir + "/KStar/hDCAPairDauVsPt", "hDCAPairDauVsPt", kTH2D, {axisDCAdau, axisPt}); + + histos.add(histodir + "/KStar/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisKStarMass}); + histos.add(histodir + "/KStar/h3dOPAngleVsMass", "h3dOPAngleVsMass", kTH3D, {{140, 0.0f, +7.0f}, axisPt, axisKStarMass}); + histos.add(histodir + "/KStar/h2dOPAngleVsPt", "h2dOPAngleVsPt", kTH2D, {{140, 0.0f, +7.0f}, axisPt}); + + if (doprocessMonteCarlo) { + + histos.add(histodir + "/MC/Photon/hV0ToCollAssoc", "hV0ToCollAssoc", kTH1D, {{2, 0.0f, 2.0f}}); + histos.add(histodir + "/MC/Photon/hPt", "hPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/Photon/hMCPt", "hMCPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/Photon/hPosTPCNSigmaEl", "hPosTPCNSigmaEl", kTH1D, {axisTPCNSigma}); + histos.add(histodir + "/MC/Photon/hNegTPCNSigmaEl", "hNegTPCNSigmaEl", kTH1D, {axisTPCNSigma}); + histos.add(histodir + "/MC/Photon/h2dPAVsPt", "h2dPAVsPt", kTH2D, {axisPA, axisPt}); + histos.add(histodir + "/MC/Photon/hPt_BadCollAssig", "hPt_BadCollAssig", kTH1D, {axisPt}); + histos.add(histodir + "/MC/Photon/h2dPAVsPt_BadCollAssig", "h2dPAVsPt_BadCollAssig", kTH2D, {axisPA, axisPt}); + + histos.add(histodir + "/MC/KShort/hV0ToCollAssoc", "hV0ToCollAssoc", kTH1D, {{2, 0.0f, 2.0f}}); + histos.add(histodir + "/MC/KShort/hPt", "hPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/KShort/hMass", "hMass", kTH1D, {axisKShortMass}); + histos.add(histodir + "/MC/KShort/hMCPt", "hMCPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/KShort/h3dTPCvsTOFNSigma_Pi", "h3dTPCvsTOFNSigma_Pi", kTH3D, {axisTPCNSigma, axisTOFNSigma, axisPt}); + + histos.add(histodir + "/MC/h2dArmenteros", "h2dArmenteros", kTH2D, {axisAPAlpha, axisAPQt}); + + histos.add(histodir + "/MC/KStar/hPt", "hPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/KStar/hMCPt", "hMCPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/KStar/hMass", "hMass", kTH1D, {axisKStarMass}); + histos.add(histodir + "/MC/KStar/hMCProcess", "hMCProcess", kTH1D, {{50, -0.5f, 49.5f}}); + histos.add(histodir + "/MC/KStar/hGenRadius", "hGenRadius", kTH1D, {axisV0PairRadius}); + histos.add(histodir + "/MC/KStar/h2dMCPtVsKShortMCPt", "h2dMCPtVsKShortMCPt", kTH2D, {axisPt, axisPt}); + histos.add(histodir + "/MC/KStar/h2dMCPtVsPhotonMCPt", "h2dMCPtVsPhotonMCPt", kTH2D, {axisPt, axisPt}); + histos.add(histodir + "/MC/KStar/h2dMCProcessVsGenRadius", "h2dMCProcessVsGenRadius", kTH2D, {{50, -0.5f, 49.5f}, axisV0PairRadius}); + histos.add(histodir + "/MC/KStar/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisKStarMass}); + histos.add(histodir + "/MC/KStar/h3dMCProcess", "h3dMCProcess", kTH3D, {{50, -0.5f, 49.5f}, axisPt, axisKStarMass}); + histos.add(histodir + "/MC/KStar/h2dOPAngleVsPt", "h2dOPAngleVsPt", kTH2D, {{140, 0.0f, +7.0f}, axisPt}); + histos.add(histodir + "/MC/KStar/h2dRadiusVspT", "h2dRadiusVspT", kTH2D, {axisV0PairRadius, axisPt}); + histos.add(histodir + "/MC/KStar/hDCAPairDauVsPt", "hDCAPairDauVsPt", kTH2D, {axisDCAdau, axisPt}); + + // 1/pT Resolution: + if (fillResoQAhistos && histodir == "BeforeSel") { + + histos.add(histodir + "/MC/Reso/h2dKShortPtResolution", "h2dKShortPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); + histos.add(histodir + "/MC/Reso/h3dKShortPtResoVsTPCCR", "h3dKShortPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); + histos.add(histodir + "/MC/Reso/h3dKShortPtResoVsTPCCR", "h3dKShortPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); + histos.add(histodir + "/MC/Reso/h2dKStarPtResolution", "h2dKStarPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); + histos.add(histodir + "/MC/Reso/h2dKStarRadiusResolution", "h2dKStarRadiusResolution", kTH2D, {axisPt, axisDeltaPt}); + } + + // For background decomposition study + if (fillBkgQAhistos) { + histos.add(histodir + "/MC/BkgStudy/h2dPtVsMassKStar_All", "h2dPtVsMassKStar_All", kTH2D, {axisPt, axisKStarMass}); + histos.add(histodir + "/MC/BkgStudy/h2dPtVsMassKStar_TrueDaughters", "h2dPtVsMassKStar_TrueDaughters", kTH2D, {axisPt, axisKStarMass}); + histos.add(histodir + "/MC/BkgStudy/h2dTrueDaughtersMatrix", "h2dTrueDaughtersMatrix", kTHnSparseD, {{10001, -5000.5f, +5000.5f}, {10001, -5000.5f, +5000.5f}}); + histos.add(histodir + "/MC/BkgStudy/h2dPtVsMassKStar_TrueGammaFakeKShort", "h2dPtVsMassKStar_TrueGammaFakeKShort", kTH2D, {axisPt, axisKStarMass}); + histos.add(histodir + "/MC/BkgStudy/h2dPtVsMassKStar_FakeGammaTrueKShort", "h2dPtVsMassKStar_FakeGammaTrueKShort", kTH2D, {axisPt, axisKStarMass}); + histos.add(histodir + "/MC/BkgStudy/h2dPtVsMassKStar_FakeDaughters", "h2dPtVsMassKStar_FakeDaughters", kTH2D, {axisPt, axisKStarMass}); + } + } + } + + // Selections + histos.add("Selection/Photon/hCandidateSel", "hCandidateSel", kTH1D, {axisCandSel}); + histos.add("Selection/KShort/hCandidateSel", "hCandidateSel", kTH1D, {axisCandSel}); + + for (size_t i = 0; i < PhotonSels.size(); ++i) { + const auto& sel = PhotonSels[i]; + + histos.add(Form("Selection/Photon/h2d%s", sel.c_str()), ("h2d" + sel).c_str(), kTH2D, {axisPt, axisPhotonMass}); + histos.get(HIST("Selection/Photon/hCandidateSel"))->GetXaxis()->SetBinLabel(i + 1, sel.c_str()); + histos.add(Form("Selection/KStar/h2dPhoton%s", sel.c_str()), ("h2dPhoton" + sel).c_str(), kTH2D, {axisPt, axisKStarMass}); + } + + for (size_t i = 0; i < KShortSels.size(); ++i) { + const auto& sel = KShortSels[i]; + + histos.add(Form("Selection/KShort/h2d%s", sel.c_str()), ("h2d" + sel).c_str(), kTH2D, {axisPt, axisKShortMass}); + histos.get(HIST("Selection/KShort/hCandidateSel"))->GetXaxis()->SetBinLabel(i + 1, sel.c_str()); + histos.add(Form("Selection/KStar/h2dKShort%s", sel.c_str()), ("h2dKShort" + sel).c_str(), kTH2D, {axisPt, axisKStarMass}); + } + } + + if (doprocessGeneratedRun3) { + + histos.add("Gen/hGenEvents", "hGenEvents", kTH2D, {{axisNch}, {2, -0.5f, +1.5f}}); + histos.get(HIST("Gen/hGenEvents"))->GetYaxis()->SetBinLabel(1, "All gen. events"); + histos.get(HIST("Gen/hGenEvents"))->GetYaxis()->SetBinLabel(2, "Gen. with at least 1 rec. events"); + + 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/hEventPVzMC", "hEventPVzMC", kTH1D, {{100, -20.0f, +20.0f}}); + histos.add("Gen/hCentralityVsPVzMC", "hCentralityVsPVzMC", kTH2D, {{101, 0.0f, 101.0f}, {100, -20.0f, +20.0f}}); + + // KStar specific + histos.add("Gen/h2dGenKStar", "h2dGenKStar", kTH2D, {axisCentrality, axisPt}); + histos.add("Gen/h2dGenKStarVsMultMC_RecoedEvt", "h2dGenKStarVsMultMC_RecoedEvt", kTH2D, {axisNch, axisPt}); + histos.add("Gen/h2dGenKStarVsMultMC", "h2dGenKStarVsMultMC", kTH2D, {axisNch, axisPt}); + } + } + + // ______________________________________________________ + // Check whether the collision passes our collision selections + // Should work with collisions, mccollisions, stracollisions and stramccollisions tables! + 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 (doPPAnalysis) { // we are in pp + 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 */); + } else { // we are in Pb-Pb + float collisionOccupancy = eventSelections.useFT0CbasedOccupancy ? collision.ft0cOccupancyInTimeRange() : collision.trackOccupancyInTimeRange(); + if (eventSelections.minOccupancy >= 0 && collisionOccupancy < eventSelections.minOccupancy) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 16 /* Below min occupancy */); + if (eventSelections.maxOccupancy >= 0 && collisionOccupancy > eventSelections.maxOccupancy) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 17 /* 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; + float centrality = doPPAnalysis ? collision.centFT0M() : collision.centFT0C(); + + if (fGetIR) { + if (interactionRate < 0) + histos.get(HIST("GeneralQA/hRunNumberNegativeIR"))->Fill(Form("%d", collision.runNumber()), 1); // This lists all run numbers without IR info! + + histos.fill(HIST("GeneralQA/hInteractionRate"), interactionRate); + histos.fill(HIST("GeneralQA/hCentralityVsInteractionRate"), centrality, interactionRate); + } + + if (eventSelections.minIR >= 0 && interactionRate < eventSelections.minIR) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 18 /* Below min IR */); + + if (eventSelections.maxIR >= 0 && interactionRate > eventSelections.maxIR) { + return false; + } + if (fillHists) + histos.fill(HIST("hEventSelection"), 19 /* Above max IR */); + + // Fill centrality histogram after event selection + if (fillHists) + histos.fill(HIST("hEventCentrality"), centrality); + + return true; + } + + // ______________________________________________________ + // Simulated processing + // Return the list of indices to the recoed collision associated to a given MC collision. + template + std::vector getListOfRecoCollIndices(TMCollisions const& mcCollisions, TCollisions const& collisions) + { + std::vector listBestCollisionIdx(mcCollisions.size()); + for (auto const& mcCollision : mcCollisions) { + auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + int biggestNContribs = -1; + int bestCollisionIndex = -1; + for (auto const& collision : groupedCollisions) { + // consider event selections in the recoed <-> gen collision association, for the denominator (or numerator) of the efficiency (or signal loss)? + if (eventSelections.useEvtSelInDenomEff) { + if (!IsEventAccepted(collision, false)) { + continue; + } + } + // Find the collision with the biggest nbr of PV contributors + // Follows what was done here: https://github.com/AliceO2Group/O2Physics/blob/master/Common/TableProducer/mcCollsExtra.cxx#L93 + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + bestCollisionIndex = collision.globalIndex(); + } + } + listBestCollisionIdx[mcCollision.globalIndex()] = bestCollisionIndex; + } + return listBestCollisionIdx; + } + + // ______________________________________________________ + // Simulated processing + // Fill generated event information (for event loss/splitting estimation) + template + void fillGeneratedEventProperties(TMCCollisions const& mcCollisions, TCollisions const& collisions) + { + std::vector listBestCollisionIdx(mcCollisions.size()); + for (auto const& mcCollision : mcCollisions) { + // Apply selections on MC collisions + if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) { + continue; + } + if (doPPAnalysis) { // we are in pp + if (eventSelections.requireINEL0 && mcCollision.multMCNParticlesEta10() < 1) { + continue; + } + + if (eventSelections.requireINEL1 && mcCollision.multMCNParticlesEta10() < 2) { + continue; + } + } + + histos.fill(HIST("Gen/hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/); + + auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + // Check if there is at least one of the reconstructed collisions associated to this MC collision + // If so, we consider it + bool atLeastOne = false; + int biggestNContribs = -1; + float centrality = 100.5f; + int nCollisions = 0; + for (auto const& collision : groupedCollisions) { + + if (!IsEventAccepted(collision, false)) { + continue; + } + + if (biggestNContribs < collision.multPVTotalContributors()) { + biggestNContribs = collision.multPVTotalContributors(); + centrality = doPPAnalysis ? collision.centFT0M() : collision.centFT0C(); + } + + nCollisions++; + atLeastOne = true; + } + + histos.fill(HIST("Gen/hCentralityVsNcoll_beforeEvSel"), centrality, groupedCollisions.size()); + histos.fill(HIST("Gen/hCentralityVsNcoll_afterEvSel"), centrality, nCollisions); + histos.fill(HIST("Gen/hCentralityVsMultMC"), centrality, mcCollision.multMCNParticlesEta05()); + histos.fill(HIST("Gen/hCentralityVsPVzMC"), centrality, mcCollision.posZ()); + histos.fill(HIST("Gen/hEventPVzMC"), mcCollision.posZ()); + + if (atLeastOne) { + histos.fill(HIST("Gen/hGenEvents"), mcCollision.multMCNParticlesEta05(), 1 /* at least 1 rec. event*/); + histos.fill(HIST("Gen/hGenEventCentrality"), centrality); + } + } + return; + } + + // ______________________________________________________ + // Simulated processing (subscribes to MC information too) + template + void analyzeGenerated(TMCCollisions const& mcCollisions, TCollisions const& collisions, TGenParticles const& genParticles) + { + fillGeneratedEventProperties(mcCollisions, collisions); + std::vector listBestCollisionIdx = getListOfRecoCollIndices(mcCollisions, collisions); + + for (auto& genParticle : genParticles) { + float centrality = 100.5f; + + // Has MC collision + if (!genParticle.has_straMCCollision()) + continue; + + // Selection on the source (generator/transport) + if (!genParticle.producedByGenerator() && mc_keepOnlyFromGenerator) + continue; + + // Select corresponding mc collision && Basic event selection + auto mcCollision = genParticle.template straMCCollision_as>(); + if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) { + continue; + } + if (doPPAnalysis) { // we are in pp + if (eventSelections.requireINEL0 && mcCollision.multMCNParticlesEta10() < 1) { + continue; + } + + if (eventSelections.requireINEL1 && mcCollision.multMCNParticlesEta10() < 2) { + continue; + } + } + + //______________________________________________________________________________ + // Generated KStar processing + if constexpr (requires { genParticle.kstarMCPt(); }) { + + float ptmc = genParticle.kstarMCPt(); + + if (listBestCollisionIdx[mcCollision.globalIndex()] > -1) { + auto collision = collisions.iteratorAt(listBestCollisionIdx[mcCollision.globalIndex()]); + centrality = doPPAnalysis ? collision.centFT0M() : collision.centFT0C(); + + if (genParticle.isKStar()) + histos.fill(HIST("Gen/h2dGenKStarVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc); + } + + if (genParticle.isKStar()) { + histos.fill(HIST("Gen/h2dGenKStar"), centrality, ptmc); + histos.fill(HIST("Gen/h2dGenKStarVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc); + } + } + } + } + + //__________________________________________ + template + int retrieveV0TrackCode(TKStarObject const& kstar) + { + + int TrkCode = 10; // 1: TPC-only, 2: TPC+Something, 3: ITS-Only, 4: ITS+TPC + Something, 10: anything else + + if (isGamma) { + if (kstar.photonPosTrackCode() == 1 && kstar.photonNegTrackCode() == 1) + TrkCode = 1; + if ((kstar.photonPosTrackCode() != 1 && kstar.photonNegTrackCode() == 1) || (kstar.photonPosTrackCode() == 1 && kstar.photonNegTrackCode() != 1)) + TrkCode = 2; + if (kstar.photonPosTrackCode() == 3 && kstar.photonNegTrackCode() == 3) + TrkCode = 3; + if (kstar.photonPosTrackCode() == 2 || kstar.photonNegTrackCode() == 2) + TrkCode = 4; + } else { + if (kstar.kshortPosTrackCode() == 1 && kstar.kshortNegTrackCode() == 1) + TrkCode = 1; + if ((kstar.kshortPosTrackCode() != 1 && kstar.kshortNegTrackCode() == 1) || (kstar.kshortPosTrackCode() == 1 && kstar.kshortNegTrackCode() != 1)) + TrkCode = 2; + if (kstar.kshortPosTrackCode() == 3 && kstar.kshortNegTrackCode() == 3) + TrkCode = 3; + if (kstar.kshortPosTrackCode() == 2 || kstar.kshortNegTrackCode() == 2) + TrkCode = 4; + } + + return TrkCode; + } + + template + void getResolution(TKStarObject const& kstar) + { + + //_______________________________________ + // Gamma MC association + if (kstar.photonPDGCode() == PDG_t::kGamma) { + if (kstar.photonmcpt() > 0) { + histos.fill(HIST("BeforeSel/MC/Reso/h3dGammaPtResoVsTPCCR"), 1.f / kstar.kshortmcpt(), 1.f / kstar.kshortPt() - 1.f / kstar.kshortmcpt(), -1 * kstar.photonNegTPCCrossedRows()); // 1/pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h3dGammaPtResoVsTPCCR"), 1.f / kstar.kshortmcpt(), 1.f / kstar.kshortPt() - 1.f / kstar.kshortmcpt(), kstar.photonPosTPCCrossedRows()); // 1/pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaPtResolution"), 1.f / kstar.photonmcpt(), 1.f / kstar.photonPt() - 1.f / kstar.photonmcpt()); // pT resolution + } + } + + //_______________________________________ + // KShort MC association + if (kstar.kshortPDGCode() == PDG_t::kK0Short) { + if (kstar.kshortmcpt() > 0) { + histos.fill(HIST("BeforeSel/MC/Reso/h2dKShortPtResolution"), 1.f / kstar.kshortmcpt(), 1.f / kstar.kshortPt() - 1.f / kstar.kshortmcpt()); // 1/pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h3dKShortPtResoVsTPCCR"), 1.f / kstar.kshortmcpt(), 1.f / kstar.kshortPt() - 1.f / kstar.kshortmcpt(), -1 * kstar.kshortNegTPCCrossedRows()); // 1/pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h3dKShortPtResoVsTPCCR"), 1.f / kstar.kshortmcpt(), 1.f / kstar.kshortPt() - 1.f / kstar.kshortmcpt(), kstar.kshortPosTPCCrossedRows()); // 1/pT resolution + } + } + + //_______________________________________ + // KStar MC association + if (kstar.isKStar()) { + histos.fill(HIST("BeforeSel/MC/Reso/h2dKStarRadiusResolution"), kstar.mcpt(), kstar.radius() - kstar.mcradius()); // pT resolution + if (kstar.mcpt() > 0) + histos.fill(HIST("BeforeSel/MC/Reso/h2dKStarPtResolution"), 1.f / kstar.mcpt(), 1.f / kstar.pt() - 1.f / kstar.mcpt()); // pT resolution + } + } + + // To save histograms for background analysis + template + void runBkgAnalysis(TKStarObject const& kstar) + { + // Check whether it is before or after selections + static constexpr std::string_view MainDir[] = {"BeforeSel", "AfterSel"}; + + bool fIsKStar = kstar.isKStar(); + int PhotonPDGCode = kstar.photonPDGCode(); + int PhotonPDGCodeMother = kstar.photonPDGCodeMother(); + int KShortPDGCode = kstar.kshortPDGCode(); + int KShortPDGCodeMother = kstar.kshortPDGCodeMother(); + float kstarpT = kstar.pt(); + float kstarMass = kstar.kstarMass(); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_All"), kstarpT, kstarMass); + + //_______________________________________ + // Real Gamma x Real KShort - but not from the same kstar! + if ((PhotonPDGCode == PDG_t::kGamma) && (KShortPDGCode == PDG_t::kK0Short) && (!fIsKStar)) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_TrueDaughters"), kstarpT, kstarMass); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dTrueDaughtersMatrix"), KShortPDGCodeMother, PhotonPDGCodeMother); + } + + //_______________________________________ + // Real Gamma x fake KShort + if ((PhotonPDGCode == PDG_t::kGamma) && (KShortPDGCode != PDG_t::kK0Short)) + histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_TrueGammaFakeKShort"), kstarpT, kstarMass); + + //_______________________________________ + // Fake Gamma x Real KShort + if ((PhotonPDGCode != PDG_t::kGamma) && ((KShortPDGCode == PDG_t::kK0Short))) + histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_FakeGammaTrueKShort"), kstarpT, kstarMass); + + //_______________________________________ + // Fake Gamma x Fake KShort + if ((PhotonPDGCode != PDG_t::kGamma) && (KShortPDGCode != PDG_t::kK0Short)) + histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_FakeDaughters"), kstarpT, kstarMass); + } + + template + void fillHistos(TKStarObject const& kstar, TCollision const& collision) + { + + // Check whether it is before or after selections + static constexpr std::string_view MainDir[] = {"BeforeSel", "AfterSel"}; + + // Get V0trackCode + int GammaTrkCode = retrieveV0TrackCode(kstar); + int KShortTrkCode = retrieveV0TrackCode(kstar); + + float photonRZLineCut = TMath::Abs(kstar.photonZconv()) * TMath::Tan(2 * TMath::ATan(TMath::Exp(-photonSelections.PhotonMaxDauEta))) - photonSelections.PhotonLineCutZ0; + float centrality = doPPAnalysis ? collision.centFT0M() : collision.centFT0C(); + //_______________________________________ + // Photon + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hTrackCode"), GammaTrkCode); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hV0Type"), kstar.photonV0Type()); + + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCANegToPV"), kstar.photonDCANegPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCAPosToPV"), kstar.photonDCAPosPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCADau"), kstar.photonDCADau()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosTPCCR"), kstar.photonPosTPCCrossedRows()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegTPCCR"), kstar.photonNegTPCCrossedRows()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosTPCNSigmaEl"), kstar.photonPosTPCNSigmaEl()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegTPCNSigmaEl"), kstar.photonNegTPCNSigmaEl()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hpT"), kstar.photonPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hY"), kstar.photonY()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosEta"), kstar.photonPosEta()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegEta"), kstar.photonNegEta()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hRadius"), kstar.photonRadius()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hZ"), kstar.photonZconv()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h2dRZCut"), kstar.photonRadius(), photonRZLineCut); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h2dRZPlane"), kstar.photonZconv(), kstar.photonRadius()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hCosPA"), kstar.photonCosPA()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPsiPair"), kstar.photonPsiPair()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPhi"), kstar.photonPhi()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dMass"), centrality, kstar.photonPt(), kstar.photonMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hMass"), kstar.photonMass()); + + //_______________________________________ + // KShorts + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hTrackCode"), KShortTrkCode); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hRadius"), kstar.kshortRadius()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hDCADau"), kstar.kshortDCADau()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hCosPA"), kstar.kshortCosPA()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hY"), kstar.kshortY()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hPosEta"), kstar.kshortPosEta()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hNegEta"), kstar.kshortNegEta()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hPosTPCCR"), kstar.kshortPosTPCCrossedRows()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hNegTPCCR"), kstar.kshortNegTPCCrossedRows()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hPosITSCls"), kstar.kshortPosITSCls()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hNegITSCls"), kstar.kshortNegITSCls()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hPosChi2PerNc"), kstar.kshortPosChi2PerNcl()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hNegChi2PerNc"), kstar.kshortNegChi2PerNcl()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hLifeTime"), kstar.kshortLifeTime()); + + //_______________________________________ + // Sigmas and KShorts + histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), kstar.photonAlpha(), kstar.photonQt()); + histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), kstar.kshortAlpha(), kstar.kshortQt()); + + if (kstar.kshortAlpha() < 1) { + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/h2dTPCvsTOFNSigma_KShortPi"), kstar.kshortPosPiTPCNSigma(), kstar.kshortPiTOFNSigma()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/h2dTPCvsTOFNSigma_KShortPi"), kstar.kshortNegPiTPCNSigma(), kstar.kshortPiTOFNSigma()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hKShortDCANegToPV"), kstar.kshortDCANegPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hKShortDCAPosToPV"), kstar.kshortDCAPosPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hKShortpT"), kstar.kshortPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/hKShortMass"), kstar.kshortMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/KShort/h3dKShortMass"), centrality, kstar.kshortPt(), kstar.kshortMass()); + + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/hMass"), kstar.kstarMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/hPt"), kstar.pt()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/hY"), kstar.kstarY()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/hRadius"), kstar.radius()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/h2dRadiusVspT"), kstar.radius(), kstar.pt()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/hDCAPairDau"), kstar.dcadaughters()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/hDCAPairDauVsPt"), kstar.dcadaughters(), kstar.pt()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/h3dMass"), centrality, kstar.pt(), kstar.kstarMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/h3dOPAngleVsMass"), kstar.opAngle(), kstar.pt(), kstar.kstarMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/KStar/h2dOPAngleVsPt"), kstar.opAngle(), kstar.pt()); + } + //_______________________________________ + // MC specific + if (doprocessMonteCarlo) { + if constexpr (requires { kstar.kshortPDGCode(); kstar.photonPDGCode(); }) { + + //_______________________________________ + // Gamma MC association + if (kstar.photonPDGCode() == PDG_t::kGamma) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hV0ToCollAssoc"), kstar.photonIsCorrectlyAssoc()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPt"), kstar.photonPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hMCPt"), kstar.photonmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPosTPCNSigmaEl"), kstar.photonPosTPCNSigmaEl()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hNegTPCNSigmaEl"), kstar.photonNegTPCNSigmaEl()); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt"), TMath::ACos(kstar.photonCosPA()), kstar.photonmcpt()); + + if (!kstar.photonIsCorrectlyAssoc()) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPt_BadCollAssig"), kstar.photonmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt_BadCollAssig"), TMath::ACos(kstar.photonCosPA()), kstar.photonmcpt()); + } + } + + //_______________________________________ + // KShort MC association + if (kstar.kshortPDGCode() == PDG_t::kK0Short) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KShort/hV0ToCollAssoc"), kstar.kshortIsCorrectlyAssoc()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KShort/hPt"), kstar.kshortPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KShort/hMCPt"), kstar.kshortmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KShort/hMass"), kstar.kshortMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KShort/h3dTPCvsTOFNSigma_Pi"), kstar.kshortPosPiTPCNSigma(), kstar.kshortPiTOFNSigma(), kstar.kshortPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KShort/h3dTPCvsTOFNSigma_Pi"), kstar.kshortNegPiTPCNSigma(), kstar.kshortPiTOFNSigma(), kstar.kshortPt()); + } + + //_______________________________________ + // KStar MC association + if (kstar.isKStar()) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/h2dArmenteros"), kstar.photonAlpha(), kstar.photonQt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/h2dArmenteros"), kstar.kshortAlpha(), kstar.kshortQt()); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/hPt"), kstar.pt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/hMCPt"), kstar.mcpt()); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h2dMCPtVsKShortMCPt"), kstar.mcpt(), kstar.kshortmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h2dMCPtVsPhotonMCPt"), kstar.mcpt(), kstar.photonmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/hMass"), kstar.kstarMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h3dMass"), centrality, kstar.mcpt(), kstar.kstarMass()); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/hMCProcess"), kstar.mcprocess()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/hGenRadius"), kstar.mcradius()); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h2dMCProcessVsGenRadius"), kstar.mcprocess(), kstar.mcradius()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h3dMCProcess"), kstar.mcprocess(), kstar.mcpt(), kstar.kstarMass()); + + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h2dOPAngleVsPt"), kstar.opAngle(), kstar.pt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/h2dRadiusVspT"), kstar.radius(), kstar.pt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/KStar/hDCAPairDauVsPt"), kstar.dcadaughters(), kstar.pt()); + } + + // For background studies: + if (fillBkgQAhistos) + runBkgAnalysis(kstar); + + //_______________________________________ + // pT resolution histos + if ((mode == 0) && fillResoQAhistos) + getResolution(kstar); + } + } + } + + template + void fillSelHistos(TKStarObject const& kstar, int PDGRequired) + { + + static constexpr std::string_view PhotonSelsLocal[] = {"NoSel", "V0Type", "DCADauToPV", + "DCADau", "DauTPCCR", "TPCNSigmaEl", "V0pT", + "Y", "V0Radius", "RZCut", "Armenteros", "CosPA", + "PsiPair", "Phi", "Mass"}; + + static constexpr std::string_view KShortSelsLocal[] = {"NoSel", "V0Radius", "DCADau", "Armenteros", + "CosPA", "Y", "TPCCR", "DauITSCls", "Lifetime", + "TPCTOFPID", "DCADauToPV", "Mass"}; + + if (PDGRequired == PDG_t::kGamma) { + if constexpr (selection_index >= 0 && selection_index < (int)std::size(PhotonSelsLocal)) { + histos.fill(HIST("Selection/Photon/hCandidateSel"), selection_index); + histos.fill(HIST("Selection/Photon/h2d") + HIST(PhotonSelsLocal[selection_index]), kstar.photonPt(), kstar.photonMass()); + histos.fill(HIST("Selection/KStar/h2dPhoton") + HIST(PhotonSelsLocal[selection_index]), kstar.pt(), kstar.kstarMass()); + } + } + + if (PDGRequired == PDG_t::kK0Short) { + if constexpr (selection_index >= 0 && selection_index < (int)std::size(KShortSelsLocal)) { + histos.fill(HIST("Selection/KShort/hCandidateSel"), selection_index); + histos.fill(HIST("Selection/KShort/h2d") + HIST(KShortSelsLocal[selection_index]), kstar.kshortPt(), kstar.kshortMass()); + histos.fill(HIST("Selection/KStar/h2dKShort") + HIST(KShortSelsLocal[selection_index]), kstar.pt(), kstar.kstarMass()); + } + } + } + + // Apply specific selections for photons + template + bool selectPhoton(TV0Object const& cand) + { + fillSelHistos<0>(cand, PDG_t::kGamma); + if (cand.photonV0Type() != photonSelections.Photonv0TypeSel && photonSelections.Photonv0TypeSel > -1) + return false; + + fillSelHistos<1>(cand, PDG_t::kGamma); + if ((TMath::Abs(cand.photonDCAPosPV()) < photonSelections.PhotonMinDCADauToPv) || (TMath::Abs(cand.photonDCANegPV()) < photonSelections.PhotonMinDCADauToPv)) + return false; + + fillSelHistos<2>(cand, PDG_t::kGamma); + if (TMath::Abs(cand.photonDCADau()) > photonSelections.PhotonMaxDCAV0Dau) + return false; + + fillSelHistos<3>(cand, PDG_t::kGamma); + if ((cand.photonPosTPCCrossedRows() < photonSelections.PhotonMinTPCCrossedRows) || (cand.photonNegTPCCrossedRows() < photonSelections.PhotonMinTPCCrossedRows)) + return false; + + fillSelHistos<4>(cand, PDG_t::kGamma); + if (((cand.photonPosTPCNSigmaEl() < photonSelections.PhotonMinTPCNSigmas) || (cand.photonPosTPCNSigmaEl() > photonSelections.PhotonMaxTPCNSigmas))) + return false; + + if (((cand.photonNegTPCNSigmaEl() < photonSelections.PhotonMinTPCNSigmas) || (cand.photonNegTPCNSigmaEl() > photonSelections.PhotonMaxTPCNSigmas))) + return false; + + fillSelHistos<5>(cand, PDG_t::kGamma); + if ((cand.photonPt() < photonSelections.PhotonMinPt) || (cand.photonPt() > photonSelections.PhotonMaxPt)) + return false; + + fillSelHistos<6>(cand, PDG_t::kGamma); + if ((TMath::Abs(cand.photonY()) > photonSelections.PhotonMaxRap) || (TMath::Abs(cand.photonPosEta()) > photonSelections.PhotonMaxDauEta) || (TMath::Abs(cand.photonNegEta()) > photonSelections.PhotonMaxDauEta)) + return false; + + fillSelHistos<7>(cand, PDG_t::kGamma); + if ((cand.photonRadius() < photonSelections.PhotonMinRadius) || (cand.photonRadius() > photonSelections.PhotonMaxRadius)) + return false; + + fillSelHistos<8>(cand, PDG_t::kGamma); + float photonRZLineCut = TMath::Abs(cand.photonZconv()) * TMath::Tan(2 * TMath::ATan(TMath::Exp(-photonSelections.PhotonMaxDauEta))) - photonSelections.PhotonLineCutZ0; + if ((TMath::Abs(cand.photonRadius()) < photonRZLineCut) || (TMath::Abs(cand.photonZconv()) > photonSelections.PhotonMaxZ)) + return false; + + fillSelHistos<9>(cand, PDG_t::kGamma); + if (cand.photonQt() > photonSelections.PhotonMaxQt) + return false; + + if (TMath::Abs(cand.photonAlpha()) > photonSelections.PhotonMaxAlpha) + return false; + + fillSelHistos<10>(cand, PDG_t::kGamma); + if (cand.photonCosPA() < photonSelections.PhotonMinV0cospa) + return false; + + fillSelHistos<11>(cand, PDG_t::kGamma); + if (TMath::Abs(cand.photonPsiPair()) > photonSelections.PhotonPsiPairMax) + return false; + + fillSelHistos<12>(cand, PDG_t::kGamma); + if ((((cand.photonPhi() > photonSelections.PhotonPhiMin1) && (cand.photonPhi() < photonSelections.PhotonPhiMax1)) || ((cand.photonPhi() > photonSelections.PhotonPhiMin2) && (cand.photonPhi() < photonSelections.PhotonPhiMax2))) && ((photonSelections.PhotonPhiMin1 != -1) && (photonSelections.PhotonPhiMax1 != -1) && (photonSelections.PhotonPhiMin2 != -1) && (photonSelections.PhotonPhiMax2 != -1))) + return false; + + fillSelHistos<13>(cand, PDG_t::kGamma); + if (TMath::Abs(cand.photonMass()) > photonSelections.PhotonMaxMass) + return false; + + fillSelHistos<14>(cand, PDG_t::kGamma); + return true; + } + + // Apply specific selections for kshortrs + template + bool selectKShort(TV0Object const& cand) + { + fillSelHistos<0>(cand, PDG_t::kK0Short); + if ((cand.kshortRadius() < kshortSelections.KShortMinv0radius) || (cand.kshortRadius() > kshortSelections.KShortMaxv0radius)) + return false; + + fillSelHistos<1>(cand, PDG_t::kK0Short); + if (TMath::Abs(cand.kshortDCADau()) > kshortSelections.KShortMaxDCAV0Dau) + return false; + + fillSelHistos<2>(cand, PDG_t::kK0Short); + if ((cand.kshortQt() < kshortSelections.KShortMinQt) || (cand.kshortQt() > kshortSelections.KShortMaxQt)) + return false; + + if ((TMath::Abs(cand.kshortAlpha()) < kshortSelections.KShortMinAlpha) || (TMath::Abs(cand.kshortAlpha()) > kshortSelections.KShortMaxAlpha)) + return false; + + fillSelHistos<3>(cand, PDG_t::kK0Short); + if (cand.kshortCosPA() < kshortSelections.KShortMinv0cospa) + return false; + + fillSelHistos<4>(cand, PDG_t::kK0Short); + if ((TMath::Abs(cand.kshortY()) > kshortSelections.KShortMaxRap) || (TMath::Abs(cand.kshortPosEta()) > kshortSelections.KShortMaxDauEta) || (TMath::Abs(cand.kshortNegEta()) > kshortSelections.KShortMaxDauEta)) + return false; + + fillSelHistos<5>(cand, PDG_t::kK0Short); + if ((cand.kshortPosTPCCrossedRows() < kshortSelections.KShortMinTPCCrossedRows) || (cand.kshortNegTPCCrossedRows() < kshortSelections.KShortMinTPCCrossedRows)) + return false; + + fillSelHistos<6>(cand, PDG_t::kK0Short); + + // check minimum number of ITS clusters + reject ITS afterburner tracks if requested + bool posIsFromAfterburner = cand.kshortPosChi2PerNcl() < 0; + bool negIsFromAfterburner = cand.kshortNegChi2PerNcl() < 0; + if (cand.kshortPosITSCls() < kshortSelections.KShortMinITSclusters && (!kshortSelections.KShortRejectPosITSafterburner || posIsFromAfterburner)) + return false; + if (cand.kshortNegITSCls() < kshortSelections.KShortMinITSclusters && (!kshortSelections.KShortRejectNegITSafterburner || negIsFromAfterburner)) + return false; + + fillSelHistos<7>(cand, PDG_t::kK0Short); + if (cand.kshortLifeTime() > kshortSelections.KShortMaxLifeTime) + return false; + + // Separating kshort selections: + fillSelHistos<8>(cand, PDG_t::kK0Short); + if (cand.kshortAlpha() < 1) { // KShort selection + // TPC Selection + if (kshortSelections.fselKShortTPCPID && (TMath::Abs(cand.kshortPosPiTPCNSigma()) > kshortSelections.KShortMaxTPCNSigmas)) + return false; + if (kshortSelections.fselKShortTPCPID && (TMath::Abs(cand.kshortNegPiTPCNSigma()) > kshortSelections.KShortMaxTPCNSigmas)) + return false; + + // // TOF Selection + // if (kshortSelections.fselKShortTOFPID && (TMath::Abs(cand.kshortPiTOFNSigma()) > kshortSelections.KShortPiMaxTOFNSigmas)) + // return false; + // if (kshortSelections.fselKShortTOFPID && (TMath::Abs(cand.lambdaPiTOFNSigma()) > kshortSelections.KShortPiMaxTOFNSigmas)) + // return false; + + // DCA Selection + fillSelHistos<9>(cand, PDG_t::kK0Short); + if ((TMath::Abs(cand.kshortDCAPosPV()) < kshortSelections.KShortMinDCAPosToPv) || (TMath::Abs(cand.kshortDCANegPV()) < kshortSelections.KShortMinDCANegToPv)) + return false; + + // Mass Selection + fillSelHistos<10>(cand, PDG_t::kK0Short); + if (TMath::Abs(cand.kshortMass() - o2::constants::physics::MassK0Short) > kshortSelections.KShortWindow) + return false; + + fillSelHistos<11>(cand, PDG_t::kK0Short); + } + return true; + } + + // Apply selections in kdyst candidates + template + bool processKStarCandidate(TKStarObject const& cand) + { + // Photon specific selections + if (!selectPhoton(cand)) + return false; + + // KShort specific selections + if (!selectKShort(cand)) + return false; + + // KStar specific selections + // Rapidity + if constexpr (requires { cand.kstarMCY(); }) { // MC + if (TMath::Abs(cand.kstarMCY()) > kstarSelections.KStarMaxRap) + return false; + } else { // Real data + if (TMath::Abs(cand.kstarY()) > kstarSelections.KStarMaxRap) + return false; + } + + // V0Pair Radius + if (cand.radius() > kstarSelections.KStarMaxRadius) + return false; + + // DCA V0Pair Daughters + if (cand.dcadaughters() > kstarSelections.KStarMaxDCADau) + return false; + + // Opening Angle + if (cand.opAngle() > kstarSelections.KStarMaxOPAngle) + return false; + + return true; + } + + // Main analysis function + template + void analyzeRecoeKStars(TCollisions const& collisions, TKStars const& fullKStars) + { + // Custom grouping + std::vector> kstargrouped(collisions.size()); + + for (const auto& kstar : fullKStars) { + kstargrouped[kstar.straCollisionId()].push_back(kstar.globalIndex()); + } + + // Collisions loop + for (const auto& coll : collisions) { + + // Event selection + if (!IsEventAccepted(coll, true)) + continue; + + // KStars loop + for (size_t i = 0; i < kstargrouped[coll.globalIndex()].size(); i++) { + auto kstar = fullKStars.rawIteratorAt(kstargrouped[coll.globalIndex()][i]); + + // if MC + if constexpr (requires { kstar.isKStar(); }) { + if (doMCAssociation && !(kstar.isKStar())) + continue; + + if (selRecoFromGenerator && !kstar.isProducedByGenerator()) + continue; + } + + // Fill histos before any selection + fillHistos<0>(kstar, coll); + + // Select kstar candidates + if (!processKStarCandidate(kstar)) + continue; + + // Fill histos after all selections + fillHistos<1>(kstar, coll); + } + } + } + + void processRealData(soa::Join const& collisions, KStars const& fullKStars) + { + analyzeRecoeKStars(collisions, fullKStars); + } + + void processMonteCarlo(soa::Join const& collisions, MCKStars const& fullKStars) + { + analyzeRecoeKStars(collisions, fullKStars); + } + + // Simulated processing in Run 3 + void processGeneratedRun3(soa::Join const& mcCollisions, soa::Join const& collisions, soa::Join const& KStarGens) + { + analyzeGenerated(mcCollisions, collisions, KStarGens); + } + + // _____________________________________________________ + PROCESS_SWITCH(k892hadronphoton, processRealData, "Do real data analysis", true); + PROCESS_SWITCH(k892hadronphoton, processMonteCarlo, "Do Monte-Carlo-based analysis", false); + PROCESS_SWITCH(k892hadronphoton, processGeneratedRun3, "process MC generated Run 3", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}