TPC IDC: optimizing outlier filtering

Outlier map:
- 1. check for obvious outlier (IDCs around median +- nSigma)
- 2. check FECs: masking FECs which are out of sync
- 3. check neighbouring outliers (in case not all outlier for larger areas are recognized)

Other changes

IDCCCDBHelper
- add function to get the mean value for IDC0

IDC Drawing
- draw plots in current canvas if one is created
- draw sector informations

IDCFactorize
- parallelize loading of large objects from file
- add function to normalize IDCZero by gain map and/or median per stack
- Moving dump to TTree functionality to CCDBHelper implementation

TPC Mapper
- mapper make function constexpr

Author: matthias-kleiner

Detectors/TPC/base/include/TPCBase/Mapper.h

@@ -405,7 +405,7 @@ class Mapper
   static constexpr unsigned short getPadsInOROC() { return mPadsInOROC; }
   static constexpr unsigned short getPadsInSector() { return mPadsInSector; }
 
-  unsigned short getNumberOfPads(const GEMstack gemStack) const
+  static constexpr unsigned short getNumberOfPads(const GEMstack gemStack)
   {
     switch (gemStack) {
       case IROCgem: {

Detectors/TPC/calibration/include/TPCCalibration/IDCCCDBHelper.h

@@ -39,7 +39,7 @@ class CalDet;
 
 /*
  Usage
- o2::tpc::IDCCCDBHelper<short> helper;
+ o2::tpc::IDCCCDBHelper<unsigned short> helper;
  // setting the IDC members manually
  helper.setIDCDelta(IDCDelta<DataT>* idcDelta);
  helper.setIDCZero(IDCZero* idcZero);
@@ -77,17 +77,8 @@ class IDCCCDBHelper
 
   /// set scaling of IDC0 to 1
   /// \param rejectOutlier do not take outlier into account
-  void setIDCZeroScale(const bool rejectOutlier = true)
-  {
-    mScaleIDC0Aside = this->scaleIDC0(Side::A, rejectOutlier);
-    mScaleIDC0Cside = this->scaleIDC0(Side::C, rejectOutlier);
-    /// check if IDC0 total is not zero, in that case no scalling is applied
-    if (mScaleIDC0Aside == 0.0 || mScaleIDC0Cside == 0.0) {
-      LOGP(error, "Please check the IDC0 total A side {} and C side {}, is zero, therefore no scaling applied!", mScaleIDC0Aside, mScaleIDC0Cside);
-      mScaleIDC0Aside = 1.0;
-      mScaleIDC0Cside = 1.0;
-    }
-  }
+  void setIDCZeroScale(const bool rejectOutlier = true);
+
   /// setting the grouping parameters
   void setGroupingParameter(IDCGroupHelperSector* helperSector, const Side side = Side::A) { mHelperSector[side] = helperSector; }
 
@@ -187,16 +178,18 @@ class IDCCCDBHelper
   TCanvas* drawFourierCoeff(TCanvas* outputCanvas, Side side, int nbins1D, float xMin1D, float xMax1D) const;
 
   /// dumping the loaded IDC0, IDC1 to a tree for debugging
+  /// \param side TPC side of the data which will be dumped
   /// \param outFileName name of the output file
-  void dumpToTree(const char* outFileName = "IDCCCDBTree.root") const;
+  void dumpToTree(const Side side, const char* outFileName = "IDCCCDBTree.root") const;
 
   /// dumping the loaded fourier coefficients to a tree
   /// \param outFileName name of the output file
   void dumpToFourierCoeffToTree(const char* outFileName = "FourierCCDBTree.root") const;
 
   /// dumping the loaded IDC0, IDC1 to a tree for debugging
+  /// \param side TPC side of the data which will be dumped
   /// \param outFileName name of the output file
-  void dumpToTreeIDCDelta(const char* outFileName = "IDCCCDBTreeDeltaIDC.root") const;
+  void dumpToTreeIDCDelta(const Side side, const char* outFileName = "IDCCCDBTreeDeltaIDC.root") const;
 
   /// convert the loaded IDC0 map to a CalDet<float>
   /// \return returns CalDet containing the IDCZero
@@ -206,10 +199,34 @@ class IDCCCDBHelper
   /// \return returns std vector of CalDets containing the IDCDelta
   std::vector<CalDet<float>> getIDCDeltaCalDet() const;
 
+  /// \return returns pointer to pad status map
+  CalDet<PadFlags>* getPadStatusMap() const { return mPadFlagsMap.get(); }
+
+  /// \return returns pointer to pad status map
+  void setPadStatusMap(const CalDet<PadFlags>& outliermap) { mPadFlagsMap = std::make_unique<CalDet<PadFlags>>(outliermap); }
+
   /// scale the stored IDC0 to 1
-  /// \param rejectOutlier do not take outlier into account
   /// \return returns the scaling factor which was used to scale the IDC0
-  float scaleIDC0(const Side side, const bool rejectOutlier = true);
+  /// \param factor to scale the IDC0s with (IDC0/=factor)
+  /// \param side TPC side of the IDCs
+  void scaleIDC0(const float factor, const Side side);
+
+  /// \returns mean of IDC0
+  /// \param side TPC side of the IDCs
+  /// \param idcZero IDCZero object for which to get the mean
+  /// \param outlierMap possible map containing the outliers which will not be taken into account (if nullptr all IDCs are taken into account)
+  static float getMeanIDC0(const Side side, const IDCZero& idcZero, const CalDet<PadFlags>* outlierMap);
+
+  /// get median of IDC0 for given stack and sector
+  /// \param idcZero IDCZero object for which to get the mean
+  /// \param sector TPC sector
+  /// \param stack stack in sector
+  static float getStackMedian(const IDCZero& idczero, const Sector sector, const GEMstack stack);
+
+  /// get medians per stack for given side
+  /// \param idcZero IDCZero object for which to get the mean
+  /// \param side TPC side of the IDCZero
+  static std::array<float, o2::tpc::GEMSTACKSPERSECTOR * o2::tpc::SECTORSPERSIDE> getStackMedian(const IDCZero& idczero, const Side side);
 
  private:
   ///\ to scale IDC0 from its total
@@ -249,7 +266,7 @@ class IDCCCDBHelper
   /// \param urow row of the ungrouped IDCs
   /// \param upad pad number of the ungrouped IDCs
   /// \param integrationInterval integration interval
-  unsigned int getUngroupedIndexGlobal(const unsigned int sector, const unsigned int region, unsigned int urow, unsigned int upad, unsigned int integrationInterval) const;
+  static unsigned int getUngroupedIndexGlobal(const unsigned int sector, const unsigned int region, unsigned int urow, unsigned int upad, unsigned int integrationInterval);
 
   ClassDefNV(IDCCCDBHelper, 3)
 };

Detectors/TPC/calibration/include/TPCCalibration/IDCContainer.h

@@ -449,7 +449,9 @@ enum class PadFlags : unsigned short {
   flagSaturatedPad = 1 << 3, ///< flag for saturated status binary 0100
   flagHighPad = 1 << 4,      ///< flag for pad with extremly high IDC value
   flagLowPad = 1 << 5,       ///< flag for pad with extremly low IDC value
-  flagSkip = 1 << 6          ///< flag for defining a pad which is just ignored during the calculation of I1 and IDCDelta
+  flagSkip = 1 << 6,         ///< flag for defining a pad which is just ignored during the calculation of I1 and IDCDelta
+  flagFEC = 1 << 7,          ///< flag for a whole masked FEC
+  flagNeighbour = 1 << 8     ///< flag if n neighbouring pads are outlier
 };
 
 inline PadFlags operator&(PadFlags a, PadFlags b) { return static_cast<PadFlags>(static_cast<int>(a) & static_cast<int>(b)); }

Detectors/TPC/calibration/include/TPCCalibration/IDCFactorization.h

@@ -75,14 +75,29 @@ class IDCFactorization : public IDCGroupHelperSector
   void factorizeIDCs(const bool norm, const bool calcDeltas);
 
   /// calculate I_0(r,\phi) = <I(r,\phi,t)>_t
+  /// \param norm normalize IDCs to pad area
   void calcIDCZero(const bool norm);
 
   /// fill I_0 values in case of dead pads,FECs etc.
   void fillIDCZeroDeadPads();
 
   /// create status map for pads which are dead or delivering extremly high values (static outliers will be mapped)
+  /// \param debug create debug output
+  void createStatusMapOutlier(const bool debug = false);
+
+  /// 1. createStatusMap, 2. checkFECs, 3. checkNeighbourOutliers
   void createStatusMap();
 
+  /// after the pad-by-pad status map has been created it can be used to mark FECs which deliver wrong values
+  /// \param maxOutliersPerFEC  if the ratio n_outliers_per_FEC / n_pads_per_FEC is larger than this value, the whole FEC is masked
+  void checkFECs(const float maxOutliersPerFEC = 0.7);
+
+  /// check for each pad where there is no outlier found the neighbouring pads for outliers. If more than n neighbours contains outliers this pad is also marked as an outlier pad
+  /// +- one pad in row direction and +- two pads in pad direction are considered
+  /// \param maxIter maximum number of iterations of the procedure
+  /// \param nOutliersNeighbours minimum number of outliers of the neighbouring pads which are required to mask a pad as outlier
+  void checkNeighbourOutliers(const int maxIter = 10, const int nOutliersNeighbours = 8);
+
   /// calculate I_1(t) = <I(r,\phi,t) / I_0(r,\phi)>_{r,\phi}
   void calcIDCOne();
 
@@ -312,13 +327,16 @@ class IDCFactorization : public IDCGroupHelperSector
 
   /// \param integrationIntervals number of integration intervals which will be dumped to the tree (-1: all integration intervalls)
   /// \param outFileName name of the output file
-  void dumpToTree(int integrationIntervals = -1, const char* outFileName = "IDCTree.root") const;
+  void dumpToTree(const Side side, const char* outFileName = "IDCTree.root");
 
   /// dumping the IDC1 to a TTree including the timestamps (the start time stamp in ms should be set with setTimeStamp())
   /// \param integrationTimeOrbits integration time in orbits
   /// \param outFileName name of the output file
   void dumpToTreeIDC1(const float integrationTimeOrbits = 12, const char* outFileName = "IDC1Tree.root") const;
 
+  /// dump IDCDelta to a TTree
+  void dumpIDCDeltaToTree(const Side side, const int chunk = 0, const char* outFileName = "IDCDeltaTree.root");
+
   /// \returns vector containing the number of integration intervals for each stored TF (dropped TFs not taken into account)
   /// \param cru cru which is used for the lookup (cru=-1: automatic cru lookup)
   std::vector<unsigned int> getIntegrationIntervalsPerTF(const int cru = -1) const;
@@ -371,6 +389,21 @@ class IDCFactorization : public IDCGroupHelperSector
   /// set the IDCOne
   void setIDCOne(const Side side, const IDCOne& idcOne) { mIDCOne[mSideIndex[side]] = idcOne; }
 
+  /// set the IDCDelta
+  void setIDCDelta(const Side side, const IDCDelta<float>& idcDelta, const int index = 0) { mIDCDelta[mSideIndex[side]][index] = idcDelta; }
+
+  /// \return returns mean of IDC0 for given side
+  float getMeanZ(const Side side) const;
+
+  /// normalize IDC0 with the set gain map
+  void normIDCZeroGain() { normIDCZero(0); }
+
+  /// normalize IDC0 with the median per stack
+  void normIDCZeroStackMedian() { normIDCZero(1); }
+
+  /// \return returns median of IDC0 per stack for all stacks
+  std::array<float, o2::tpc::GEMSTACKS> getStackMedian() const;
+
  private:
   const unsigned int mTimeFrames{};                                 ///< number of timeframes which are stored
   const unsigned int mTimeFramesDeltaIDC{};                         ///< number of timeframes of which Delta IDCs are stored
@@ -400,7 +433,7 @@ class IDCFactorization : public IDCGroupHelperSector
   void drawIDCZeroHelper(const bool type, const Sector sector, const std::string filename, const float minZ, const float maxZ) const;
 
   /// get time frame and index of integrationInterval in the TF
-  void getTF(const unsigned int region, unsigned int integrationInterval, unsigned int& timeFrame, unsigned int& interval) const;
+  void getTF(unsigned int integrationInterval, unsigned int& timeFrame, unsigned int& interval) const;
 
   /// \returns chunk from timeframe
   unsigned int getChunk(const unsigned int timeframe) const;
@@ -411,6 +444,9 @@ class IDCFactorization : public IDCGroupHelperSector
   /// helper function for drawing
   void drawPadFlagMap(const bool type, const Sector sector, const std::string filename, const PadFlags flag) const;
 
+  /// normalize IDC0
+  void normIDCZero(const int type);
+
   ClassDefNV(IDCFactorization, 2)
 };
 

Detectors/TPC/calibration/include/TPCCalibration/IDCGroupingParameter.h

@@ -45,10 +45,9 @@ struct ParameterIDCGroup : public o2::conf::ConfigurableParamHelper<ParameterIDC
   unsigned int groupPadsSectorEdges{0};                                                 ///< decoded number of pads at the sector edges which are grouped differently. First digit specifies the EdgePadGroupingMethod  (example: 0: no pads are grouped, 110: first two pads are not grouped, 3210: first pad is not grouped, second + third pads are grouped, fourth + fifth + sixth pads are grouped)
   AveragingMethod method = AveragingMethod::FAST;                                       ///< method which is used for averaging
   float sigma = 3.f;                                                                    ///< sigma cut which can be used during the grouping for outlier filtering
-  float minIDC0Median = 6;                                                              ///< this value is used for identifying outliers (pads with high IDC0 values): "accepted IDC 0 values > median_IDC0 + stdDev * minIDC0Median"
-  float maxIDC0Median = 6;                                                              ///< this value is used for identifying outliers (pads with high IDC0 values): "accepted IDC 0 values < median_IDC0 + stdDev * maxIDC0Median"
-  int minIDC0Val = 1;                                                                   ///< minimum accepted IDC0 value (should be larger than 0 in case of real data)
-  int maxIDC0Val = 150;                                                                 ///< maximum accepted IDC0 value (max IDC0 for 3MHz pp ~30)
+  float minIDC0Median = 8;                                                              ///< this value is used for identifying outliers (pads with low IDC0 values): "accepted IDC 0 values > median_IDC0 + stdDev * minIDC0Median"
+  float maxIDC0Median = 15;                                                             ///< this value is used for identifying outliers (pads with high IDC0 values): "accepted IDC 0 values < median_IDC0 + stdDev * maxIDC0Median"
+  float minmaxIDC0MedianCentreEdge = 20;                                                ///< this value is used for identifying outliers at the centre of the cross and the edges of the sector i.e. where the gain is significantly lower (cross) and sometimes higher at the edges"
 
   /// Helper function for setting the groupimg parameters from a string (can be "X": parameters in all regions are "X" or can be "1,2,3,4,5,6,7,8,9,10" for setting individual regions)
   /// \param sgroupPads string for grouping parameter in pad direction

Detectors/TPC/calibration/include/TPCCalibration/RobustAverage.h

@@ -71,7 +71,8 @@ class RobustAverage
 
   /// returns the filtered average value
   /// \param sigma maximum accepted standard deviation: sigma*stdev
-  float getFilteredAverage(const float sigma = 3);
+  ///\param interQuartileRange number of points in inner quartile to consider
+  std::pair<float, float> getFilteredAverage(const float sigma = 3, const float interQuartileRange = 0.9);
 
   /// \return returns mean of stored values
   float getMean() const { return mValues.empty() ? 0 : getMean(mValues.begin(), mValues.end()); }
@@ -83,7 +84,10 @@ class RobustAverage
   float getWeightedMean() const { return getWeightedMean(mValues.begin(), mValues.end(), mWeights.begin(), mWeights.end()); }
 
   /// \return returns standard deviation of stored values
-  float getStdDev() { return getStdDev(getMean()); }
+  float getStdDev() { return getStdDev(getMean(), mValues.begin(), mValues.end()); }
+
+  /// \return returns stored values
+  const auto& getValues() { return mValues; }
 
   /// values which will be averaged and filtered
   void print() const;
@@ -99,13 +103,13 @@ class RobustAverage
 
   /// \return returns standard deviation of stored values
   /// \param mean mean of stored values
-  float getStdDev(const float mean);
+  float getStdDev(const float mean, std::vector<float>::const_iterator begin, std::vector<float>::const_iterator end);
 
   /// performing outlier filtering of the stored values by defining range of included values in terms of standard deviation
   /// \param mean mean of the stored values
   /// \param stdev standard deviation of the values
   /// \param sigma maximum accepted standard deviation: sigma*stdev
-  float getFilteredMean(const float mean, const float stdev, const float sigma);
+  float getFilteredMean(const float mean, const float stdev, const float sigma) const;
 };
 
 } // namespace o2::tpc

Detectors/TPC/calibration/src/IDCAverageGroupHelper.cxx

@@ -20,7 +20,7 @@ float o2::tpc::IDCAverageGroupHelper<o2::tpc::IDCAverageGroupCRU>::getGroupedIDC
   switch (paramIDCGroup.method) {
     case o2::tpc::AveragingMethod::SLOW:
     default:
-      return mRobustAverage[mThreadNum].getFilteredAverage(paramIDCGroup.sigma);
+      return mRobustAverage[mThreadNum].getFilteredAverage(paramIDCGroup.sigma).first;
       break;
     case o2::tpc::AveragingMethod::FAST:
       return withweights ? mRobustAverage[mThreadNum].getWeightedMean() : mRobustAverage[mThreadNum].getMean();
@@ -80,7 +80,7 @@ float o2::tpc::IDCAverageGroupHelper<o2::tpc::IDCAverageGroupTPC>::getGroupedIDC
   switch (paramIDCGroup.method) {
     case o2::tpc::AveragingMethod::SLOW:
     default:
-      return mRobustAverage[mThreadNum].getFilteredAverage(paramIDCGroup.sigma);
+      return mRobustAverage[mThreadNum].getFilteredAverage(paramIDCGroup.sigma).first;
       break;
     case o2::tpc::AveragingMethod::FAST:
       return withweights ? mRobustAverage[mThreadNum].getWeightedMean() : mRobustAverage[mThreadNum].getMean();