CbmStsTimeCal.cxx

Go to the documentation of this file.

/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Dario Ramirez [committer] */
 
#include "CbmStsTimeCal.h"
 
#include "CbmBmonDigi.h"
#include "CbmFsdDigi.h"
#include "CbmMuchDigi.h"
#include "CbmMvdDigi.h"
#include "CbmRichDigi.h"
#include "CbmStsDigi.h"
#include "CbmTofDigi.h"
#include "CbmTrdDigi.h"
 
#include <TColor.h>
#include <TDirectory.h>
#include <TF1.h>
#include <TFitResult.h>
#include <TLine.h>
#include <TPaveText.h>
#include <TStyle.h>
 
#include <cassert>
#include <ctime>
#include <fstream>
#include <iostream>
#include <typeinfo>
 
#include <fmt/core.h>
#include <yaml-cpp/emitter.h>
#include <yaml-cpp/emittermanip.h>
#include <yaml-cpp/node/node.h>
 
CbmStsTimeCal::CbmStsTimeCal(ECbmModuleId ref_sys, double t_min, double t_max)
  : fTimeWindowMin(t_min)
  , fTimeWindowMax(t_max)
  , fRefSystem(ref_sys)
{
}
 
CbmStsTimeCal::CbmStsTimeCal(int run_id, ECbmModuleId ref_sys, double t_min, double t_max)
{
  fRunId         = run_id;
  fTimeWindowMin = t_min;
  fTimeWindowMax = t_max;
  fRefSystem     = ref_sys;
}
 
 
InitStatus CbmStsTimeCal::Init()
{
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman != nullptr) {
    fDigiManager = CbmDigiManager::Instance();
    fDigiManager->Init();
 
    fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
 
    if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) LOG(fatal) << GetName() << ": No StsDigi branch in input!";
    if (!fDigiManager->IsPresent(fRefSystem)) LOG(fatal) << GetName() << ": No " << fRefSystem << " branch in input!";
 
    LoadWalkFromFile();
 
    fOutputPath = fRunId > 0 ? Form("%d", fRunId) : ".";
    if (gSystem->AccessPathName(fOutputPath.c_str(), kFileExists)) {
      if (system(("mkdir -p " + fOutputPath).c_str())) {  // Check output folder
        LOG(warning) << "Could not create output directory\n Setting output path at current location:\n";
      }
    }
    fReportOffset = std::make_unique<TFile>(Form("%s/cbm_sts_time_cal_offset.root", fOutputPath.c_str()), "RECREATE");
    fReportFile   = std::make_unique<TFile>(Form("%s/cbm_sts_time_cal_report.root", fOutputPath.c_str()), "RECREATE");
    fReportFit    = std::make_unique<TFile>(Form("%s/cbm_sts_time_cal_fits.root", fOutputPath.c_str()), "RECREATE");
 
    return kSUCCESS;
  }
  return kERROR;
}
 
void CbmStsTimeCal::InitTimeWalkMap(int32_t address)
{
  for (unsigned int asic_idx = 0; asic_idx < 16; asic_idx++) {
    fWalkMapRaw[address][asic_idx] = std::vector<double>(32, 0);
    auto loaded_par                = fWalkMap.Get(address, asic_idx);
    if (loaded_par.size() == 31) {
      for (int adc = 1; adc <= 31; adc++) {
        fWalkMapRaw[address][asic_idx][adc] = loaded_par[adc - 1];
      }
    }
  }
}
 
void CbmStsTimeCal::SetWalkFile(std::string f_name)
{
  LOG(debug) << Form("Setting user define time calibration file: %s", f_name.c_str());
  fParFile = f_name;
}
 
void CbmStsTimeCal::LoadWalkFromFile()
{
  if (!fParFile.length()) return;
 
  if (TString(fParFile.c_str()).EndsWith(".yaml") || TString(fParFile.c_str()).EndsWith(".yml")) {
    LOG(info) << Form("Loading time calibration from parameter file: %s", fParFile.c_str());
    fWalkMap          = cbm::algo::yaml::ReadFromFile<cbm::algo::sts::WalkMap>(fParFile);
    fGlobalTimeOffset = fWalkMap.GetSystemTimeOffset();
    LOG(info) << "Current system offset: " << fGlobalTimeOffset << " ns";
  }
}
 
void CbmStsTimeCal::BookHistograms(int32_t address)
{
  LOG(debug) << Form("Booking histograms for module: 0x%x", address);
 
  std::string h_name, h_title;
  int nb_time_bins = (std::abs(fTimeWindowMax - fTimeWindowMin) + kStsClock) / kStsClock;
 
  int unit               = CbmStsAddress::GetElementId(address, kStsUnit);
  int ladd               = CbmStsAddress::GetElementId(address, kStsLadder);
  int modu               = CbmStsAddress::GetElementId(address, kStsModule);
  std::string smart_name = Form("U%d_L%d_M%d", unit, ladd, modu);
 
  for (int asic_idx = 0; asic_idx < 16; asic_idx++) {
    h_name  = Form("0x%x_%02d", address, asic_idx);
    h_title = Form("%s_%02d", smart_name.c_str(), asic_idx);
 
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_title.c_str(), nb_time_bins, fTimeWindowMin - 0.5 * kStsClock,
                             fTimeWindowMax + 0.5 * kStsClock, 31, 1, 32);
 
    fH2D[h_name]->GetXaxis()->SetTitle(
      Form("t_{%s} - t_{StsDigi} [ns]", std::string(cbm::util::ToString(fRefSystem)).c_str()));
    fH2D[h_name]->GetYaxis()->SetTitle("Charge_{StsDigi} [ADC]");
  }
 
  h_name       = Form("0x%x_dt_vs_channel", address);
  h_title      = smart_name;
  fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_title.c_str(), 2048, 0, 2048, nb_time_bins,
                                        fTimeWindowMin - 0.5 * kStsClock, fTimeWindowMax + 0.5 * kStsClock);
 
  fH2D[h_name]->GetYaxis()->SetTitle(
    Form("t_{%s} - t_{StsDigi} [ns]", std::string(cbm::util::ToString(fRefSystem)).c_str()));
  fH2D[h_name]->GetXaxis()->SetTitle("Channel_{StsDigi}");
 
  fAddressBook.insert(address);
}
 
void CbmStsTimeCal::Exec(Option_t*)
{
  LOG(info) << "Running CbmStsTimeCal - Entry " << entry_;
 
  if (fRefSystem == ECbmModuleId::kBmon) CheckTiming<CbmBmonDigi>();
  if (fRefSystem == ECbmModuleId::kTof) CheckTiming<CbmTofDigi>();
  if (fRefSystem == ECbmModuleId::kRich) CheckTiming<CbmRichDigi>();
 
  entry_++;
}
 
template<class Digi>
void CbmStsTimeCal::CheckTiming()
{
  auto sts_digis_ = fDigiManager->GetArray<CbmStsDigi>();
  auto ref_digis_ = fDigiManager->GetArray<Digi>();
 
  size_t nb_sts_digis = fDigiManager->GetNofDigis(ECbmModuleId::kSts);
  size_t nb_ref_digis = fDigiManager->GetNofDigis(Digi::GetSystem());
 
  LOG(info) << nb_sts_digis << "\t" << nb_ref_digis;
 
  // Loops over Ref Detector
  for (size_t ref_digi_idx = 0; ref_digi_idx < nb_ref_digis - 1; ref_digi_idx++) {
    const Digi* ref_digi = &ref_digis_[ref_digi_idx];
    double ref_digi_time = ref_digi->GetTime();
 
    if (fAnalysisCuts != nullptr
        && !fAnalysisCuts->Check(CbmCutId::kBmonDigiSide, CbmTofAddress::GetChannelSide(ref_digi->GetAddress())))
      continue;
 
    double first_digi_in_window = ref_digi_time + fTimeWindowMin;
    double last_digi_in_window  = ref_digi_time + fTimeWindowMax;
    // Find first_digi_in_window
    int lo = 0, hi = nb_sts_digis - 1;
    while (lo < hi) {
      int m = (lo + hi) / 2;
      if (sts_digis_[m].GetTime() >= first_digi_in_window) {
        hi = m;
      }
      else {
        lo = m + 1;
      }
    }
    int lower_hitB = lo;
 
    // Find last_digi_in_window
    lo = 0, hi = nb_sts_digis - 1;
    while (lo < hi) {
      int m = (lo + hi + 1) / 2;
      if (sts_digis_[m].GetTime() <= last_digi_in_window) {
        lo = m;
      }
      else {
        hi = m - 1;
      }
    }
    int upper_hitB = lo;
    if (lower_hitB > upper_hitB) {
      LOG(debug) << "Not time match found for the current hit within the time window";
      continue;
    }
 
    for (int sts_digi_idx = lower_hitB; sts_digi_idx < upper_hitB; sts_digi_idx++) {
      const CbmStsDigi* sts_digi = &sts_digis_[sts_digi_idx];
      int32_t sts_digi_addr      = sts_digi->GetAddress();
 
      if (!fAddressBook.count(sts_digi_addr)) {
        BookHistograms(sts_digi_addr);
        InitTimeWalkMap(sts_digi_addr);
      }
 
      int32_t sts_digi_chan = sts_digi->GetChannel();
      int32_t sts_digi_char = sts_digi->GetCharge();
      int32_t asic_idx      = sts_digi_chan / 128;
      double sts_digi_time  = sts_digi->GetTime();
 
      fH2D[Form("0x%x_%02d", sts_digi_addr, asic_idx)]->Fill(ref_digi_time - sts_digi_time, sts_digi_char);
      fH2D[Form("0x%x_dt_vs_channel", sts_digi_addr)]->Fill(sts_digi_chan, ref_digi_time - sts_digi_time);
    }  // end sts loop
  }    // end ref loop
}
 
 
struct FitStsTime {
  std::unique_ptr<TH1> h_x;
  std::unique_ptr<TF1> f_x;
  std::unique_ptr<TF1> f_peak;
  std::unique_ptr<TF1> f_bkgd;
  TFitResultPtr fit_ptr;
  bool use_fit_result{false};
  double time_offset{0};
 
  TCanvas* Draw()
  {
    if (h_x == nullptr) {
      return nullptr;
    }
    int pad_x = 1024;
    int pad_y = 1024;
 
    std::unique_ptr<TCanvas> canvas = std::make_unique<TCanvas>("_c0", "", pad_x, pad_y);
    canvas->cd();
 
    double line_width     = 0.06;
    double pave0_top_edge = 0.95;
    double pave0_bot_edge = pave0_top_edge - line_width * 3;
    double pave1_top_edge = pave0_bot_edge;
    double pave1_bot_edge = pave1_top_edge - line_width * 3;
 
    h_x->SetLineStyle(1);
    h_x->SetLineWidth((0.002 * pad_y));
    h_x->SetLineColor(kBlack);
    h_x->DrawClone("histo");
 
    if (use_fit_result) {
      f_x->SetTitle("Fitted function");
      f_x->SetLineStyle(1);
      f_x->SetLineWidth((0.002 * pad_y));
      f_x->SetLineColor(kRed);
      f_x->DrawClone("same");
    }
 
    if (use_fit_result) {
      f_peak->SetTitle("Peak");
      f_peak->SetLineStyle(kDotted);
      f_peak->SetLineWidth((0.002 * pad_y));
      f_peak->SetLineColor(kBlue);
      f_peak->DrawClone("same");
      std::unique_ptr<TPaveText> pave = std::make_unique<TPaveText>(0.16, pave0_bot_edge, 0.45, pave0_top_edge, "NDC");
      pave->SetBorderSize(1);
      pave->SetFillColor(0);
      pave->SetTextAlign(12);
      pave->AddText(Form("p0 = %.3f", f_peak->GetParameter(0)))->SetTextColor(kBlue);
      pave->AddText(Form("p1 = %.3f", f_peak->GetParameter(1)))->SetTextColor(kBlue);
      pave->AddText(Form("p2 = %.3f", f_peak->GetParameter(2)))->SetTextColor(kBlue);
      pave->DrawClone();
    }
 
    if (use_fit_result) {
      f_bkgd->SetTitle("Background");
      f_bkgd->SetLineStyle(kDotted);
      f_bkgd->SetLineWidth((0.002 * pad_y));
      f_bkgd->SetLineColor(kOrange);
      f_bkgd->DrawClone("same");
      std::unique_ptr<TPaveText> pave = std::make_unique<TPaveText>(0.16, pave1_bot_edge, 0.45, pave1_top_edge, "NDC");
      pave->SetBorderSize(1);
      pave->SetFillColor(0);
      pave->SetTextAlign(12);
      pave->AddText(Form("p0 = %.3f", f_bkgd->GetParameter(0)))->SetTextColor(kOrange);
      pave->AddText(Form("p1 = %.3f", f_bkgd->GetParameter(1)))->SetTextColor(kOrange);
      pave->AddText(Form("p2 = %.3f", f_bkgd->GetParameter(2)))->SetTextColor(kOrange);
      pave->DrawClone();
    }
 
    TLine offset(time_offset, h_x->GetMinimum(), time_offset, h_x->GetMaximum());
    offset.SetLineWidth(0.004 * pad_y);
    offset.SetLineColor(kGray);
    offset.SetLineStyle(kDotted);
    offset.DrawClone();
 
    // Fit convergence
    std::unique_ptr<TPaveText> fit_pave = std::make_unique<TPaveText>(0.65, 0.75, 0.95, 0.95, "NDC");
    fit_pave->SetBorderSize(1);
    fit_pave->SetFillColor(0);
    fit_pave->SetTextAlign(12);
    fit_pave->AddText(Form("NDF   = %d", fit_ptr->Ndf()));
    fit_pave->AddText(Form("Chi2  = %.3f", fit_ptr->Chi2()));
    fit_pave->AddText(Form("PRob  = %E", fit_ptr->Prob()));
    fit_pave->DrawClone();
 
    return (TCanvas*) canvas->Clone();
  }
 
  FitStsTime(TH1D* h)
  {
    if (h == nullptr) {
      return;
    }
    h_x = std::make_unique<TH1D>(*h);
    h_x->Scale(1. / h_x->Integral());
 
    // Fit configuration
    double time_resolution = 4.8;
    double h_max           = h_x->GetBinCenter(h_x->GetMaximumBin());
    double fit_x_min       = h_max - 10. * time_resolution;
    double fit_x_max       = h_max + 10. * time_resolution;
 
    f_x    = std::make_unique<TF1>("f_x", "[0]*TMath::Gaus(x, [1], [2]) + [3] + [4]*x + [5]*x*x", fit_x_min, fit_x_max);
    f_peak = std::make_unique<TF1>("f_peak", "[0]*TMath::Gaus(x, [1], [2])", fit_x_min, fit_x_max);
    f_bkgd = std::make_unique<TF1>("f_bkgd", "[0] + [1]*x + [2]*x*x", fit_x_min, fit_x_max);
 
    // Initial values
    f_x->SetParameter(0, 0.5);
    f_x->SetParameter(1, h_max);
    f_x->SetParameter(2, time_resolution);
    f_x->SetParameter(5, -1e-4);
 
    // Limits
    f_x->SetParLimits(0, 0.0, 1.0);
    f_x->SetParLimits(1, h_max - 2. * time_resolution, h_max + 2. * time_resolution);
    f_x->SetParLimits(2, 3.125, 25);
    f_x->SetParLimits(5, -1, 0.0);
 
 
    fit_ptr = h_x->Fit("f_x", "SQN0", "", fit_x_min, fit_x_max);
 
    f_x->SetParameters(fit_ptr->GetParams());
    f_peak->SetParameter(0, fit_ptr->Parameter(0));
    f_peak->SetParameter(1, fit_ptr->Parameter(1));
    f_peak->SetParameter(2, fit_ptr->Parameter(2));
    f_bkgd->SetParameter(0, fit_ptr->Parameter(3));
    f_bkgd->SetParameter(1, fit_ptr->Parameter(4));
    f_bkgd->SetParameter(2, fit_ptr->Parameter(5));
 
    bool par_out_of_limits =
      fit_ptr->Parameter(1) < h_max - 2. * time_resolution || h_max + 2. * time_resolution < fit_ptr->Parameter(1);
 
    use_fit_result = !fit_ptr && fit_ptr->IsValid() && !par_out_of_limits;
 
    if (use_fit_result) {
      time_offset = f_peak->GetParameter(1);
    }
    else {
      time_offset = h_max;
    }
  }
};
 
double CbmStsTimeCal::FindModuleOffset(int32_t address)
{
 
  double module_offset = 0;
  TH2D* dt_vs_chn      = fH2D[Form("0x%x_dt_vs_channel", address)].get();
  if (dt_vs_chn != nullptr) {
    TH1D* p_y = (TH1D*) dt_vs_chn->ProjectionY();
 
    double peak = p_y->GetBinCenter(p_y->GetMaximumBin());
 
    FitStsTime offset_fit(p_y);
    TCanvas* drawing = offset_fit.Draw();
    if (drawing != nullptr) {
      fReportOffset->cd();
      drawing->Write(Form("0x%x_offset", address));
    }
 
    if (offset_fit.use_fit_result) {
      module_offset = offset_fit.time_offset;
    }
    else {
      module_offset = peak;
    }
  }
  return module_offset;
}
 
void CbmStsTimeCal::CheckTimeWalk()
{
  /* This function can be optimized by parallelizing the loop over modules or ASIC idx, depending on the ratio
                * For mCBM case, there are usually less than 16 modules so it makes more sense to parallelize the ASIC idx loop
                * but for the full CBM case, the amount of modules greatly exceeds the number of ASIC per modules (16)
                */
  for (int32_t module : fAddressBook) {  // loop over modules
    double module_offset = FindModuleOffset(module);
    for (unsigned int asic_idx = 0; asic_idx < 16; asic_idx++) {  // loop over ASICs
 
      double adc_channel[32];
      double adc_channel_err[32];
      double time_offset[32];
      double time_offset_err[32];
      int n_pts = 0;
 
      std::string h_name = Form("0x%x_%02d", module, asic_idx);
      if (!fH2D.count(h_name)) return;
 
      auto h = fH2D[h_name].get();
 
      for (int adc = 1; adc <= 31; adc++) {
        TH1D* h_x = h->ProjectionX(Form("%s_dt", h_name.c_str()), adc, adc);
 
        FitStsTime offset_fit(h_x);
        TCanvas* drawing = offset_fit.Draw();
        if (drawing != nullptr) {
          std::string fit_folder = Form("0x%x/%02d", module, asic_idx);
          TDirectory* dir        = (TDirectory*) fReportFit->Get(fit_folder.c_str());
          if (!dir) {
            dir = fReportFit->mkdir(fit_folder.c_str());
          }
          dir->cd();
          drawing->Write(Form("0x%x_%02d_%02d", module, asic_idx, adc));
        }
 
        if (offset_fit.use_fit_result) {
          time_offset[n_pts]     = offset_fit.fit_ptr->Parameter(1);
          time_offset_err[n_pts] = offset_fit.fit_ptr->Parameter(2);
        }
        else {
          /* If the fitting for individual ADC of the ASIC was unsuccessful
           * the time offset is taken as the one estimated for the whole module
           */
          time_offset[n_pts]     = module_offset;
          time_offset_err[n_pts] = h_x->GetRMS();
        }
        adc_channel[n_pts]     = adc;
        adc_channel_err[n_pts] = 0;
 
        // Update time walk map
        fWalkMapRaw[module][asic_idx][adc] += time_offset[n_pts];
 
        n_pts++;
 
      }  // ---- end ADC loop ----
 
      fG1D[h_name] = std::make_unique<TGraphErrors>(n_pts, adc_channel, time_offset, adc_channel_err, time_offset_err);
      fG1D[h_name]->SetTitle(
        Form("%s : StsDigi_0x%x time off-set", std::string(cbm::util::ToString(fRefSystem)).c_str(), module));
      fG1D[h_name]->GetXaxis()->SetTitle("Charge_{StsDigi} [ADC]");
      fG1D[h_name]->GetYaxis()->SetTitle(
        Form("t_{%s} - t_{StsDigi} [ns]", std::string(cbm::util::ToString(fRefSystem)).c_str()));
      fG1D[h_name]->SetMinimum(fTimeWindowMin);
      fG1D[h_name]->SetMaximum(fTimeWindowMax);
 
    }  // ---- end module loop ----
  }    // ---- end ASICs loop ----
}
 
double CbmStsTimeCal::FindGlobalOffset()
{
  double setup_offset  = 0;
  double min_sum_sigma = 999;
  for (int32_t module : fAddressBook) {  // module loop
    double module_avg_sigma  = 0;
    double module_avg_offset = 0;
    int n_of_values          = 0;
 
    for (unsigned int asic_idx = 0; asic_idx < 16; asic_idx++) {  // asic loop
      std::string h_name = Form("0x%x_%02d", module, asic_idx);
      auto g             = fG1D[h_name].get();
 
      int n_of_pto    = g->GetN();
      int lower_bound = n_of_pto > 5 ? n_of_pto - 5 : 0;
 
      // Calculate the average value of the sigma for the last ADC
      for (int pto = n_of_pto; pto > lower_bound - 5; pto--) {
        double offset = g->GetPointY(pto);
        double sigma  = g->GetErrorY(pto);
        if (sigma) {
          module_avg_sigma += sigma;
          module_avg_offset += offset;
          n_of_values++;
        }
      }
    }  // end asic loop
 
    module_avg_sigma  = n_of_values ? module_avg_sigma / n_of_values : -1;
    module_avg_offset = n_of_values ? module_avg_offset / n_of_values : -1;
 
    if (module_avg_sigma && module_avg_sigma < min_sum_sigma) {
      min_sum_sigma = module_avg_sigma;
      setup_offset  = module_avg_offset;
    }
  }  // end module loop
 
  return setup_offset;
}
 
void CbmStsTimeCal::WriteYAML()
{
  std::string o_yaml_file = fOutputPath + "/StsTimeCalibration.yaml";
  LOG(info) << "Writing parameter file to:" << o_yaml_file;
 
  YAML::Emitter walk_map;
  walk_map << YAML::BeginMap;
  walk_map << YAML::Key << "timeOffset";
  walk_map << YAML::Value << (int) fGlobalTimeOffset;
  walk_map << YAML::Key << "WalkMap";
  walk_map << YAML::Value << YAML::BeginMap;
  for (const auto& [module, asics] : fWalkMapRaw) {
    walk_map << YAML::Key << fmt::format("0x{:x}", module);
    walk_map << YAML::Value << YAML::BeginSeq;
    for (const auto& [asic, pars] : asics) {
      // Begin a flow sequence
      walk_map << YAML::Flow << YAML::BeginSeq;
      for (const auto& value : pars) {
        walk_map << fmt::format("{:+.3f}", value);  // Add '+' for positive numbers
      }
      // End the flow sequence
      walk_map << YAML::EndSeq;
    }
    walk_map << YAML::EndSeq;
  }
  walk_map << YAML::EndMap << YAML::EndMap;
  assert(walk_map.good());
 
  std::ofstream walk_map_out(o_yaml_file);
  assert(walk_map_out.is_open());
 
  walk_map_out << walk_map.c_str();
  walk_map_out.close();
}
 
void CbmStsTimeCal::Finish()
{
  CheckTimeWalk();
 
  // resize ADC offset vector to 31 channels from 0 [0-30`]
  for (auto& [module, asics] : fWalkMapRaw) {   // Loop over modules
    for (auto& [asic_idx, asic_par] : asics) {  // Loop over ASICs
      for (int adc = 1; adc <= 31; adc++) {     // loop over ADC channels 31 channels [1 - 31]
        asic_par[adc - 1] = asic_par[adc];
      }
      asic_par.pop_back();
    }
  }
 
  WriteYAML();
 
  DrawResults();
  SaveToFile();
 
 
  fH1D.clear();
  fH2D.clear();
}
 
 
void CbmStsTimeCal::DrawResults()
{
  std::unique_ptr<TCanvas> adc_vs_dt_1d =
    std::make_unique<TCanvas>("adc_vs_dt_1d", "adc_vs_dt_1d", 10, 10, 4 * 1024, 4 * 1024);
  adc_vs_dt_1d->Divide(4, 4);
 
  std::unique_ptr<TCanvas> adc_vs_dt_2d =
    std::make_unique<TCanvas>("adc_vs_dt_2d", "adc_vs_dt_2d", 10, 10, 4 * 1024, 4 * 1024);
  adc_vs_dt_2d->Divide(4, 4);
 
  std::unique_ptr<TCanvas> dt_vs_chn = std::make_unique<TCanvas>("dt_vs_chn", "dt_vs_chn", 10, 10, 1024, 1024);
 
  fReportFile->cd();
  for (int32_t sensor : fAddressBook) {
    for (int asic_idx = 0; asic_idx < 16; asic_idx++) {
      adc_vs_dt_1d->cd(asic_idx + 1);
      if (fG1D[Form("0x%x_%02d", sensor, asic_idx)] != nullptr) {
        fG1D[Form("0x%x_%02d", sensor, asic_idx)]->Draw("APL");
      }
 
      adc_vs_dt_2d->cd(asic_idx + 1);
      if (fH2D[Form("0x%x_%02d", sensor, asic_idx)] != nullptr) {
        fH2D[Form("0x%x_%02d", sensor, asic_idx)]->Draw("colz");
      }
    }
 
    dt_vs_chn->cd();
    if (fH2D[Form("0x%x_dt_vs_channel", sensor)] != nullptr) {
      fH2D[Form("0x%x_dt_vs_channel", sensor)]->Draw("colz");
    }
 
    adc_vs_dt_1d->Write(Form("adc_vs_dt_1d_0x%x", sensor));
    adc_vs_dt_2d->Write(Form("adc_vs_dt_2d_0x%x", sensor));
    dt_vs_chn->Write(Form("dt_vs_chn_0x%x", sensor));
  }
}