TPC: Multiple unrelated updates

Detectors/TPC/base/CMakeLists.txt

@@ -38,6 +38,7 @@ o2_add_library(TPCBase
                        src/IonTailSettings.cxx
                        src/FEEConfig.cxx
                        src/DeadChannelMapCreator.cxx
+                       src/CommonModeCorrection.cxx
                PUBLIC_LINK_LIBRARIES Vc::Vc Boost::boost O2::DataFormatsTPC
                                      O2::DetectorsRaw O2::CCDB FairRoot::Base)
 
@@ -70,6 +71,7 @@ o2_target_root_dictionary(TPCBase
                                   include/TPCBase/IonTailSettings.h
                                   include/TPCBase/FEEConfig.h
                                   include/TPCBase/DeadChannelMapCreator.h
+                                  include/TPCBase/CommonModeCorrection.h
                                   include/TPCBase/CDBTypes.h)
 o2_add_test(Base
             COMPONENT_NAME tpc

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

@@ -0,0 +1,246 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file CommonModeCorrection.h
+/// \brief  Calculate the common mode correction factor
+/// \author Jens Wiechula, Jens.Wiechula@ikf.uni-frankfurt.de
+
+#ifndef AliceO2_TPC_CommonModeCorrection_H_
+#define AliceO2_TPC_CommonModeCorrection_H_
+
+#include <gsl/span>
+#include <string_view>
+#include <vector>
+
+#include "DataFormatsTPC/Digit.h"
+#include "TPCBase/FEEConfig.h"
+
+namespace o2::tpc
+{
+
+/// Class to calculate the common mode correction
+///
+/// Calculation of the common mode correction, based on the algorithm propsed by Marian Ivanov
+/// The calculation is done for one single CRU and time bin
+class CommonModeCorrection
+{
+ public:
+  struct CMdata {
+    std::vector<float> adcValues;
+    std::vector<float> cmKValues;
+    std::vector<float> pedestals;
+
+    void resize(size_t newSize)
+    {
+      adcValues.resize(newSize);
+      cmKValues.resize(newSize);
+      pedestals.resize(newSize);
+    }
+
+    void clear()
+    {
+      adcValues.clear();
+      cmKValues.clear();
+      pedestals.clear();
+    }
+  };
+
+  struct CMInfo {
+    float cmValue{};        ///< common mode value from pseudo code
+    float cmValueStd{};     ///< std dev of common mode values from pseudo code
+    float cmValueCRU{};     ///< common mode value from firmware, if available
+    float sumPos{};         ///< sum of positive signals > mSumPosThreshold
+    float sumNeg{};         ///< sum of negative signals <= mSumPosThreshold, corrected for k-factor
+    uint16_t nPadsUsed{};   ///< number of pads used for CM calculation
+    uint16_t nNeg{};        ///< number of pads used for sumNeg
+    uint16_t nOccupancy{};  ///< number of CM corrected pads larger than mOccupancyThreshold
+    uint16_t nSaturation{}; ///< number of pads in saturation
+  };
+
+  struct CMDebug {
+    std::vector<uint8_t> nPadsOk{};
+    std::vector<uint16_t> adcDist{};
+  };
+
+  using CalPadMapType = std::unordered_map<std::string, CalPad>;
+
+  /// Calculation of the common mode value
+  ///
+  /// \param value pad-by-pad charge values
+  /// \param cmKValues corresponding pad-by-pad common mode k-factors
+  /// \param pedestals corresponding pad-by-pad pedestals
+  /// \param
+  CMInfo getCommonMode(gsl::span<const float> values, gsl::span<const float> cmKValues, gsl::span<const float> pedestals, CMDebug* cmDebug = nullptr) const;
+  CMInfo getCommonMode(const std::vector<float>& values, const std::vector<float>& cmKValues, const std::vector<float>& pedestals) const { return getCommonMode(gsl::span(values), gsl::span(cmKValues), gsl::span(pedestals)); }
+
+  CMInfo getCommonMode(const CMdata& cmData) const { return getCommonMode(std::span(cmData.adcValues), std::span(cmData.cmKValues), std::span(cmData.pedestals)); }
+
+  void setNPadsCompRandom(int n) { mNPadsCompRamdom = n; }
+  int getNPadsCompRandom() const { return mNPadsCompRamdom; }
+
+  void setNPadsCompMin(int n) { mNPadsCompMin = n; }
+  int getNPadsCompMin() const { return mNPadsCompMin; }
+
+  /// Minimum number of pads required in the CM calculation to be used for digit correction
+  void setNPadsMinCM(int n) { mNPadsMinCM = n; }
+  int getNPadsMinCM() const { return mNPadsMinCM; }
+
+  void setQEmpty(float q) { mQEmpty = q; }
+  float getQEmpty() const { return mQEmpty; }
+
+  void setQComp(float q) { mQComp = q; }
+  float getQComp() const { return mQComp; }
+
+  /// The mQComp will be set to (cm - mQCompScaleThreshold) * mQCompScale, if cm > mQCompScaleThreshold
+  void setQCompScaleThreshold(float q) { mQCompScaleThreshold = q; }
+  float getQCompScaleThreshold() const { return mQCompScaleThreshold; }
+
+  /// The mQComp will be set to (cm - mQCompScaleThreshold) * mQCompScale, if cm > mQCompScaleThreshold
+  void setQCompScale(float q) { mQCompScale = q; }
+  float getQCompScale() const { return mQCompScale; }
+
+  /// Threshold above which a signal is considered for sumPos, if debug information is used
+  void setSumPosThreshold(float threshold) { mSumPosThreshold = threshold; }
+  float getSumPosThreshold() const { return mSumPosThreshold; }
+
+  /// Threshold above which a signal is considered for the occupancy
+  void setOccupancyThreshold(float threshold) { mOccupancyThreshold = threshold; }
+  float getOccupancyThreshold() const { return mOccupancyThreshold; }
+
+  /// Pad maps loaded from FEEConfig
+  void setPadMaps(CalPadMapType& padMaps) { mPadMaps = padMaps; }
+
+  /// load a CalPad from file and add it to the local mPadMaps
+  /// \param fileName input file name
+  /// \param nameInFile name of the CalPad object in the file
+  /// \param namePadMap name under which to store the object in the mPadMaps, if empty use the same as nameInFile
+  void loadCalPad(std::string_view fileName, std::string_view nameInFile, std::string_view namePadMap = "");
+
+  /// load CMkValues from file, assuming it is stored under the name "CMkValues
+  void loadCMkValues(std::string_view fileName) { loadCalPad(fileName, "CMkValues"); }
+
+  /// load Pedestals from file, assuming it is stored under the name "Pedestals
+  void loadPedestals(std::string_view fileName) { loadCalPad(fileName, "Pedestals"); }
+
+  /// Custom setting of CalPad, overwriting what was set in mPadMaps
+  void setCalPad(const CalPad& calPad, std::string_view name) { mPadMaps[name.data()] = calPad; }
+
+  /// cmk value
+  float getCMkValue(int sector, int row, int pad) { return mPadMaps["CMkValues"].getValue(sector, row, pad); }
+
+  /// pedestal value
+  float getPedestalValue(int sector, int row, int pad) { return mPadMaps["Pedestals"].getValue(sector, row, pad); }
+
+  /// load the Pad maps from CCDB
+  void
+    loadDefaultPadMaps(FEEConfig::Tags feeTag = FEEConfig::Tags::Physics30sigma);
+
+  CMdata collectCMdata(const std::vector<Digit>& digits, int cru, int timeBin);
+
+  int getCommonMode(std::vector<Digit>& digits, std::vector<std::vector<CMInfo>>& cmValues, bool negativeOnly = false, bool hasInjectedCMValue = false, std::vector<std::vector<CMDebug>>* cmDebug = nullptr, int minTimeBin = -1, int maxTimeBin = -1) const;
+
+  /// corret digits for common mode
+  /// \param cmValues will contain CM information for each CRU and time bin
+  /// \param negativeOnly only correct negative common mode signals
+  /// \return maximum
+  int correctDigits(std::vector<Digit>& digits, std::vector<std::vector<CMInfo>>& cmValues, bool negativeOnly = false, bool hasInjectedCMValue = false, std::vector<std::vector<CMDebug>>* cmDebug = nullptr, int minTimeBin = -1, int maxTimeBin = -1) const;
+
+  void correctDigits(std::string_view digiFileIn, Long64_t maxEntries = -1, std::string_view digitFileOut = "tpcdigit_cmcorr.root", std::string_view cmFileOut = "CommonModeValues.root", bool negativeOnly = false, int nThreads = 1, bool writeOnlyCM = false, bool writeDebug = false, bool hasInjectedCMValue = false, int minTimeBin = -1, int maxTimeBin = -1);
+
+  void limitKFactorPrecision(bool limit = true) { mLimitKFactor = limit; }
+  void limitPedestalPrecision(bool limit = true) { mLimitPedestal = limit; }
+
+  /// set the number of threads used for CM calculation
+  /// \param nThreads number of threads
+  static void setNThreads(const int nThreads) { sNThreads = nThreads; }
+
+  /// \return returns the number of threads used for decoding
+  static int getNThreads() { return sNThreads; }
+
+  /// add artificial common mode, only works when using the 'correctDigits' function
+  void addCommonMode(float cm) { mArtificialCM = cm; }
+
+  void setCorrectOutputForPedestal(bool corret = true) { mCorrectOutputForPedestal = corret; }
+  bool getCorrectOutputForPedestal() const { return mCorrectOutputForPedestal; }
+
+  /// Add zeros for pads without signal
+  void setAddSubthreshold(bool addSubthreshold) { mSubthreshold = addSubthreshold; }
+  bool getAddSubthreshold() const { return mSubthreshold; }
+
+  static float decodeInjectedCMValue(float lower, float upper);
+
+ private:
+  inline static int sNThreads{1};        ///< Number of parallel threads for the CM calculation
+  int mNPadsCompRamdom{10};              ///< Number of random pads to compare with to check if the present pad is empty
+  int mNPadsCompMin{7};                  ///< Minimum number of neighbouring pads with q close to present pad to define this as empty
+  int mNPadsMinCM{0};                    ///< Minimum number of pads required in the CM calculation to be used for digit correction
+  float mQEmpty{2};                      ///< Threshold to enter check for empty pad
+  float mQComp{1};                       ///< Threshold for comparison with random pads
+  float mQCompScaleThreshold{0};         ///< Charge threshold from which on to increase mQComp
+  float mQCompScale{0};                  ///< Slope with which to increase mQComp if below mQCompScaleThreshold
+  float mSumPosThreshold{2};             ///< calculate sumPos > mSumPosThreshold, sumNeg M<= mSumPosThreshold
+  float mOccupancyThreshold{3};          ///< calculate number of pads > mQCompScaleThreshold after CM correction
+  bool mLimitKFactor{false};             ///< Limit the k-factor precision to 2I6F
+  bool mLimitPedestal{false};            ///< Limit the preestal precision to 10I2F
+  int mSubthreshold{0};                  ///< Add data for pads without signal. 1 = add zeros; 2 = add random noise
+  float mArtificialCM{};                 ///< artificial common mode signals
+  bool mCorrectOutputForPedestal{false}; ///< correct the writte out ADC for the pedestal value
+
+  CalPadMapType mPadMaps; ///< Pad-by-pad CRU configuration values (Pedestal, Noise, ITF + CM parameters)
+
+  struct pos {
+    int row;
+    int pad;
+  };
+
+  // positions of lower words per CRU in sector
+  const std::array<pos, 10> mCMInjectIDLower{
+    // row0 pad0 row1 pad1
+    pos{0, 2},
+    pos{20, 1},
+    pos{32, 2},
+    pos{51, 1},
+    pos{63, 1},
+    pos{84, 1},
+    pos{97, 1},
+    pos{116, 2},
+    pos{127, 2},
+    pos{142, 0},
+  };
+
+  // positions of upper words per CRU in sector
+  const std::array<pos, 10> mCMInjectIDUpper{
+    // row0 pad0 row1 pad1
+    pos{0, 3},
+    pos{20, 3},
+    pos{32, 3},
+    pos{51, 3},
+    pos{63, 2},
+    pos{84, 4},
+    pos{97, 2},
+    pos{115, 5},
+    pos{127, 3},
+    pos{142, 4},
+  };
+
+  /// Return the value stored in mPadMaps["calibName"]
+  /// \param calibName name of calibraion in mPadMaps
+  /// \param cru CRU number
+  /// \param pad Pad number within the CRU
+  float getCalPadValue(const std::string calibName, int icru, int pad) const;
+
+  bool padMapExists(const std::string& calibName);
+
+  ClassDefNV(CommonModeCorrection, 2);
+};
+
+} // namespace o2::tpc
+#endif

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

@@ -396,6 +396,14 @@ class Mapper
 
   bool isOutOfSector(GlobalPosition3D posEle, const Sector& sector, const float margin = 0.f) const;
 
+  static bool isEdgePad(int rowInSector, int padInRow);
+  static bool isFirstOrLastRowInStack(int rowInSector);
+  static bool isBelowSpacerCross(int rowInSector, int padInRow);
+  static bool isHighCouplingPad(int rowInSector, int padInRow)
+  {
+    return isEdgePad(rowInSector, padInRow) || isFirstOrLastRowInStack(rowInSector) || isBelowSpacerCross(rowInSector, padInRow);
+  }
+
   static constexpr unsigned short getNumberOfIROCs() { return 36; }
   static constexpr unsigned short getNumberOfOROCs() { return 36; }
   static constexpr unsigned short getPadsInIROC() { return mPadsInIROC; }
@@ -523,6 +531,7 @@ class Mapper
   static constexpr unsigned int GLOBALPADOFFSET[NREGIONS]{0, 1200, 2400, 3840, 5280, 6720, 8160, 9760, 11360, 12960};                                                 ///< offset of number of pads for region
   static constexpr unsigned int ROWSPERREGION[NREGIONS]{17, 15, 16, 15, 18, 16, 16, 14, 13, 12};                                                                      ///< number of pad rows for region
   static constexpr unsigned int ROWOFFSET[NREGIONS]{0, 17, 32, 48, 63, 81, 97, 113, 127, 140};                                                                        ///< offset to calculate local row from global row
+  static constexpr unsigned int ROWOFFSETSTACK[4]{0, 63, 97, 127};                                                                                                    ///< offset to calculate local row from global row
   static constexpr float REGIONAREA[NREGIONS]{374.4f, 378.f, 453.6f, 470.88f, 864.f, 864.f, 1167.36f, 1128.96f, 1449.6f, 1456.8f};                                    ///< volume of each region in cm^2
   static constexpr float INVPADAREA[NREGIONS]{1 / 0.312f, 1 / 0.315f, 1 / 0.315f, 1 / 0.327f, 1 / 0.6f, 1 / 0.6f, 1 / 0.7296f, 1 / 0.7056f, 1 / 0.906f, 1 / 0.9105f}; ///< inverse size of the pad area padwidth*padLength
   static constexpr unsigned REGION[PADROWS] = {
@@ -542,7 +551,7 @@ class Mapper
     {0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4},             // region 7
     {0, 0, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5},                // region 8
     {0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 5}                    // region 9
-  };                                                        ///< additional pads per row compared to first row
+  }; ///< additional pads per row compared to first row
   const inline static std::vector<unsigned int> OFFSETCRULOCAL[NREGIONS]{
     {0, 66, 132, 198, 266, 334, 402, 472, 542, 612, 684, 756, 828, 902, 976, 1050, 1124},         // region 0
     {0, 76, 152, 228, 306, 384, 462, 542, 622, 702, 784, 866, 948, 1032, 1116},                   // region 1
@@ -554,7 +563,7 @@ class Mapper
     {0, 110, 220, 332, 444, 556, 670, 784, 898, 1014, 1130, 1246, 1364, 1482},                    // region 7
     {0, 118, 236, 356, 476, 598, 720, 844, 968, 1092, 1218, 1344, 1472},                          // region 8
     {0, 128, 258, 388, 520, 652, 784, 918, 1052, 1188, 1324, 1462}                                // region 9
-  };                                                                                              ///< row offset in cru for given local pad row
+  }; ///< row offset in cru for given local pad row
   const inline static std::vector<unsigned int> PADSPERROW[NREGIONS]{
     {66, 66, 66, 68, 68, 68, 70, 70, 70, 72, 72, 72, 74, 74, 74, 74, 76},      // region 0
     {76, 76, 76, 78, 78, 78, 80, 80, 80, 82, 82, 82, 84, 84, 84},              // region 1
@@ -566,7 +575,7 @@ class Mapper
     {110, 110, 112, 112, 112, 114, 114, 114, 116, 116, 116, 118, 118, 118},    // region 7
     {118, 118, 120, 120, 122, 122, 124, 124, 124, 126, 126, 128, 128},         // region 8
     {128, 130, 130, 132, 132, 132, 134, 134, 136, 136, 138, 138}               // region 9
-  };                                                                           ///< number of pads per row in region
+  }; ///< number of pads per row in region
   static constexpr unsigned int OFFSETCRUGLOBAL[PADROWS]{
     0, 66, 132, 198, 266, 334, 402, 472, 542, 612, 684, 756, 828, 902, 976, 1050, 1124,         // region 0
     0, 76, 152, 228, 306, 384, 462, 542, 622, 702, 784, 866, 948, 1032, 1116,                   // region 1
@@ -578,7 +587,7 @@ class Mapper
     0, 110, 220, 332, 444, 556, 670, 784, 898, 1014, 1130, 1246, 1364, 1482,                    // region 7
     0, 118, 236, 356, 476, 598, 720, 844, 968, 1092, 1218, 1344, 1472,                          // region 8
     0, 128, 258, 388, 520, 652, 784, 918, 1052, 1188, 1324, 1462                                // region 9
-  };                                                                                            ///< row offset in cru for given global pad row
+  }; ///< row offset in cru for given global pad row
 
   static constexpr unsigned int LinksPerRegionPerEndpoint[NREGIONS][NENDPOINTS]{
     {8, 7},   // region 0
@@ -591,7 +600,7 @@ class Mapper
     {10, 10}, // region 7
     {10, 10}, // region 8
     {10, 10}, // region 9
-  };          ///< number of links per region per end point
+  }; ///< number of links per region per end point
 
  private:
   Mapper(const std::string& mappingDir);