Skip to content

Commit

Permalink
Fix and extension
Browse files Browse the repository at this point in the history
  • Loading branch information
@wiechula
wiechula committed Sep 27, 2024
1 parent ed02730 commit 9890a03
Show file tree
Hide file tree
Showing 2 changed files with 122 additions and 32 deletions.
32 changes: 26 additions & 6 deletions Detectors/TPC/base/include/TPCBase/CommonModeCorrection.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,9 +48,15 @@ class CommonModeCorrection

struct CMInfo {
float cmValue{};
float cmValueStd{};
int nPadsUsed{};
};

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
Expand All @@ -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<const float> values, gsl::span<const float> cmKValues, gsl::span<const float> pedestals) const;
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)); }
Expand All @@ -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; }

Expand Down Expand Up @@ -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<Digit>& digits, std::vector<std::vector<CMInfo>>& cmValues, bool negativeOnly = false) const;
int correctDigits(std::vector<Digit>& digits, std::vector<std::vector<CMInfo>>& cmValues, bool negativeOnly = 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);
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; }
Expand All @@ -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)

Expand Down
122 changes: 96 additions & 26 deletions Detectors/TPC/base/src/CommonModeCorrection.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
/// \file CommonModeCorrection.cxx
/// \brief Calculate the common mode correction factor

#include <random>
// #include <random>
#include <algorithm>
#include <thread>
#include <mutex>
Expand All @@ -28,7 +28,7 @@

using namespace o2::tpc;
using namespace o2::tpc::cru_calib_helpers;
CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span<const float> values, gsl::span<const float> cmKValues, gsl::span<const float> pedestals) const
CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span<const float> values, gsl::span<const float> cmKValues, gsl::span<const float> pedestals, CMDebug* cmDebug) const
{
if (values.size() == 0) {
return CMInfo{};
Expand All @@ -41,6 +41,12 @@ CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span<const
static math_utils::RandomRing random(math_utils::RandomRing<>::RandomType::Flat);
std::vector<float> 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];
Expand All @@ -51,6 +57,12 @@ CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span<const
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) {
Expand All @@ -61,23 +73,41 @@ CommonModeCorrection::CMInfo CommonModeCorrection::getCommonMode(gsl::span<const
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;
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)
Expand Down Expand Up @@ -126,61 +156,90 @@ CommonModeCorrection::CMdata CommonModeCorrection::collectCMdata(const std::vect
return data;
}

int CommonModeCorrection::correctDigits(std::vector<Digit>& digits, std::vector<std::vector<CMInfo>>& cmValues, bool negativeOnly) const
int CommonModeCorrection::correctDigits(std::vector<Digit>& digits, std::vector<std::vector<CMInfo>>& cmValues, bool negativeOnly, std::vector<std::vector<CMDebug>>* 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");
Expand Down Expand Up @@ -217,7 +276,12 @@ void CommonModeCorrection::correctDigits(std::string_view digiFileIn, Long64_t m
LOGP(info, "Processing entry {}/{}", iTF + 1, nEntries);

std::vector<std::vector<CMInfo>> cmValues; // CRU * timeBin
std::vector<std::vector<CMDebug>> cmDebug; // CRU * timeBin

cmValues.resize(CRU::MaxCRU);
if (writeDebug) {
cmDebug.resize(CRU::MaxCRU);
}
int maxTimeBin = -1;

auto worker = [&](int iTread) {
Expand All @@ -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};
Expand All @@ -247,22 +311,28 @@ 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) {
pcstream << "cm"
<< "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();
Expand Down

0 comments on commit 9890a03

Please sign in to comment.