Extra tolerances for P/S-vertexing

Author: shahor02

Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h

@@ -70,7 +70,7 @@ struct PVertexerParams : public o2::conf::ConfigurableParamHelper<PVertexerParam
   int minNContributorsForIRcut = 4;     ///< do not apply IR cut to vertices below IR tagging efficiency threshold
   float maxTError = 0.2;                ///< use min of vertex time error or this for nsigma evaluation
   float minTError = 0.003;              ///< don't use error smaller than that (~BC/2/minNContributorsForFT0cut)
-  float nSigmaTimeCut = 4.;             ///< eliminate vertex if there is no FT0 signal within this cut
+  float nSigmaTimeCut = 4.;             ///< eliminate vertex if there is no FT0 or BC signal within this cut
   float timeBiasMS = 0;                 ///< relative bias in ms to add to TPCITS-based time stamp
   //
   // stopping condition params
@@ -81,6 +81,11 @@ struct PVertexerParams : public o2::conf::ConfigurableParamHelper<PVertexerParam
   int maxNScaleSlowConvergence = 3; ///< max number of weak scaling decrease iterations
   bool useTimeInChi2 = true;        ///< use track-vertex time difference in chi2 calculation
 
+  // track vertex time-wise association
+  float nSigmaTimeTrack = 4.;       ///< define track time bracker as +- this number of sigmas of its time resolution
+  float timeMarginTrackTime = 0.5;  ///< additive marginal error in \mus to track time bracket
+  float timeMarginVertexTime = 0.0; ///< additive marginal error to \mus vertex time bracket
+
   O2ParamDef(PVertexerParams, "pvertexer");
 };
 

Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h

@@ -63,12 +63,14 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper<SVertexerParam
   float maxTglCasc = 2.;    // maximum tgLambda of cascade
 
   // 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
   //
   // 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
 

Detectors/Vertexing/include/DetectorsVertexing/VertexTrackMatcher.h

@@ -59,7 +59,6 @@ class VertexTrackMatcher
  private:
   void updateTimeDependentParams();
   void extractTracks(const o2::globaltracking::RecoContainer& data, const std::unordered_map<GIndex, bool>& vcont);
-
   std::vector<TrackTBracket> mTBrackets;
   float mITSROFrameLengthMUS = 0;       ///< ITS RO frame in mus
   float mMFTROFrameLengthMUS = 0;       ///< MFT RO frame in mus

Detectors/Vertexing/src/PVertexer.cxx

@@ -765,7 +765,7 @@ bool PVertexer::setCompatibleIR(PVertex& vtx)
   // assign compatible IRs accounting for the bunch filling scheme
   const auto& vtxT = vtx.getTimeStamp();
   o2::InteractionRecord irMin(mStartIR), irMax(mStartIR);
-  auto rangeT = mPVParams->nSigmaTimeCut * std::max(mPVParams->minTError, std::min(mPVParams->maxTError, vtxT.getTimeStampError()));
+  auto rangeT = std::max(mPVParams->minTError, mPVParams->nSigmaTimeCut * std::min(mPVParams->maxTError, vtxT.getTimeStampError())) + mPVParams->timeMarginVertexTime;
   float t = vtxT.getTimeStamp() + mPVParams->timeBiasMS;
   if (t > rangeT) {
     irMin += o2::InteractionRecord(1.e3 * (t - rangeT));
@@ -777,11 +777,13 @@ bool PVertexer::setCompatibleIR(PVertex& vtx)
   irMax++; // to account for rounding
   // restrict using bunch filling
   int bc = mClosestBunchAbove[irMin.bc];
+  LOG(DEBUG) << "irMin.bc = " << irMin.bc << " bcAbove = " << bc;
   if (bc < irMin.bc) {
     irMin.orbit++;
   }
   irMin.bc = bc;
   bc = mClosestBunchBelow[irMax.bc];
+  LOG(DEBUG) << "irMax.bc = " << irMax.bc << " bcBelow = " << bc;
   if (bc > irMax.bc) {
     if (irMax.orbit == 0) {
       return false;
@@ -791,6 +793,9 @@ bool PVertexer::setCompatibleIR(PVertex& vtx)
   irMax.bc = bc;
   vtx.setIRMin(irMin);
   vtx.setIRMax(irMax);
+  if (irMin > irMax) {
+    LOG(DEBUG) << "Reject VTX " << vtx.asString() << " trange = " << rangeT;
+  }
   return irMax >= irMin;
 }
 

Detectors/Vertexing/src/SVertexer.cxx

@@ -49,11 +49,13 @@ void SVertexer::process(const o2::globaltracking::RecoContainer& recoData) // ac
     auto& seedP = mTracksPool[POS][itp];
     int firstN = mVtxFirstTrack[NEG][seedP.vBracket.getMin()];
     if (firstN < 0) {
+      LOG(DEBUG) << "No partner is found for pos.track " << itp << " out of " << ntrP;
       continue;
     }
     for (int itn = firstN; itn < ntrN; itn++) { // start from the 1st negative track of lowest-ID vertex of positive
       auto& seedN = mTracksPool[NEG][itn];
       if (seedN.vBracket > seedP.vBracket) { // all vertices compatible with seedN are in future wrt that of seedP
+        LOG(DEBUG) << "Brackets do not match";
         break;
       }
 #ifdef WITH_OPENMP
@@ -231,6 +233,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
   float rv0 = std::sqrt(r2v0), drv0P = rv0 - seedP.minR, drv0N = rv0 - seedN.minR;
   if (drv0P > mSVParams->causalityRTolerance || drv0P < -mSVParams->maxV0ToProngsRDiff ||
       drv0N > mSVParams->causalityRTolerance || drv0N < -mSVParams->maxV0ToProngsRDiff) {
+    LOG(DEBUG) << "RejCausality " << drv0P << " " << drv0N;
     return false;
   }
 
@@ -250,9 +253,11 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
   std::array<float, 3> pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]};
   float pt2V0 = pV0[0] * pV0[0] + pV0[1] * pV0[1], prodXYv0 = dxv0 * pV0[0] + dyv0 * pV0[1], tDCAXY = prodXYv0 / pt2V0;
   if (pt2V0 < mMinPt2V0) { // pt cut
+    LOG(DEBUG) << "RejPt2 " << pt2V0;
     return false;
   }
   if (pV0[2] * pV0[2] / pt2V0 > mMaxTgl2V0) { // tgLambda cut
+    LOG(DEBUG) << "RejTgL " << pV0[2] * pV0[2] / pt2V0;
     return false;
   }
   float p2V0 = pt2V0 + pV0[2] * pV0[2], ptV0 = std::sqrt(pt2V0);
@@ -266,17 +271,19 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
       goodHyp = hypCheckStatus[ipid] = true;
     }
   }
-  if (!goodHyp) {
+  if (!goodHyp && mSVParams->checkV0Hypothesis) {
+    LOG(DEBUG) << "RejHypo";
     return false;
   }
 
-  bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (hypCheckStatus[HypV0::Lambda] || hypCheckStatus[HypV0::AntiLambda]);
+  bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (hypCheckStatus[HypV0::Lambda] || hypCheckStatus[HypV0::AntiLambda]));
   bool rejectIfNotCascade = false;
   float dcaX = dxv0 - pV0[0] * tDCAXY, dcaY = dyv0 - pV0[1] * tDCAXY, dca2 = dcaX * dcaX + dcaY * dcaY;
   float cosPAXY = prodXYv0 / std::sqrt(r2v0 * pt2V0);
 
   if (checkForCascade) { // use loser cuts for cascade v0 candidates
     if (dca2 > mMaxDCAXY2ToMeanVertexV0Casc || cosPAXY < mSVParams->minCosPAXYMeanVertexCascV0) {
+      LOG(DEBUG) << "Rej for cascade DCAXY2: " << dca2 << " << cosPAXY: " << cosPAXY;
       return false;
     }
   }
@@ -298,6 +305,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
     float dx = v0XYZ[0] - pv.getX(), dy = v0XYZ[1] - pv.getY(), dz = v0XYZ[2] - pv.getZ(), prodXYZv0 = dx * pV0[0] + dy * pV0[1] + dz * pV0[2];
     float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0);
     if (cosPA < bestCosPA) {
+      LOG(DEBUG) << "Rej. cosPA: " << cosPA;
       continue;
     }
     if (!added) {
@@ -318,10 +326,10 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
   // check cascades
   if (checkForCascade) {
     int nCascAdded = 0;
-    if (hypCheckStatus[HypV0::Lambda]) {
+    if (hypCheckStatus[HypV0::Lambda] || !mSVParams->checkCascadeHypothesis) {
       nCascAdded += checkCascades(rv0, pV0, p2V0, iN, NEG, ithread);
     }
-    if (hypCheckStatus[HypV0::AntiLambda]) {
+    if (hypCheckStatus[HypV0::AntiLambda] || !mSVParams->checkCascadeHypothesis) {
       nCascAdded += checkCascades(rv0, pV0, p2V0, iP, POS, ithread);
     }
     if (!nCascAdded && rejectIfNotCascade) { // v0 would be accepted only if it creates a cascade

Detectors/Vertexing/src/VertexTrackMatcher.cxx

@@ -86,8 +86,8 @@ void VertexTrackMatcher::process(const o2::globaltracking::RecoContainer& recoDa
     }
     const auto& vtx = vertices[iv];
     const auto& vto = vtxOrdBrack.emplace_back(VtxTBracket{
-      {float((vtx.getIRMin().differenceInBC(recoData.startIR) - 0.5f) * o2::constants::lhc::LHCBunchSpacingMUS),
-       float((vtx.getIRMax().differenceInBC(recoData.startIR) + 0.5f) * o2::constants::lhc::LHCBunchSpacingMUS)},
+      {float((vtx.getIRMin().differenceInBC(recoData.startIR) - 0.5f) * o2::constants::lhc::LHCBunchSpacingMUS - mPVParams->timeMarginVertexTime),
+       float((vtx.getIRMax().differenceInBC(recoData.startIR) + 0.5f) * o2::constants::lhc::LHCBunchSpacingMUS + mPVParams->timeMarginVertexTime)},
       iv});
     if (vto.tBracket.delta() > maxVtxSpan) {
       maxVtxSpan = vto.tBracket.delta();
@@ -196,9 +196,12 @@ void VertexTrackMatcher::extractTracks(const o2::globaltracking::RecoContainer&
       terr *= this->mMFTROFrameLengthMUS;
     } else if constexpr (isGlobalFwdTrack<decltype(_tr)>()) {
       t0 = _tr.getTimeMUS().getTimeStamp();
-      terr = _tr.getTimeMUS().getTimeStampError();
+      terr = _tr.getTimeMUS().getTimeStampError() * mPVParams->nSigmaTimeTrack; // gaussian errors must be scaled by requested n-sigma
+    } else {
+      terr *= mPVParams->nSigmaTimeTrack; // gaussian errors must be scaled by requested n-sigma
     }
     // for all other tracks the time is in \mus with gaussian error
+    terr += mPVParams->timeMarginTrackTime;
     mTBrackets.emplace_back(TrackTBracket{{t0 - terr, t0 + terr}, _origID});
     return true;
   };