Add mass hypothesis for H4L, He4L, He5L and charge hypothesis for dec…

…ays3Body

Fix the description of decays3Body

Fix the charge of hyperhelium4

Set the ID of hyperhelium4 as the last one

Enable H4L and He4L

Add mass hypothesis of He5L

Use the hypothesis of daughter's charge for mass check

Add new calcMass function and code format fix

Remove space

Parameters fix

Move the momentum calculation of candidate outside

Only modify the charge hypothesis for bachelors

Small updates

Please consider the following formatting changes

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<float, 3>& xyz, const std::array<float, 3>& pxyz, const std::array<float, 6>& covxyz, const Track& tr0, const Track& tr1, const Track& tr2);
+  Decay3Body(const std::array<float, 3>& xyz, const std::array<float, 3>& pxyz, const std::array<float, 6>& 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]; }

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!");
 

DataFormats/Reconstruction/src/Decay3Body.cxx

@@ -13,7 +13,7 @@
 
 using namespace o2::dataformats;
 
-Decay3Body::Decay3Body(PID pid, const std::array<float, 3>& xyz, const std::array<float, 3>& pxyz, const std::array<float, 6>& covxyz, const Track& tr0, const Track& tr1, const Track& tr2)
+Decay3Body::Decay3Body(const std::array<float, 3>& xyz, const std::array<float, 3>& pxyz, const std::array<float, 6>& covxyz, const Track& tr0, const Track& tr1, const Track& tr2, o2::track::PID pid)
   : mProngs{tr0, tr1, tr2}
 {
   std::array<float, 21> cov{}, cov1{}, cov2{};

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
   {

Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h

@@ -87,6 +87,8 @@ class SVertexer
   enum Hyp3body {
     H3L3body,
     AntiH3L3body,
+    H4L3body,
+    AntiH4L3body,
     He4L3body,
     AntiHe4L3body,
     He5L3body,

Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h

@@ -138,7 +138,9 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper<SVertexerParam
   // cuts on different 3 body PID params
   bool check3bodyHypothesis = true;
   float pidCutsH3L3body[SVertex3Hypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5};  // H3L -> 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");
 };

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,60 +1096,68 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s
     tr0.getPxPyPzGlo(p0);
     tr1.getPxPyPzGlo(p1);
     tr2.getPxPyPzGlo(p2);
-    std::array<float, 3> 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<float, 3> 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.) ||
         std::abs(dca.getY()) > mSVParams->maxDCAXY3Body || std::abs(dca.getZ()) > mSVParams->maxDCAZ3Body) {
       continue;
     }
     if (mSVParams->createFull3Bodies) {
-      candidate3B.setCosPA(bestCosPA);
+      candidate3B.setCosPA(vtxCosPA);
       candidate3B.setDCA(fitter3body.getChi2AtPCACandidate());
       m3bodyTmp[ithread].push_back(candidate3B);
     }

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