svertexer: adaptive mass windows (#12136)

* svertexer: adaptive mass windows

* One more step: Lucia's changes

* Tighter windows when doing decay chains

* Please consider the following formatting changes

* Latest bugfixes

* Please consider the following formatting changes

---------

Co-authored-by: David Dobrigkeit Chinellato <david.dobrigkeit.chinellato.cern.ch>
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: 2 people

Detectors/Vertexing/include/DetectorsVertexing/SVertexHypothesis.h

@@ -30,15 +30,20 @@ class SVertexHypothesis
 
  public:
   using PID = o2::track::PID;
-  enum PIDParams { SigmaM,  // sigma of mass res at 0 pt
-                   NSigmaM, // number of sigmas of mass res
-                   MarginM, // additive safety margin in mass cut
-                   CPt };   // pT dependence of mass resolution parameterized as mSigma*(1+mC1*pt);
+  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 = 4;
+  static constexpr int NPIDParams = 9;
 
-  void set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float cpt, float bz = 0.f);
-  void set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz = 0.f);
+  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); }
   float getMassPosProng() const { return PID::getMass(mPIDPosProng); }
@@ -57,14 +62,44 @@ 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 { return mPars[NSigmaM] * getSigma(pt) + mPars[MarginM]; }
+  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
+      float sigmaV0Cascade = getSigmaV0Cascade(pt);
+      if (sigmaV0Cascade > maxSigma) { // insuring that at low pt one gets reasonable width as the parametrisation function may explode to unphysical values
+        return mPars[idxNsigma] * maxSigma + mPars[idxMargin];
+      } else {
+        return mPars[idxNsigma] * sigmaV0Cascade + mPars[idxMargin];
+      }
+    } else if (mPIDV0 == PID::K0 || mPIDV0 == PID::Lambda) { // case for V0s, AntiLambda is included in PID::Lambda
+      return mPars[idxNsigma] * getSigmaV0Cascade(pt) + mPars[idxMargin];
+    } else {
+      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); }
@@ -76,8 +111,9 @@ class SVertexHypothesis
 
  public: // to be deleted
   std::array<float, NPIDParams> mPars{};
+  float maxSigma;
 
-  ClassDefNV(SVertexHypothesis, 1);
+  ClassDefNV(SVertexHypothesis, 2);
 };
 
 class SVertex3Hypothesis

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};    // Photon
-  float pidCutsK0[SVertexHypothesis::NPIDParams] = {0.003, 20, 0.07, 0.5};        // K0
-  float pidCutsLambda[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.07, 0.5};    // Lambda
-  float pidCutsHTriton[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5};  // HyperTriton
-  float pidCutsHhydrog4[SVertexHypothesis::NPIDParams] = {0.0025, 14, 0.07, 0.5}; // 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., 5.0, 0.0, 2.84798e-03, 9.84206e-04, 3.31951e-03, 2.39438};     // K0
+  float pidCutsLambda[SVertexHypothesis::NPIDParams] = {0., 20, 0., 5.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.001, 20, 0.07, 0.5};    // XiMinus
-  float pidCutsOmegaMinus[SVertexHypothesis::NPIDParams] = {0.001, 20, 0.07, 0.5}; // 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,21 +17,29 @@
 
 using namespace o2::vertexing;
 
-void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, float sig, float nSig, float margin, float cpt, float bz)
+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 maxSigmaInput)
 {
   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;
+  mPars[CPt3] = cpt3;
+  maxSigma = maxSigmaInput;
   float absBz{std::abs(bz)};
-  mPars[CPt] = absBz > 1e-3 ? cpt * 5.0066791 / absBz : 0.; // assume that pT dependent sigma is linear with B
+  if (cpt3 < 1)
+    mPars[CPt] = absBz > 1e-3 ? cpt * 5.0066791 / absBz : 0.; // assume that pT dependent sigma is linear with B; case for HyperTriton and Hyperhydrog4
 }
 
-void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz)
+void SVertexHypothesis::set(PID v0, PID ppos, PID pneg, const float pars[NPIDParams], float bz, float maxSigmaInput)
 {
-  set(v0, ppos, pneg, pars[SigmaM], pars[NSigmaM], pars[MarginM], pars[CPt], bz);
+  set(v0, ppos, pneg, pars[SigmaM], pars[NSigmaM], pars[MarginM], pars[NSigmaTightM], pars[MarginTightM], pars[CPt], pars[CPt1], pars[CPt2], pars[CPt3], bz, maxSigmaInput);
 }
 
 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

@@ -283,8 +283,8 @@ void SVertexer::updateTimeDependentParams()
   mV0Hyps[HypV0::AntiHyperTriton].set(PID::HyperTriton, PID::Pion, PID::Helium3, mSVParams->pidCutsHTriton, bz);
   mV0Hyps[HypV0::Hyperhydrog4].set(PID::Hyperhydrog4, PID::Alpha, PID::Pion, mSVParams->pidCutsHhydrog4, bz);
   mV0Hyps[HypV0::AntiHyperhydrog4].set(PID::Hyperhydrog4, PID::Pion, PID::Alpha, mSVParams->pidCutsHhydrog4, bz);
-  mCascHyps[HypCascade::XiMinus].set(PID::XiMinus, PID::Lambda, PID::Pion, mSVParams->pidCutsXiMinus, bz);
-  mCascHyps[HypCascade::OmegaMinus].set(PID::OmegaMinus, PID::Lambda, PID::Kaon, mSVParams->pidCutsOmegaMinus, bz);
+  mCascHyps[HypCascade::XiMinus].set(PID::XiMinus, PID::Lambda, PID::Pion, mSVParams->pidCutsXiMinus, bz, mSVParams->maximalCascadeWidth);
+  mCascHyps[HypCascade::OmegaMinus].set(PID::OmegaMinus, PID::Lambda, PID::Kaon, mSVParams->pidCutsOmegaMinus, bz, mSVParams->maximalCascadeWidth);
 
   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);
@@ -567,6 +567,12 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
       goodHyp = hypCheckStatus[ipid] = true;
     }
   }
+  // check tight lambda mass only
+  bool goodLamForCascade = false, goodALamForCascade = false;
+  if (mV0Hyps[Lambda].checkTight(p2Pos, p2Neg, p2V0, ptV0))
+    goodLamForCascade = true;
+  if (mV0Hyps[AntiLambda].checkTight(p2Pos, p2Neg, p2V0, ptV0))
+    goodALamForCascade = true;
 
   // apply mass selections for 3-body decay
   bool good3bodyV0Hyp = false;
@@ -581,7 +587,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 || (goodLamForCascade || goodALamForCascade));
   bool rejectIfNotCascade = false;
 
   if (!goodHyp && mSVParams->checkV0Hypothesis) {
@@ -673,10 +679,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 (goodLamForCascade || !mSVParams->checkCascadeHypothesis) {
       nV0Used += checkCascades(v0Idxnew, v0new, rv0, pV0, p2V0, iN, NEG, vlist, ithread);
     }
-    if (hypCheckStatus[HypV0::AntiLambda] || !mSVParams->checkCascadeHypothesis) {
+    if (goodALamForCascade || !mSVParams->checkCascadeHypothesis) {
       nV0Used += checkCascades(v0Idxnew, v0new, rv0, pV0, p2V0, iP, POS, vlist, ithread);
     }
   }