From fcc6aba5c4e00eaff713eb6c6a43f4cf6a3d75ad Mon Sep 17 00:00:00 2001 From: wiechula Date: Tue, 24 Sep 2024 11:18:43 +0200 Subject: [PATCH 1/8] Add option to rescale IT fraction for OROCs --- Detectors/TPC/calibration/macro/prepareITFiles.C | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Detectors/TPC/calibration/macro/prepareITFiles.C b/Detectors/TPC/calibration/macro/prepareITFiles.C index ff58d57583342..eac0355e0ddfd 100644 --- a/Detectors/TPC/calibration/macro/prepareITFiles.C +++ b/Detectors/TPC/calibration/macro/prepareITFiles.C @@ -31,7 +31,7 @@ using namespace o2::tpc::cru_calib_helpers; using namespace o2::tpc; -void prepareCMFiles(const std::string_view itDataFile, std::string outputDir = "./") +void prepareITFiles(const std::string_view itDataFile, std::string outputDir = "./", float orocFractionScale = 1.f) { const auto& mapper = Mapper::instance(); @@ -102,6 +102,10 @@ void prepareCMFiles(const std::string_view itDataFile, std::string outputDir = " float fractionVal = rocFraction.getValue(ipad); float expLambdaVal = rocExpLambda.getValue(ipad); + if (roc.isOROC()) { + fractionVal *= orocFractionScale; + } + if ((fractionVal <= 0) || (fractionVal > 0.6)) { LOGP(error, "Too fraction value in ROC {:2}, CRU {:3}, fec in CRU: {:2}, SAMPA: {}, channel: {:2}: {:.4f}, setting value to roc mean {}", iroc, cruID, fecInPartition, sampa, sampaChannel, fractionVal, meanFraction); fractionVal = meanFraction; From 34ca379b2c57fddcf2febf5f35b28df411f4383e Mon Sep 17 00:00:00 2001 From: wiechula Date: Fri, 6 Sep 2024 12:58:27 +0200 Subject: [PATCH 2/8] Common mode calculation --- Detectors/TPC/base/CMakeLists.txt | 2 + .../include/TPCBase/CommonModeCorrection.h | 246 ++++++++ .../TPC/base/src/CommonModeCorrection.cxx | 568 ++++++++++++++++++ Detectors/TPC/base/src/TPCBaseLinkDef.h | 3 + 4 files changed, 819 insertions(+) create mode 100644 Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h create mode 100644 Detectors/TPC/base/src/CommonModeCorrection.cxx diff --git a/Detectors/TPC/base/CMakeLists.txt b/Detectors/TPC/base/CMakeLists.txt index c13fec6f03ab7..d4c1bc4602d54 100644 --- a/Detectors/TPC/base/CMakeLists.txt +++ b/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 diff --git a/Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h b/Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h new file mode 100644 index 0000000000000..a222327d2b434 --- /dev/null +++ b/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 +#include +#include + +#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 adcValues; + std::vector cmKValues; + std::vector 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 nPadsOk{}; + std::vector adcDist{}; + }; + + using CalPadMapType = std::unordered_map; + + /// 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 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)); } + + 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& digits, int cru, int timeBin); + + int getCommonMode(std::vector& digits, std::vector>& cmValues, bool negativeOnly = false, bool hasInjectedCMValue = false, std::vector>* 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& digits, std::vector>& cmValues, bool negativeOnly = false, bool hasInjectedCMValue = 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, 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 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 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 diff --git a/Detectors/TPC/base/src/CommonModeCorrection.cxx b/Detectors/TPC/base/src/CommonModeCorrection.cxx new file mode 100644 index 0000000000000..729fb408eb204 --- /dev/null +++ b/Detectors/TPC/base/src/CommonModeCorrection.cxx @@ -0,0 +1,568 @@ +// 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.cxx +/// \brief Calculate the common mode correction factor + +// #include +#include +#include +#include +#include "CCDB/CcdbApi.h" +#include "TPCBase/CommonModeCorrection.h" +#include "TPCBase/Mapper.h" +#include "TPCBase/Utils.h" +#include "TPCBase/CRUCalibHelpers.h" +#include "TChain.h" +#include "TROOT.h" +#include "TFile.h" +#include "MathUtils/RandomRing.h" +#include "CommonUtils/TreeStreamRedirector.h" + +using namespace o2::tpc; +using namespace o2::tpc::cru_calib_helpers; +CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span values, gsl::span cmKValues, gsl::span pedestals, CMDebug* cmDebug) const +{ + if (values.size() == 0) { + return CMInfo{}; + } + // sanity check + if (values.size() != cmKValues.size() || values.size() != pedestals.size()) { + LOGP(error, "vector sizes of input values, cmKValues and pedestals don't match: {}, {}, {}", values.size(), cmKValues.size(), pedestals.size()); + return CMInfo{}; + } + static math_utils::RandomRing random(math_utils::RandomRing<>::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]; + const float adcPadRaw = values[iPad]; + const float adcPad = adcPadRaw - pedestal; + const float adcPadNorm = (kCM > 0) ? adcPad / kCM : 0; + + if (adcPadRaw > 1023.7) { + ++cmInfo.nSaturation; + } + + if (adcPad > mQEmpty) { + continue; + } + + float qCompAdd = 0; + if ((mQCompScaleThreshold < 0) && (adcPadNorm < mQCompScaleThreshold)) { + qCompAdd = (mQCompScaleThreshold - adcPadNorm) * mQCompScale; + LOGP(info, "Setting qCompAdd to {} for {}", qCompAdd, adcPadNorm); + } + + int nPadsOK = 0; + + for (int iRnd = 0; iRnd < mNPadsCompRamdom; ++iRnd) { + int padRnd = 0; + do { + padRnd = int(random.getNextValue() * (values.size() - 1)); + } while (padRnd == iPad); + const float kCMRnd = mLimitKFactor ? fixedSizeToFloat<6>(floatToFixedSize<8, 6>(cmKValues[padRnd])) : cmKValues[padRnd]; + const float pedestalRnd = mLimitPedestal ? fixedSizeToFloat(floatToFixedSize(pedestals[padRnd])) : pedestals[padRnd]; + const float adcPadRnd = values[padRnd] - pedestalRnd; + const float adcPadRndNorm = (kCMRnd > 0) ? adcPadRnd / kCMRnd : 0; + const float adcDist = std::abs(adcPadNorm - adcPadRndNorm); + 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 = 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)); + } + cmInfo.cmValue = commonMode; + cmInfo.cmValueStd = commonModeStd; + cmInfo.nPadsUsed = entriesCM; + + 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]; + const float adcPadRaw = values[iPad]; + const float adcPad = adcPadRaw - pedestal; + const float adcPadNorm = (kCM > 0) ? adcPad / kCM : 0; + const float adcPadCorr = adcPad - kCM * commonMode; + + if (adcPadCorr > mSumPosThreshold) { + cmInfo.sumPos += adcPadCorr; + } else { + cmInfo.sumNeg += adcPadNorm; + ++cmInfo.nNeg; + } + + if (mOccupancyThreshold > 0) { + if (adcPadCorr > mOccupancyThreshold) { + ++cmInfo.nOccupancy; + } + } + } + + return cmInfo; +} + +void CommonModeCorrection::loadDefaultPadMaps(FEEConfig::Tags tag) +{ + o2::ccdb::CcdbApi cdbApi; + cdbApi.init("http://alice-ccdb.cern.ch"); + const auto feeConfig = cdbApi.retrieveFromTFileAny("TPC/Config/FEE", {}, long(tag)); + if (!feeConfig) { + LOGP(error, "Could not retrieve pad maps"); + return; + } + mPadMaps = feeConfig->padMaps; + delete feeConfig; +} + +CommonModeCorrection::CMdata CommonModeCorrection::collectCMdata(const std::vector& digits, int cru, int timeBin) +{ + + CMdata data; + if (!padMapExists("CMkValues") || padMapExists("Pedestals")) { + return data; + } + + for (const auto& digit : digits) { + if (digit.getTimeStamp() < timeBin) { + continue; + } + + if (digit.getTimeStamp() > timeBin) { + break; + } + + if (digit.getCRU() < cru) { + continue; + } + + if (digit.getCRU() > cru) { + break; + } + + const auto sector = CRU(digit.getCRU()).sector(); + data.adcValues.emplace_back(digit.getChargeFloat()); + data.cmKValues.emplace_back(mPadMaps["CMkValues"].getValue(sector, digit.getRow(), digit.getPad())); + data.pedestals.emplace_back(mPadMaps["Pedestals"].getValue(sector, digit.getRow(), digit.getPad())); + } + return data; +} + +int CommonModeCorrection::getCommonMode(std::vector& digits, std::vector>& cmValues, bool negativeOnly, bool hasInjectedCMValue, std::vector>* cmDebug, int minTimeBin, int maxTimeBin) const +{ + // calculation common mode values + 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 decoding of the injected common mode signals + float cmInjectedLower{}; + float cmInjectedUpper{}; + + for (size_t iDigit = 0; iDigit < digits.size(); ++iDigit) { + auto& digit = digits[iDigit]; + 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(lastTimeBin + 500); + if (cmDebug) { + (*cmDebug)[lastCRU].resize(lastTimeBin + 500); + } + } + if (mSubthreshold > 0) { + const size_t nPadsCRU = Mapper::PADSPERREGION[lastCRU % 10]; + const auto dataSize = data.adcValues.size(); + if (dataSize < nPadsCRU) { + data.resize(nPadsCRU); + if (mSubthreshold == 2) { + for (size_t i = dataSize; i < nPadsCRU; ++i) { + data.adcValues[i] = gRandom->Gaus(); + } + } + } + } + cmValuesCRU[lastTimeBin] = getCommonMode(data.adcValues, data.cmKValues, data.pedestals, cmDebug ? &((*cmDebug)[lastCRU][lastTimeBin]) : nullptr); + if (hasInjectedCMValue) { + cmValuesCRU[lastTimeBin].cmValueCRU = decodeInjectedCMValue(cmInjectedLower, cmInjectedUpper); + } + // LOGP(info, "processing CRU {}, timeBin {}, CM = {}", lastCRU, lastTimeBin, cmValuesCRU[lastTimeBin].cmValue); + + data.clear(); + } + const auto sector = CRU(digit.getCRU()).sector(); + 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); + } + lastCRU = digit.getCRU(); + lastTimeBin = timeBin; + maxTimeBinProcessed = std::max(lastTimeBin, maxTimeBinProcessed); + + bool isInjectedCMPad = false; + if (hasInjectedCMValue) { + const auto posLow = mCMInjectIDLower[lastCRU % 10]; + const auto posUpper = mCMInjectIDUpper[lastCRU % 10]; + const auto row = digit.getRow(); + const auto pad = digit.getPad(); + if (row == posLow.row) { + if (pad == posLow.pad) { + cmInjectedLower = digit.getChargeFloat(); + isInjectedCMPad = true; + // LOGP(info, "setting lower CM value cru {}, row {}, pad {}: {:012b}", digit.getCRU(), row, pad, floatToFixedSize(digit.getChargeFloat())); + } + } + if (row == posUpper.row) { + if (pad == posUpper.pad) { + cmInjectedUpper = digit.getChargeFloat(); + isInjectedCMPad = true; + // LOGP(info, "setting upper CM value cru {}, row {}, pad {}: {:012b}", digit.getCRU(), row, pad, floatToFixedSize(digit.getChargeFloat())); + if (cmInjectedUpper == 0) { + LOGP(info, "cm upper = 0 cru {}, row {}, pad {}", digit.getCRU(), row, pad); + } + } + } + } + + if (!isInjectedCMPad) { + data.adcValues.emplace_back(charge); + data.cmKValues.emplace_back(cmkValue); + data.pedestals.emplace_back(pedestal); + } + } + { + auto& cmValuesCRU = cmValues[lastCRU]; + if (cmValuesCRU.size() <= lastTimeBin) { + cmValuesCRU.resize(lastTimeBin + 500); + if (cmDebug) { + (*cmDebug)[lastCRU].resize(lastTimeBin + 500); + } + } + 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); + + if (hasInjectedCMValue) { + cmValuesCRU[lastTimeBin].cmValueCRU = decodeInjectedCMValue(cmInjectedLower, cmInjectedUpper); + } + + data.clear(); + } + return maxTimeBinProcessed; +} + +int CommonModeCorrection::correctDigits(std::vector& digits, std::vector>& cmValues, bool negativeOnly, bool hasInjectedCMValue, std::vector>* cmDebug, int minTimeBin, int maxTimeBin) const +{ + const auto maxTimeBinProcessed = getCommonMode(digits, cmValues, negativeOnly, hasInjectedCMValue, cmDebug, minTimeBin, maxTimeBin); + const auto& cmkValues = mPadMaps.at("CMkValues"); + const auto& pedestals = mPadMaps.at("Pedestals"); + // ===| 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; + const auto cmNPads = cmValues[digit.getCRU()][digit.getTimeStamp()].nPadsUsed; + if ((!negativeOnly || cmValue < 0) && (cmNPads > mNPadsMinCM)) { + digit.setCharge(digit.getCharge() - cmValue * cmKValue); + if (mCorrectOutputForPedestal) { + const auto sector = CRU(digit.getCRU()).sector(); + const auto pedestal = pedestals.getValue(sector, digit.getRow(), digit.getPad()); + digit.setCharge(digit.getChargeFloat() - pedestal); + } + } + } + + 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, bool writeDebug, bool hasInjectedCMValue, int minTimeBin, int maxTimeBin) +{ + ROOT::EnableThreadSafety(); + + TChain* tree = o2::tpc::utils::buildChain(fmt::format("ls {}", digiFileIn), "o2sim", "o2sim"); + Long64_t nEntries = tree->GetEntries(); + if (maxEntries > 0) { + nEntries = std::min(nEntries, maxEntries); + } + + if (mPadMaps.find("Pedestals") == mPadMaps.end()) { + LOGP(info, "Using empty pedestals"); + mPadMaps["Pedestals"] = CalPad("Pedestals"); + } + + std::unique_ptr fOut; + std::unique_ptr tOut; + if (!writeOnlyCM) { + fOut.reset(TFile::Open(digitFileOut.data(), "RECREATE")); + fOut->SetCompressionLevel(5); // zstd default level + fOut->SetCompressionAlgorithm(5); // zstd + tOut = std::make_unique("o2sim", "o2sim"); + } + + std::array*, 36> digitizedSignal; + std::array outBranches{}; + for (size_t iSec = 0; iSec < digitizedSignal.size(); ++iSec) { + digitizedSignal[iSec] = nullptr; + tree->SetBranchAddress(Form("TPCDigit_%zu", iSec), &digitizedSignal[iSec]); + if (tOut) { + outBranches[iSec] = tOut->Branch(Form("TPCDigit_%zu", iSec), &digitizedSignal[iSec]); + } + } + + o2::utils::TreeStreamRedirector pcstream(cmFileOut.data(), "recreate"); + pcstream.GetFile()->SetCompressionAlgorithm(5); + pcstream.GetFile()->SetCompressionLevel(5); + + for (Long64_t iTF = 0; iTF < nEntries; ++iTF) { + tree->GetEntry(iTF); + 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 maxTimeBinSeen = -1; + + auto worker = [&](int iTread) { + // for (size_t iSector = 0; iSector < 36; ++iSector) { + for (size_t iSector = iTread; iSector < 36; iSector += nThreads) { + LOGP(info, "Processing entry {}/{}, starting sector {}", iTF + 1, nEntries, iSector); + auto digits = digitizedSignal[iSector]; + int maxTimeBinSector = 0; + if (digits && (digits->size() > 0)) { + maxTimeBinSector = correctDigits(*digits, cmValues, negativeOnly, hasInjectedCMValue, writeDebug ? &cmDebug : nullptr, minTimeBin, maxTimeBin); + } + { + static std::mutex maxMutex; + std::lock_guard lock{maxMutex}; + maxTimeBinSeen = std::max(maxTimeBinSeen, maxTimeBinSector); + if (outBranches[iSector]) { + outBranches[iSector]->Fill(); + LOGP(info, "Filling branch for sector {}", iSector); + } + } + } + }; + + std::vector threads(nThreads); + + for (int i = 0; i < threads.size(); i++) { + threads[i] = std::thread(worker, i); + } + + // wait for the threads to finish + for (auto& th : threads) { + th.join(); + } + + size_t maxTimeCRU = 0; + for (int iCRU = 0; iCRU < cmValues.size(); ++iCRU) { + maxTimeCRU = std::max(maxTimeCRU, cmValues[iCRU].size()); + } + const int maxTBCRU = std::min(maxTimeBinSeen, int(maxTimeCRU)); + + for (int iTimeBin = 0; iTimeBin < maxTBCRU; ++iTimeBin) { + + std::vector cm(CRU::MaxCRU); + std::vector cmD(CRU::MaxCRU); + std::vector sumPosStack(36 * 4); + std::vector nPosStack(36 * 4); + std::vector nSaturationStack(36 * 4); + std::vector sumPosStackCRU(CRU::MaxCRU); + std::vector sumPosStackCRUCorr(CRU::MaxCRU); + std::vector nSaturationStackCRU(CRU::MaxCRU); + + for (int iCRU = 0; iCRU < cmValues.size(); ++iCRU) { + if (cmValues[iCRU].size() == 0) { + continue; + } + cm[iCRU] = cmValues[iCRU][iTimeBin]; + if (writeDebug) { + cmD[iCRU] = cmDebug[iCRU][iTimeBin]; + } + const CRU cru(iCRU); + const StackID stackID{cru.sector(), cru.gemStack()}; + const auto index = stackID.getIndex(); + sumPosStack[index] += cm[iCRU].sumPos; + nPosStack[index] += (Mapper::PADSPERREGION[cru.region()] - cm[iCRU].nNeg); + nSaturationStack[index] += cm[iCRU].nSaturation; + } + + for (int iCRU = 0; iCRU < cmValues.size(); ++iCRU) { + if (cmValues[iCRU].size() == 0) { + continue; + } + const CRU cru(iCRU); + const StackID stackID{cru.sector(), cru.gemStack()}; + const auto index = stackID.getIndex(); + sumPosStackCRU[iCRU] = sumPosStack[index]; + sumPosStackCRUCorr[iCRU] = sumPosStack[index] - nPosStack[index] * cm[iCRU].cmValue; + nSaturationStackCRU[iCRU] = nSaturationStack[index]; + } + + pcstream << "cm" + << "iTF=" << iTF + << "iTimeBin=" << iTimeBin + << "cmInfo=" << cm + << "sumPosStack=" << sumPosStackCRU + << "sumPosStackCorr=" << sumPosStackCRUCorr + << "nSaturationStack=" << nSaturationStackCRU; + + if (writeDebug) { + pcstream << "cm" + << "cmDebug=" << cmD; + } + + pcstream << "cm" + << "\n"; + } + + // if (tOut) { + // tOut->Fill(); + // } + } + + pcstream.Close(); + if (fOut && tOut) { + tOut->SetEntries(nEntries); + fOut->cd(); + tOut->Write(); + tOut.reset(); + fOut->Close(); + } +} + +float CommonModeCorrection::decodeInjectedCMValue(float lower, float upper) +{ + // CRU row0 pad0 row1 pad1 + // 0 0 2 0 3 + // 1 20 1 20 3 + // 2 32 2 32 3 + // 3 51 1 51 3 + // 4 62 1 62 2 + // 5 84 1 84 4 + // 6 97 1 97 2 + // 7 116 2 115 5 + // 8 127 2 127 3 + // 9 142 0 142 4 + // + // CM Value encoding: + // Kanal 0 : Bit 11 ... 8 = 0x8. Bit 7..0 CM-Werte Bits 7...0 + // Kanal 1 : Bit 11.. 9 = "100". Bit 8 = CM Positive, Bits 6..0 = CM-Wert Bits 14..8 + const int ilower = floatToFixedSize(lower); + const int iupper = floatToFixedSize(upper); + if (!(ilower & 0x800) || !(iupper & 0x800)) { + LOGP(error, "Not a CM word: lower: {:012b} upper: {:012b}", ilower, iupper); + return 0; + } + const int fixedSizeCM = ((iupper & 0x7F) << 8) + (ilower & 0xFF); + const float floatCM = fixedSizeToFloat<8>(fixedSizeCM); + + // bit 8 of upper word is the sign 1 = positive + return (iupper & 0x100) ? floatCM : -floatCM; +} + +float CommonModeCorrection::getCalPadValue(const std::string calibName, int icru, int pad) const +{ + if (mPadMaps.find(calibName) == mPadMaps.end()) { + LOGP(error, "{} not set, cannot be used", calibName); + return 0; + } + const auto& calPad = mPadMaps.at(calibName); + const CRU cru(icru); + const int roc = cru.roc(); + const int padOffset = (cru.isIROC()) ? Mapper::GLOBALPADOFFSET[cru.region()] : Mapper::GLOBALPADOFFSET[cru.region()] - Mapper::GLOBALPADOFFSET[4]; + + const auto& calArray = calPad.getCalArray(roc); + + return calArray.getValue(padOffset + pad); +} + +bool CommonModeCorrection::padMapExists(const std::string& calibName) +{ + if (mPadMaps.find(calibName) == mPadMaps.end()) { + LOGP(error, "{} not in mPadMaps", calibName); + return false; + } + return true; +} + +void CommonModeCorrection::loadCalPad(std::string_view fileName, std::string_view nameInFile, std::string_view namePadMap) +{ + if (fileName.size() == 0) { + return; + } + + auto pads = o2::tpc::utils::readCalPads(fileName, nameInFile); + if ((pads.size() == 0) || (pads.at(0) == nullptr)) { + LOGP(error, "Could not load object {} from file {}", nameInFile, fileName); + return; + } + + if (namePadMap.size() == 0) { + namePadMap = nameInFile; + } + + mPadMaps[namePadMap.data()] = *pads[0]; +} diff --git a/Detectors/TPC/base/src/TPCBaseLinkDef.h b/Detectors/TPC/base/src/TPCBaseLinkDef.h index 33b6cf2c03392..60924db3953e2 100644 --- a/Detectors/TPC/base/src/TPCBaseLinkDef.h +++ b/Detectors/TPC/base/src/TPCBaseLinkDef.h @@ -66,6 +66,9 @@ #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::tpc::IonTailSettings> + ; #pragma link C++ class o2::tpc::FEEConfig + ; #pragma link C++ class o2::tpc::CRUConfig + ; +#pragma link C++ class o2::tpc::CommonModeCorrection + ; +#pragma link C++ class std::vector < o2::tpc::CommonModeCorrection::CMInfo> + ; +#pragma link C++ class std::vector < o2::tpc::CommonModeCorrection::CMDebug> + ; #pragma link C++ namespace o2::tpc::utils; #pragma link C++ function o2::tpc::utils::tokenize(const std::string_view, const std::string_view); From 9c26a3d5d71b60508d944c7917f117d8fcf0eb65 Mon Sep 17 00:00:00 2001 From: wiechula Date: Fri, 1 Nov 2024 16:27:20 +0100 Subject: [PATCH 3/8] Allow setting maxZ2X via config params --- .../SpacePoints/SpacePointsCalibConfParam.h | 39 ++++++++++--------- .../include/SpacePoints/TrackResiduals.h | 1 - .../SpacePoints/src/TrackResiduals.cxx | 20 +++++----- 3 files changed, 31 insertions(+), 29 deletions(-) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h index 2465cbf512d2b..9a4d7c1474287 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h @@ -29,35 +29,35 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelpermaxZ2X; mIsInitialized = true; + + if (doBinning) { + // initialize binning + initBinning(); + } + LOG(info) << "Initialization complete"; } @@ -182,10 +184,10 @@ void TrackResiduals::initBinning() } // // Z/X binning - mDZ2XI = mNZ2XBins / sMaxZ2X; + mDZ2XI = mNZ2XBins / mMaxZ2X; mDZ2X = 1.0f / mDZ2XI; // for uniform case only if (mUniformBins[VoxZ]) { - LOGF(info, "Z/X-binning is uniform with %i bins from 0 to %f", mNZ2XBins, sMaxZ2X); + LOGF(info, "Z/X-binning is uniform with %i bins from 0 to %f", mNZ2XBins, mMaxZ2X); for (int iz = 0; iz < mNZ2XBins; ++iz) { mZ2XBinsDH.push_back(.5f * mDZ2X); mZ2XBinsDI.push_back(mDZ2XI); @@ -265,7 +267,7 @@ int TrackResiduals::getRowID(float x) const bool TrackResiduals::findVoxelBin(int secID, float x, float y, float z, std::array& bvox) const { // Z/X bin - if (fabs(z / x) > sMaxZ2X) { + if (fabs(z / x) > mMaxZ2X) { return false; } int bz = getZ2XBinExact(secID < SECTORSPERSIDE ? z / x : -z / x); @@ -601,7 +603,7 @@ int TrackResiduals::validateVoxels(int iSec) resVox.flags |= Masked; } } // loop over Z - } // loop over Y/X + } // loop over Y/X mValidFracXBins[iSec][ix] = static_cast(cntValid) / (mNY2XBins * mNZ2XBins); LOGP(debug, "Sector {}: xBin {} has {} % of voxels valid. Total masked due to fit: {} ,and sigma: {}", iSec, ix, mValidFracXBins[iSec][ix] * 100., cntMaskedFit, cntMaskedSigma); From 196c6af84598a852289efa7bf7717a1367e4830b Mon Sep 17 00:00:00 2001 From: wiechula Date: Fri, 8 Nov 2024 22:19:40 +0100 Subject: [PATCH 4/8] Add functions for specific pad selections * edge pads * stack boundary rows * cross region --- Detectors/TPC/base/include/TPCBase/Mapper.h | 19 ++++++++--- Detectors/TPC/base/src/Mapper.cxx | 36 +++++++++++++++++++++ 2 files changed, 50 insertions(+), 5 deletions(-) diff --git a/Detectors/TPC/base/include/TPCBase/Mapper.h b/Detectors/TPC/base/include/TPCBase/Mapper.h index cee3d76db85b2..f2ff425675df6 100644 --- a/Detectors/TPC/base/include/TPCBase/Mapper.h +++ b/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 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 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); diff --git a/Detectors/TPC/base/src/Mapper.cxx b/Detectors/TPC/base/src/Mapper.cxx index 56ce283178da0..2796d488f014d 100644 --- a/Detectors/TPC/base/src/Mapper.cxx +++ b/Detectors/TPC/base/src/Mapper.cxx @@ -298,5 +298,41 @@ void Mapper::setTraceLengths(std::string_view inputFile, std::vector& len } } +bool Mapper::isEdgePad(int rowInSector, int padInRow) +{ + const auto& mapper = instance(); + return (padInRow == 0) || (padInRow == mapper.getNumberOfPadsInRowSector(rowInSector) - 1); +} + +bool Mapper::isFirstOrLastRowInStack(int rowInSector) +{ + if (rowInSector == 0 || rowInSector == PADROWS - 1) { + return true; + } + + const auto& mapper = instance(); + for (int i = 1; i < 4; ++i) { + if (rowInSector == ROWOFFSETSTACK[i] || rowInSector == ROWOFFSETSTACK[i] - 1) { + return true; + } + } + return false; +} + +bool Mapper::isBelowSpacerCross(int rowInSector, int padInRow) +{ + static std::vector ROWSBELOWCROSS{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + if (ROWSBELOWCROSS[rowInSector]) { + return true; + } + + const auto& mapper = instance(); + const auto padCenter = mapper.getNumberOfPadsInRowSector(rowInSector) / 2; + if (padInRow == padCenter || padInRow == padCenter - 1) { + return true; + } + return false; +} + } // namespace tpc } // namespace o2 From 945c08cf304acd38b473aa570e36703504de95bb Mon Sep 17 00:00:00 2001 From: wiechula Date: Wed, 16 Oct 2024 09:13:55 +0200 Subject: [PATCH 5/8] Add possibility to limit CM k-value in high coupling regions * Fix conversetion of fixed point to float --- .../TPC/calibration/macro/prepareCMFiles.C | 48 +++++++++++++++---- 1 file changed, 40 insertions(+), 8 deletions(-) diff --git a/Detectors/TPC/calibration/macro/prepareCMFiles.C b/Detectors/TPC/calibration/macro/prepareCMFiles.C index dc14bc61aa793..08880ccbe4862 100644 --- a/Detectors/TPC/calibration/macro/prepareCMFiles.C +++ b/Detectors/TPC/calibration/macro/prepareCMFiles.C @@ -28,7 +28,10 @@ using namespace o2::tpc::cru_calib_helpers; using namespace o2::tpc; -void prepareCMFiles(const std::string_view pulserFile, std::string outputDir = "./") +/// \param limitHighCouplingPads if > 0 limit pads in the high coupling region to this value +/// \param replaceHighCouplingPads if > 0 replace pads in the high coupling region by this value (take preceedence over limitHighCouplingPads) +/// \param maxValue if > 0 limit to this maximum value +void prepareCMFiles(const std::string_view pulserFile, std::string outputDir = "./", float limitHighCouplingPads = 0, float replaceHighCouplingPads = 0, float maxValue = 0) { constexpr uint32_t DataBits = 8; constexpr uint32_t FractionalBits = 6; @@ -115,6 +118,7 @@ void prepareCMFiles(const std::string_view pulserFile, std::string outputDir = " const int fecInPartition = fecInfo.getIndex() - partInfo.getSectorFECOffset(); const int dataWrapperID = fecInPartition >= fecOffset; const int globalLinkID = (fecInPartition % fecOffset) + dataWrapperID * 12; + const auto& padPos = mapper.padPos(globalPad); float pulserVal = rocPulserQtot.getValue(ipad); @@ -128,6 +132,20 @@ void prepareCMFiles(const std::string_view pulserFile, std::string outputDir = " pulserVal = MaxVal; } + if (replaceHighCouplingPads > 0) { + if (Mapper::isHighCouplingPad(padPos.getRow(), padPos.getPad())) { + pulserVal = replaceHighCouplingPads; + } + } else if (limitHighCouplingPads > 0) { + if (Mapper::isHighCouplingPad(padPos.getRow(), padPos.getPad())) { + pulserVal = std::min(pulserVal, limitHighCouplingPads); + } + } + + if (maxValue > 0) { + pulserVal = std::min(pulserVal, maxValue); + } + const int hwChannel = getHWChannel(sampa, sampaChannel, region % 2); // for debugging // printf("%4d %4d %4d %4d %4d: %u\n", cru.number(), globalLinkID, hwChannel, fecInfo.getSampaChip(), fecInfo.getSampaChannel(), getADCValue(pedestal)); @@ -143,22 +161,36 @@ void prepareCMFiles(const std::string_view pulserFile, std::string outputDir = " const bool onlyFilled = false; // ===| k-Values full float precision |=== - const auto outFileFloatTxt = (outputDir + "/commonMode_K_values_float.txt"); - const auto outFileFloatRoot = (outputDir + "/commonMode_K_values_float.root"); + string nameAdd; + if (replaceHighCouplingPads > 0) { + nameAdd = fmt::format(".replaceHC_{:.2}", replaceHighCouplingPads); + } else if (limitHighCouplingPads > 0) { + nameAdd = fmt::format(".limitHC_{:.2}", limitHighCouplingPads); + } + + if (maxValue > 0) { + nameAdd += fmt::format(".maxValue_{:.2}", maxValue); + } + + string outNameBase = "commonMode_K_values" + nameAdd; + string outNameInvBase = "commonMode_inv_K_values" + nameAdd; + + const auto outFileFloatTxt = (outputDir + "/" + outNameBase + "_float.txt"); + const auto outFileFloatRoot = (outputDir + "/" + outNameBase + "_float.root"); writeValues(outFileFloatTxt, commonModeKValuesFloat, onlyFilled); - getCalPad(outFileFloatTxt, outFileFloatRoot, "CMkValues"); + getCalPad<0>(outFileFloatTxt, outFileFloatRoot, "CMkValues"); // ===| k-Values limited precision 2I6F |=== - const auto outFileTxt = (outputDir + "/commonMode_K_values.txt"); - const auto outFileRoot = (outputDir + "/commonMode_K_values.root"); + const auto outFileTxt = (outputDir + "/" + outNameBase + ".txt"); + const auto outFileRoot = (outputDir + "/" + outNameBase + ".root"); writeValues(outFileTxt, commonModeKValues, onlyFilled); getCalPad(outFileTxt, outFileRoot, "CMkValues"); // ===| inverse k-Values limited precision 2I6F |=== - const auto outFileInvTxt = (outputDir + "/commonMode_inv_K_values.txt"); - const auto outFileInvRoot = (outputDir + "/commonMode_inv_K_values.root"); + const auto outFileInvTxt = (outputDir + "/" + outNameInvBase + ".txt"); + const auto outFileInvRoot = (outputDir + "/" + outNameInvBase + ".root"); writeValues(outFileInvTxt, commonModeInvKValues, onlyFilled); getCalPad(outFileInvTxt, outFileInvRoot, "InvCMkValues"); From 8359622b805f79dfd1e816768c1f4196291078e7 Mon Sep 17 00:00:00 2001 From: wiechula Date: Fri, 13 Dec 2024 12:54:03 +0100 Subject: [PATCH 6/8] Add PadFlags treatment, add protection --- Detectors/TPC/base/src/Painter.cxx | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/Detectors/TPC/base/src/Painter.cxx b/Detectors/TPC/base/src/Painter.cxx index 863547a666611..9f143d3fa45ce 100644 --- a/Detectors/TPC/base/src/Painter.cxx +++ b/Detectors/TPC/base/src/Painter.cxx @@ -334,9 +334,9 @@ TCanvas* painter::draw(const CalDet& calDet, int nbins1D, float xMin1D, float const GlobalPosition2D pos = mapper.getPadCentre(PadROCPos(roc, irow, ipad)); const int bin = hist2D->FindBin(pos.X(), pos.Y()); if (!hist2D->GetBinContent(bin)) { - hist2D->SetBinContent(bin, val); + hist2D->SetBinContent(bin, double(val)); } - hist1D->Fill(val); + hist1D->Fill(double(val)); } } } @@ -430,7 +430,7 @@ void painter::fillHistogram2D(TH2& h2D, const CalDet& calDet, Side side) const GlobalPosition2D pos = mapper.getPadCentre(PadROCPos(roc, irow, ipad)); const int bin = h2D.FindBin(pos.X(), pos.Y()); if (!h2D.GetBinContent(bin)) { - h2D.SetBinContent(bin, val); + h2D.SetBinContent(bin, double(val)); } } } @@ -454,7 +454,7 @@ void painter::fillHistogram2D(TH2& h2D, const CalArray& calArray) const GlobalPadNumber pad = mapper.getPadNumber(padSubset, position, irow, ipad); const auto val = calArray.getValue(pad); const int cpad = ipad - padsInRow / 2; - h2D.Fill(irow, cpad, val); + h2D.Fill(irow, cpad, double(val)); } } } @@ -523,6 +523,17 @@ std::enable_if_t::value, bool> hasData(const CalArray& ca return cal.getSum() > T{0}; } +template +std::enable_if_t::value, bool> hasData(const CalArray& cal) +{ + for (const auto v : cal.getData()) { + if (int(v) > 0) { + return true; + } + } + return false; +} + template std::vector painter::makeSummaryCanvases(const CalDet& calDet, int nbins1D, float xMin1D, float xMax1D, bool onlyFilled, std::vector* outputCanvases) { @@ -589,7 +600,7 @@ std::vector painter::makeSummaryCanvases(const CalDet& calDet, int // ===| 1D histogram |=== auto h1D = new TH1F(fmt::format("h1_{}_{:02d}", calName, iroc).data(), fmt::format("{} distribution ROC {:02d} ({});ADC value", calName, iroc, getROCTitle(iroc)).data(), nbins1D, xMin1D, xMax1D); for (const auto& val : roc.getData()) { - h1D->Fill(val); + h1D->Fill(double(val)); } // ===| 2D histogram |=== @@ -1342,6 +1353,9 @@ void painter::adjustPalette(TH1* h, float x2ndc, float tickLength) gPad->Modified(); gPad->Update(); auto palette = (TPaletteAxis*)h->GetListOfFunctions()->FindObject("palette"); + if (!palette) { + return; + } palette->SetX2NDC(x2ndc); auto ax = h->GetZaxis(); ax->SetTickLength(tickLength); @@ -1425,6 +1439,12 @@ template TCanvas* painter::draw(const CalArray& calArray); template TH2* painter::getHistogram2D(const CalDet& calDet, Side side); template TH2* painter::getHistogram2D(const CalArray& calArray); +template TCanvas* painter::draw(const CalDet& calDet, int, float, float, TCanvas*); +template std::vector painter::makeSummaryCanvases(const CalDet& calDet, int, float, float, bool, std::vector*); +template TCanvas* painter::draw(const CalArray& calArray); +template TH2* painter::getHistogram2D(const CalDet& calDet, Side side); +template TH2* painter::getHistogram2D(const CalArray& calArray); + template TCanvas* painter::draw(const CalDet& calDet, int, float, float, TCanvas*); template std::vector painter::makeSummaryCanvases(const CalDet& calDet, int, float, float, bool, std::vector*); template TCanvas* painter::draw(const CalArray& calArray); From 952031360672aafcc384a3d2e7c542652bec9773 Mon Sep 17 00:00:00 2001 From: wiechula Date: Fri, 13 Dec 2024 12:55:17 +0100 Subject: [PATCH 7/8] Add possibility to store canvases in single PDF --- Detectors/TPC/base/src/Utils.cxx | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/Detectors/TPC/base/src/Utils.cxx b/Detectors/TPC/base/src/Utils.cxx index 8879e8ab342d2..d8a420ea4f03c 100644 --- a/Detectors/TPC/base/src/Utils.cxx +++ b/Detectors/TPC/base/src/Utils.cxx @@ -133,7 +133,7 @@ void utils::addFECInfo() h->SetTitle(title.data()); } -void utils::saveCanvases(TObjArray& arr, std::string_view outDir, std::string_view types, std::string_view rootFileName, std::string nameAdd) +void utils::saveCanvases(TObjArray& arr, std::string_view outDir, std::string_view types, std::string_view singleOutFileName, std::string nameAdd) { if (types.size()) { for (auto c : arr) { @@ -141,21 +141,39 @@ void utils::saveCanvases(TObjArray& arr, std::string_view outDir, std::string_vi } } - if (rootFileName.size()) { - std::unique_ptr outFile(TFile::Open(fmt::format("{}/{}", outDir, rootFileName).data(), "recreate")); - arr.Write(arr.GetName(), TObject::kSingleKey); - outFile->Close(); + if (singleOutFileName.size()) { + const auto outFileNames = o2::utils::Str::tokenize(singleOutFileName.data(), ','); + for (const auto& outFileName : outFileNames) { + auto fileName = fmt::format("{}/{}", outDir, outFileName); + if (o2::utils::Str::endsWith(outFileName, ".root")) { + std::unique_ptr outFile(TFile::Open(fileName.data(), "recreate")); + arr.Write(arr.GetName(), TObject::kSingleKey); + outFile->Close(); + } else if (o2::utils::Str::endsWith(outFileName, ".pdf")) { + const auto nCanv = arr.GetEntries(); + for (int i = 0; i < nCanv; ++i) { + auto fileName2 = fileName; + if (i == 0) { + fileName2 += "("; + } else if (i == nCanv - 1) { + fileName2 += ")"; + } + auto c = static_cast(arr.UncheckedAt(i)); + c->Print(fileName2.data(), fmt::format("Title:{}", c->GetTitle()).data()); + } + } + } } } -void utils::saveCanvases(std::vector& canvases, std::string_view outDir, std::string_view types, std::string_view rootFileName, std::string nameAdd) +void utils::saveCanvases(std::vector& canvases, std::string_view outDir, std::string_view types, std::string_view singleOutFileName, std::string nameAdd) { TObjArray arr; for (auto c : canvases) { arr.Add(c); } - saveCanvases(arr, outDir, types, rootFileName, nameAdd); + saveCanvases(arr, outDir, types, singleOutFileName, nameAdd); } void utils::saveCanvas(TCanvas& c, std::string_view outDir, std::string_view types, std::string nameAdd) From 95aa1888479d7a9d88824648a557d705f59ca6d2 Mon Sep 17 00:00:00 2001 From: wiechula Date: Thu, 19 Dec 2024 16:03:50 +0100 Subject: [PATCH 8/8] Add treatment of clustom dE/dx file and disabling dE/dx input --- Detectors/TPC/workflow/src/CalibdEdxSpec.cxx | 36 ++++++++++++++++--- .../TPC/workflow/src/CalibratordEdxSpec.cxx | 34 +++++++++++++++--- 2 files changed, 61 insertions(+), 9 deletions(-) diff --git a/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx b/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx index a32a4a1bb3089..97b69156a2a6d 100644 --- a/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx +++ b/Detectors/TPC/workflow/src/CalibdEdxSpec.cxx @@ -18,17 +18,18 @@ // o2 includes #include "CCDB/CcdbApi.h" #include "CCDB/CcdbObjectInfo.h" -#include "CommonUtils/NameConf.h" +// #include "CommonUtils/NameConf.h" #include "DataFormatsTPC/TrackTPC.h" -#include "DataFormatsParameters/GRPObject.h" +// #include "DataFormatsParameters/GRPObject.h" #include "DetectorsCalibration/Utils.h" #include "Framework/Task.h" #include "Framework/DataProcessorSpec.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/CCDBParamSpec.h" +#include "GPUO2InterfaceConfigurableParam.h" #include "TPCCalibration/CalibdEdx.h" #include "TPCWorkflow/ProcessingHelpers.h" -#include "TPCBase/CDBInterface.h" +#include "TPCBase/CDBTypes.h" #include "TPCBase/Utils.h" #include "DetectorsBase/GRPGeomHelper.h" @@ -68,6 +69,29 @@ class CalibdEdxDevice : public Task mCalib->set2DFitThreshold(minEntries2D); mCalib->setElectronCut(fitThreshold, fitPasses, fitThresholdLowFactor); mCalib->setMaterialType(mMatType); + + mCustomdEdxFileName = o2::gpu::GPUConfigurableParamGPUSettingsO2::Instance().dEdxCorrFile; + mDisableTimeGain = o2::gpu::GPUConfigurableParamGPUSettingsO2::Instance().dEdxDisableResidualGain; + + if (mDisableTimeGain) { + LOGP(info, "TimeGain correction was disabled via GPU_global.dEdxDisableResidualGain=1"); + } + + if (!mDisableTimeGain && !mCustomdEdxFileName.empty()) { + std::unique_ptr fdEdxCustom(TFile::Open(mCustomdEdxFileName.data())); + if (!fdEdxCustom || !fdEdxCustom->IsOpen() || fdEdxCustom->IsZombie()) { + LOGP(error, "Could not open custom TimeGain file {}", mCustomdEdxFileName); + } else { + const auto timeGain = fdEdxCustom->Get("CalibdEdxCorrection"); + if (!timeGain) { + LOGP(error, "Could not load 'CalibdEdxCorrection' from file {}", mCustomdEdxFileName); + } else { + const auto meanParamTot = timeGain->getMeanParams(ChargeType::Tot); + LOGP(info, "Loaded custom TimeGain from file {} with {} dimensions and mean qTot Params {}", mCustomdEdxFileName, timeGain->getDims(), utils::elementsToString(meanParamTot)); + mCalib->setCalibrationInput(*timeGain); + } + } + } } void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final @@ -75,7 +99,7 @@ class CalibdEdxDevice : public Task if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { return; } - if (matcher == ConcreteDataMatcher("TPC", "TIMEGAIN", 0)) { + if ((mDisableTimeGain == 0) && mCustomdEdxFileName.empty() && (matcher == ConcreteDataMatcher("TPC", "TIMEGAIN", 0))) { mCalib->setCalibrationInput(*(o2::tpc::CalibdEdxCorrection*)obj); const auto meanParamTot = mCalib->getCalibrationInput().getMeanParams(ChargeType::Tot); LOGP(info, "Updating TimeGain with {} dimensions and mean qTot Params {}", mCalib->getCalibrationInput().getDims(), utils::elementsToString(meanParamTot)); @@ -143,7 +167,9 @@ class CalibdEdxDevice : public Task uint64_t mRunNumber{0}; ///< processed run number uint64_t mTimeStampStart{0}; ///< time stamp for first TF for CCDB output std::unique_ptr mCalib; - bool mMakeGaussianFits{true}; ///< make gaussian fits or take the mean + bool mMakeGaussianFits{true}; ///< make gaussian fits or take the mean + bool mDisableTimeGain{false}; ///< if time gain is disabled via GPU_global.dEdxDisableResidualGain=1 + std::string mCustomdEdxFileName{}; ///< name of the custom dE/dx file configured via GPU_global.dEdxCorrFile }; DataProcessorSpec getCalibdEdxSpec(const o2::base::Propagator::MatCorrType matType) diff --git a/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx b/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx index 6e477084d992c..ce45356aa28c8 100644 --- a/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx +++ b/Detectors/TPC/workflow/src/CalibratordEdxSpec.cxx @@ -21,18 +21,19 @@ // o2 includes #include "CCDB/CcdbApi.h" #include "CCDB/CcdbObjectInfo.h" -#include "CommonUtils/NameConf.h" +// #include "CommonUtils/NameConf.h" #include "DataFormatsTPC/TrackTPC.h" -#include "DataFormatsParameters/GRPObject.h" +// #include "DataFormatsParameters/GRPObject.h" #include "DetectorsCalibration/Utils.h" #include "Framework/Task.h" #include "Framework/DataProcessorSpec.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/CCDBParamSpec.h" +#include "GPUO2InterfaceConfigurableParam.h" #include "TPCCalibration/CalibratordEdx.h" #include "TPCWorkflow/ProcessingHelpers.h" #include "DetectorsBase/GRPGeomHelper.h" -#include "TPCBase/CDBInterface.h" +#include "TPCBase/CDBTypes.h" #include "TPCBase/Utils.h" using namespace o2::framework; @@ -85,6 +86,29 @@ class CalibratordEdxDevice : public Task mCalibrator->setTrackDebug(trackDebug); mCalibrator->setMakeGaussianFits(makeGaussianFits); + mCustomdEdxFileName = o2::gpu::GPUConfigurableParamGPUSettingsO2::Instance().dEdxCorrFile; + mDisableTimeGain = o2::gpu::GPUConfigurableParamGPUSettingsO2::Instance().dEdxDisableResidualGain; + + if (mDisableTimeGain) { + LOGP(info, "TimeGain correction was disabled via GPU_global.dEdxDisableResidualGain=1"); + } + + if (!mDisableTimeGain && !mCustomdEdxFileName.empty()) { + std::unique_ptr fdEdxCustom(TFile::Open(mCustomdEdxFileName.data())); + if (!fdEdxCustom || !fdEdxCustom->IsOpen() || fdEdxCustom->IsZombie()) { + LOGP(error, "Could not open custom TimeGain file {}", mCustomdEdxFileName); + } else { + const auto timeGain = fdEdxCustom->Get("CalibdEdxCorrection"); + if (!timeGain) { + LOGP(error, "Could not load 'CalibdEdxCorrection' from file {}", mCustomdEdxFileName); + } else { + mTimeGain = *timeGain; + const auto meanParamTot = mTimeGain.getMeanParams(ChargeType::Tot); + LOGP(info, "Loaded custom TimeGain from file {} with {} dimensions and mean qTot Params {}", mCustomdEdxFileName, mTimeGain.getDims(), utils::elementsToString(meanParamTot)); + } + } + } + if (dumpData) { const auto dumpDataName = ic.options().get("file-dump-name"); mCalibrator->enableDebugOutput(dumpDataName); @@ -96,7 +120,7 @@ class CalibratordEdxDevice : public Task if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { return; } - if (matcher == ConcreteDataMatcher("TPC", "TIMEGAIN", 0)) { + if ((mDisableTimeGain == 0) && mCustomdEdxFileName.empty() && (matcher == ConcreteDataMatcher("TPC", "TIMEGAIN", 0))) { mTimeGain = *(o2::tpc::CalibdEdxCorrection*)obj; const auto meanParamTot = mTimeGain.getMeanParams(ChargeType::Tot); LOGP(info, "Updating TimeGain with {} dimensions and mean qTot Params {}", mTimeGain.getDims(), utils::elementsToString(meanParamTot)); @@ -181,6 +205,8 @@ class CalibratordEdxDevice : public Task uint32_t mRunNumber{0}; ///< processed run number long mCalibIntervalExtensionMS{0}; ///< Extension of the calibration interval end in ms o2::tpc::CalibdEdxCorrection mTimeGain{}; ///< currently valid TimeGain + bool mDisableTimeGain{false}; ///< if time gain is disabled via GPU_global.dEdxDisableResidualGain=1 + std::string mCustomdEdxFileName{}; ///< name of the custom dE/dx file configured via GPU_global.dEdxCorrFile }; DataProcessorSpec getCalibratordEdxSpec(const o2::base::Propagator::MatCorrType matType)