EMCALChannelCalibrator.h

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.
 
 
#ifndef EMCAL_CHANNEL_CALIBRATOR_H_
#define EMCAL_CHANNEL_CALIBRATOR_H_
 
#include "EMCALCalibration/EMCALTimeCalibData.h"
#include "EMCALCalibration/EMCALChannelData.h"
#include "EMCALCalibration/EMCALCalibExtractor.h"
#include "EMCALCalibration/EMCALCalibParams.h"
#include "EMCALCalib/GainCalibrationFactors.h"
#include "DetectorsCalibration/TimeSlotCalibration.h"
#include "DetectorsCalibration/TimeSlot.h"
#include "DataFormatsEMCAL/Cell.h"
#include "EMCALBase/Geometry.h"
#include "CCDB/CcdbObjectInfo.h"
#include "EMCALCalib/CalibDB.h"
#include "DataFormatsParameters/GRPECSObject.h"
 
#include "Framework/Logger.h"
#include "CommonUtils/MemFileHelper.h"
#include "CCDB/CcdbApi.h"
#include "DetectorsCalibration/Utils.h"
#include <boost/histogram.hpp>
#include <boost/histogram/ostream.hpp>
#include <boost/format.hpp>
 
#include <array>
#include <boost/histogram.hpp>
#include <fstream>
 
namespace o2
{
namespace emcal
{
template <typename DataInput, typename DataOutput>
class EMCALChannelCalibrator : public o2::calibration::TimeSlotCalibration<DataInput>
{
  using TFType = o2::calibration::TFType;
  using Slot = o2::calibration::TimeSlot<DataInput>;
  using Cell = o2::emcal::Cell;
  using CcdbObjectInfo = o2::ccdb::CcdbObjectInfo;
  using CcdbObjectInfoVector = std::vector<CcdbObjectInfo>;
 
 public:
  EMCALChannelCalibrator(int nb = 1000, float r = 0.35) : mNBins(nb), mRange(r)
  {
    this->setSaveDirectory(EMCALCalibParams::Instance().filePathSave);
  };
 
  ~EMCALChannelCalibrator() final = default;
 
  bool hasEnoughData(const Slot& slot) const final;
  void initOutput() final;
  void finalizeSlot(Slot& slot) final;
  o2::calibration::TimeSlot<DataInput>& emplaceNewSlot(bool front, TFType tstart, TFType tend) final;
 
  bool saveLastSlotData(TFile& fl) final;
 
  bool adoptSavedData(const o2::calibration::TimeSlotMetaData& metadata, TFile& fl) final;
 
  void setIsTest(bool isTest) { mTest = isTest; }
  bool isTest() const { return mTest; }
 
  const CcdbObjectInfoVector& getInfoVector() const { return mInfoVector; }
  const std::vector<DataOutput>& getOutputVector() const { return mCalibObjectVector; }
 
  bool getSaveAtEOR() const { return mSaveAtEOR; }
  void setSaveAtEOR(bool tmp) { mSaveAtEOR = tmp; }
 
  bool getLoadAtSOR() const { return mLoadAtSOR; }
  void setLoadAtSOR(bool tmp) { mLoadAtSOR = tmp; }
 
  std::shared_ptr<EMCALCalibExtractor> getCalibExtractor() { return mCalibrator; } // return shared pointer!
  void SetCalibExtractor(std::shared_ptr<EMCALCalibExtractor> extr) { mCalibrator = extr; };
 
  bool setGainCalibrationFactors(o2::emcal::GainCalibrationFactors* gainCalibFactors);
 
  void setFillNr(int fn) { mFillNr = fn; }
 
  void setRunType(o2::parameters::GRPECSObject::RunType rt) { mRunType = rt; }
 
  void setCurrTSInHours(int ts) { mStartTSCalib = ts; }
 
 private:
  int mNBins = 0;                                                                               
  float mRange = 0.;                                                                            
  bool mTest = false;                                                                           
  bool mSaveAtEOR = false;                                                                      
  bool mLoadAtSOR = false;                                                                      
  o2::parameters::GRPECSObject::RunType mRunType = o2::parameters::GRPECSObject::RunType::NONE; 
  int mFillNr = 0;                                                                              
  int mStartTSCalib = -1;                                                                       
  std::shared_ptr<EMCALCalibExtractor> mCalibrator;
 
  // output
  CcdbObjectInfoVector mInfoVector; // vector of CCDB Infos , each element is filled with the CCDB description of the accompanying TimeSlewing object
  std::vector<DataOutput> mCalibObjectVector;
 
  ClassDefOverride(EMCALChannelCalibrator, 1);
};
 
//_____________________________________________
template <typename DataInput, typename DataOutput>
void EMCALChannelCalibrator<DataInput, DataOutput>::initOutput()
{
  mInfoVector.clear();
  mCalibObjectVector.clear();
  std::string nameFile = "tcp.root";
  if constexpr (std::is_same<DataInput, o2::emcal::EMCALChannelData>::value) {
    nameFile = "bcm.root";
  }
 
  // this->setSaveFileName(nameFile);
  // mNEvents = 0;
  return;
}
 
//_____________________________________________
template <typename DataInput, typename DataOutput>
bool EMCALChannelCalibrator<DataInput, DataOutput>::hasEnoughData(const o2::calibration::TimeSlot<DataInput>& slot) const
{
  // in case of the end of run, we pretend to have enough data to trigger the saving of the calib objects to load them in the next run
  if (mSaveAtEOR) {
    return true;
  }
  const DataInput* c = slot.getContainer();
  return (mTest ? true : c->hasEnoughData());
}
 
//_____________________________________________
template <typename DataInput, typename DataOutput>
void EMCALChannelCalibrator<DataInput, DataOutput>::finalizeSlot(o2::calibration::TimeSlot<DataInput>& slot)
{
 
  // Extract results for the single slot
  DataInput* c = slot.getContainer();
  LOG(info) << "Finalize slot " << slot.getTFStart() << " <= TF <= " << slot.getTFEnd();
  // check if slot contains a minimum amount of data
  if (c->getNEvents() < EMCALCalibParams::Instance().minNEventsSaveSlot) {
    LOG(info) << "Slot only contains " << c->getNEvents() << " events. Not saving this slot. " << EMCALCalibParams::Instance().minNEventsSaveSlot << " required";
    return;
  }
 
  if constexpr (std::is_same<DataInput, o2::emcal::EMCALChannelData>::value) {
    if (c->getNEvents() < EMCALCalibParams::Instance().minNEvents_bc) {
      LOG(info) << "saving the slot with " << c->getNEvents() << " events. " << EMCALCalibParams::Instance().minNEvents_tc << " events needed for calibration";
      this->saveLastSlot();
      return;
    }
  } else if constexpr (std::is_same<DataInput, o2::emcal::EMCALTimeCalibData>::value) {
    if (c->getNEvents() < EMCALCalibParams::Instance().minNEvents_tc) {
      LOG(info) << "saving the slot with " << c->getNEvents() << " events. " << EMCALCalibParams::Instance().minNEvents_tc << " events needed for calibration";
      this->saveLastSlot();
      return;
    }
  }
 
  // get the start timestamp of the slot. Either use the manually set start ts if available, or just the start time of the slot.
  long tsStart = slot.getStartTimeMS();
  if (EMCALCalibParams::Instance().useStaticStartTimeSlot && slot.getStaticStartTimeMS() > 0) {
    tsStart = slot.getStaticStartTimeMS();
    LOG(info) << "Adjusting the start TS of the slot from " << slot.getStartTimeMS() << " to " << slot.getStaticStartTimeMS();
  }
 
  std::map<std::string, std::string> md;
  if constexpr (std::is_same<DataInput, o2::emcal::EMCALChannelData>::value) {
    LOG(debug) << "Launching the calibration.";
    o2::emcal::BadChannelMap bcm = mCalibrator->calibrateBadChannels(c->getHisto(), c->getHistoTime());
    LOG(debug) << "Done with the calibraiton";
    // for the CCDB entry
    auto clName = o2::utils::MemFileHelper::getClassName(bcm);
    auto flName = o2::ccdb::CcdbApi::generateFileName(clName);
    mInfoVector.emplace_back(CalibDB::getCDBPathBadChannelMap(), clName, flName, md, tsStart, slot.getEndTimeMS() + EMCALCalibParams::Instance().endTimeMargin, true);
    mCalibObjectVector.push_back(bcm);
 
    if ((EMCALCalibParams::Instance().localRootFilePath).find(".root") != std::string::npos) {
      std::ifstream ffile(EMCALCalibParams::Instance().localRootFilePath.c_str());
 
      TFile fLocalStorage((EMCALCalibParams::Instance().localRootFilePath).c_str(), ffile.good() == true ? "update" : "recreate");
      fLocalStorage.cd();
 
      TH2F* histBCMap = (TH2F*)bcm.getHistogramRepresentation();
      std::string nameBCHist = "BadChannels_" + std::to_string(tsStart);
      histBCMap->Write(nameBCHist.c_str(), TObject::kOverwrite);
 
      TH2F hCalibHist = o2::utils::TH2FFromBoost(c->getHisto());
      std::string nameBCInputHist = "EnergyVsCellID_" + std::to_string(tsStart);
      hCalibHist.Write(nameBCInputHist.c_str(), TObject::kOverwrite);
 
      TH2F hCalibHistTime = o2::utils::TH2FFromBoost(c->getHistoTime());
      std::string nameBCInputHistTime = "TimeVsCellID_" + std::to_string(tsStart);
      hCalibHistTime.Write(nameBCInputHistTime.c_str(), TObject::kOverwrite);
 
      fLocalStorage.Close();
    }
  } else if constexpr (std::is_same<DataInput, o2::emcal::EMCALTimeCalibData>::value) {
    o2::emcal::TimeCalibrationParams tcd = mCalibrator->calibrateTime(c->getHisto(), EMCALCalibParams::Instance().minTimeForFit_tc, EMCALCalibParams::Instance().maxTimeForFit_tc, EMCALCalibParams::Instance().restrictFitRangeToMax_tc);
 
    // for the CCDB entry
    auto clName = o2::utils::MemFileHelper::getClassName(slot);
    auto flName = o2::ccdb::CcdbApi::generateFileName(clName);
 
    // prepareCCDBobjectInfo
    mInfoVector.emplace_back(CalibDB::getCDBPathTimeCalibrationParams(), clName, flName, md, tsStart, slot.getEndTimeMS() + EMCALCalibParams::Instance().endTimeMargin, true);
    mCalibObjectVector.push_back(tcd);
 
    if ((EMCALCalibParams::Instance().localRootFilePath).find(".root") != std::string::npos) {
      std::ifstream ffile(EMCALCalibParams::Instance().localRootFilePath.c_str());
 
      TFile fLocalStorage((EMCALCalibParams::Instance().localRootFilePath).c_str(), ffile.good() == true ? "update" : "recreate");
      fLocalStorage.cd();
      TH1F* histTCparams = (TH1F*)tcd.getHistogramRepresentation(false); // high gain calibration
      std::string nameTCHist = "TCParams_HG_" + std::to_string(tsStart);
      histTCparams->Write(nameTCHist.c_str(), TObject::kOverwrite);
 
      TH1F* histTCparams_LG = (TH1F*)tcd.getHistogramRepresentation(true); // low gain calibration
      std::string nameTCHist_LG = "TCParams_LG_" + std::to_string(tsStart);
      histTCparams_LG->Write(nameTCHist_LG.c_str(), TObject::kOverwrite);
 
      TH2F hCalibHist = o2::utils::TH2FFromBoost(c->getHisto());
      std::string nameTCInputHist = "TimeVsCellID_" + std::to_string(tsStart);
      hCalibHist.Write(nameTCInputHist.c_str(), TObject::kOverwrite);
      fLocalStorage.Close();
    }
  }
}
 
template <typename DataInput, typename DataOutput>
o2::calibration::TimeSlot<DataInput>& EMCALChannelCalibrator<DataInput, DataOutput>::emplaceNewSlot(bool front, TFType tstart, TFType tend)
{
  auto& cont = o2::calibration::TimeSlotCalibration<DataInput>::getSlots();
  auto& slot = front ? cont.emplace_front(tstart, tend) : cont.emplace_back(tstart, tend);
  slot.setContainer(std::make_unique<DataInput>());
  return slot;
}
 
template <typename DataInput, typename DataOutput>
bool EMCALChannelCalibrator<DataInput, DataOutput>::saveLastSlotData(TFile& fl)
{
  LOG(info) << "EMC calib histos are saved in " << fl.GetName();
  // we only have 1 slot
  auto& cont = o2::calibration::TimeSlotCalibration<DataInput>::getSlots();
  auto& slot = cont.at(0);
  DataInput* c = slot.getContainer();
 
  // get the start timestamp of the slot. Either use the manually set start ts if available, or just the start time of the slot.
  long tsStart = slot.getStartTimeMS();
  if (EMCALCalibParams::Instance().useStaticStartTimeSlot && slot.getStaticStartTimeMS() > 0) {
    tsStart = slot.getStaticStartTimeMS();
  }
  // timestamp in hours and seconds
  int timeNowHour = static_cast<int>(tsStart / o2::ccdb::CcdbObjectInfo::HOUR);
  int timeNowSec = static_cast<int>(tsStart / o2::ccdb::CcdbObjectInfo::SECOND);
 
  // create global event properties histogram. Same for both calibration types
  TH1I hGlobalProperties("hGlobalProperties", "hGlobalProperties", 4, -0.5, 3.5);
  hGlobalProperties.GetXaxis()->SetBinLabel(1, "Fill nr.");
  hGlobalProperties.GetXaxis()->SetBinLabel(2, "run type");
  hGlobalProperties.GetXaxis()->SetBinLabel(3, "ts in hours");
  hGlobalProperties.GetXaxis()->SetBinLabel(4, "ts in seconds");
  hGlobalProperties.SetBinContent(1, mFillNr);
  hGlobalProperties.SetBinContent(2, mRunType);
  hGlobalProperties.SetBinContent(3, timeNowHour);
  hGlobalProperties.SetBinContent(4, timeNowSec);
 
  if constexpr (std::is_same<DataInput, o2::emcal::EMCALChannelData>::value) {
    auto hist = c->getHisto();
    auto histTime = c->getHistoTime();
 
    TH2F hEnergy = o2::utils::TH2FFromBoost(hist);
    TH2F hTime = o2::utils::TH2FFromBoost(histTime, "histTime");
    TH1D hNEvents("hNEvents", "hNEvents", 1, 0, 1);
    hNEvents.SetBinContent(1, c->getNEvents());
 
    fl.cd();
    hEnergy.Write("EnergyVsCellID");
    hTime.Write("TimeVsCellID");
    hNEvents.Write("NEvents");
    hGlobalProperties.Write("GlobalProperties");
 
  } else if constexpr (std::is_same<DataInput, o2::emcal::EMCALTimeCalibData>::value) {
    auto histTime = c->getHisto();
    TH2F hTime = o2::utils::TH2FFromBoost(histTime);
    TH1D hNEvents("hNEvents", "hNEvents", 1, 0, 1);
    hNEvents.SetBinContent(1, c->getNEvents());
 
    fl.cd();
    hTime.Write("TimeVsCellID");
    hNEvents.Write("NEvents");
    hGlobalProperties.Write("GlobalProperties");
  }
 
  return true;
}
 
template <typename DataInput, typename DataOutput>
bool EMCALChannelCalibrator<DataInput, DataOutput>::adoptSavedData(const o2::calibration::TimeSlotMetaData& metadata, TFile& fl)
{
  LOG(info) << "Loading data from previous run";
 
  if (!this->getSavedSlotAllowed() || !this->getLoadAtSOR())
    return true;
 
  auto& cont = o2::calibration::TimeSlotCalibration<DataInput>::getSlots();
 
  if (cont.size() == 0) {
    LOG(warning) << "cont.size() is 0, calibration objects from previous run cannot be loaded...";
    return true;
  }
  auto& slot = cont.at(0);
  DataInput* c = slot.getContainer();
 
  // check run type and fill
  TH1I* hGlobalProperties = (TH1I*)fl.Get("GlobalProperties");
  if (!hGlobalProperties) {
    LOG(error) << "GlobalProperties histogram not found. Will not load previous calibration histograms";
  } else {
    int fillNr = hGlobalProperties->GetBinContent(1);
    int runType = hGlobalProperties->GetBinContent(2);
    int tsOld = hGlobalProperties->GetBinContent(3);
    int tsOldSec = (hGlobalProperties->GetNbinsX() > 3) ? hGlobalProperties->GetBinContent(4) : 0;                                                 // Protection as 4th bin was only added later
    int tsDiff = (mStartTSCalib > 0 ? mStartTSCalib : static_cast<int>(o2::ccdb::getCurrentTimestamp() / o2::ccdb::CcdbObjectInfo::HOUR)) - tsOld; // get current timestamp if mStartTSCalib is not set
    LOG(debug) << "tsOld " << tsOld << "   tsNow " << (mStartTSCalib > 0 ? mStartTSCalib : static_cast<int>(o2::ccdb::getCurrentTimestamp() / o2::ccdb::CcdbObjectInfo::HOUR)) << "   tsDiff " << tsDiff;
 
    if (EMCALCalibParams::Instance().requireSameRunType && runType != static_cast<int>(mRunType)) {
      LOG(info) << "adoptSavedData: Same run type required but run types differ: " << runType << " != " << static_cast<int>(mRunType);
      return false;
    }
    if (EMCALCalibParams::Instance().requireSameFill && fillNr != mFillNr) {
      LOG(info) << "adoptSavedData: Same fill nr. required but fills differ: " << fillNr << " != " << mFillNr;
      return false;
    }
    if (EMCALCalibParams::Instance().tsDiffMax > 0 && (EMCALCalibParams::Instance().tsDiffMax < tsDiff || tsDiff < 0)) {
      LOG(info) << "adoptSavedData: Maximum difference in ts is: " << EMCALCalibParams::Instance().tsDiffMax << " but " << tsDiff << " is given";
      return false;
    }
    if (EMCALCalibParams::Instance().useStaticStartTimeSlot && tsOldSec > 0) {
      slot.setStaticStartTimeMS(static_cast<long>(tsOldSec * o2::ccdb::CcdbObjectInfo::SECOND));
      LOG(info) << "adoptSavedData: Setting the start timestamp to " << static_cast<long>(tsOldSec * o2::ccdb::CcdbObjectInfo::SECOND);
    }
  }
 
  if constexpr (std::is_same<DataInput, o2::emcal::EMCALChannelData>::value) {
    TH2D* hEnergy = (TH2D*)fl.Get("EnergyVsCellID");
    TH2D* hTime = (TH2D*)fl.Get("TimeVsCellID");
    if (!hEnergy || !hTime) {
      return false;
    }
    auto hEnergyBoost = o2::utils::boostHistoFromRoot_2D<boostHisto2d>(hEnergy);
    auto hTimeBoost = o2::utils::boostHistoFromRoot_2D<boostHisto2d>(hTime);
 
    c->setHisto(hEnergyBoost);
    c->setHistoTime(hTimeBoost);
 
  } else if constexpr (std::is_same<DataInput, o2::emcal::EMCALTimeCalibData>::value) {
    TH2D* hTime = (TH2D*)fl.Get("TimeVsCellID");
    if (!hTime) {
      return false;
    }
    auto hTimeBoost = o2::utils::boostHistoFromRoot_2D<boostHisto2d>(hTime);
 
    c->setHisto(hTimeBoost);
  }
  TH1D* hEvents = (TH1D*)fl.Get("NEvents");
  if (!hEvents) {
    return false;
  }
  c->setNEvents(hEvents->GetBinContent(1));
  LOG(info) << "successfully loaded histograms with " << hEvents->GetBinContent(1) << " events";
 
  return true;
}
 
template <typename DataInput, typename DataOutput>
bool EMCALChannelCalibrator<DataInput, DataOutput>::setGainCalibrationFactors(o2::emcal::GainCalibrationFactors* gainCalibFactors)
{
 
  auto& cont = o2::calibration::TimeSlotCalibration<DataInput>::getSlots();
  if (cont.size() == 0) {
    return false; // time slot object not yet there
  }
 
  auto& slot = cont.at(0);
  DataInput* c = slot.getContainer();
  c->setGainCalibFactors(gainCalibFactors);
 
  return true;
}
 
} // end namespace emcal
} // end namespace o2
 
#endif /*EMCAL_CHANNEL_CALIBRATOR_H_ */