CbmStsChannelQA.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 "CbmStsChannelQA.h"
 
#include "CbmStsDigi.h"
#include "TArrow.h"
#include "TBox.h"
#include "TPaveLabel.h"
#include "TPaveText.h"
#include "TText.h"
 
#include <TParameter.h>
 
#include <ctime>
#include <iostream>
#include <typeinfo>
 
CbmStsChannelQA::CbmStsChannelQA(int dead_threshold) : fActiveMinEntries(dead_threshold) {}
CbmStsChannelQA::CbmStsChannelQA(int dead_threshold, double noise_threshold,
                                 std::optional<std::pair<size_t, size_t>> spill_on_off_threshold)
  : fActiveMinEntries(dead_threshold)
  , fSBThreshold(noise_threshold)
  , fSpillThresholds(spill_on_off_threshold)
{
  assert(noise_threshold > 0);
  if (fSpillThresholds.has_value()) {
    assert(fSpillThresholds->first <= fSpillThresholds->second);
    fSpillSections = {":ramp", ":spill_on", ":spill_off"};
  }
}
 
InitStatus CbmStsChannelQA::Init()
{
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman != nullptr) {
    fDigiManager = CbmDigiManager::Instance();
    fDigiManager->Init();
 
    if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) LOG(fatal) << GetName() << ": No StsDigi branch in input!";
 
    LOG(info) << "CbmStsChannelQA configuration:";
    LOG(info) << "Report lvl:                     " << fReportLvl;
    LOG(info) << "Active channel minimum entries: " << fActiveMinEntries;
    LOG(info) << "Noisy channel S/B thresholds:    " << fSBThreshold;
    if (fSpillThresholds.has_value()) {
      LOG(info) << Form("Spill OFF: RefDet_DigiRate < %lu", fSpillThresholds->first);
      LOG(info) << Form("Spill  ON: RefDet_DigiRate > %lu", fSpillThresholds->second);
    }
    else {
      LOG(info) << "No Spill section defined. Using all time slices";
    }
 
    return kSUCCESS;
  }
  return kERROR;
}
 
 
void CbmStsChannelQA::BookHistograms(int32_t address)
{
  LOG(debug) << Form("Booking histograms for module_addr: 0x%x", address);
 
  double dt_max  = 120;  // HARDCODE
  auto t_binning = cbm_sts_utils::HBinning{uint32_t(dt_max / 3.125), 0, dt_max};
 
  std::string h_name;
  for (auto modifier : fSpillSections) {
    h_name = Form("0x%x_charge_vs_channel%s", address, modifier);
    LOG(debug) << h_name;
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), 2048, 0, 2048, 31, 1, 32);
    fH2D[h_name]->GetXaxis()->SetTitle("Channel_{StsDigi}");
    fH2D[h_name]->GetYaxis()->SetTitle("Charge_{StsDigi} [ADC]");
 
    LOG(debug) << h_name;
    h_name       = Form("0x%x_dt_vs_charge%s", address, modifier);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), 31, 1, 32, t_binning.n_of_bins,
                                          t_binning.x_min, t_binning.x_max);
    fH2D[h_name]->GetYaxis()->SetTitle("Charge_{StsDigi} [ADC]");
    fH2D[h_name]->GetYaxis()->SetTitle("Time difference [ns]");
 
    LOG(debug) << h_name;
    h_name       = Form("0x%x_dt_vs_channel%s", address, modifier);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), 2048, 0, 2048, t_binning.n_of_bins,
                                          t_binning.x_min, t_binning.x_max);
    fH2D[h_name]->GetXaxis()->SetTitle("Channel_{StsDigi}");
    fH2D[h_name]->GetYaxis()->SetTitle("Time difference [ns]");
  }
 
  fAddressBook.insert(address);
}
 
void CbmStsChannelQA::ResetLastTime()
{
  for (auto& [address, time] : fLastDigiTime) {
    for (int chn = 0; chn < 2048; chn++) {
      time[chn] = -1;
    }
  }
}
 
void CbmStsChannelQA::Exec(Option_t*)
{
  LOG(info) << "Running CbmStsChannelQA ...";
 
  // Get Spill status
  size_t nb_ref_digis = fDigiManager->GetNofDigis(ECbmModuleId::kBmon);

  int fSpillStatus = 0;
  if (fSpillThresholds.has_value()) {
    fSpillStatus = nb_ref_digis <= fSpillThresholds->first ? -1 : (nb_ref_digis <= fSpillThresholds->second ? 0 : 1);
  }
 
  std::string str_spill;
  switch (fSpillStatus) {
    case +1: {
      fNbTsSpillOn++;
      str_spill = fSpillSections[1];
      break;
    }
    case -1: {
      fNbTsSpillOff++;
      str_spill = fSpillSections[2];
      break;
    }
    default: {
      str_spill = fSpillSections[0];
      break;
    }
  }
 
  LOG(info) << Form("TS %d\t-\t Spill section: %s", entry_, str_spill.c_str());
 
 
  auto sts_digis_     = fDigiManager->GetArray<CbmStsDigi>();
  size_t nb_sts_digis = fDigiManager->GetNofDigis(ECbmModuleId::kSts);
  for (size_t sts_digi_idx = 0; sts_digi_idx < nb_sts_digis; sts_digi_idx++) {  // sts digi loop
    const CbmStsDigi* sts_digi = &sts_digis_[sts_digi_idx];
    int32_t sts_digi_addr      = sts_digi->GetAddress();
    int32_t sts_digi_chan      = sts_digi->GetChannel();
    int32_t sts_digi_char      = sts_digi->GetCharge();
    double sts_digi_time       = sts_digi->GetTime();
 
    if (!fAddressBook.count(sts_digi_addr)) {
      BookHistograms(sts_digi_addr);
 
      for (int chn = 0; chn < 2048; chn++) {
        fLastDigiTime[sts_digi_addr][chn] = -1;
      }
    }
 
    fH2D[Form("0x%x_charge_vs_channel%s", sts_digi_addr, str_spill.c_str())]->Fill(sts_digi_chan, sts_digi_char);
 
    if (fLastDigiTime[sts_digi_addr][sts_digi_chan] > 0) {
      fH2D[Form("0x%x_dt_vs_charge%s", sts_digi_addr, str_spill.c_str())]->Fill(
        sts_digi_char, sts_digi_time - fLastDigiTime[sts_digi_addr][sts_digi_chan]);
      fH2D[Form("0x%x_dt_vs_channel%s", sts_digi_addr, str_spill.c_str())]->Fill(
        sts_digi_chan, sts_digi_time - fLastDigiTime[sts_digi_addr][sts_digi_chan]);
    }
 
    fLastDigiTime[sts_digi_addr][sts_digi_chan] = sts_digi_time;
 
  }  // end sts digi loop
 
  ResetLastTime();
  entry_++;
}
 
void CbmStsChannelQA::CheckDeadChannels()
{
  const char* mod = fSpillThresholds.has_value() ? fSpillSections[1] : fSpillSections[0];
  for (auto& module_addr : fAddressBook) {
    auto h = (TH1D*) fH2D[Form("0x%x_charge_vs_channel%s", module_addr, mod)]->ProjectionX();
    for (int chn = 0; chn < 2048; chn++) {
      double entries = h->GetBinContent(chn + 1);
 
      if (entries < fActiveMinEntries) {
        fDeadChannelList[module_addr].push_back(chn);
      }
    }
  }
}
 
void CbmStsChannelQA::CheckNoisyChannels()
{
  for (auto& module_addr : fAddressBook) {
    TH1D* h_sgn = fH2D[Form("0x%x_charge_vs_channel:spill_on", module_addr)]->ProjectionX();
    TH1D* h_bkg = fH2D[Form("0x%x_charge_vs_channel:spill_off", module_addr)]->ProjectionX();
 
    h_sgn->Scale(1. / fNbTsSpillOn);
    h_bkg->Scale(1. / fNbTsSpillOff);
    h_sgn->Add(h_bkg, -1);
 
    for (int chn = 0; chn < 2048; chn++) {
      double sgn = h_sgn->GetBinContent(chn + 1);
      double bkg = h_bkg->GetBinContent(chn + 1);
 
      if (bkg == 0) continue;
 
      double sb_ratio = sgn / bkg;
 
      if (sb_ratio < fSBThreshold) {
        LOG(debug) << Form("Noisy channel 0x%x: %d, \tS/B: %0.4f", module_addr, chn, sb_ratio);
        fNoisyChannelList[module_addr].push_back(chn);
      }
    }
  }
}
 
void CbmStsChannelQA::GenerateReport()
{
  std::string f_name                 = "cbm_sts_channel_qa_report.root";
  std::unique_ptr<TFile> o_root_file = std::make_unique<TFile>(f_name.c_str(), "RECREATE");
  int pad_size_x                     = 1000;
  int pad_size_y                     = 1000;
  for (auto& module_addr : fAddressBook) {
    std::string c_name         = Form("0x%x_channel", module_addr);
    std::unique_ptr<TCanvas> c = std::make_unique<TCanvas>(c_name.c_str(), c_name.c_str(), pad_size_x, pad_size_y);
    TH1D* h_sgn                = fH2D[Form("0x%x_charge_vs_channel:spill_on", module_addr)]->ProjectionX();
    TH1D* h_bkg                = fH2D[Form("0x%x_charge_vs_channel:spill_off", module_addr)]->ProjectionX();
 
    LOG(debug) << Form("Building report for 0x%x", module_addr);
    h_sgn->Scale(1. / fNbTsSpillOn);
    h_bkg->Scale(1. / fNbTsSpillOff);
 
    h_sgn->SetLineColor(kRed);
    h_bkg->SetLineColor(kBlack);
 
    h_sgn->SetTitle("");
    h_bkg->SetTitle("");
 
    h_sgn->SetFillColorAlpha(kRed, 0.2);
    h_bkg->SetFillColorAlpha(kBlack, 0.2);
 
    double h_sgn_max = h_sgn->GetMaximum();
    double h_bkg_max = h_bkg->GetMaximum();
    double max       = std::max(h_sgn_max, h_bkg_max);
 
    h_sgn->SetMaximum(1.05 * max);
    h_bkg->SetMaximum(1.05 * max);
 
    c->cd(1);
    gPad->SetLogy(1);
    gPad->SetLeftMargin(0.15);
    gPad->SetBottomMargin(0.12);
    gPad->SetTopMargin(0.10);
    gPad->SetRightMargin(0.12);
 
    h_sgn->Draw("histo");
    h_bkg->Draw("histo same");
 
    // Draw a line for bad channels
    if (fReportLvl > 1) {
      LOG(debug) << Form("Drawing %d noisy channel markers: %lu", module_addr, fNoisyChannelList[module_addr].size());
      for (auto& chn : fNoisyChannelList[module_addr]) {
        TBox chn_box = TBox(chn + 0.2, 0.99 * max, chn + 0.8, 1.01 * max);
        chn_box.SetLineColor(kBlue);
        chn_box.SetFillColor(kBlue);
        chn_box.DrawClone();
      }
    }
 
    c->Write(Form("0x%x_channel:spill_on_vs_off.png", module_addr));
  }
}
 
void CbmStsChannelQA::Finish()
{
  CheckDeadChannels();
  if (fSpillThresholds.has_value()) CheckNoisyChannels();
 
  if (fReportLvl) GenerateReport();
 
  SaveToFile();
 
  fH1D.clear();
  fH2D.clear();
 
  // Write normalization factor to file
  std::unique_ptr<TFile> out_file = std::make_unique<TFile>("cbm_sts_channel_qa.root", "UPDATE");
  TParameter<int>("fNbTsSpillOn", fNbTsSpillOn).Write();
  TParameter<int>("fNbTsSpillOff", fNbTsSpillOff).Write();
 
  std::ofstream dead_info("sts_dead_channels.info");
  std::ofstream dead_par("sts_dead_channels.par");
  dead_info << Form("# Module\t Number of dead channels\n");
  int n_of_dead_chn = 0;
  for (auto& [module_addr, list] : fDeadChannelList) {
    dead_info << Form("0x%x\t%lu\n", module_addr, list.size());
    for (int& dead_chn : list) {
      dead_par << module_addr << "\t" << dead_chn << std::endl;
    }
    n_of_dead_chn += list.size();
  }
  dead_info << Form("# Total Number of dead channels: %d\n", n_of_dead_chn);
 
  std::ofstream noisy_info("sts_noisy_channels.info");
  std::ofstream noisy_par("sts_noisy_channels.par");
  noisy_info << Form("# Module\t Number of noisy channels\n");
  int n_of_noisy_chn = 0;
  for (auto& [module_addr, list] : fNoisyChannelList) {
    noisy_info << Form("0x%x\t%lu\n", module_addr, list.size());
    for (int& noisy_chn : list) {
      noisy_par << module_addr << "\t" << noisy_chn << std::endl;
    }
    n_of_noisy_chn += list.size();
  }
  noisy_info << Form("# Total Number of noisy channels: %d\n", n_of_noisy_chn);
}