CalibPulser.cxx

Go to the documentation of this file.

// 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.
 
 
#include <string>
#include <vector>
#include <algorithm>
 
#include "TObjArray.h"
#include "TFile.h"
 
#include "TPCBase/ROC.h"
#include "MathUtils/fit.h"
#include "TPCCalibration/CalibPulser.h"
#include "TPCCalibration/CalibPulserParam.h"
 
using namespace o2::tpc;
using o2::math_utils::getStatisticsData;
using o2::math_utils::StatisticsData;
 
CalibPulser::CalibPulser(PadSubset padSubset)
  : mNbinsT0{200},
    mXminT0{-2},
    mXmaxT0{2},
    mNbinsQtot{200},
    mXminQtot{50},
    mXmaxQtot{700},
    mNbinsWidth{100},
    mXminWidth{0.1},
    mXmaxWidth{5.1},
    mFirstTimeBin{10},
    mLastTimeBin{490},
    mADCMin{5},
    mADCMax{1023},
    mNumberOfADCs{mADCMax - mADCMin + 1},
    mPeakIntMinus{2},
    mPeakIntPlus{2},
    mMinimumQtot{20},
    mMinimumQmax{10},
    mMaxTimeBinRange{5},
    mPedestal{nullptr},
    mNoise{nullptr},
    mPulserData{},
    mT0Histograms{},
    mWidthHistograms{},
    mQtotHistograms{}
{
  mT0Histograms.resize(ROC::MaxROC);
  mWidthHistograms.resize(ROC::MaxROC);
  mQtotHistograms.resize(ROC::MaxROC);
 
  // TODO: until automatic T0 calibration is done we use the same time range
  //       as for the time bin selection
  mXminT0 = mFirstTimeBin;
  mXmaxT0 = mLastTimeBin;
 
  mCalDets["T0"] = CalPad("T0");
  mCalDets["Width"] = CalPad("Width");
  mCalDets["Qtot"] = CalPad("Qtot");
}
 
//______________________________________________________________________________
void CalibPulser::init()
{
  const auto& param = CalibPulserParam::Instance();
 
  mNbinsT0 = param.NbinsT0;
  mXminT0 = param.XminT0;
  mXmaxT0 = param.XmaxT0;
  mNbinsQtot = param.NbinsQtot;
  mXminQtot = param.XminQtot;
  mXmaxQtot = param.XmaxQtot;
  mNbinsWidth = param.NbinsWidth;
  mXminWidth = param.XminWidth;
  mXmaxWidth = param.XmaxWidth;
  mFirstTimeBin = param.FirstTimeBin;
  mLastTimeBin = param.LastTimeBin;
  mADCMin = param.ADCMin;
  mADCMax = param.ADCMax;
  mNumberOfADCs = mADCMax - mADCMin + 1;
  mPeakIntMinus = param.PeakIntMinus;
  mPeakIntPlus = param.PeakIntPlus;
  mMinimumQtot = param.MinimumQtot;
  mMinimumQmax = param.MinimumQmax;
  mMaxTimeBinRange = param.MaxTimeBinRange;
}
 
//______________________________________________________________________________
Int_t CalibPulser::updateROC(const Int_t roc, const Int_t row, const Int_t pad,
                             const Int_t timeBin, Float_t signal)
{
  // ===| range checks |========================================================
  if (timeBin < mFirstTimeBin || timeBin > mLastTimeBin) {
    return 0;
  }
 
  // ---| pedestal subtraction |---
  if (mPedestal) {
    signal -= mPedestal->getValue(ROC(roc), row, pad);
  }
 
  if (signal < mADCMin || signal > mADCMax) {
    return 0;
  }
 
  // ===| temporary calibration data |==========================================
  const PadROCPos padROCPos(roc, row, pad);
  VectorType& adcData = mPulserData[padROCPos];
 
  if (!adcData.size()) {
    // accept first and last time bin, so difference +1
    adcData.resize(mLastTimeBin - mFirstTimeBin + 1);
    // printf("new pad pos\n");
  }
 
  adcData[timeBin - mFirstTimeBin] = signal;
  // printf("%2d, %3d, %3d, %3d: %.2f\n", roc, row, pad, timeBin, signal);
 
  // ---| entries per time bin |---
  auto& timeBinEntries = mTimeBinEntries[roc];
  if (!timeBinEntries.size()) {
    timeBinEntries.resize(mLastTimeBin - mFirstTimeBin + 1);
  }
  ++timeBinEntries[timeBin - mFirstTimeBin];
 
  return 1;
}
 
//______________________________________________________________________________
void CalibPulser::endReader()
{
  // loop over all signals of all pads filled for the present raw reader
  // for (auto& keyValue : mPulserData) {
  // const auto& padROCPos = keyValue.first;
  // const auto& adcData = keyValue.second;
  //
  auto timeRangeROCs = getTimeRangeROCs();
 
  for (const auto& [padROCPos, adcData] : mPulserData) {
    const auto currentChannel = mMapper.getPadNumberInROC(padROCPos);
 
    // std::cout << (int)padROCPos.getROC().getRoc() << ", " << padROCPos.getRow() << ", " << padROCPos.getPad() << std::endl;
    // for (auto& val: adcData) std::cout << val << ", ";
    // std::cout << std::endl;
 
    const auto& elemPair = timeRangeROCs[int(padROCPos.getROC())];
    const auto data = processPadData(padROCPos, adcData, elemPair);
    // std::cout << data.mT0+mFirstTimeBin << " : " << data.mQtot << " : " << data.mWidth << "\n";
 
    // fill histograms
    if (data.isValid()) {
      getHistoT0(padROCPos.getROC(), true)->Fill(data.mT0 + mFirstTimeBin, currentChannel);
      getHistoQtot(padROCPos.getROC(), true)->Fill(data.mQtot, currentChannel);
      getHistoSigma(padROCPos.getROC(), true)->Fill(data.mWidth, currentChannel);
    }
  }
 
  // reset the adc data to free space
  mPulserData.clear();
  for (auto& v : mTimeBinEntries) {
    std::fill(v.begin(), v.end(), 0);
  }
}
 
//______________________________________________________________________________
TH2S* CalibPulser::getHistogram(ROC roc, CalibPulser::PtrVectorType& rocVector,
                                int nbins, float min, float max,
                                std::string_view type, bool create /*=kFALSE*/)
{
  TH2S* vec = rocVector[roc].get();
  if (vec || !create) {
    return vec;
  }
 
  const size_t nChannels = mMapper.getNumberOfPads(roc);
  rocVector[roc] = std::make_unique<TH2S>(Form("hCalib%s%02d", type.data(), roc.getRoc()),
                                          Form("%s calibration histogram ROC %02d", type.data(), roc.getRoc()),
                                          nbins, min, max,
                                          nChannels, 0, nChannels);
  // printf("new histogram %s for ROC %2d\n", type.data(), int(roc));
  return rocVector[roc].get();
}
 
//______________________________________________________________________________
CalibPulser::PulserData CalibPulser::processPadData(const PadROCPos& padROCPos, const CalibPulser::VectorType& adcData, const CalibPulser::ElemPair& range)
{
  // data to return
  PulserData data;
 
  // find time bin with maximum number of entreis
 
  const auto vectorSize = adcData.size();
  // const auto maxElement = std::max_element(std::begin(adcData), std::end(adcData));
  const auto maxElement = std::max_element(std::begin(adcData) + range.first, std::begin(adcData) + range.last);
 
  if (*maxElement < mMinimumQmax) {
    return data;
  }
 
  const auto maxPosition = std::distance(std::begin(adcData), maxElement);
 
  double weightedSum = 0.;
  double weightedSum2 = 0.;
  double chargeSum = 0.;
  for (int t = maxPosition - mPeakIntMinus; t <= maxPosition + mPeakIntPlus; ++t) {
    const auto signal = adcData[t];
    // check time bounds
    if (t < 0 || t >= vectorSize) {
      continue;
    }
    weightedSum += signal * (t + 0.5); // +0.5 to get the center of the time bin
    weightedSum2 += signal * (t + 0.5) * (t + 0.5);
    chargeSum += signal;
  }
 
  if (chargeSum > mMinimumQtot) {
    weightedSum /= chargeSum;
    weightedSum2 = std::sqrt(std::abs(weightedSum2 / chargeSum - weightedSum * weightedSum));
    // L1 phase correction?
    data.mT0 = weightedSum;
    data.mWidth = weightedSum2;
    data.mQtot = chargeSum;
  }
 
  return data;
}
 
//______________________________________________________________________________
std::array<CalibPulser::ElemPair, ROC::MaxROC> CalibPulser::getTimeRangeROCs()
{
  std::array<ElemPair, ROC::MaxROC> maxTimeBin;
 
  for (size_t iROC = 0; iROC < mTimeBinEntries.size(); ++iROC) {
    const auto& timeBinEntries = mTimeBinEntries[iROC];
    const auto size = timeBinEntries.size();
    const auto maxElement = std::max_element(std::begin(timeBinEntries), std::end(timeBinEntries));
    const auto maxPosition = std::distance(std::begin(timeBinEntries), maxElement);
 
    maxTimeBin[iROC].first = (maxPosition > mMaxTimeBinRange) ? maxPosition - mMaxTimeBinRange : 0;
    maxTimeBin[iROC].last = std::min(size_t(maxPosition + mMaxTimeBinRange), size - 1);
  }
 
  return maxTimeBin;
}
 
//______________________________________________________________________________
void CalibPulser::resetData()
{
  mPulserData.clear();
  for (auto& v : mTimeBinEntries) {
    std::fill(v.begin(), v.end(), 0);
  }
 
  std::vector<PtrVectorType*> v{&mT0Histograms, &mWidthHistograms, &mQtotHistograms};
  for (auto histArray : v) {
    for (auto& histPtr : *histArray) {
      auto ptr = histPtr.get();
      if (ptr) {
        ptr->Reset();
      }
    }
  }
}
 
//______________________________________________________________________________
void CalibPulser::analyse()
{
  for (ROC roc; !roc.looped(); ++roc) {
    auto histT0 = mT0Histograms.at(roc).get();
    auto histWidth = mWidthHistograms.at(roc).get();
    auto histQtot = mQtotHistograms.at(roc).get();
 
    if (!histT0 || !histWidth || !histQtot) {
      continue;
    }
 
    // array pointer
    const auto arrT0 = histT0->GetArray();
    const auto arrWidth = histWidth->GetArray();
    const auto arrQtot = histQtot->GetArray();
 
    const auto numberOfPads = mMapper.getNumberOfPads(roc);
    for (auto iChannel = 0; iChannel < numberOfPads; ++iChannel) {
      const int offsetT0 = (mNbinsT0 + 2) * (iChannel + 1) + 1;
      const int offsetQtot = (mNbinsQtot + 2) * (iChannel + 1) + 1;
      const int offsetWidth = (mNbinsWidth + 2) * (iChannel + 1) + 1;
 
      StatisticsData dataT0 = getStatisticsData(arrT0 + offsetT0, mNbinsT0, mXminT0, mXmaxT0);
      StatisticsData dataWidth = getStatisticsData(arrWidth + offsetWidth, mNbinsWidth, mXminWidth, mXmaxWidth);
      StatisticsData dataQtot = getStatisticsData(arrQtot + offsetQtot, mNbinsQtot, mXminQtot, mXmaxQtot);
 
      mCalDets["T0"].getCalArray(roc).setValue(iChannel, dataT0.mCOG);
      mCalDets["Width"].getCalArray(roc).setValue(iChannel, dataWidth.mCOG);
      mCalDets["Qtot"].getCalArray(roc).setValue(iChannel, dataQtot.mCOG);
    }
  }
}
 
//______________________________________________________________________________
void CalibPulser::dumpToFile(const std::string filename, uint32_t type /* = 0*/)
{
  auto f = std::unique_ptr<TFile>(TFile::Open(filename.c_str(), "recreate"));
  if (type == 0) {
    f->WriteObject(&mCalDets["T0"], "T0");
    f->WriteObject(&mCalDets["Width"], "Width");
    f->WriteObject(&mCalDets["Qtot"], "Qtot");
 
    if (mDebugLevel) {
      printf("dump debug info\n");
      // temporary arrays for writing the objects
      TObjArray vT0;
      for (auto& val : mT0Histograms) {
        vT0.Add(val.get());
      }
      TObjArray vWidth;
      for (auto& val : mWidthHistograms) {
        vWidth.Add(val.get());
      }
      TObjArray vQtot;
      for (auto& val : mQtotHistograms) {
        vQtot.Add(val.get());
      }
 
      vT0.Write("T0Histograms", TObject::kSingleKey);
      vWidth.Write("WidthHistograms", TObject::kSingleKey);
      vQtot.Write("QtotHistograms", TObject::kSingleKey);
    }
  } else if (type == 1) {
    f->WriteObject(this, "CalibPulser");
  }
 
  f->Close();
}