Use DBSCAN for Time-Z clustering in PVertexing + debris reduction

Author: shahor02

Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h

@@ -80,7 +80,7 @@ class PVertexer
   bool getValidateWithFT0() const { return mValidateWithFT0; }
 
   auto& getTracksPool() const { return mTracksPool; }
-  auto& getTimeClusters() const { return mTimesClusters; }
+  auto& getTimeZClusters() const { return mTimeZClusters; }
   auto& getSortedTrackIndices() const { return mSortedTrackID; }
 
   auto& getMeanVertex() const { return mMeanVertex; }
@@ -92,7 +92,6 @@ class PVertexer
 
   float estimateScale2()
   {
-    float minrange = std::min(mPVParams->zHistoBinSize, mPVParams->minZSeedRange);
     auto sc = mPVParams->zHistoBinSize * mPVParams->zHistoBinSize * mTukey2I / (mStatZErr.getMean() * mStatZErr.getMean());
     return sc;
   }
@@ -108,11 +107,12 @@ class PVertexer
   void applyConstraint(VertexSeed& vtxSeed) const;
   bool upscaleSigma(VertexSeed& vtxSeed) const;
   void createTracksPool(gsl::span<const o2d::TrackTPCITS> tracksITSTPC);
-  void clusterizeTimeBruteForce(float margin = 0.1, float cut = 25);
-  void clusterizeTime(float binSize = 0.1, float maxTDist = 0.6);
   int findVertices(const VertexingInput& input, std::vector<PVertex>& vertices, std::vector<int>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs);
   std::pair<int, int> getBestFT0Trigger(const PVertex& vtx, gsl::span<const o2::ft0::RecPoints> ft0Data, int& currEntry) const;
 
+  int dbscan_RangeQuery(int idxs, std::vector<int>& cand, const std::vector<int>& status);
+  void dbscan_clusterize();
+
   o2::BunchFilling mBunchFilling;
   std::array<int16_t, o2::constants::lhc::LHCMaxBunches> mClosestBunchAbove; // closest filled bunch from above
   std::array<int16_t, o2::constants::lhc::LHCMaxBunches> mClosestBunchBelow; // closest filled bunch from below
@@ -123,7 +123,9 @@ class PVertexer
   o2::math_utils::StatAccumulator mStatTErr;
   std::vector<TrackVF> mTracksPool;        ///< tracks in internal representation used for vertexing
   std::vector<int> mSortedTrackID;         ///< indices of tracks sorted in time
-  std::vector<TimeCluster> mTimesClusters; ///< set of time clusters
+  std::vector<TimeZCluster> mTimeZClusters; ///< set of time clusters
+  std::vector<int> mClusterTrackIDs;        ///< IDs of tracks making the clusters
+
   float mBz = 0.;                          ///< mag.field at beam line
   bool mValidateWithFT0 = false;           ///< require vertex validation with FT0 (if available)
 

Detectors/Vertexing/include/DetectorsVertexing/PVertexerHelpers.h

@@ -116,6 +116,23 @@ struct TrackVF {
     return z + tgL * (vx * cosAlp + vy * sinAlp - x);
   }
 
+  // weighted distance^2 to other track (accounting for own errors only)
+  float getDist2(const TrackVF& o) const
+  {
+    auto dtnorm2 = (timeEst.getTimeStamp() - o.timeEst.getTimeStamp()) / timeEst.getTimeStampError();
+    auto dz = z - o.z;
+    return dtnorm2 + dz * dz * sig2ZI;
+  }
+
+  // weighted distance^2 to other track (accounting for both track errors errors only)
+  float getDist2Sym(const TrackVF& o) const
+  {
+    auto dt = timeEst.getTimeStamp() - o.timeEst.getTimeStamp();
+    auto dz = z - o.z;
+    float dte2 = o.timeEst.getTimeStampError() * o.timeEst.getTimeStampError() + timeEst.getTimeStampError() * timeEst.getTimeStampError();
+    return dt / dte2 + dz * dz / (1. / sig2ZI + 1. / o.sig2ZI);
+  }
+
   float getResiduals(const PVertex& vtx, float& dy, float& dz) const
   {
     // get residuals (Y and Z DCA in track frame) and calculate chi2
@@ -139,7 +156,6 @@ struct TrackVF {
 };
 
 struct VertexingInput {
-  gsl::span<TrackVF> tracks;
   gsl::span<int> idRange;
   TimeEst timeEst{0, -1.}; // negative error means don't use time info
   float scaleSigma2 = 10;
@@ -227,7 +243,7 @@ struct SeedHisto {
   }
 };
 
-struct TimeCluster {
+struct TimeZCluster {
   TimeEst timeEst;
   int first = -1;
   int last = -1;
@@ -275,7 +291,7 @@ struct TimeCluster {
     }
   }
 
-  void merge(TimeCluster& c)
+  void merge(TimeZCluster& c)
   {
     if (c.first < last) {
       first = c.first;

Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h

@@ -25,8 +25,19 @@ namespace vertexing
 struct PVertexerParams : public o2::conf::ConfigurableParamHelper<PVertexerParams> {
   static constexpr float kDefTukey = 5.0f; ///< def.value for tukey constant
 
+  // DBSCAN clustering settings
+  float dbscanMaxDist2 = 9.;   ///< distance^2 cut (eps^2).
+  float dbscanDeltaT = 10.;    ///< abs. time difference cut
+  float dbscanAdaptCoef = 0.1; ///< adapt dbscan minPts for each cluster as minPts=max(minPts, currentSize*dbscanAdaptCoef).
+
+  int maxVerticesPerCluster = 1; ///< max vertices per time-z cluster to look for
+  int maxTrialsPerCluster = 3;   ///< max unsucessful trials for vertex search per vertex
+
+  // track selection
+  float dcaTolerance = 1.3; ///< consider tracks within this abs DCA to mean vertex
+  float pullIniCut = 9;     ///< cut on pull (n^2 sigma) on dca to mean vertex
+
   // parameters
-  float minZSeedRange = 0.5;         ///< min proximity of Zseed which should be covered by 1st iteration
   float zHistoRange = 20.;           ///< +-range of the Zseed histo
   float zHistoBinSize = 0.5;         ///< size of the Zseed histo bin
   float tukey = kDefTukey;           ///< 1./[Tukey parameter]^2