diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/Decay3Body.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/Decay3Body.h index eb7229252b586..f7f18976fa508 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/Decay3Body.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/Decay3Body.h @@ -29,7 +29,7 @@ class Decay3Body : public o2::track::TrackParCov /// TO BE DONE: extend to gener using PID = o2::track::PID; Decay3Body() = default; - Decay3Body(PID pid, const std::array& xyz, const std::array& pxyz, const std::array& covxyz, const Track& tr0, const Track& tr1, const Track& tr2); + Decay3Body(const std::array& xyz, const std::array& pxyz, const std::array& covxyz, const Track& tr0, const Track& tr1, const Track& tr2, o2::track::PID pid = o2::track::PID::HyperTriton); const Track& getProng(int i) const { return mProngs[i]; } Track& getProng(int i) { return mProngs[i]; } diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h index e01ce253156a7..ce70e69aa6ddd 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h @@ -33,19 +33,19 @@ namespace o2cp = o2::constants::physics; namespace pid_constants // GPUs currently cannot have static constexpr array members { typedef uint8_t ID; -static constexpr ID NIDsTot = 17; +static constexpr ID NIDsTot = 19; #if !defined(GPUCA_GPUCODE_DEVICE) || defined(GPUCA_GPU_DEBUG_PRINT) GPUconstexpr() const char* sNames[NIDsTot + 1] = ///< defined particle names {"Electron", "Muon", "Pion", "Kaon", "Proton", "Deuteron", "Triton", "He3", "Alpha", - "Pion0", "Photon", "K0", "Lambda", "HyperTriton", "Hyperhydrog4", "XiMinus", "OmegaMinus", nullptr}; + "Pion0", "Photon", "K0", "Lambda", "HyperTriton", "Hyperhydrog4", "XiMinus", "OmegaMinus", "HyperHelium4", "HyperHelium5", nullptr}; #endif GPUconstexpr() const float sMasses[NIDsTot] = ///< defined particle masses {o2cp::MassElectron, o2cp::MassMuon, o2cp::MassPionCharged, o2cp::MassKaonCharged, o2cp::MassProton, o2cp::MassDeuteron, o2cp::MassTriton, o2cp::MassHelium3, o2cp::MassAlpha, o2cp::MassPionNeutral, o2cp::MassPhoton, - o2cp::MassKaonNeutral, o2cp::MassLambda, o2cp::MassHyperTriton, o2cp::MassHyperhydrog4, o2cp::MassXiMinus, o2cp::MassOmegaMinus}; + o2cp::MassKaonNeutral, o2cp::MassLambda, o2cp::MassHyperTriton, o2cp::MassHyperhydrog4, o2cp::MassXiMinus, o2cp::MassOmegaMinus, o2cp::MassHyperHelium4, o2cp::MassHyperHelium5}; GPUconstexpr() const float sMasses2[NIDsTot] = ///< defined particle masses^2 {o2cp::MassElectron * o2cp::MassElectron, @@ -64,7 +64,9 @@ GPUconstexpr() const float sMasses2[NIDsTot] = ///< defined particle masses^2 o2cp::MassHyperTriton* o2cp::MassHyperTriton, o2cp::MassHyperhydrog4* o2cp::MassHyperhydrog4, o2cp::MassXiMinus* o2cp::MassXiMinus, - o2cp::MassOmegaMinus* o2cp::MassOmegaMinus}; + o2cp::MassOmegaMinus* o2cp::MassOmegaMinus, + o2cp::MassHyperHelium4* o2cp::MassHyperHelium4, + o2cp::MassHyperHelium5* o2cp::MassHyperHelium5}; GPUconstexpr() const float sMasses2Z[NIDsTot] = ///< defined particle masses / Z {o2cp::MassElectron, o2cp::MassMuon, @@ -73,12 +75,14 @@ GPUconstexpr() const float sMasses2Z[NIDsTot] = ///< defined particle masses / Z o2cp::MassTriton, o2cp::MassHelium3 / 2., o2cp::MassAlpha / 2., 0, 0, 0, 0, o2cp::MassHyperTriton, o2cp::MassHyperhydrog4, - o2cp::MassXiMinus, o2cp::MassOmegaMinus}; + o2cp::MassXiMinus, o2cp::MassOmegaMinus, + o2cp::MassHyperHelium4 / 2., o2cp::MassHyperHelium5 / 2.}; GPUconstexpr() const int sCharges[NIDsTot] = ///< defined particle charges {1, 1, 1, 1, 1, 1, 1, 2, 2, 0, 0, 0, 0, 1, 1, - 1, 1}; + 1, 1, + 2, 2}; } // namespace pid_constants class PID @@ -110,8 +114,10 @@ class PID static constexpr ID Hyperhydrog4 = 14; static constexpr ID XiMinus = 15; static constexpr ID OmegaMinus = 16; + static constexpr ID HyperHelium4 = 17; + static constexpr ID HyperHelium5 = 18; static constexpr ID FirstExt = PI0; - static constexpr ID LastExt = OmegaMinus; + static constexpr ID LastExt = HyperHelium5; static constexpr ID NIDsTot = pid_constants::NIDsTot; ///< total number of defined IDs static_assert(NIDsTot == LastExt + 1, "Incorrect NIDsTot, please update!"); diff --git a/DataFormats/Reconstruction/src/Decay3Body.cxx b/DataFormats/Reconstruction/src/Decay3Body.cxx index 46f5c9447c2fb..aa071cea675cd 100644 --- a/DataFormats/Reconstruction/src/Decay3Body.cxx +++ b/DataFormats/Reconstruction/src/Decay3Body.cxx @@ -13,7 +13,7 @@ using namespace o2::dataformats; -Decay3Body::Decay3Body(PID pid, const std::array& xyz, const std::array& pxyz, const std::array& covxyz, const Track& tr0, const Track& tr1, const Track& tr2) +Decay3Body::Decay3Body(const std::array& xyz, const std::array& pxyz, const std::array& covxyz, const Track& tr0, const Track& tr1, const Track& tr2, o2::track::PID pid) : mProngs{tr0, tr1, tr2} { std::array cov{}, cov1{}, cov2{}; diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h index 0510290fb5ae8..1450e0c15e98c 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h @@ -131,10 +131,14 @@ class SVertex3Hypothesis void set(PID v0, PID ppos, PID pneg, PID pbach, float sig, float nSig, float margin, float cpt, float bz = 0.f); void set(PID v0, PID ppos, PID pneg, PID pbach, const float pars[NPIDParams], float bz = 0.f); + PID getPIDHyp() const { return mPIDV0; } float getMassV0Hyp() const { return PID::getMass(mPIDV0); } float getMassPosProng() const { return PID::getMass(mPIDPosProng); } float getMassNegProng() const { return PID::getMass(mPIDNegProng); } float getMassBachProng() const { return PID::getMass(mPIDBachProng); } + float getChargePosProng() const { return PID::getCharge(mPIDPosProng); } + float getChargeNegProng() const { return PID::getCharge(mPIDNegProng); } + float getChargeBachProng() const { return PID::getCharge(mPIDBachProng); } float calcMass2(float p2Pos, float p2Neg, float p2Bach, float p2Tot) const { diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h index 9d0823a04e625..b933363bb352d 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h @@ -87,6 +87,8 @@ class SVertexer enum Hyp3body { H3L3body, AntiH3L3body, + H4L3body, + AntiH4L3body, He4L3body, AntiHe4L3body, He5L3body, diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h index b1ce6c1dfd4ad..1a21cf1d89393 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h @@ -138,7 +138,9 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper d p pi- + float pidCutsH4L3body[SVertex3Hypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5}; // H4L -> t p pi- float pidCutsHe4L3body[SVertex3Hypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5}; // He4L -> He3 p pi- + float pidCutsHe5L3body[SVertex3Hypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5}; // He5L -> He4 p pi- O2ParamDef(SVertexerParams, "svertexer"); }; diff --git a/Detectors/Vertexing/src/SVertexer.cxx b/Detectors/Vertexing/src/SVertexer.cxx index 78e0ce309d6ba..1d48bcceb0097 100644 --- a/Detectors/Vertexing/src/SVertexer.cxx +++ b/Detectors/Vertexing/src/SVertexer.cxx @@ -298,6 +298,12 @@ void SVertexer::updateTimeDependentParams() m3bodyHyps[Hyp3body::H3L3body].set(PID::HyperTriton, PID::Proton, PID::Pion, PID::Deuteron, mSVParams->pidCutsH3L3body, bz); m3bodyHyps[Hyp3body::AntiH3L3body].set(PID::HyperTriton, PID::Pion, PID::Proton, PID::Deuteron, mSVParams->pidCutsH3L3body, bz); + m3bodyHyps[Hyp3body::H4L3body].set(PID::Hyperhydrog4, PID::Proton, PID::Pion, PID::Triton, mSVParams->pidCutsH4L3body, bz); + m3bodyHyps[Hyp3body::AntiH4L3body].set(PID::Hyperhydrog4, PID::Pion, PID::Proton, PID::Triton, mSVParams->pidCutsH4L3body, bz); + m3bodyHyps[Hyp3body::He4L3body].set(PID::HyperHelium4, PID::Proton, PID::Pion, PID::Helium3, mSVParams->pidCutsHe4L3body, bz); + m3bodyHyps[Hyp3body::AntiHe4L3body].set(PID::HyperHelium4, PID::Pion, PID::Proton, PID::Helium3, mSVParams->pidCutsHe4L3body, bz); + m3bodyHyps[Hyp3body::He5L3body].set(PID::HyperHelium5, PID::Proton, PID::Pion, PID::Alpha, mSVParams->pidCutsHe5L3body, bz); + m3bodyHyps[Hyp3body::AntiHe5L3body].set(PID::HyperHelium5, PID::Pion, PID::Proton, PID::Alpha, mSVParams->pidCutsHe5L3body, bz); for (auto& ft : mFitterV0) { ft.setBz(bz); @@ -1090,52 +1096,60 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s tr0.getPxPyPzGlo(p0); tr1.getPxPyPzGlo(p1); tr2.getPxPyPzGlo(p2); - std::array p3B = {p0[0] + p1[0] + p2[0], p0[1] + p1[1] + p2[1], p0[2] + p1[2] + p2[2]}; - float pt2candidate = p3B[0] * p3B[0] + p3B[1] * p3B[1], p2candidate = pt2candidate + p3B[2] * p3B[2]; - if (pt2candidate < mMinPt23Body) { // pt cut - continue; - } - if (p3B[2] * p3B[2] / pt2candidate > mMaxTgl23Body) { // tgLambda cut - continue; - } - - // compute primary vertex and cosPA of the 3-body decay - auto bestCosPA = mSVParams->minCosPA3body; + bool goodHyp = false; + o2::track::PID pidHyp = o2::track::PID::Electron; // Update if goodHyp is true auto decay3bodyVtxID = -1; + auto vtxCosPA = -1; + + std::array pbach = {0, 0, 0}, p3B = {0, 0, 0}; // Update during the check of invariant mass + for (int ipid = 0; ipid < NHyp3body; ipid++) { + // check mass based on hypothesis of charge of bachelor (pos and neg expected to be proton/pion) + float bachChargeFactor = m3bodyHyps[ipid].getChargeBachProng() / tr2.getAbsCharge(); + pbach = {bachChargeFactor * p2[0], bachChargeFactor * p2[1], bachChargeFactor * p2[2]}; + p3B = {p0[0] + p1[0] + pbach[0], p0[1] + p1[1] + pbach[1], p0[2] + p1[2] + pbach[2]}; + float sqP0 = p0[0] * p0[0] + p0[1] * p0[1] + p0[2] * p0[2], sqP1 = p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2], sqPBach = pbach[0] * pbach[0] + pbach[1] * pbach[1] + pbach[2] * pbach[2]; + float pt2Candidate = p3B[0] * p3B[0] + p3B[1] * p3B[1], p2Candidate = pt2Candidate + p3B[2] * p3B[2]; + float ptCandidate = std::sqrt(pt2Candidate); + if (m3bodyHyps[ipid].check(sqP0, sqP1, sqPBach, p2Candidate, ptCandidate)) { + if (pt2Candidate < mMinPt23Body) { // pt cut + continue; + } + if (p3B[2] * p3B[2] > pt2Candidate * mMaxTgl23Body) { // tgLambda cut + continue; + } - for (int iv = decay3bodyVlist.getMin(); iv <= decay3bodyVlist.getMax(); iv++) { - const auto& pv = mPVertices[iv]; - // check cos of pointing angle - float dx = vertexXYZ[0] - pv.getX(), dy = vertexXYZ[1] - pv.getY(), dz = vertexXYZ[2] - pv.getZ(), prodXYZ3body = dx * p3B[0] + dy * p3B[1] + dz * p3B[2]; - float cosPA = prodXYZ3body / std::sqrt((dx * dx + dy * dy + dz * dz) * p2candidate); - if (cosPA < bestCosPA) { - LOG(debug) << "Rej. cosPA: " << cosPA; - continue; - } - decay3bodyVtxID = iv; - bestCosPA = cosPA; - } - if (decay3bodyVtxID == -1) { - LOG(debug) << "3-body decay not compatible with any vertex"; - continue; - } - - const auto& decay3bodyPv = mPVertices[decay3bodyVtxID]; - float sqP0 = p0[0] * p0[0] + p0[1] * p0[1] + p0[2] * p0[2], sqP1 = p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2], sqP2 = p2[0] * p2[0] + p2[1] * p2[1] + p2[2] * p2[2]; - float pt3B = std::sqrt(pt2candidate); + // compute primary vertex and cosPA of the 3-body decay + auto bestCosPA = mSVParams->minCosPA3body; + for (int iv = decay3bodyVlist.getMin(); iv <= decay3bodyVlist.getMax(); iv++) { + const auto& pv = mPVertices[iv]; + // check cos of pointing angle + float dx = vertexXYZ[0] - pv.getX(), dy = vertexXYZ[1] - pv.getY(), dz = vertexXYZ[2] - pv.getZ(), prodXYZ3body = dx * p3B[0] + dy * p3B[1] + dz * p3B[2]; + float cosPA = prodXYZ3body / std::sqrt((dx * dx + dy * dy + dz * dz) * p2Candidate); + if (cosPA < bestCosPA) { + LOG(debug) << "Rej. cosPA: " << cosPA; + continue; + } + decay3bodyVtxID = iv; + bestCosPA = cosPA; + } + if (decay3bodyVtxID == -1) { + LOG(debug) << "3-body decay not compatible with any vertex"; + continue; + } - bool goodHyp = false; - for (int ipid = 0; ipid < 2; ipid++) { // TODO: expand this loop to cover all the 3body cases if (m3bodyHyps[ipid].check(sqP0, sqP1, sqP2, sqPtot, pt3B)) - if (m3bodyHyps[ipid].check(sqP0, sqP1, sqP2, p2candidate, pt3B)) { goodHyp = true; + pidHyp = m3bodyHyps[ipid].getPIDHyp(); + vtxCosPA = bestCosPA; break; } } if (!goodHyp) { continue; } - Decay3Body candidate3B(PID::HyperTriton, vertexXYZ, p3B, fitter3body.calcPCACovMatrixFlat(cand3B), tr0, tr1, tr2); + + const auto& decay3bodyPv = mPVertices[decay3bodyVtxID]; + Decay3Body candidate3B(vertexXYZ, p3B, fitter3body.calcPCACovMatrixFlat(cand3B), tr0, tr1, tr2, pidHyp); o2::track::TrackParCov trc = candidate3B; o2::dataformats::DCA dca; if (!trc.propagateToDCA(decay3bodyPv, fitter3body.getBz(), &dca, 5.) || @@ -1143,7 +1157,7 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s continue; } if (mSVParams->createFull3Bodies) { - candidate3B.setCosPA(bestCosPA); + candidate3B.setCosPA(vtxCosPA); candidate3B.setDCA(fitter3body.getChi2AtPCACandidate()); m3bodyTmp[ithread].push_back(candidate3B); } diff --git a/Framework/Core/include/Framework/AnalysisDataModel.h b/Framework/Core/include/Framework/AnalysisDataModel.h index 0833c6e4aae65..9c939412cd941 100644 --- a/Framework/Core/include/Framework/AnalysisDataModel.h +++ b/Framework/Core/include/Framework/AnalysisDataModel.h @@ -1480,7 +1480,7 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Track2, track2, int, Tracks, "_2"); //! Track 2 in DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! Collision index } // namespace decay3body -DECLARE_SOA_TABLE(Decay3Bodys, "AOD", "DECAY3BODY", //! Run 2 cascade table +DECLARE_SOA_TABLE(Decay3Bodys, "AOD", "DECAY3BODY", //! 3-body decay table o2::soa::Index<>, decay3body::CollisionId, decay3body::Track0Id, decay3body::Track1Id, decay3body::Track2Id); using Decay3Bodys = Decay3Bodys; //! this defines the current default version