Tighter windows when doing decay chains

Author: David Dobrigkeit Chinellato

Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h

@@ -33,14 +33,16 @@ class SVertexHypothesis
   enum PIDParams { SigmaM,  // sigma of mass res at 0 pt
                    NSigmaM, // number of sigmas of mass res
                    MarginM, // additive safety margin in mass cut
+                   NSigmaTightM, // number of sigmas of mass res when doing tight cut around mass (V0s used in cascades)
+                   MarginTightM, // additive safety margin in mass cut when doing tight cut around mass (V0s used in cascades)
                    CPt,     // pT dependence of mass resolution parameterized as mSigma*(1+mC1*pt);
                    CPt1,
                    CPt2,
                    CPt3 }; // pT dependence of mass resolution of Cascade parameterized as CPt+CPt1*pt +CPt2*TMath::Exp(-CPt3*pt);
 
-  static constexpr int NPIDParams = 7;
+  static constexpr int NPIDParams = 9;
 
-  void set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float cpt, float cpt1, float cpt2, float cpt3, float bz = 0.f, float maxSigma = 0.01);
+  void set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float nSigTight, float marginTight, float cpt, float cpt1, float cpt2, float cpt3, float bz = 0.f, float maxSigma = 0.01);
   void set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz = 0.f, float maxSigma = 0.01);
 
   float getMassV0Hyp() const { return PID::getMass(mPIDV0); }
@@ -60,28 +62,43 @@ class SVertexHypothesis
   { // check if given mass and pt is matching to hypothesis
     return check(calcMass(p2Pos, p2Neg, p2V0), ptV0);
   }
-
   bool check(float mass, float pt) const
   { // check if given mass and pt is matching to hypothesis
     return std::abs(mass - getMassV0Hyp()) < getMargin(pt);
   }
 
+  bool checkTight(float p2Pos, float p2Neg, float p2V0, float ptV0) const
+  { // check if given mass and pt is matching to hypothesis
+    return checkTight(calcMass(p2Pos, p2Neg, p2V0), ptV0);
+  }
+  bool checkTight(float mass, float pt) const
+  { // check if given mass and pt is matching to hypothesis
+    return std::abs(mass - getMassV0Hyp()) < getMarginTight(pt);
+  }
+
   float getSigmaV0Cascade(float pt) const { return mPars[CPt] + mPars[CPt1] * pt + mPars[CPt2] * std::exp(-mPars[CPt3] * pt); }
   float getSigma(float pt) const { return mPars[SigmaM] * (1.f + mPars[CPt] * pt); }
-  float getMargin(float pt) const
+  float getMargin(float pt, bool tight = false) const
   {
+    int idxNsigma = NSigmaM;
+    int idxMargin = MarginM;
+    if(tight){ // move to indices for tight variables in case asked to do so (tighter peak cuts for decay chains)
+      idxNsigma = NSigmaTightM;
+      idxMargin = MarginTightM;  
+    }
     if (mPIDV0 == PID::XiMinus || mPIDV0 == PID::OmegaMinus) { // case for cascades, antiparticles included
       if (getSigmaV0Cascade(pt) > maxSigma) {                  // insuring that at low pt one gets reasonable width as the parametrisation function may explode to unphysical values
-        return mPars[NSigmaM] * maxSigma + mPars[MarginM];
+        return mPars[idxNsigma] * maxSigma + mPars[idxMargin];
       } else {
-        return mPars[NSigmaM] * getSigmaV0Cascade(pt) + mPars[MarginM];
+        return mPars[idxNsigma] * getSigmaV0Cascade(pt) + mPars[idxMargin];
       }
     } else if (mPIDV0 == PID::K0 || mPIDV0 == PID::Lambda) { // case for V0s, AntiLambda is included in PID::Lambda
-      return mPars[NSigmaM] * getSigmaV0Cascade(pt) + mPars[MarginM];
+      return mPars[idxNsigma] * getSigmaV0Cascade(pt) + mPars[idxMargin];
     } else {
-      return mPars[NSigmaM] * getSigma(pt) + mPars[MarginM]; // case for HyperTriton and Hyperhydrog4
+      return mPars[idxNsigma] * getSigma(pt) + mPars[idxMargin]; // case for HyperTriton and Hyperhydrog4
     }
   }
+  float getMarginTight(float pt) const { return getMargin(pt, true); }
 
  private:
   float getMass2PosProng() const { return PID::getMass2(mPIDPosProng); }

Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h

@@ -96,16 +96,17 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper<SVertexerParam
 
   // cuts on different V0 PID params
   bool checkV0Hypothesis = true;
-  float pidCutsPhoton[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.60, 0.0, 0.0, 0.0, 0.0};    // Photon
-  float pidCutsK0[SVertexHypothesis::NPIDParams] = {0., 20, 0., 2.84798e-03, 9.84206e-04, 3.31951e-03, 2.39438};        // K0
-  float pidCutsLambda[SVertexHypothesis::NPIDParams] = {0., 20, 0., 1.09004e-03, 2.62291e-04, 8.93179e-03, 2.83121};    // Lambda
-  float pidCutsHTriton[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5, 0.0, 0.0, 0.0};  // HyperTriton
-  float pidCutsHhydrog4[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5, 0.0, 0.0, 0.0}; // Hyperhydrog4 - Need to update
+  float pidCutsPhoton[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.60, 20, 0.60, 0.0, 0.0, 0.0, 0.0};                        // Photon
+  float pidCutsK0[SVertexHypothesis::NPIDParams] = {0., 20, 0., 4.0, 0.0, 2.84798e-03, 9.84206e-04, 3.31951e-03, 2.39438};     // K0
+  float pidCutsLambda[SVertexHypothesis::NPIDParams] = {0., 20, 0., 4.0, 0.0, 1.09004e-03, 2.62291e-04, 8.93179e-03, 2.83121}; // Lambda
+  float pidCutsHTriton[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 14, 0.0, 0.5, 0.0, 0.0, 0.0};                      // HyperTriton
+  float pidCutsHhydrog4[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 14, 0.0, 0.5, 0.0, 0.0, 0.0};                     // Hyperhydrog4 - Need to update
   //
   // cuts on different Cascade PID params
   bool checkCascadeHypothesis = true;
-  float pidCutsXiMinus[SVertexHypothesis::NPIDParams] = {0.0, 10, 0.0, 1.56315e-03, 2.23279e-04, 2.75136e-02, 3.309};    // XiMinus
-  float pidCutsOmegaMinus[SVertexHypothesis::NPIDParams] = {0.0, 10, 0.0, 1.43572e-03, 6.94416e-04, 2.13534e+05, 1.48889e+01}; // OmegaMinus
+  float pidCutsXiMinus[SVertexHypothesis::NPIDParams] = {0.0, 10, 0.0, 4.0, 0.0, 1.56315e-03, 2.23279e-04, 2.75136e-02, 3.309};          // XiMinus
+  float pidCutsOmegaMinus[SVertexHypothesis::NPIDParams] = {0.0, 10, 0.0, 4.0, 0.0, 1.43572e-03, 6.94416e-04, 2.13534e+05, 1.48889e+01}; // OmegaMinus
+  float maximalCascadeWidth = 0.006;
   //
   // cuts on different 3 body PID params
   bool check3bodyHypothesis = true;

Detectors/Vertexing/src/SVertexHypothesis.cxx

@@ -17,14 +17,16 @@
 
 using namespace o2::vertexing;
 
-void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float cpt, float cpt1, float cpt2, float cpt3, float bz, float maxSigma)
+void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float nSigTight, float marginTight, float cpt, float cpt1, float cpt2, float cpt3, float bz, float maxSigma)
 {
   mPIDV0 = v0;
   mPIDPosProng = ppos;
   mPIDNegProng = pneg;
   mPars[SigmaM] = sig;
   mPars[NSigmaM] = nSig;
   mPars[MarginM] = margin;
+  mPars[NSigmaTightM] = nSigTight;
+  mPars[MarginTightM] = marginTight;
   mPars[CPt] = cpt;
   mPars[CPt1] = cpt1;
   mPars[CPt2] = cpt2;
@@ -37,7 +39,7 @@ void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, float sig, float nSig, f
 
 void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz, float maxSigma)
 {
-  set(v0, ppos, pneg, pars[SigmaM], pars[NSigmaM], pars[MarginM], pars[CPt], pars[CPt1], pars[CPt2], pars[CPt3], bz, maxSigma);
+  set(v0, ppos, pneg, pars[SigmaM], pars[NSigmaM], pars[MarginM], pars[NSigmaTightM], pars[MarginTightM], pars[CPt], pars[CPt1], pars[CPt2], pars[CPt3], bz, maxSigma);
 }
 
 void SVertex3Hypothesis::set(PID v0, PID ppos, PID pneg, PID pbach, float sig, float nSig, float margin, float cpt, float bz)

Detectors/Vertexing/src/SVertexer.cxx

@@ -561,11 +561,14 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
   float p2Pos = pP[0] * pP[0] + pP[1] * pP[1] + pP[2] * pP[2], p2Neg = pN[0] * pN[0] + pN[1] * pN[1] + pN[2] * pN[2];
 
   bool goodHyp = false;
-  std::array<bool, NHypV0> hypCheckStatus{};
+  std::array<bool, NHypV0> hypCheckStatus{}, hypCheckTightStatus{};
   for (int ipid = 0; ipid < NHypV0; ipid++) {
     if (mV0Hyps[ipid].check(p2Pos, p2Neg, p2V0, ptV0)) {
       goodHyp = hypCheckStatus[ipid] = true;
     }
+    if (mV0Hyps[ipid].checkTight(p2Pos, p2Neg, p2V0, ptV0)) {
+      hypCheckTightStatus[ipid] = true;
+    }
   }
 
   // apply mass selections for 3-body decay
@@ -581,7 +584,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
   // we want to reconstruct the 3 body decay of hypernuclei starting from the V0 of a proton and a pion (e.g. H3L->d + (p + pi-), or He4L->He3 + (p + pi-)))
   bool checkFor3BodyDecays = mEnable3BodyDecays && (!mSVParams->checkV0Hypothesis || good3bodyV0Hyp) && (pt2V0 > 0.5);
   bool rejectAfter3BodyCheck = false; // To reject v0s which can be 3-body decay candidates but not cascade or v0
-  bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (hypCheckStatus[HypV0::Lambda] || hypCheckStatus[HypV0::AntiLambda]));
+  bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (hypCheckTightStatus[HypV0::Lambda] || hypCheckTightStatus[HypV0::AntiLambda]));
   bool rejectIfNotCascade = false;
 
   if (!goodHyp && mSVParams->checkV0Hypothesis) {
@@ -673,10 +676,10 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
   // check cascades
   int nCascIni = mCascadesIdxTmp[ithread].size(), nV0Used = 0; // number of times this particular v0 (with assigned PV) was used (not counting using its clones with other PV)
   if (checkForCascade) {
-    if (hypCheckStatus[HypV0::Lambda] || !mSVParams->checkCascadeHypothesis) {
+    if (hypCheckTightStatus[HypV0::Lambda] || !mSVParams->checkCascadeHypothesis) {
       nV0Used += checkCascades(v0Idxnew, v0new, rv0, pV0, p2V0, iN, NEG, vlist, ithread);
     }
-    if (hypCheckStatus[HypV0::AntiLambda] || !mSVParams->checkCascadeHypothesis) {
+    if (hypCheckTightStatus[HypV0::AntiLambda] || !mSVParams->checkCascadeHypothesis) {
       nV0Used += checkCascades(v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread);
     }
   }