Extra cleanup with t,z MAD estimators

DataFormats/Reconstruction/include/ReconstructionDataFormats/PrimaryVertex.h

@@ -37,6 +37,11 @@ class PrimaryVertex : public Vertex<TimeStampWithError<float, float>>
   void setIR(const InteractionRecord& ir) { mIRMin = mIRMax = ir; }
   bool hasUniqueIR() const { return !mIRMin.isDummy() && (mIRMin == mIRMax); }
 
+  float getTMAD() const { return mTMAD; }
+  void setTMAD(float v) { mTMAD = v; }
+  float getZMAD() const { return mZMAD; }
+  void setZMAD(float v) { mZMAD = v; }
+
 #ifndef GPUCA_ALIGPUCODE
   void print() const;
   std::string asString() const;
@@ -45,8 +50,10 @@ class PrimaryVertex : public Vertex<TimeStampWithError<float, float>>
  protected:
   InteractionRecord mIRMin{}; ///< by default not assigned!
   InteractionRecord mIRMax{}; ///< by default not assigned!
+  float mTMAD = -1;           ///< MAD estimator for Tsigma
+  float mZMAD = -1;           ///< MAD estimator for Zsigma
 
-  ClassDefNV(PrimaryVertex, 1);
+  ClassDefNV(PrimaryVertex, 2);
 };
 
 #ifndef GPUCA_ALIGPUCODE

Detectors/GlobalTrackingWorkflow/src/PrimaryVertexingSpec.cxx

@@ -81,6 +81,7 @@ void PrimaryVertexingSpec::init(InitContext& ic)
     throw std::runtime_error(fmt::format("directory {} for raw data dumps does not exist", dumpDir));
   }
   mVertexer.setPoolDumpDirectory(dumpDir);
+  mVertexer.setTrackSources(mTrackSrc);
 }
 
 void PrimaryVertexingSpec::run(ProcessingContext& pc)
@@ -102,14 +103,19 @@ void PrimaryVertexingSpec::run(ProcessingContext& pc)
     std::vector<o2d::GlobalTrackID> gids;
     auto maxTrackTimeError = PVertexerParams::Instance().maxTimeErrorMUS;
     auto trackMaxX = PVertexerParams::Instance().trackMaxX;
+    auto minIBHits = PVertexerParams::Instance().minIBHits;
     auto halfROFITS = 0.5 * mITSROFrameLengthMUS + mITSROFBiasMUS;
     auto hw2ErrITS = 2.f / std::sqrt(12.f) * mITSROFrameLengthMUS; // conversion from half-width to error for ITS
 
-    auto creator = [maxTrackTimeError, hw2ErrITS, halfROFITS, trackMaxX, &tracks, &gids](auto& _tr, GTrackID _origID, float t0, float terr) {
+    auto creator = [maxTrackTimeError, hw2ErrITS, halfROFITS, trackMaxX, minIBHits, &tracks, &gids, &recoData](auto& _tr, GTrackID _origID, float t0, float terr) {
       if constexpr (isBarrelTrack<decltype(_tr)>()) {
         if (!_origID.includesDet(DetID::ITS) || _tr.getX() > trackMaxX) {
           return true; // just in case this selection was not done on RecoContainer filling level
         }
+        auto itsID = recoData.getITSContributorGID(_origID);
+        if (!itsID.isSourceSet() || o2::math_utils::numberOfBitsSet(recoData.getITSTrack(itsID).getPattern() & 7) < minIBHits) {
+          return true;
+        }
         if constexpr (isITSTrack<decltype(_tr)>()) {
           t0 += halfROFITS;  // ITS time is supplied in \mus as beginning of ROF
           terr *= hw2ErrITS; // error is supplied as a half-ROF duration, convert to \mus
@@ -129,13 +135,17 @@ void PrimaryVertexingSpec::run(ProcessingContext& pc)
       recoData.fillTrackMCLabels(gids, tracksMCInfo);
     }
     mVertexer.setStartIR(recoData.startIR);
-    std::vector<o2::InteractionRecord> ft0Data;
+    static std::vector<InteractionCandidate> ft0Data;
     if (mValidateWithIR) { // select BCs for validation
+      ft0Data.clear();
       const o2::ft0::InteractionTag& ft0Params = o2::ft0::InteractionTag::Instance();
       auto ft0all = recoData.getFT0RecPoints();
       for (const auto& ftRP : ft0all) {
         if (ft0Params.isSelected(ftRP)) {
-          ft0Data.push_back(ftRP.getInteractionRecord());
+          ft0Data.emplace_back(InteractionCandidate{ftRP.getInteractionRecord(),
+                                                    float(ftRP.getInteractionRecord().differenceInBC(recoData.startIR) * o2::constants::lhc::LHCBunchSpacingMUS),
+                                                    float(ftRP.getTrigger().getAmplA() + ftRP.getTrigger().getAmplC()),
+                                                    GTrackID::FT0});
         }
       }
     }
@@ -151,9 +161,9 @@ void PrimaryVertexingSpec::run(ProcessingContext& pc)
   }
 
   mTimer.Stop();
-  LOGP(info, "Found {} PVs, Time CPU/Real:{:.3f}/{:.3f} (DBScan: {:.4f}, Finder:{:.4f}, Rej.Debris:{:.4f}, Reattach:{:.4f}) | {} trials for {} TZ-clusters, max.trials: {}, Slowest TZ-cluster: {} ms of mult {}",
+  LOGP(info, "Found {} PVs, Time CPU/Real:{:.3f}/{:.3f} (DBScan: {:.4f}, Finder:{:.4f}, MADSel:{:.4f}, Rej.Debris:{:.4f}, Reattach:{:.4f}) | {} trials for {} TZ-clusters, max.trials: {}, Slowest TZ-cluster: {} ms of mult {}",
        vertices.size(), mTimer.CpuTime() - timeCPU0, mTimer.RealTime() - timeReal0,
-       mVertexer.getTimeDBScan().CpuTime(), mVertexer.getTimeVertexing().CpuTime(), mVertexer.getTimeDebris().CpuTime(), mVertexer.getTimeReAttach().CpuTime(),
+       mVertexer.getTimeDBScan().CpuTime(), mVertexer.getTimeVertexing().CpuTime(), mVertexer.getTimeMADSel().CpuTime(), mVertexer.getTimeDebris().CpuTime(), mVertexer.getTimeReAttach().CpuTime(),
        mVertexer.getTotTrials(), mVertexer.getNTZClusters(), mVertexer.getMaxTrialsPerCluster(),
        mVertexer.getLongestClusterTimeMS(), mVertexer.getLongestClusterMult());
 }

Detectors/Vertexing/include/DetectorsVertexing/PVertexer.h

@@ -31,6 +31,7 @@
 #include "DetectorsVertexing/PVertexerParams.h"
 #include "ReconstructionDataFormats/GlobalTrackID.h"
 #include "DataFormatsCalibration/MeanVertexObject.h"
+#include "ITSMFTBase/DPLAlpideParam.h"
 #include "gsl/span"
 #include <numeric>
 #include <TTree.h>
@@ -59,7 +60,11 @@ class PVertexer
   void init();
   void end();
   template <typename TR>
-  int process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> bcData,
+  int process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<InteractionCandidate> intCand,
+              std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
+              const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx);
+  template <typename TR>
+  int process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> intRec,
               std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
               const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx);
 
@@ -82,6 +87,7 @@ class PVertexer
   void setBz(float bz) { mBz = bz; }
   void setValidateWithIR(bool v) { mValidateWithIR = v; }
   bool getValidateWithIR() const { return mValidateWithIR; }
+  void setTrackSources(GTrackID::mask_t s);
 
   auto& getTracksPool() const { return mTracksPool; }
   auto& getTimeZClusters() const { return mTimeZClusters; }
@@ -117,6 +123,7 @@ class PVertexer
   TStopwatch& getTimeDBScan() { return mTimeDBScan; }
   TStopwatch& getTimeVertexing() { return mTimeVertexing; }
   TStopwatch& getTimeDebris() { return mTimeDebris; }
+  TStopwatch& getTimeMADSel() { return mTimeMADSel; }
   TStopwatch& getTimeReAttach() { return mTimeReAttach; }
 
   void setPoolDumpDirectory(const std::string& d) { mPoolDumpDirectory = d; }
@@ -127,7 +134,7 @@ class PVertexer
   static constexpr int DBS_UNDEF = -2, DBS_NOISE = -1, DBS_INCHECK = -10;
 
   SeedHistoTZ buildHistoTZ(const VertexingInput& input);
-  int runVertexing(gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> bcData,
+  int runVertexing(gsl::span<o2d::GlobalTrackID> gids, const gsl::span<InteractionCandidate> intCand,
                    std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
                    gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx);
   void createMCLabels(gsl::span<const o2::MCCompLabel> lblTracks, const std::vector<uint32_t>& trackIDs, const std::vector<V2TRef>& v2tRefs, std::vector<o2::MCEventLabel>& lblVtx);
@@ -144,14 +151,17 @@ class PVertexer
   bool upscaleSigma(VertexSeed& vtxSeed) const;
   bool relateTrackToMeanVertex(o2::track::TrackParCov& trc, float vtxErr2);
   bool relateTrackToVertex(o2::track::TrackParCov& trc, const o2d::VertexBase& vtxSeed) const;
+  void applyMADSelection(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const std::vector<V2TRef>& v2tRefs, const std::vector<uint32_t>& trackIDs);
+  void applITSOnlyFractionCut(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const std::vector<V2TRef>& v2tRefs, const std::vector<uint32_t>& trackIDs);
+  void applInteractionValidation(std::vector<PVertex>& vertices, std::vector<int>& timeSort, const gsl::span<InteractionCandidate> intCand, int minContrib);
 
   template <typename TR>
   void createTracksPool(const TR& tracks, gsl::span<const o2d::GlobalTrackID> gids);
 
   int findVertices(const VertexingInput& input, std::vector<PVertex>& vertices, std::vector<uint32_t>& trackIDs, std::vector<V2TRef>& v2tRefs);
   void reAttach(std::vector<PVertex>& vertices, std::vector<int>& timeSort, std::vector<uint32_t>& trackIDs, std::vector<V2TRef>& v2tRefs);
 
-  std::pair<int, int> getBestIR(const PVertex& vtx, const gsl::span<o2::InteractionRecord> bcData, int& currEntry) const;
+  std::pair<int, int> getBestIR(const PVertex& vtx, const gsl::span<InteractionCandidate> intCand, int& currEntry) const;
 
   int dbscan_RangeQuery(int idxs, std::vector<int>& cand, std::vector<int>& status);
   void dbscan_clusterize();
@@ -173,7 +183,7 @@ class PVertexer
   float mBz = 0.;                           ///< mag.field at beam line
   float mDBScanDeltaT = 0.;                 ///< deltaT cut for DBScan check
   float mDBSMaxZ2InvCorePoint = 0;          ///< inverse of max sigZ^2 of the track which can be core point in the DBScan
-  bool mValidateWithIR = false;             ///< require vertex validation with InteractionRecords (if available)
+  bool mValidateWithIR = false;             ///< require vertex validation with InteractionCandidates (if available)
 
   o2::InteractionRecord mStartIR{0, 0}; ///< IR corresponding to the start of the TF
 
@@ -183,6 +193,8 @@ class PVertexer
   //
 
   ///========== Parameters to be set externally, e.g. from CCDB ====================
+  GTrackID::mask_t mTrackSrc{};
+  std::vector<int> mSrcVec{};
   const PVertexerParams* mPVParams = nullptr;
   float mTukey2I = 0;                        ///< 1./[Tukey parameter]^2
   static constexpr float kDefTukey = 5.0f;   ///< def.value for tukey constant
@@ -192,12 +204,18 @@ class PVertexer
   size_t mNTZClustersIni = 0;
   size_t mTotTrials = 0;
   size_t mMaxTrialPerCluster = 0;
+  float mMaxTDiffDebris = 0;      ///< when reducing debris, don't consider vertices separated by time > this value in \mus
+  float mMaxTDiffDebrisExtra = 0; ///< when reducing debris, don't consider vertices separated by time > this value in \mus (optional additional cut
+  float mMaxTDiffDebrisFiducial = 0;
+  float mMaxZDiffDebrisFiducial = 0;
+  float mMaxMultRatDebrisFiducial = 0;
   long mLongestClusterTimeMS = 0;
   int mLongestClusterMult = 0;
   bool mPoolDumpProduced = false;
   TStopwatch mTimeDBScan;
   TStopwatch mTimeVertexing;
   TStopwatch mTimeDebris;
+  TStopwatch mTimeMADSel;
   TStopwatch mTimeReAttach;
   std::string mPoolDumpDirectory{};
 #ifdef _PV_DEBUG_TREE_
@@ -263,6 +281,7 @@ void PVertexer::createTracksPool(const TR& tracks, gsl::span<const o2d::GlobalTr
 
   // check all containers
   float vtxErr2 = 0.5 * (mMeanVertex.getSigmaX2() + mMeanVertex.getSigmaY2()) + mPVParams->meanVertexExtraErrSelection * mPVParams->meanVertexExtraErrSelection;
+  const auto& alpParams = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance();
 
   for (uint32_t i = 0; i < ntGlo; i++) {
     int id = sortedTrackID[i];
@@ -275,6 +294,18 @@ void PVertexer::createTracksPool(const TR& tracks, gsl::span<const o2d::GlobalTr
     auto& tvf = mTracksPool.emplace_back(trc, tracks[id].getTimeMUS(), id, gids[id], mPVParams->addTimeSigma2, mPVParams->addZSigma2);
     if (tvf.wghHisto < 0) {
       mTracksPool.pop_back(); // discard bad track
+      continue;
+    }
+    if (gids[id].getSource() == GTrackID::ITSTPC) { // if the track was adjusted to ITS ROF boundary, flag it
+      float bcf = tvf.timeEst.getTimeStamp() / o2::constants::lhc::LHCBunchSpacingMUS + o2::constants::lhc::LHCMaxBunches;
+      int bcWrtROF = int(bcf - alpParams.roFrameBiasInBC) % alpParams.roFrameLengthInBC;
+      if (bcWrtROF == 0) {
+        float dbc = bcf - (int(bcf / alpParams.roFrameBiasInBC)) * alpParams.roFrameBiasInBC;
+        if (std::abs(dbc) < 1e-6) {
+          tvf.setITSTPCAdjusted();
+          LOGP(debug, "Adjusted t={} -> bcf={} dbc = {}", tvf.timeEst.getTimeStamp(), bcf, dbc);
+        }
+      }
     }
   }
 
@@ -288,12 +319,28 @@ void PVertexer::createTracksPool(const TR& tracks, gsl::span<const o2d::GlobalTr
 
 //___________________________________________________________________
 template <typename TR>
-int PVertexer::process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> bcData,
+int PVertexer::process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<InteractionCandidate> intCand,
+                       std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
+                       const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx)
+{
+  createTracksPool(tracks, gids);
+  return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx);
+}
+
+//___________________________________________________________________
+template <typename TR>
+int PVertexer::process(const TR& tracks, const gsl::span<o2d::GlobalTrackID> gids, const gsl::span<o2::InteractionRecord> intRec,
                        std::vector<PVertex>& vertices, std::vector<o2d::VtxTrackIndex>& vertexTrackIDs, std::vector<V2TRef>& v2tRefs,
                        const gsl::span<const o2::MCCompLabel> lblTracks, std::vector<o2::MCEventLabel>& lblVtx)
 {
   createTracksPool(tracks, gids);
-  return runVertexing(gids, bcData, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx);
+  static std::vector<InteractionCandidate> intCand;
+  intCand.clear();
+  intCand.reserve(intRec.size());
+  for (const auto& ir : intRec) {
+    intCand.emplace_back(InteractionCandidate{ir, float(ir.differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingMUS), 0, 0});
+  }
+  return runVertexing(gids, intCand, vertices, vertexTrackIDs, v2tRefs, lblTracks, lblVtx);
 }
 
 //___________________________________________________________________

Detectors/Vertexing/include/DetectorsVertexing/PVertexerHelpers.h

@@ -105,6 +105,8 @@ struct TrackVF {
   enum { kUsed,
          kNoVtx = -1,
          kDiscarded = kNoVtx - 1 };
+  enum { kITSTPCAdjust = 0x1,
+         kDummyHBin = 0xffff };
   float x;      ///< reference X
   float y;      ///< Y at X
   float z;      ///< Z at X
@@ -119,11 +121,14 @@ struct TrackVF {
   TimeEst timeEst;
   float wgh = 0.; ///< track weight wrt current vertex seed
   float wghHisto = 0.; // weight based on track errors, used for histogramming
-  int entry;      ///< track entry in the input vector
-  int32_t bin = -1; // seeds histo bin
-  GTrackID gid{};
+  int entry;           ///< track entry in the input vector
   int vtxID = kNoVtx; ///< assigned vertex
+  GTrackID gid{};
+  uint16_t bin = kDummyHBin; // seeds histo bin
+  uint16_t flags = 0;
   //
+  void setITSTPCAdjusted() { flags &= kITSTPCAdjust; }
+  bool isITSTPCAdjusted() const { return flags & kITSTPCAdjust; }
   bool canAssign() const { return wgh > 0. && vtxID == kNoVtx; }
   bool canUse() const { return vtxID == kNoVtx; }
   bool canUse(float zmin, float zmax) const
@@ -215,7 +220,7 @@ struct TrackVF {
 struct SeedHistoTZ : public o2::dataformats::FlatHisto2D_f {
   using o2::dataformats::FlatHisto2D<float>::FlatHisto2D;
 
-  int fillAndFlagBin(float x, float y, float w)
+  uint16_t fillAndFlagBin(float x, float y, float w)
   {
     uint32_t bin = getBin(x, y);
     if (isValidBin(bin)) {
@@ -224,9 +229,9 @@ struct SeedHistoTZ : public o2::dataformats::FlatHisto2D_f {
       }
       fillBin(bin, w);
       nEntries++;
-      return bin;
+      return uint16_t(bin);
     }
-    return -1;
+    return 0xffff;
   }
 
   void clear()
@@ -268,6 +273,13 @@ struct TrackVFDump {
   ClassDefNV(TrackVFDump, 1);
 };
 
+struct InteractionCandidate : public o2::InteractionRecord {
+  float time = 0;
+  float amplitude = 0;
+  uint32_t flag = 0; // origin, etc.
+  InteractionCandidate() = default;
+};
+
 } // namespace vertexing
 } // namespace o2
 

Detectors/Vertexing/include/DetectorsVertexing/PVertexerParams.h

@@ -45,6 +45,7 @@ struct PVertexerParams : public o2::conf::ConfigurableParamHelper<PVertexerParam
   float pullIniCut = 9;     ///< cut on pull (n^2 sigma) on dca to mean vertex
   float maxTimeErrorMUS = 10.0; ///< max time error in ms of the track to account
   float trackMaxX = 5.;         ///< lowest updtate point must be below this X
+  int minIBHits = 1.;           ///< min number of IB hits (set to 1 for bwd compatibility, but better use 2)
 
   // histogramming and its weigths params
   float histoBinZSize = 0.05;       ///< size of the seedTZ histo bin Z
@@ -66,15 +67,35 @@ struct PVertexerParams : public o2::conf::ConfigurableParamHelper<PVertexerParam
   // cleanup
   bool applyDebrisReduction = true;        ///< apply algorithm reducing split vertices
   bool applyReattachment = true;           ///< refit vertices reattaching tracks to closest found vertex
+
   float timeMarginReattach = 1.;           ///< safety marging for track time vs vertex time difference during reattachment
-  float maxTDiffDebris = 7.0;              ///< when reducing debris, don't consider vertices separated by time > this value in \mus
+
+  float maxTDiffDebris = 7.0;              ///< when reducing debris, don't consider vertices separated by time > this value in \mus if >0, if <0: mult factor to ITS ROF
   float maxZDiffDebris = 1.0;              ///< don't consider vertices separated by Z > this value in cm
   float maxMultRatDebris = 0.05;           ///< don't consider vertices with multiplicity ratio above this
   float maxChi2TZDebris = 2000.;           ///< don't consider vertices with mutual chi2 exceeding this (for pp should be ~10)
   float addTimeSigma2Debris = 0.05 * 0.05; ///< increment time error^2 by this amount when calculating vertex-to-vertex chi2
   float addZSigma2Debris = 0.005 * 0.005;  ///< increment z error^2 by this amount when calculating vertex-to-vertex chi2
 
+  // extra debris reduction cut, ignored if maxTDiffDebrisExtra == 0. The maxTDiffDebrisExtra must not exceed maxTDiffDebris
+  float maxTDiffDebrisExtra = 0;                ///< when reducing debris, don't consider vertices separated by time > this value in \mus if >0, if <0: mult factor to ITS ROF
+  float maxZDiffDebrisExtra = 1.0;              ///< don't consider vertices separated by Z > this value in cm
+  float maxMultRatDebrisExtra = 0.05;           ///< don't consider vertices with multiplicity ratio above this
+  float maxChi2TZDebrisExtra = 2000.;           ///< don't consider vertices with mutual chi2 exceeding this (for pp should be ~10)
+  float addTimeSigma2DebrisExtra = 0.05 * 0.05; ///< increment time error^2 by this amount when calculating vertex-to-vertex chi2
+  float addZSigma2DebrisExtra = 0.005 * 0.005;  ///< increment z error^2 by this amount when calculating vertex-to-vertex chi2
+
+  float maxITSOnlyFraction = 1.0; ///< max ITS-only tracks fraction to accept, recommended value for PbPb = 0.85
+  float minITSOnlyFraction = 0.0; ///< min ITS-only tracks fraction to accept
+  //
+  // MAD clean up (median abs. deviation in T and Z)
+  float maxTMAD = -1.; ///< max accepted tMAD, not tMAD cleanup if negative
+  float minTMAD = -1.; ///< min accepted tMAD
+  float maxZMAD = -1.; ///< max accepted zMAD, not zMAD cleanup if negative
+  float minZMAD = -1.; ///< min accepted zMAD
+
   // validation with externally provided InteractionRecords (e.g. from FT0)
+  int minNContributorsForIRcutIni = -1; ///< min multiplicity to reject before cleanup if no matching interaction found (if >= 0)
   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)