From 9890a03c47d07624734a55e9dd680d187a3218d7 Mon Sep 17 00:00:00 2001 From: wiechula Date: Thu, 19 Sep 2024 13:22:12 +0200 Subject: [PATCH] Fix and extension --- .../include/TPCBase/CommonModeCorrection.h | 32 ++++- .../TPC/base/src/CommonModeCorrection.cxx | 122 ++++++++++++++---- 2 files changed, 122 insertions(+), 32 deletions(-) diff --git a/Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h b/Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h index 416a55877e2ec..286897fdf8312 100644 --- a/Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h +++ b/Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h @@ -48,9 +48,15 @@ class CommonModeCorrection struct CMInfo { float cmValue{}; + float cmValueStd{}; int nPadsUsed{}; }; + struct CMDebug { + std::vector nPadsOk{}; + std::vector adcDist{}; + }; + using CalPadMapType = std::unordered_map; /// Calculation of the common mode value @@ -59,7 +65,7 @@ class CommonModeCorrection /// \param cmKValues corresponding pad-by-pad common mode k-factors /// \param pedestals corresponding pad-by-pad pedestals /// \param - CMInfo getCommonMode(gsl::span values, gsl::span cmKValues, gsl::span pedestals) const; + CMInfo getCommonMode(gsl::span values, gsl::span cmKValues, gsl::span pedestals, CMDebug* cmDebug = nullptr) const; CMInfo getCommonMode(const std::vector& values, const std::vector& cmKValues, const std::vector& 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)); } @@ -76,6 +82,14 @@ class CommonModeCorrection 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; } + /// Pad maps loaded from FEEConfig void setPadMaps(CalPadMapType& padMaps) { mPadMaps = padMaps; } @@ -103,9 +117,9 @@ class CommonModeCorrection /// \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& digits, std::vector>& cmValues, bool negativeOnly = false) const; + int correctDigits(std::vector& digits, std::vector>& cmValues, bool negativeOnly = false, std::vector>* 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); + 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, int minTimeBin = -1, int maxTimeBin = -1); void limitKFactorPrecision(bool limit = true) { mLimitKFactor = limit; } void limitPedestalPrecision(bool limit = true) { mLimitPedestal = limit; } @@ -117,14 +131,20 @@ class CommonModeCorrection /// \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; } + 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 + 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 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 bool mLimitKFactor{false}; ///< Limit the k-factor precision to 2I6F bool mLimitPedestal{false}; ///< Limit the preestal precision to 10I2F + float mArtificialCM{}; ///< artificial common mode signals CalPadMapType mPadMaps; ///< Pad-by-pad CRU configuration values (Pedestal, Noise, ITF + CM parameters) diff --git a/Detectors/TPC/base/src/CommonModeCorrection.cxx b/Detectors/TPC/base/src/CommonModeCorrection.cxx index 2efecf396d796..ffb039bc6e789 100644 --- a/Detectors/TPC/base/src/CommonModeCorrection.cxx +++ b/Detectors/TPC/base/src/CommonModeCorrection.cxx @@ -12,7 +12,7 @@ /// \file CommonModeCorrection.cxx /// \brief Calculate the common mode correction factor -#include +// #include #include #include #include @@ -28,7 +28,7 @@ using namespace o2::tpc; using namespace o2::tpc::cru_calib_helpers; -CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span values, gsl::span cmKValues, gsl::span pedestals) const +CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span values, gsl::span cmKValues, gsl::span pedestals, CMDebug* cmDebug) const { if (values.size() == 0) { return CMInfo{}; @@ -41,6 +41,12 @@ CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span::RandomType::Flat); std::vector adcCM; //< ADC values used for common mode calculation + CMInfo cmInfo; + if (cmDebug) { + cmDebug->nPadsOk.resize(mNPadsCompRamdom + 1); + cmDebug->adcDist.resize(10); + } + for (size_t iPad = 0; iPad < values.size(); ++iPad) { const float kCM = mLimitKFactor ? fixedSizeToFloat<6>(floatToFixedSize<8, 6>(cmKValues[iPad])) : cmKValues[iPad]; const float pedestal = mLimitPedestal ? fixedSizeToFloat(floatToFixedSize(pedestals[iPad])) : pedestals[iPad]; @@ -51,6 +57,12 @@ CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span(floatToFixedSize<8, 6>(cmKValues[padRnd])) : cmKValues[padRnd]; const float pedestalRnd = mLimitPedestal ? fixedSizeToFloat(floatToFixedSize(pedestals[padRnd])) : pedestals[padRnd]; const float adcPadRnd = values[padRnd] - pedestalRnd; - if (std::abs(adcPadNorm - adcPadRnd / kCMRnd) < mQComp) { + const float adcDist = std::abs(adcPadNorm - adcPadRnd / kCMRnd); + if (cmDebug) { + const size_t distPos = std::min(cmDebug->adcDist.size() - 1, size_t(adcDist / 0.5)); + ++cmDebug->adcDist[distPos]; + } + if (adcDist < mQComp) { ++nPadsOK; } } + if (cmDebug) { + ++cmDebug->nPadsOk[nPadsOK]; + } + if (nPadsOK >= mNPadsCompMin) { adcCM.emplace_back(adcPadNorm); } } const int entriesCM = int(adcCM.size()); - float commonMode = std::accumulate(adcCM.begin(), adcCM.end(), 0.f); + float commonMode = 0; // std::accumulate(adcCM.begin(), adcCM.end(), 0.f); + float commonModeStd = 0; if (entriesCM > 0) { + std::for_each(adcCM.begin(), adcCM.end(), [&commonMode, &commonModeStd](const auto val) { + commonMode += val; + commonModeStd += val * val; + }); commonMode /= float(entriesCM); + commonModeStd = std::sqrt(std::abs(commonModeStd / entriesCM - commonMode * commonMode)); } - return CMInfo{commonMode, entriesCM}; + cmInfo.cmValue = commonMode; + cmInfo.cmValueStd = commonModeStd; + cmInfo.nPadsUsed = entriesCM; + return cmInfo; } void CommonModeCorrection::loadDefaultPadMaps(FEEConfig::Tags tag) @@ -126,61 +156,90 @@ CommonModeCorrection::CMdata CommonModeCorrection::collectCMdata(const std::vect return data; } -int CommonModeCorrection::correctDigits(std::vector& digits, std::vector>& cmValues, bool negativeOnly) const +int CommonModeCorrection::correctDigits(std::vector& digits, std::vector>& cmValues, bool negativeOnly, std::vector>* cmDebug, int minTimeBin, int maxTimeBin) const { // calculation common mode values - int maxTimeBin = -1; + int maxTimeBinProcessed = -1; int lastCRU = -1; int lastTimeBin = -1; CMdata data; const auto& cmkValues = mPadMaps.at("CMkValues"); const auto& pedestals = mPadMaps.at("Pedestals"); + bool doArtificialCM = std::abs(mArtificialCM) > 0; + for (size_t iDigit = 0; iDigit < digits.size(); ++iDigit) { auto& digit = digits[iDigit]; - if (((lastCRU > 0) && (digit.getCRU() != lastCRU)) || ((lastTimeBin > 0) && (digit.getTimeStamp() != lastTimeBin))) { + const auto timeBin = digit.getTimeStamp(); + if ((minTimeBin > -1) && (timeBin < minTimeBin)) { + continue; + } + if ((maxTimeBin > -1) && (timeBin > maxTimeBin)) { + continue; + } + if ((lastCRU > -1) && ((digit.getCRU() != lastCRU) || (digit.getTimeStamp() != lastTimeBin))) { auto& cmValuesCRU = cmValues[lastCRU]; if (cmValuesCRU.size() <= lastTimeBin) { cmValuesCRU.resize(cmValuesCRU.size() + 500); + if (cmDebug) { + (*cmDebug)[lastCRU].resize((*cmDebug)[lastCRU].size() + 500); + } } - cmValuesCRU[lastTimeBin] = getCommonMode(data.adcValues, data.cmKValues, data.pedestals); - // LOGP(info, "processing CRU {}, timeBin {}, CM = {}", lastCRU, lastTimeBin, cmVal); + cmValuesCRU[lastTimeBin] = getCommonMode(data.adcValues, data.cmKValues, data.pedestals, cmDebug ? &((*cmDebug)[lastCRU][lastTimeBin]) : nullptr); + // LOGP(info, "processing CRU {}, timeBin {}, CM = {}", lastCRU, lastTimeBin, cmValuesCRU[lastTimeBin].cmValue); data.clear(); } const auto sector = CRU(digit.getCRU()).sector(); - data.adcValues.emplace_back(digit.getChargeFloat()); - data.cmKValues.emplace_back(cmkValues.getValue(sector, digit.getRow(), digit.getPad())); - data.pedestals.emplace_back(pedestals.getValue(sector, digit.getRow(), digit.getPad())); + const auto cmkValue = cmkValues.getValue(sector, digit.getRow(), digit.getPad()); + const auto pedestal = pedestals.getValue(sector, digit.getRow(), digit.getPad()); + float charge = digit.getChargeFloat(); + if (doArtificialCM) { + charge = std::clamp(charge + mArtificialCM * cmkValue, 0.f, 1023.f); + } + data.adcValues.emplace_back(charge); + data.cmKValues.emplace_back(cmkValue); + data.pedestals.emplace_back(pedestal); lastCRU = digit.getCRU(); - lastTimeBin = digit.getTimeStamp(); - maxTimeBin = std::max(lastTimeBin, maxTimeBin); + lastTimeBin = timeBin; + maxTimeBinProcessed = std::max(lastTimeBin, maxTimeBinProcessed); } { auto& cmValuesCRU = cmValues[lastCRU]; if (cmValuesCRU.size() <= lastTimeBin) { cmValuesCRU.resize(cmValuesCRU.size() + 500); + if (cmDebug) { + (*cmDebug)[lastCRU].resize((*cmDebug)[lastCRU].size() + 500); + } } - cmValuesCRU[lastTimeBin] = getCommonMode(data.adcValues, data.cmKValues, data.pedestals); - // LOGP(info, "processing CRU {}, timeBin {}, CM = {}", lastCRU, lastTimeBin, cmVal); + cmValuesCRU[lastTimeBin] = getCommonMode(data.adcValues, data.cmKValues, data.pedestals, cmDebug ? &((*cmDebug)[lastCRU][lastTimeBin]) : nullptr); + // LOGP(info, "processing CRU {}, timeBin {}, CM = {}", lastCRU, lastTimeBin, cmValuesCRU[lastTimeBin].cmValue); data.clear(); } // apply correction for (auto& digit : digits) { + const auto timeBin = digit.getTimeStamp(); + if ((minTimeBin > -1) && (timeBin < minTimeBin)) { + continue; + } + if ((maxTimeBin > -1) && (timeBin > maxTimeBin)) { + continue; + } const auto sector = CRU(digit.getCRU()).sector(); const auto cmKValue = cmkValues.getValue(sector, digit.getRow(), digit.getPad()); + // LOGP(info, "correcting value for CRU {}, time bin {}", digit.getCRU(), digit.getTimeStamp()); const auto cmValue = cmValues[digit.getCRU()][digit.getTimeStamp()].cmValue; if (!negativeOnly || cmValue < 0) { digit.setCharge(digit.getCharge() - cmValue * cmKValue); } } - return maxTimeBin; + return maxTimeBinProcessed; } -void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t maxEntries, std::string_view digitFileOut, std::string_view cmFileOut, bool negativeOnly, int nThreads, bool writeOnlyCM) +void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t maxEntries, std::string_view digitFileOut, std::string_view cmFileOut, bool negativeOnly, int nThreads, bool writeOnlyCM, bool writeDebug, int minTimeBin, int maxTimeBin) { TChain* tree = o2::tpc::utils::buildChain(fmt::format("ls {}", digiFileIn), "o2sim", "o2sim"); @@ -217,7 +276,12 @@ void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t m LOGP(info, "Processing entry {}/{}", iTF + 1, nEntries); std::vector> cmValues; // CRU * timeBin + std::vector> cmDebug; // CRU * timeBin + cmValues.resize(CRU::MaxCRU); + if (writeDebug) { + cmDebug.resize(CRU::MaxCRU); + } int maxTimeBin = -1; auto worker = [&](int iTread) { @@ -227,7 +291,7 @@ void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t m if (!digits || (digits->size() == 0)) { continue; } - const int maxTimeBinSector = correctDigits(*digits, cmValues, negativeOnly); + const int maxTimeBinSector = correctDigits(*digits, cmValues, negativeOnly, writeDebug ? &cmDebug : nullptr, minTimeBin, maxTimeBin); { static std::mutex maxMutex; std::lock_guard lock{maxMutex}; @@ -247,10 +311,6 @@ void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t m th.join(); } - if (tOut) { - tOut->Fill(); - } - for (int iCRU = 0; iCRU < cmValues.size(); ++iCRU) { int maxTBCRU = std::min(maxTimeBin, int(cmValues[iCRU].size())); for (int iTimeBin = 0; iTimeBin < maxTBCRU; ++iTimeBin) { @@ -258,11 +318,21 @@ void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t m << "iTF=" << iTF << "iCRU=" << iCRU << "iTimeBin=" << iTimeBin - << "cm=" << cmValues[iCRU][iTimeBin].cmValue - << "cmEntries=" << cmValues[iCRU][iTimeBin].nPadsUsed + << "cmInfo=" << cmValues[iCRU][iTimeBin]; + if (writeDebug) { + pcstream << "cm" + << "cmDebug=" << cmDebug[iCRU][iTimeBin]; + } + // << "cm=" << cmValues[iCRU][iTimeBin].cmValue + // << "cmEntries=" << cmValues[iCRU][iTimeBin].nPadsUsed + pcstream << "cm" << "\n"; } } + + if (tOut) { + tOut->Fill(); + } } pcstream.Close();