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 24, 2024
1 parent 2bbaa14 commit 0e4888c
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 27 deletions.
22 changes: 16 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 Down Expand Up @@ -103,9 +109,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) 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);

void limitKFactorPrecision(bool limit = true) { mLimitKFactor = limit; }
void limitPedestalPrecision(bool limit = true) { mLimitPedestal = limit; }
Expand All @@ -117,14 +123,18 @@ 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
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
90 changes: 69 additions & 21 deletions Detectors/TPC/base/src/CommonModeCorrection.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 @@ -61,23 +67,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,7 +150,7 @@ 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) const
{
// calculation common mode values
int maxTimeBin = -1;
Expand All @@ -138,20 +162,29 @@ int CommonModeCorrection::correctDigits(std::vector<Digit>& digits, std::vector<

for (size_t iDigit = 0; iDigit < digits.size(); ++iDigit) {
auto& digit = digits[iDigit];
if (((lastCRU > 0) && (digit.getCRU() != lastCRU)) || ((lastTimeBin > 0) && (digit.getTimeStamp() != lastTimeBin))) {
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 (mArtificialCM > 0) {
charge = std::clamp(charge + mArtificialCM, 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();
Expand All @@ -161,16 +194,20 @@ int CommonModeCorrection::correctDigits(std::vector<Digit>& digits, std::vector<
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 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);
Expand All @@ -180,7 +217,7 @@ int CommonModeCorrection::correctDigits(std::vector<Digit>& digits, std::vector<
return maxTimeBin;
}

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)
{

TChain* tree = o2::tpc::utils::buildChain(fmt::format("ls {}", digiFileIn), "o2sim", "o2sim");
Expand Down Expand Up @@ -217,7 +254,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 +269,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);
{
static std::mutex maxMutex;
std::lock_guard lock{maxMutex};
Expand All @@ -247,22 +289,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 0e4888c

Please sign in to comment.