AliceO2Group · shahor02 · Nov 13, 2020 · Oct 29, 2020 · Nov 7, 2020
@@ -82,7 +82,7 @@ class ClusterPattern
  /// Sets the whole bitmask: the number of rows, the number of columns and the pattern
  void setPattern(const unsigned char bitmask[kExtendedPatternBytes]);
  /// Static: Compute pattern's COG position. Returns the number of fired pixels
-  static int getCOG(int nRow, int nCol, const unsigned char patt[MaxPatternBytes], float& xCOG, float& zCOG);
+  static int getCOG(int rowSpan, int colSpan, const unsigned char patt[MaxPatternBytes], float& xCOG, float& zCOG);
  /// Compute pattern's COG position. Returns the number of fired pixels
  int getCOG(float& xCOG, float& zCOG) const;


@@ -87,11 +87,49 @@ class NoiseMap
    }
    return n;
  }
+  int dumpAboveProbThreshold(float p = 1e-7) const
+  {
+    return dumpAboveThreshold(p * mNumOfStrobes);
+  }
+
+  void applyProbThreshold(float t, long int n)
+  {
+    // Remove from the maps all pixels with the firing probability below the threshold
+    mProbThreshold = t;
+    mNumOfStrobes = n;
+    for (auto& map : mNoisyPixels) {
+      for (auto it = map.begin(); it != map.end();) {
+        float prob = float(it->second) / mNumOfStrobes;
+        if (prob < mProbThreshold) {
+          it = map.erase(it);
+        } else {
+          ++it;
+        }
+      }
+    }
+  }
+  float getProbThreshold() const { return mProbThreshold; }
+  long int getNumOfStrobes() const { return mNumOfStrobes; }
+
+  bool isNoisy(int chip, int row, int col) const
+  {
+    if (chip > mNoisyPixels.size()) {
+      return false;
+    }
+    auto key = row * 1024 + col;
+    const auto keyIt = mNoisyPixels[chip].find(key);
+    if (keyIt != mNoisyPixels[chip].end()) {
+      return true;
+    }
+    return false;
+  }

 private:
  std::vector<std::map<int, int>> mNoisyPixels; ///< Internal noise map representation
+  long int mNumOfStrobes = 0;                   ///< Accumulated number of ALPIDE strobes
+  float mProbThreshold = 0;                     ///< Probability threshold for noisy pixels

-  ClassDefNV(NoiseMap, 1);
+  ClassDefNV(NoiseMap, 2);
 };
 } // namespace itsmft
 } // namespace o2

@@ -157,7 +157,7 @@ class TopologyDictionary
  ///Returns the local position of a compact cluster
  math_utils::Point3D<float> getClusterCoordinates(const CompCluster& cl) const;
  ///Returns the local position of a compact cluster
-  static math_utils::Point3D<float> getClusterCoordinates(const CompCluster& cl, const ClusterPattern& patt);
+  static math_utils::Point3D<float> getClusterCoordinates(const CompCluster& cl, const ClusterPattern& patt, bool isGroup = true);

  friend BuildTopologyDictionary;
  friend LookUp;

@@ -135,12 +135,18 @@ math_utils::Point3D<float> TopologyDictionary::getClusterCoordinates(const CompC
  return locCl;
 }

-math_utils::Point3D<float> TopologyDictionary::getClusterCoordinates(const CompCluster& cl, const ClusterPattern& patt)
+math_utils::Point3D<float> TopologyDictionary::getClusterCoordinates(const CompCluster& cl, const ClusterPattern& patt, bool isGroup)
 {
+  auto refRow = cl.getRow();
+  auto refCol = cl.getCol();
  float xCOG = 0, zCOG = 0;
  patt.getCOG(xCOG, zCOG);
+  if (isGroup) {
+    refRow -= round(xCOG);
+    refCol -= round(zCOG);
+  }
  math_utils::Point3D<float> locCl;
-  o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(cl.getRow() - round(xCOG) + xCOG, cl.getCol() - round(zCOG) + zCOG, locCl);
+  o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(refRow + xCOG, refCol + zCOG, locCl);
  return locCl;
 }


@@ -13,7 +13,7 @@

 #endif

-void MakeNoiseMapFromClusters(std::string input = "o2clus_its.root", std::string output = "noise.root", int threshold = 3)
+void MakeNoiseMapFromClusters(std::string input = "o2clus_its.root", std::string output = "noise.root", bool only1pix = false, float probT = 3e-6)
 {
  TFile in(input.data());
  if (!in.IsOpen()) {
@@ -43,33 +43,67 @@ void MakeNoiseMapFromClusters(std::string input = "o2clus_its.root", std::string

  o2::itsmft::NoiseMap noiseMap(24120);

-  int n1pix = 0;
+  long int nStrobes = 0;
  auto nevents = clusTree->GetEntries();
  for (int n = 0; n < nevents; n++) {
    clusTree->GetEntry(n);
    auto pattIt = patternsPtr->cbegin();
    for (const auto& rof : *rofVec) {
+      nStrobes++;
      auto clustersInFrame = rof.getROFData(*clusters);
      for (const auto& c : clustersInFrame) {
        if (c.getPatternID() != o2::itsmft::CompCluster::InvalidPatternID)
          continue;
+
        o2::itsmft::ClusterPattern patt(pattIt);
-        if (patt.getRowSpan() != 1)
-          continue;
-        if (patt.getColumnSpan() != 1)
-          continue;
+
        auto id = c.getSensorID();
        auto row = c.getRow();
        auto col = c.getCol();
-        noiseMap.increaseNoiseCount(id, row, col);
-        n1pix++;
+        auto colSpan = patt.getColumnSpan();
+        auto rowSpan = patt.getRowSpan();
+
+        if ((rowSpan == 1) && (colSpan == 1)) {
+          noiseMap.increaseNoiseCount(id, row, col);
+          continue;
+        }
+
+        if (only1pix)
+          continue;
+
+        auto nBits = rowSpan * colSpan;
+        int ic = 0, ir = 0;
+        for (unsigned int i = 2; i < patt.getUsedBytes() + 2; i++) {
+          unsigned char tempChar = patt.getByte(i);
+          int s = 128; // 0b10000000
+          while (s > 0) {
+            if ((tempChar & s) != 0) {
+              noiseMap.increaseNoiseCount(id, row + ir, col + ic);
+            }
+            ic++;
+            s >>= 1;
+            if ((ir + 1) * ic == nBits) {
+              break;
+            }
+            if (ic == colSpan) {
+              ic = 0;
+              ir++;
+            }
+          }
+          if ((ir + 1) * ic == nBits) {
+            break;
+          }
+        }
      }
    }
  }

-  int ncalib = noiseMap.dumpAboveThreshold(threshold);
-  std::cout << "Threshold: " << threshold << "  number of pixels: " << ncalib << '\n';
-  std::cout << "Number of 1-pixel clusters: " << n1pix << '\n';
+  noiseMap.applyProbThreshold(probT, nStrobes);
+
+  int fired = probT * nStrobes;
+  int ncalib = noiseMap.dumpAboveThreshold(fired);
+  std::cout << "Probalibity threshold: " << probT
+            << "  number of pixels: " << ncalib << '\n';

  TFile out(output.data(), "new");
  if (!out.IsOpen()) {

@@ -208,7 +208,7 @@ void CheckTopologies(std::string clusfile = "o2clus_its.root",
            auto locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local
            auto locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart());
            locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z()));
-            const auto locC = o2::itsmft::TopologyDictionary::getClusterCoordinates(cluster, pattern);
+            const auto locC = o2::itsmft::TopologyDictionary::getClusterCoordinates(cluster, pattern, false);
            dX = locH.X() - locC.X();
            dZ = locH.Z() - locC.Z();
          } else {

@@ -516,7 +516,7 @@ void CookedTracker::process(gsl::span<const o2::itsmft::CompClusterExt> const& c
      }
    } else {
      o2::itsmft::ClusterPattern patt(pattIt);
-      locXYZ = dict.getClusterCoordinates(comp, patt);
+      locXYZ = dict.getClusterCoordinates(comp, patt, false);
    }
    auto sensorID = comp.getSensorID();
    // Inverse transformation to the local --> tracking

@@ -73,7 +73,7 @@ void ioutils::convertCompactClusters(gsl::span<const itsmft::CompClusterExt> clu
      }
    } else {
      o2::itsmft::ClusterPattern patt(pattIt);
-      locXYZ = dict.getClusterCoordinates(c, patt);
+      locXYZ = dict.getClusterCoordinates(c, patt, false);
    }
    auto& cl3d = output.emplace_back(c.getSensorID(), geom->getMatrixT2L(c.getSensorID()) ^ locXYZ); // local --> tracking
    cl3d.setErrors(sigmaY2, sigmaZ2, sigmaYZ);
@@ -172,7 +172,7 @@ void ioutils::loadEventData(ROframe& event, gsl::span<const itsmft::CompClusterE
      }
    } else {
      o2::itsmft::ClusterPattern patt(pattIt);
-      locXYZ = dict.getClusterCoordinates(c, patt);
+      locXYZ = dict.getClusterCoordinates(c, patt, false);
    }
    auto sensorID = c.getSensorID();
    // Inverse transformation to the local --> tracking
@@ -221,7 +221,7 @@ int ioutils::loadROFrameData(const o2::itsmft::ROFRecord& rof, ROframe& event, g
      }
    } else {
      o2::itsmft::ClusterPattern patt(pattIt);
-      locXYZ = dict.getClusterCoordinates(c, patt);
+      locXYZ = dict.getClusterCoordinates(c, patt, false);
    }
    auto sensorID = c.getSensorID();
    // Inverse transformation to the local --> tracking

@@ -115,7 +115,6 @@ class AlpideCoder
  static bool isEmptyChip(uint8_t b) { return (b & CHIPEMPTY) == CHIPEMPTY; }

  static void setNoisyPixels(const NoiseMap* noise) { mNoisyPixels = noise; }
-  static void setNoiseThreshold(int t) { mNoiseThreshold = t; }

  /// decode alpide data for the next non-empty chip from the buffer
  template <class T, typename CG>
@@ -323,21 +322,14 @@ class AlpideCoder
  //
  void print() const;
  void reset();
-  //
-  template <class T>
-  static int getChipID(T& buffer)
-  {
-    uint8_t id = 0;
-    return (buffer.current(id) && isChipHeaderOrEmpty(id)) ? (id & AlpideCoder::MaskChipID) : -1;
-  }

 private:
  /// Output a non-noisy fired pixel
  static void addHit(ChipPixelData& chipData, short row, short col)
  {
    if (mNoisyPixels) {
      auto chipID = chipData.getChipID();
-      if (mNoisyPixels->getNoiseLevel(chipID, row, col) > mNoiseThreshold) {
+      if (mNoisyPixels->isNoisy(chipID, row, col)) {
        return;
      }
    }
@@ -423,13 +415,12 @@ class AlpideCoder
  //

  static const NoiseMap* mNoisyPixels;
-  static int mNoiseThreshold;

  // cluster map used for the ENCODING only
  std::vector<int> mFirstInRow;     //! entry of 1st pixel of each non-empty row in the mPix2Encode
  std::vector<PixLink> mPix2Encode; //! pool of links: fired pixel + index of the next one in the row
  //
-  ClassDefNV(AlpideCoder, 2);
+  ClassDefNV(AlpideCoder, 3);
 };

 } // namespace itsmft

@@ -18,7 +18,6 @@
 using namespace o2::itsmft;

 const NoiseMap* AlpideCoder::mNoisyPixels = nullptr;
-int AlpideCoder::mNoiseThreshold = 3;

 //_____________________________________
 void AlpideCoder::print() const

@@ -309,10 +309,13 @@ void Clusterer::ClustererThread::streamCluster(const std::vector<PixelData>& pix
  }
  uint16_t pattID = (isHuge || parent->mPattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(rowSpanW, colSpanW, patt);
  if (pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID)) {
-    float xCOG = 0., zCOG = 0.;
-    ClusterPattern::getCOG(rowSpanW, colSpanW, patt, xCOG, zCOG);
-    rowMin += round(xCOG);
-    colMin += round(zCOG);
+    if (pattID != CompCluster::InvalidPatternID) {
+      //For groupped topologies, the reference pixel is the COG pixel
+      float xCOG = 0., zCOG = 0.;
+      ClusterPattern::getCOG(rowSpanW, colSpanW, patt, xCOG, zCOG);
+      rowMin += round(xCOG);
+      colMin += round(zCOG);
+    }
    if (patternsPtr) {
      patternsPtr->emplace_back((unsigned char)rowSpanW);
      patternsPtr->emplace_back((unsigned char)colSpanW);