packages/mvtx/MvtxNoiseMap.h

0001 /**
0002  * @file mvtx/MvtxHitPruner.h
0003  * @author Yasser Corrales Morales
0004  * @date April 2024
0005  * @brief  Definition of the MVTX NoiseMap
0006  */
0007
0008 #ifndef MVTX_NOISEMAP_H
0009 #define MVTX_NOISEMAP_H
0010
0011 #include "Rtypes.h"
0012
0013 #include <iostream>
0014 #include <climits>
0015 #include <cassert>
0016 #include <cmath>
0017 #include <vector>
0018 #include <map>
0019
0020 //class CompClusterExt;
0021
0022 /**
0023  * @brief NoiseMap class for the MVTX
0024  */
0025
0026 class MvtxNoiseMap
0027 {
0028  public:
0029   /// Constructor, initializing values for position, charge and readout frame
0030   MvtxNoiseMap(std::vector<std::map<int, int>>& noise) { mNoisyPixels.swap(noise); }
0031
0032   /// Constructor
0033   MvtxNoiseMap() = default;
0034   /// Constructor
0035   MvtxNoiseMap(int nchips)
0036   {
0037     mNoisyPixels.assign(nchips, std::map<int, int>());
0038   }
0039   /// Destructor
0040   ~MvtxNoiseMap() = default;
0041
0042   /// Get the noise level for this pixels
0043   float getNoiseLevel(int chip, int row, int col) const
0044   {
0045     assert(chip < (int)mNoisyPixels.size());
0046     const auto keyIt = mNoisyPixels[chip].find(getKey(row, col));
0047     if (keyIt != mNoisyPixels[chip].end()) {
0048       return keyIt->second;
0049     }
0050     return 0;
0051   }
0052
0053   void increaseNoiseCount(int chip, int row, int col)
0054   {
0055     assert(chip < (int)mNoisyPixels.size());
0056     mNoisyPixels[chip][getKey(row, col)]++;
0057   }
0058
0059   void increaseNoiseCount(int chip, const std::vector<int>& rowcolKey)
0060   {
0061     assert(chip < (int)mNoisyPixels.size());
0062     auto& ch = mNoisyPixels[chip];
0063     for (const auto k : rowcolKey) {
0064       ch[k]++;
0065     }
0066   }
0067
0068   int dumpAboveThreshold(int t = 3) const
0069   {
0070     int n = 0;
0071     auto chipID = mNoisyPixels.size();
0072     while (chipID--) {
0073       const auto& map = mNoisyPixels[chipID];
0074       for (const auto& pair : map) {
0075         if (pair.second <= t) {
0076           continue;
0077         }
0078         n++;
0079         auto key = pair.first;
0080         auto row = key2Row(key);
0081         auto col = key2Col(key);
0082         std::cout << "chip, row, col, noise: " << chipID << ' ' << row << ' ' << col << ' ' << pair.second << '\n';
0083       }
0084     }
0085     return n;
0086   }
0087   int dumpAboveProbThreshold(float p = 1e-7) const
0088   {
0089     return dumpAboveThreshold(p * mNumOfStrobes);
0090   }
0091
0092   void applyProbThreshold(float t, long int n, float relErr = 0.2f, int minChipID = 0, int maxChipID = 24119)
0093   {
0094     // Remove from the maps all pixels with the firing probability below the threshold
0095     // Apply the cut only for chips between minChipID and maxChipID (included)
0096     if (n < 1) {
0097       std::cerr << "Cannot apply threshold with " << n <<" ROFs scanned" << std::endl;
0098       return;
0099     }
0100     mProbThreshold = t;
0101     mNumOfStrobes = n;
0102     float minFiredForErr = 0.f;
0103     if (relErr > 0) {
0104       minFiredForErr = relErr * relErr - 1. / n;
0105       if (minFiredForErr <= 0.f) {
0106         std::cerr << "Noise threshold " << t << " with relative error " << relErr << " is not reachable with " << n << " ROFs processed, mask all permanently fired pixels" << std::endl;
0107         minFiredForErr = n;
0108       } else {
0109         minFiredForErr = 1. / minFiredForErr;
0110       }
0111     }
0112     int minFired = std::ceil(std::max(t * mNumOfStrobes, minFiredForErr)); // min number of fired pixels exceeding requested threshold
0113     auto req = getMinROFs(t, relErr);
0114     if (n < req) {
0115       mProbThreshold = float(minFired) / n;
0116       std::cerr << "Requested relative error " << relErr << " with prob.threshold " << t << " needs > " << req << " ROFs, " << n << " provided: pixels with noise >" << mProbThreshold << " will be masked" << std::endl;
0117     }
0118
0119     int currChipID = 0;
0120     for (auto& map : mNoisyPixels) {
0121       if (currChipID < minChipID || currChipID > maxChipID) { // chipID range
0122         currChipID++;
0123         continue;
0124       }
0125       for (auto it = map.begin(); it != map.end();) {
0126         if (it->second < minFired) {
0127           it = map.erase(it);
0128         } else {
0129           ++it;
0130         }
0131       }
0132       currChipID++;
0133     }
0134   }
0135   float getProbThreshold() const { return mProbThreshold; }
0136   long int getNumOfStrobes() const { return mNumOfStrobes; }
0137
0138   bool isNoisy(int chip, int row, int col) const
0139   {
0140     assert(chip < (int)mNoisyPixels.size());
0141     return (mNoisyPixels[chip].find(getKey(row, col)) != mNoisyPixels[chip].end());
0142   }
0143
0144   bool isNoisyOrFullyMasked(int chip, int row, int col) const
0145   {
0146     assert(chip < (int)mNoisyPixels.size());
0147     return isNoisy(chip, row, col) || isFullChipMasked(chip);
0148   }
0149
0150   bool isNoisy(int chip) const
0151   {
0152     assert(chip < (int)mNoisyPixels.size());
0153     return !mNoisyPixels[chip].empty();
0154   }
0155
0156   // Methods required by the calibration framework
0157   void print();
0158 //  void fill(const gsl::span<const CompClusterExt> data);
0159 //  void merge(const MvtxNoiseMap* prev) {}
0160
0161   const std::map<int, int>* getChipMap(int chip) const { return chip < (int)mNoisyPixels.size() ? &mNoisyPixels[chip] : nullptr; }
0162
0163   std::map<int, int>& getChip(int chip) { return mNoisyPixels[chip]; }
0164   const std::map<int, int>& getChip(int chip) const { return mNoisyPixels[chip]; }
0165
0166   void maskFullChip(int chip, bool cleanNoisyPixels = false)
0167   {
0168     if (cleanNoisyPixels) {
0169       resetChip(chip);
0170     }
0171     increaseNoiseCount(chip, -1, -1);
0172   }
0173
0174   bool isFullChipMasked(int chip) const
0175   {
0176     return isNoisy(chip, -1, -1);
0177   }
0178
0179   void resetChip(int chip)
0180   {
0181     assert(chip < (int)mNoisyPixels.size());
0182     mNoisyPixels[chip].clear();
0183   }
0184
0185   static long getMinROFs(float t, float relErr)
0186   {
0187     // calculate min number of ROFs needed to reach threshold t with relative error relErr
0188     relErr = relErr >= 0.f ? relErr : 0.1;
0189     t = t >= 0.f ? t : 1e-6;
0190     return std::ceil((1. + 1. / t) / (relErr * relErr));
0191   }
0192
0193   size_t size() const { return mNoisyPixels.size(); }
0194   void setNumOfStrobes(long n) { mNumOfStrobes = n; }
0195   void addStrobes(long n) { mNumOfStrobes += n; }
0196   long getNumberOfStrobes() const { return mNumOfStrobes; }
0197   static int getKey(int row, int col) { return (row << SHIFT) + col; }
0198   static int key2Row(int key) { return key >> SHIFT; }
0199   static int key2Col(int key) { return key & MASK; }
0200
0201  private:
0202   static constexpr int SHIFT = 10, MASK = (0x1 << SHIFT) - 1;
0203   std::vector<std::map<int, int>> mNoisyPixels; ///< Internal noise map representation
0204   long int mNumOfStrobes = 0;                   ///< Accumulated number of ALPIDE strobes
0205   float mProbThreshold = 0;                     ///< Probability threshold for noisy pixels
0206
0207   ClassDefNV(MvtxNoiseMap, 1);
0208 };
0209
0210 #endif /* MVTX_NOISEMAP_H */