IDCCCDBHelper.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 "TPCCalibration/IDCCCDBHelper.h"
#include "TPCCalibration/IDCDrawHelper.h"
#include "TPCCalibration/IDCGroupHelperSector.h"
#include "TPCCalibration/IDCContainer.h"
#include "TPCCalibration/RobustAverage.h"
 
#include "TPCBase/CalDet.h"
#include "TPCBase/Mapper.h"
#include "CommonUtils/TreeStreamRedirector.h"
#include "TPCBase/Painter.h"
 
#include "TStyle.h"
#include "TLine.h"
#include "TCanvas.h"
#include "TH2Poly.h"
#include "TGraph.h"
#include "TText.h"
#include "TPCCalibration/IDCFactorization.h"
#include <map>
 
#include "TProfile.h"
 
template <typename DataT>
unsigned int o2::tpc::IDCCCDBHelper<DataT>::getNIntegrationIntervalsIDCDelta(const o2::tpc::Side side) const
{
  return (mIDCDelta[side] && mHelperSector[side]) ? mIDCDelta[side]->getNIDCs() / (mHelperSector[side]->getNIDCsPerSector() * SECTORSPERSIDE) : 0;
}
 
template <typename DataT>
unsigned int o2::tpc::IDCCCDBHelper<DataT>::getNIntegrationIntervalsIDCOne(const o2::tpc::Side side) const
{
  return mIDCOne[side] ? mIDCOne[side]->getNIDCs() : 0;
}
 
template <typename DataT>
float o2::tpc::IDCCCDBHelper<DataT>::getIDCZeroVal(const unsigned int sector, const unsigned int region, unsigned int urow, unsigned int upad) const
{
  const auto side = Sector(sector).side();
  return !mIDCZero[side] ? -1 : (mIDCZero[side]->getNIDC0() == Mapper::getNumberOfPadsPerSide()) ? mIDCZero[side]->getValueIDCZero(getUngroupedIndexGlobal(sector, region, urow, upad, 0))
                                                                                                 : mIDCZero[side]->getValueIDCZero(mHelperSector[side]->getIndexUngrouped(sector, region, urow, upad, 0));
}
 
template <typename DataT>
float o2::tpc::IDCCCDBHelper<DataT>::getIDCDeltaVal(const unsigned int sector, const unsigned int region, unsigned int urow, unsigned int upad, unsigned int integrationInterval) const
{
  const auto side = Sector(sector).side();
  return (!mIDCDelta[side] || !mHelperSector[side]) ? -1 : mIDCDelta[side]->getValue(mHelperSector[side]->getIndexUngrouped(sector, region, urow, upad, integrationInterval));
}
 
template <typename DataT>
float o2::tpc::IDCCCDBHelper<DataT>::getIDCOneVal(const o2::tpc::Side side, const unsigned int integrationInterval) const
{
  return !mIDCOne[side] ? -1 : mIDCOne[side]->getValueIDCOne(integrationInterval);
}
 
template <typename DataT>
float o2::tpc::IDCCCDBHelper<DataT>::getIDCVal(const unsigned int sector, const unsigned int region, unsigned int urow, unsigned int upad, unsigned int integrationInterval) const
{
  return (getIDCDeltaVal(sector, region, urow, upad, integrationInterval) + 1.f) * getIDCZeroVal(sector, region, urow, upad) * getIDCOneVal(Sector(sector).side(), integrationInterval);
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::drawIDCZeroHelper(const bool type, const o2::tpc::Sector sector, const std::string filename, const float minZ, const float maxZ) const
{
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc = [this](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
    return this->getIDCZeroVal(sector, region, irow, pad);
  };
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
  const std::string zAxisTitle = IDCDrawHelper::getZAxisTitle(IDCType::IDCZero);
  type ? IDCDrawHelper::drawSide(drawFun, sector.side(), zAxisTitle, filename, minZ, maxZ) : IDCDrawHelper::drawSector(drawFun, 0, Mapper::NREGIONS, sector, zAxisTitle, filename, minZ, maxZ);
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::drawIDCDeltaHelper(const bool type, const Sector sector, const unsigned int integrationInterval, const std::string filename, const float minZ, const float maxZ) const
{
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc = [this, integrationInterval](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
    return this->getIDCDeltaVal(sector, region, irow, pad, integrationInterval);
  };
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
  const std::string zAxisTitle = IDCDrawHelper::getZAxisTitle(IDCType::IDCDelta);
  type ? IDCDrawHelper::drawSide(drawFun, sector.side(), zAxisTitle, filename, minZ, maxZ) : IDCDrawHelper::drawSector(drawFun, 0, Mapper::NREGIONS, sector, zAxisTitle, filename, minZ, maxZ);
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::drawIDCHelper(const bool type, const Sector sector, const unsigned int integrationInterval, const std::string filename, const float minZ, const float maxZ) const
{
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc = [this, integrationInterval](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
    return this->getIDCVal(sector, region, irow, pad, integrationInterval);
  };
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
  const std::string zAxisTitle = IDCDrawHelper::getZAxisTitle(IDCType::IDC);
  type ? IDCDrawHelper::drawSide(drawFun, sector.side(), zAxisTitle, filename, minZ, maxZ) : IDCDrawHelper::drawSector(drawFun, 0, Mapper::NREGIONS, sector, zAxisTitle, filename, minZ, maxZ);
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::drawPadFlagMap(const bool type, const Sector sector, const std::string filename, const PadFlags flag) const
{
  if (!mPadFlagsMap) {
    LOGP(info, "Status map not set returning");
    return;
  }
 
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc = [this, flag](const unsigned int sector, const unsigned int region, const unsigned int row, const unsigned int pad) {
    const unsigned int padInRegion = Mapper::OFFSETCRULOCAL[region][row] + pad;
    const auto flagDraw = mPadFlagsMap->getCalArray(region + sector * Mapper::NREGIONS).getValue(padInRegion);
    if ((flagDraw & flag) == flag) {
      return 1;
    } else {
      return 0;
    }
  };
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
  const std::string zAxisTitle = "status flag";
  type ? IDCDrawHelper::drawSide(drawFun, sector.side(), zAxisTitle, filename) : IDCDrawHelper::drawSector(drawFun, 0, Mapper::NREGIONS, sector, zAxisTitle, filename);
}
 
template <typename DataT>
unsigned int o2::tpc::IDCCCDBHelper<DataT>::getUngroupedIndexGlobal(const unsigned int sector, const unsigned int region, unsigned int urow, unsigned int upad, unsigned int integrationInterval)
{
  return IDCGroupHelperSector::getUngroupedIndexGlobal(sector, region, urow, upad, integrationInterval);
}
 
template <typename DataT>
TCanvas* o2::tpc::IDCCCDBHelper<DataT>::drawIDCZeroCanvas(TCanvas* outputCanvas, std::string_view type, int nbins1D, float xMin1D, float xMax1D, int integrationInterval) const
{
  TCanvas* canv = nullptr;
 
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc;
 
  if (type == "IDC0") {
    idcFunc = [this](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
      return this->getIDCZeroVal(sector, region, irow, pad);
    };
  } else if (type == "IDCDelta") {
    idcFunc = [this, integrationInterval](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
      return this->getIDCDeltaVal(sector, region, irow, pad, integrationInterval);
    };
    if (integrationInterval <= 0) {
      integrationInterval = std::min(getNIntegrationIntervalsIDCDelta(Side::A), getNIntegrationIntervalsIDCDelta(Side::C));
    }
  } else if (type == "IDC") {
    idcFunc = [this, integrationInterval](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
      return this->getIDCVal(sector, region, irow, pad, integrationInterval);
    };
    if (integrationInterval <= 0) {
      integrationInterval = std::min(getNIntegrationIntervalsIDCDelta(Side::A), getNIntegrationIntervalsIDCDelta(Side::C));
    }
  } else {
    LOGP(error, "Please provide a valid IDC data type. 'IDC0', 'IDCDelta', or 'IDC'.");
    return canv;
  }
 
  if (outputCanvas) {
    canv = outputCanvas;
  } else {
    canv = new TCanvas(fmt::format("c_sides_{}", type.data()).data(), fmt::format("{}", type.data()).data(), 1000, 1000);
  }
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
  const std::string zAxisTitle = IDCDrawHelper::getZAxisTitle(IDCType::IDCZero);
 
  auto hAside2D = IDCDrawHelper::drawSide(drawFun, Side::A, zAxisTitle);
  auto hCside2D = IDCDrawHelper::drawSide(drawFun, Side::C, zAxisTitle);
  hAside2D->SetTitle(fmt::format("{} (A-Side)", type.data()).data());
  hCside2D->SetTitle(fmt::format("{} (C-Side)", type.data()).data());
  hAside2D->SetMinimum(xMin1D);
  hAside2D->SetMaximum(xMax1D);
  hCside2D->SetMinimum(xMin1D);
  hCside2D->SetMaximum(xMax1D);
 
  auto hAside1D = IDCDrawHelper::drawSide(drawFun, fmt::format("{}", type.data()).data(), Side::A, nbins1D, xMin1D, xMax1D);
  auto hCside1D = IDCDrawHelper::drawSide(drawFun, fmt::format("{}", type.data()).data(), Side::C, nbins1D, xMin1D, xMax1D);
 
  canv->Divide(2, 2);
  canv->cd(1);
  hAside2D->Draw("colz");
  hAside2D->SetStats(0);
  hAside2D->SetTitleOffset(1.05, "XY");
  hAside2D->SetTitleSize(0.05, "XY");
  o2::tpc::painter::drawSectorsXY(Side::A);
 
  canv->cd(2);
  hCside2D->Draw("colz");
  hCside2D->SetStats(0);
  hCside2D->SetTitleOffset(1.05, "XY");
  hCside2D->SetTitleSize(0.05, "XY");
  o2::tpc::painter::drawSectorsXY(Side::C);
 
  canv->cd(3);
  hAside1D->Draw();
 
  canv->cd(4);
  hCside1D->Draw();
 
  hAside2D->SetBit(TObject::kCanDelete);
  hCside2D->SetBit(TObject::kCanDelete);
  hAside1D->SetBit(TObject::kCanDelete);
  hCside1D->SetBit(TObject::kCanDelete);
 
  return canv;
}
 
template <typename DataT>
TCanvas* o2::tpc::IDCCCDBHelper<DataT>::drawIDCZeroScale(TCanvas* outputCanvas, bool rejectOutlier) const
{
  TCanvas* canv = nullptr;
 
  if (outputCanvas) {
    canv = outputCanvas;
  } else {
    canv = new TCanvas("c_sides_IDC0_scale", "IDC0", 1000, 1000);
  }
  canv->cd();
  double xsides[] = {0, 1};
  double yTotIDC0[] = {mScaleIDC0Aside, mScaleIDC0Cside};
  auto gIDC0Scale = new TGraph(2, xsides, yTotIDC0);
  gIDC0Scale->SetName("g_IDC0ScaleFactor");
  gIDC0Scale->SetTitle("Scaling Factor (IDC0);Sides;IDC0 Total (arb. unit)");
  gIDC0Scale->SetMarkerColor(kBlue);
  gIDC0Scale->SetMarkerStyle(21);
  gIDC0Scale->SetMarkerSize(1);
  gIDC0Scale->GetXaxis()->SetLabelColor(0);
  gIDC0Scale->GetXaxis()->CenterTitle();
 
  gIDC0Scale->Draw("ap");
  // Draw labels on the x axis
  TText* t = new TText();
  t->SetTextSize(0.035);
  const char* labels[2] = {"A", "C"};
  for (Int_t i = 0; i < 2; i++) {
    t->DrawText(xsides[i], 0.5, labels[i]);
  }
  gIDC0Scale->GetXaxis()->SetLimits(-0.5, 1.5);
  canv->Update();
  canv->Modified();
  gIDC0Scale->SetBit(TObject::kCanDelete);
  t->SetBit(TObject::kCanDelete);
  return canv;
}
 
template <typename DataT>
TCanvas* o2::tpc::IDCCCDBHelper<DataT>::drawIDCZeroRadialProfile(TCanvas* outputCanvas, int nbinsY, float yMin, float yMax) const
{
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc = [this](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
    return this->getIDCZeroVal(sector, region, irow, pad);
  };
 
  TCanvas* canv = nullptr;
 
  if (outputCanvas) {
    canv = outputCanvas;
  } else {
    canv = new TCanvas("c_sides_IDC0_radialProfile", "IDC0", 1000, 1000);
  }
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
 
  const auto radialBinning = o2::tpc::painter::getRowBinningCM();
 
  auto hAside2D = new TH2F("h_IDC0_radialProfile_Aside", "IDC0: Radial profile (A-Side)", radialBinning.size() - 1, radialBinning.data(), nbinsY, yMin, yMax);
  hAside2D->GetXaxis()->SetTitle("x (cm)");
  hAside2D->GetYaxis()->SetTitle("IDC0 (ADC)");
  hAside2D->SetTitleOffset(1.05, "XY");
  hAside2D->SetTitleSize(0.05, "XY");
  hAside2D->SetStats(0);
 
  auto hCside2D = new TH2F("h_IDC0_radialProfile_Cside", "IDC0: Radial profile (C-Side)", radialBinning.size() - 1, radialBinning.data(), nbinsY, yMin, yMax);
  hCside2D->GetXaxis()->SetTitle("x (cm)");
  hCside2D->GetYaxis()->SetTitle("IDC0 (ADC)");
  hCside2D->SetTitleOffset(1.05, "XY");
  hCside2D->SetTitleSize(0.05, "XY");
  hCside2D->SetStats(0);
 
  IDCDrawHelper::drawRadialProfile(drawFun, *hAside2D, Side::A);
  IDCDrawHelper::drawRadialProfile(drawFun, *hCside2D, Side::C);
 
  canv->Divide(1, 2);
  canv->cd(1);
  hAside2D->Draw("colz");
  hAside2D->ProfileX("profile_ASide", 1, -1, "d,same");
  // profA->Draw("d,same");
  hAside2D->SetStats(0);
 
  canv->cd(2);
  hCside2D->Draw("colz");
  hCside2D->ProfileX("profile_CSide", 1, -1, "d,same");
  // profC->Draw("d,same");
  hAside2D->SetStats(0);
 
  hAside2D->SetBit(TObject::kCanDelete);
  hCside2D->SetBit(TObject::kCanDelete);
 
  return canv;
}
 
template <typename DataT>
TCanvas* o2::tpc::IDCCCDBHelper<DataT>::drawIDCZeroStackCanvas(TCanvas* outputCanvas, Side side, std::string_view type, int nbins1D, float xMin1D, float xMax1D, int integrationInterval) const
{
  TCanvas* canv = nullptr;
 
  std::function<float(const unsigned int, const unsigned int, const unsigned int, const unsigned int)> idcFunc;
 
  if (type == "IDC0") {
    idcFunc = [this](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
      return this->getIDCZeroVal(sector, region, irow, pad);
    };
  } else if (type == "IDCDelta") {
    idcFunc = [this, integrationInterval](const unsigned int sector, const unsigned int region, const unsigned int irow, const unsigned int pad) {
      return this->getIDCDeltaVal(sector, region, irow, pad, integrationInterval);
    };
    if (integrationInterval <= 0) {
      integrationInterval = std::min(getNIntegrationIntervalsIDCDelta(Side::A), getNIntegrationIntervalsIDCDelta(Side::C));
    }
  } else {
    LOGP(error, "Please provide a valid IDC data type. 'IDC0' or 'IDCDelta'.");
    return canv;
  }
 
  if (outputCanvas) {
    canv = outputCanvas;
  } else {
    canv = new TCanvas(fmt::format("c_GEMStacks_{}_1D_{}Side", type.data(), (side == Side::A) ? "A" : "C").data(), fmt::format("{} value 1D distribution for each GEM stack {}-Side", type.data(), (side == Side::A) ? "A" : "C").data(), 1000, 1000);
  }
 
  IDCDrawHelper::IDCDraw drawFun;
  drawFun.mIDCFunc = idcFunc;
  IDCDrawHelper::drawIDCZeroStackCanvas(drawFun, side, type, nbins1D, xMin1D, xMax1D, *canv, integrationInterval);
 
  return canv;
}
 
template <typename DataT>
TCanvas* o2::tpc::IDCCCDBHelper<DataT>::drawIDCOneCanvas(TCanvas* outputCanvas, int nbins1D, float xMin1D, float xMax1D, int integrationIntervals) const
{
  TCanvas* canv = nullptr;
 
  if (outputCanvas) {
    canv = outputCanvas;
  } else {
    canv = new TCanvas("c_sides_IDC1_1D", "IDC1 1D distribution for each side", 1000, 1000);
  }
 
  auto hAside1D = new TH1F("h_IDC1_1D_ASide", "IDC1 distribution over integration intervals A-Side", nbins1D, xMin1D, xMax1D);
  auto hCside1D = new TH1F("h_IDC1_1D_CSide", "IDC1 distribution over integration intervals C-Side", nbins1D, xMin1D, xMax1D);
 
  hAside1D->GetXaxis()->SetTitle("IDC1");
  hAside1D->SetTitleOffset(1.05, "XY");
  hAside1D->SetTitleSize(0.05, "XY");
  hCside1D->GetXaxis()->SetTitle("IDC1");
  hCside1D->SetTitleOffset(1.05, "XY");
  hCside1D->SetTitleSize(0.05, "XY");
 
  if (integrationIntervals <= 0) {
    integrationIntervals = std::min(mIDCOne[Side::A]->getNIDCs(), mIDCOne[Side::C]->getNIDCs());
  }
 
  for (unsigned int integrationInterval = 0; integrationInterval < integrationIntervals; ++integrationInterval) {
    hAside1D->Fill(getIDCOneVal(Side::A, integrationInterval));
    hCside1D->Fill(getIDCOneVal(Side::C, integrationInterval));
  }
 
  canv->Divide(1, 2);
  canv->cd(1);
  hAside1D->Draw();
  canv->cd(2);
  hCside1D->Draw();
 
  hAside1D->SetBit(TObject::kCanDelete);
  hCside1D->SetBit(TObject::kCanDelete);
 
  return canv;
}
 
template <typename DataT>
TCanvas* o2::tpc::IDCCCDBHelper<DataT>::drawFourierCoeff(TCanvas* outputCanvas, Side side, int nbins1D, float xMin1D, float xMax1D) const
{
  TCanvas* canv = nullptr;
 
  if (outputCanvas) {
    canv = outputCanvas;
  } else {
    canv = new TCanvas(fmt::format("c_FourierCoefficients_1D_{}Side", (side == Side::A) ? "A" : "C").data(), fmt::format("1D distributions of Fourier Coefficients ({}-Side)", (side == Side::A) ? "A" : "C").data(), 1000, 1000);
  }
 
  std::vector<TH1F*> histos;
 
  for (int i = 0; i < mFourierCoeff[side]->getNCoefficientsPerTF(); i++) {
    histos.emplace_back(new TH1F(fmt::format("h_FourierCoeff{}_{}Side", i, (side == Side::A) ? "A" : "C").data(), fmt::format("1D distribution of Fourier Coefficient {} ({}-Side)", i, (side == Side::A) ? "A" : "C").data(), nbins1D, xMin1D, xMax1D));
    histos.back()->GetXaxis()->SetTitle(fmt::format("Fourier Coefficient {}", i).data());
    histos.back()->SetBit(TObject::kCanDelete);
  }
 
  const auto& coeffs = mFourierCoeff[side]->getFourierCoefficients();
  const auto nCoeffPerTF = mFourierCoeff[side]->getNCoefficientsPerTF();
 
  for (int i = 0; i < mFourierCoeff[side]->getNCoefficients(); i++) {
    histos.at(i % nCoeffPerTF)->Fill(coeffs.at(i));
  }
 
  canv->DivideSquare(mFourierCoeff[side]->getNCoefficientsPerTF());
 
  size_t pad = 1;
 
  for (const auto& hist : histos) {
    canv->cd(pad);
    hist->SetTitleOffset(1.05, "XY");
    hist->SetTitleSize(0.05, "XY");
    hist->Draw();
    pad++;
  }
 
  return canv;
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::dumpToTree(const Side side, const char* outFileName) const
{
  const Mapper& mapper = Mapper::instance();
  o2::utils::TreeStreamRedirector pcstream(outFileName, "RECREATE");
  pcstream.GetFile()->cd();
 
  const unsigned int nIDCsSector = Mapper::getPadsInSector() * Mapper::NSECTORS;
  std::vector<int> vRow(nIDCsSector);
  std::vector<int> vPad(nIDCsSector);
  std::vector<float> vXPos(nIDCsSector);
  std::vector<float> vYPos(nIDCsSector);
  std::vector<float> vGlobalXPos(nIDCsSector);
  std::vector<float> vGlobalYPos(nIDCsSector);
  std::vector<float> idcsZero(nIDCsSector);
  std::vector<unsigned int> sectorv(nIDCsSector);
  std::vector<int> outlier(nIDCsSector);
  std::vector<float> stackMedian(nIDCsSector);
 
  if (!mIDCZero[side]) {
    LOGP(info, "IDC0 not set. returning");
    return;
  }
  auto stacksMedianTmp = getStackMedian(*mIDCZero[side], side);
  const std::array<int, Mapper::NREGIONS> stackInSector{0, 0, 0, 0, 1, 1, 2, 2, 3, 3};
  unsigned int index = 0;
  const unsigned int secStart = (side == Side::A) ? 0 : SECTORSPERSIDE;
  const unsigned int secEnd = (side == Side::A) ? SECTORSPERSIDE : Mapper::NSECTORS;
  for (unsigned int sector = secStart; sector < secEnd; ++sector) {
    for (unsigned int region = 0; region < Mapper::NREGIONS; ++region) {
      for (unsigned int irow = 0; irow < Mapper::ROWSPERREGION[region]; ++irow) {
        for (unsigned int ipad = 0; ipad < Mapper::PADSPERROW[region][irow]; ++ipad) {
          const auto padNum = Mapper::getGlobalPadNumber(irow, ipad, region);
          const auto padTmp = (sector < SECTORSPERSIDE) ? ipad : (Mapper::PADSPERROW[region][irow] - ipad - 1); // C-Side is mirrored
          const auto& padPosLocal = mapper.padPos(padNum);
          vRow[index] = padPosLocal.getRow();
          vPad[index] = padPosLocal.getPad();
          vXPos[index] = mapper.getPadCentre(padPosLocal).X();
          vYPos[index] = mapper.getPadCentre(padPosLocal).Y();
          const GlobalPosition2D globalPos = mapper.LocalToGlobal(LocalPosition2D(vXPos[index], -vYPos[index]), Sector(sector));
          vGlobalXPos[index] = globalPos.X();
          vGlobalYPos[index] = globalPos.Y();
          idcsZero[index] = getIDCZeroVal(sector, region, irow, padTmp);
          if (mPadFlagsMap) {
            const unsigned int padInRegion = Mapper::OFFSETCRULOCAL[region][irow] + padTmp;
            outlier[index] = static_cast<int>(mPadFlagsMap->getCalArray(region + sector * Mapper::NREGIONS).getValue(padInRegion));
          }
          stackMedian[index] = stacksMedianTmp[(sector - secStart) * GEMSTACKSPERSECTOR + stackInSector[region]];
          sectorv[index++] = sector;
        }
      }
    }
  }
  std::vector<float> idcOne = !mIDCOne[side] ? std::vector<float>() : mIDCOne[side]->mIDCOne;
 
  pcstream << "tree"
           << "IDC0.=" << idcsZero
           << "IDC1=" << idcOne
           << "pad.=" << vPad
           << "row.=" << vRow
           << "lx.=" << vXPos
           << "ly.=" << vYPos
           << "gx.=" << vGlobalXPos
           << "gy.=" << vGlobalYPos
           << "outlier.=" << outlier
           << "sector.=" << sectorv
           << "stacksMedian=" << stackMedian
           << "\n";
  pcstream.Close();
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::dumpToTreeIDCDelta(const Side side, const char* outFileName) const
{
  o2::utils::TreeStreamRedirector pcstream(outFileName, "RECREATE");
  pcstream.GetFile()->cd();
  const unsigned int nIDCsSector = Mapper::getPadsInSector() * SECTORSPERSIDE;
  std::vector<float> idcsDelta(nIDCsSector);
  for (int integrationInterval = 0; integrationInterval < getNIntegrationIntervalsIDCDelta(side); ++integrationInterval) {
    unsigned int index = 0;
    const unsigned int secStart = (side == Side::A) ? 0 : SECTORSPERSIDE;
    const unsigned int secEnd = (side == Side::A) ? SECTORSPERSIDE : Mapper::NSECTORS;
    for (unsigned int sector = secStart; sector < secEnd; ++sector) {
      for (unsigned int region = 0; region < Mapper::NREGIONS; ++region) {
        for (unsigned int irow = 0; irow < Mapper::ROWSPERREGION[region]; ++irow) {
          for (unsigned int ipad = 0; ipad < Mapper::PADSPERROW[region][irow]; ++ipad) {
            const auto padTmp = (sector < SECTORSPERSIDE) ? ipad : (Mapper::PADSPERROW[region][irow] - ipad - 1); // C-Side is mirrored
            idcsDelta[index++] = getIDCDeltaVal(sector, region, irow, padTmp, integrationInterval);
          }
        }
      }
    }
    float idcOne = getIDCOneVal(side, integrationInterval);
 
    pcstream << "tree"
             << "IDC1=" << idcOne
             << "IDCDelta.=" << idcsDelta
             << "\n";
  }
  pcstream.Close();
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::dumpToFourierCoeffToTree(const char* outFileName) const
{
  o2::utils::TreeStreamRedirector pcstream(outFileName, "RECREATE");
  pcstream.GetFile()->cd();
 
  for (int iside = 0; iside < SIDES; ++iside) {
    const Side side = (iside == 0) ? Side::A : Side::C;
 
    if (!mFourierCoeff[side]) {
      continue;
    }
 
    const int nTFs = mFourierCoeff[side]->getNCoefficients() / mFourierCoeff[side]->getNCoefficientsPerTF();
    for (int iTF = 0; iTF < nTFs; ++iTF) {
      std::vector<float> coeff;
      std::vector<int> ind;
      int coeffPerTF = mFourierCoeff[side]->getNCoefficientsPerTF();
      for (int i = 0; i < coeffPerTF; ++i) {
        const int index = mFourierCoeff[side]->getIndex(iTF, i);
        coeff.emplace_back((*mFourierCoeff[side])(index));
        ind.emplace_back(i);
      }
 
      pcstream << "tree"
               << "iTF=" << iTF
               << "index=" << ind
               << "coeffPerTF=" << coeffPerTF
               << "coeff.=" << coeff
               << "side=" << iside
               << "\n";
    }
  }
  pcstream.Close();
}
 
template <typename DataT>
o2::tpc::CalDet<float> o2::tpc::IDCCCDBHelper<DataT>::getIDCZeroCalDet() const
{
  CalDet<float> calIDC0("IDC0");
  for (unsigned int cru = 0; cru < CRU::MaxCRU; ++cru) {
    const o2::tpc::CRU cruTmp(cru);
    const Side side = cruTmp.side();
    const int region = cruTmp.region();
    const int sector = cruTmp.sector();
    for (unsigned int lrow = 0; lrow < Mapper::ROWSPERREGION[region]; ++lrow) {
      const unsigned int integrationInterval = 0;
      for (unsigned int pad = 0; pad < Mapper::PADSPERROW[region][lrow]; ++pad) {
        const auto idcZero = getIDCZeroVal(sector, region, lrow, pad);
        calIDC0.setValue(sector, Mapper::ROWOFFSET[region] + lrow, pad, idcZero);
      }
    }
  }
  return calIDC0;
}
 
template <typename DataT>
std::vector<o2::tpc::CalDet<float>> o2::tpc::IDCCCDBHelper<DataT>::getIDCDeltaCalDet() const
{
  const unsigned int nIntervalsA = getNIntegrationIntervalsIDCDelta(Side::A);
  const unsigned int nIntervalsC = getNIntegrationIntervalsIDCDelta(Side::C);
  if (nIntervalsA != nIntervalsC) {
    LOGP(info, "Number of integration interval for A Side {} unequal for C side {}", nIntervalsA, nIntervalsC);
  }
 
  std::vector<CalPad> calIDCDelta(std::max(nIntervalsA, nIntervalsC));
  for (auto& calpadIDC : calIDCDelta) {
    calpadIDC = CalPad("IDCDelta", PadSubset::ROC);
  }
 
  for (unsigned int cru = 0; cru < CRU::MaxCRU; ++cru) {
    const o2::tpc::CRU cruTmp(cru);
    const Side side = cruTmp.side();
    const int region = cruTmp.region();
    const int sector = cruTmp.sector();
 
    const unsigned int nIntervals = getNIntegrationIntervalsIDCDelta(side);
    for (unsigned int integrationInterval = 0; integrationInterval < nIntervals; ++integrationInterval) {
      for (unsigned int lrow = 0; lrow < Mapper::ROWSPERREGION[region]; ++lrow) {
        for (unsigned int pad = 0; pad < Mapper::PADSPERROW[region][lrow]; ++pad) {
          const auto idcdelta = getIDCDeltaVal(sector, region, lrow, pad, integrationInterval);
          calIDCDelta[integrationInterval].setValue(sector, Mapper::ROWOFFSET[region] + lrow, pad, idcdelta);
        }
      }
    }
  }
  return calIDCDelta;
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::createOutlierMap()
{
  if (!mIDCZero[Side::A] && !mIDCZero[Side::C]) {
    LOGP(info, "IDC0 not set. returning");
  }
 
  int nCRUS = 0;
  if (mIDCZero[Side::A]) {
    nCRUS += CRU::MaxCRU / 2;
  }
 
  if (mIDCZero[Side::C]) {
    nCRUS += CRU::MaxCRU / 2;
  }
 
  std::vector<uint32_t> crus(nCRUS);
  std::iota(crus.begin(), crus.end(), mIDCZero[Side::A] ? 0 : (CRU::MaxCRU / 2));
 
  IDCFactorization idc(1, 1, crus);
  if (mIDCZero[Side::A]) {
    idc.setIDCZero(Side::A, *mIDCZero[Side::A]);
  }
  if (mIDCZero[Side::C]) {
    idc.setIDCZero(Side::C, *mIDCZero[Side::C]);
  }
  idc.createStatusMap();
  mPadFlagsMap = idc.getPadStatusMap();
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::scaleIDC0(const float factor, const Side side)
{
  if (!mIDCZero[side]) {
    LOGP(info, "IDC0 not set. returning");
    return;
  }
 
  *mIDCZero[side] /= factor;
}
 
template <typename DataT>
void o2::tpc::IDCCCDBHelper<DataT>::setIDCZeroScale(const bool rejectOutlier)
{
  if (rejectOutlier) {
    createOutlierMap();
  }
 
  if (mIDCZero[Side::A]) {
    mScaleIDC0Aside = getMeanIDC0(Side::A, *mIDCZero[Side::A], rejectOutlier ? mPadFlagsMap.get() : nullptr);
    scaleIDC0(mScaleIDC0Aside, Side::A);
  } else {
    mScaleIDC0Aside = 0;
  }
 
  if (mIDCZero[Side::C]) {
    mScaleIDC0Cside = getMeanIDC0(Side::C, *mIDCZero[Side::C], rejectOutlier ? mPadFlagsMap.get() : nullptr);
    scaleIDC0(mScaleIDC0Cside, Side::C);
  } else {
    mScaleIDC0Cside = 0;
  }
 
  if (mScaleIDC0Aside == 0.0 || mScaleIDC0Cside == 0.0) {
    LOGP(error, "Please check the IDC0 total A side {} and C side {}, is zero, therefore no scaling applied!", mScaleIDC0Aside, mScaleIDC0Cside);
    mScaleIDC0Aside = 1.0;
    mScaleIDC0Cside = 1.0;
  }
}
 
template <typename DataT>
float o2::tpc::IDCCCDBHelper<DataT>::getMeanIDC0(const Side side, const IDCZero& idcZero, const CalDet<PadFlags>* outlierMap)
{
  // 1. calculate mean per row
  std::vector<double> idc0Sum(Mapper::PADROWS);
  std::vector<unsigned int> idcChannelsCount(Mapper::PADROWS);
  const unsigned int firstCRU = (side == Side::A) ? 0 : CRU::MaxCRU / 2;
  for (unsigned int cru = firstCRU; cru < firstCRU + CRU::MaxCRU / 2; ++cru) {
    const o2::tpc::CRU cruTmp(cru);
    const int region = cruTmp.region();
    const int sector = cruTmp.sector();
    for (unsigned int lrow = 0; lrow < Mapper::ROWSPERREGION[region]; ++lrow) {
      const unsigned int glRow = lrow + Mapper::ROWOFFSET[region];
      for (unsigned int pad = 0; pad < Mapper::PADSPERROW[region][lrow]; ++pad) {
        const unsigned int padInRegion = Mapper::OFFSETCRULOCAL[region][lrow] + pad;
        if (outlierMap) {
          const o2::tpc::PadFlags flag = outlierMap->getCalArray(cru).getValue(padInRegion);
          if ((flag & PadFlags::flagSkip) == PadFlags::flagSkip) {
            continue;
          }
        }
        const auto index = getUngroupedIndexGlobal(sector, region, lrow, pad, 0);
        idc0Sum[glRow] += idcZero.getValueIDCZero(index);
        ++idcChannelsCount[glRow];
      }
    }
  }
 
  // 2. calculate mean of means
  double idc0SumGlobal = 0;
  unsigned int idc0CountGlobal = 0;
  for (unsigned int i = 0; i < Mapper::PADROWS; ++i) {
    if (idcChannelsCount[i]) {
      idc0SumGlobal += idc0Sum[i] / idcChannelsCount[i];
      ++idc0CountGlobal;
    }
  }
 
  if (idc0CountGlobal == 0) {
    return -1;
  }
 
  return idc0SumGlobal /= idc0CountGlobal;
}
 
template <typename DataT>
float o2::tpc::IDCCCDBHelper<DataT>::getStackMedian(const IDCZero& idczero, const Sector sector, const GEMstack stack)
{
  int cruStart = sector * CRU::CRUperSector;
  if (stack == GEMstack::OROC1gem) {
    cruStart += CRU::CRUperIROC;
  } else if (stack == GEMstack::OROC2gem) {
    cruStart += CRU::CRUperIROC + CRU::CRUperPartition;
  } else if (stack == GEMstack::OROC3gem) {
    cruStart += CRU::CRUperIROC + 2 * CRU::CRUperPartition;
  }
  int cruEnd = (stack == GEMstack::IROCgem) ? (cruStart + CRU::CRUperIROC) : (cruStart + CRU::CRUperPartition);
 
  RobustAverage average;
  average.reserve(Mapper::getNumberOfPads(stack));
  for (int cru = cruStart; cru < cruEnd; ++cru) {
    const o2::tpc::CRU cruTmp(cru);
    const int region = cruTmp.region();
    for (unsigned int lrow = 0; lrow < Mapper::ROWSPERREGION[region]; ++lrow) {
      const unsigned int integrationInterval = 0;
      for (unsigned int pad = 0; pad < Mapper::PADSPERROW[region][lrow]; ++pad) {
        const auto index = getUngroupedIndexGlobal(sector, region, lrow, pad, integrationInterval);
        const float idcZero = idczero.mIDCZero[index];
        if ((idcZero != -1) && (idcZero != 0)) {
          average.addValue(idcZero);
        }
      }
    }
  }
  return average.getMedian();
}
 
template <typename DataT>
std::array<float, o2::tpc::GEMSTACKSPERSECTOR * o2::tpc::SECTORSPERSIDE> o2::tpc::IDCCCDBHelper<DataT>::getStackMedian(const IDCZero& idczero, const Side side)
{
  const int stacksPerSide = o2::tpc::GEMSTACKSPERSECTOR * o2::tpc::SECTORSPERSIDE;
  std::array<float, stacksPerSide> median;
  unsigned int sectorStart = (side == Side::A) ? 0 : o2::tpc::SECTORSPERSIDE;
  unsigned int sectorEnd = (side == Side::A) ? o2::tpc::SECTORSPERSIDE : Mapper::NSECTORS;
  for (unsigned int sector = sectorStart; sector < sectorEnd; ++sector) {
    for (int stackInSector = 0; stackInSector < GEMSTACKSPERSECTOR; ++stackInSector) {
      median[(sector * GEMSTACKSPERSECTOR + stackInSector) % stacksPerSide] = IDCCCDBHelper<>::getStackMedian(idczero, Sector(sector), static_cast<GEMstack>(stackInSector));
    }
  }
  return median;
}
 
template <typename DataT>
std::pair<int, int> o2::tpc::IDCCCDBHelper<DataT>::getNOutliers() const
{
  std::pair<int, int> nOutliers{};
  if (!mPadFlagsMap) {
    LOGP(info, "Status map not set returning");
    return nOutliers;
  }
  for (CRU cru; !cru.looped(); ++cru) {
    const auto& calArray = mPadFlagsMap->getCalArray(cru);
    for (auto flag : calArray.getData()) {
      if ((flag & PadFlags::flagSkip) == PadFlags::flagSkip) {
        (cru.side() == Side::A) ? ++nOutliers.first : ++nOutliers.second;
      }
    }
  }
  return nOutliers;
}
 
template class o2::tpc::IDCCCDBHelper<float>;
template class o2::tpc::IDCCCDBHelper<unsigned short>;
template class o2::tpc::IDCCCDBHelper<unsigned char>;