DQMHistAnalysisCDCEpics.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#include <dqm/analysis/modules/DQMHistAnalysisCDCEpics.h>
#include <cdc/geometry/CDCGeometryPar.h>
 
#include <TLatex.h>
 
using namespace std;
using namespace Belle2;
 
//-----------------------------------------------------------------
//                 Register module
//-----------------------------------------------------------------
REG_MODULE(DQMHistAnalysisCDCEpics);
 
DQMHistAnalysisCDCEpicsModule::DQMHistAnalysisCDCEpicsModule()
  : DQMHistAnalysisModule()
{
  addParam("HistDirectory", m_name_dir, "CDC Dir of DQM Histogram", std::string("CDC"));
  addParam("RefDirectory", m_name_refdir, "Reference histogram dir", std::string("ref/CDC/default"));
  addParam("PvPrefix", m_name_pvpfx, "PV Prefix and Name", std::string("CDC:"));
  addParam("RefFilePhi", m_fname_refphi, "Reference histogram file name", std::string("CDCDQM_PhiRef.root"));
  addParam("HistLayerADC", m_hname_ladc, "Layer ADC Histogram Name", std::string("hADCLayer"));
  addParam("HistBoardADC", m_hname_badc, "Board ADC Histogram Name", std::string("hADCBoard"));
  addParam("HistBoardTDC", m_hname_btdc, "Board TDC Histogram Name", std::string("hTDC"));
  addParam("HistPhiIndex", m_hname_idxphi, "Phi Index Histogram Name", std::string("hPhiIndex"));
  addParam("HistPhiEff", m_hname_effphi, "Phi Eff Histogram Name", std::string("hPhiEff"));
  addParam("HistHitsPhi", m_hname_hitsphi, "Phi Hits Histogram Name", std::string("hPhiNCDC"));
  addParam("MinEvt", m_minevt, "Min events for intra-run point", 1000);
  addParam("HistTrackingWireEff", m_histoTrackingWireEff, "Wire Eff Histogram Name", std::string("hTrackingWireEff"));
  addParam("DoTH2PolyTrackingWireEff", m_doTH2PolyTrackingWireEff,
           "If true, creates TH2Poly instead of TH2F for TrackingWireEff Histos", m_doTH2PolyTrackingWireEff);
  addParam("FirstEffBoundary", m_firstEffBoundary, "The first boundary of the efficiency range", m_firstEffBoundary);
  addParam("SecondEffBoundary", m_secondEffBoundary, "The second boundary of the efficiency range", m_secondEffBoundary);
  for (int i = 0; i < 300; i++) {
    m_hists_bADC[i] = nullptr;
    m_hists_bTDC[i] = nullptr;
    //layer
    if (i < 56)m_hists_lADC[i] = nullptr;
  }
 
  B2DEBUG(20, "DQMHistAnalysisCDCEpics: Constructor done.");
}
 
DQMHistAnalysisCDCEpicsModule::~DQMHistAnalysisCDCEpicsModule()
{
}
 
void DQMHistAnalysisCDCEpicsModule::initialize()
{
 
  gROOT->cd();
  c_histmd_ladc = new TCanvas("CDC/c_histmd_ladc", "c_histmd_ladc", 500, 400);
  m_histmd_ladc = new TH1F("CDC/histmd_ladc", "m_histmd_ladc", 56, 0, 56);
  m_histmd_ladc->SetTitle("ADC Medians vs Layers (SL-lines); CDC Layer index; ADC medians");
 
  c_hist_adc = new TCanvas("CDC/c_hist_adc", "c_hist_adc", 500, 400);
  m_hist_adc = new TH1F("CDC/hist_adc", "m_hist_adc", 300, 0, 300);
  m_hist_adc->SetTitle("ADC Medians; CDC board index; ADC medians");
 
  c_hist_tdc = new TCanvas("CDC/c_hist_tdc", "c_hist_tdc", 500, 400);
  m_hist_tdc = new TH1F("CDC/hist_tdc", "m_hist_tdc", 300, 0, 300);
  m_hist_tdc->SetTitle("TDC Medians; CDC board index; TDC medians");
 
  //array of various phi histograms
  for (int ic = 0; ic < 8; ic++) {
    c_hist_skimphi[ic] = new TCanvas(Form("CDC/c_hist_skimphi_c%d", ic), Form("hist_skimphi_c%d", ic), 500, 400);
  }
 
  c_hist_crphi = new TCanvas("CDC/c_hist_crphi", "c_hist_crphi", 500, 400);
  c_hist_hitsphi = new TCanvas("CDC/c_hist_hitsphi", "c_hist_hitsphi", 500, 400);
 
  //CR alarm reference
  if (m_fname_refphi != "") {
    m_fileRefPhi = TFile::Open(m_fname_refphi.data(), "READ");
    if (m_fileRefPhi && m_fileRefPhi->IsOpen()) {
      B2INFO("DQMHistAnalysisCDCEpics: reference (" << m_fname_refphi << ") found OK");
      m_histref_phiindex = (TH2F*)m_fileRefPhi->Get((m_name_refdir + "/hPhiIndex").data());
      if (!m_histref_phiindex)B2INFO("\t .. but (histogram) not found");
      else B2INFO("\t ..and (cdcdqm_phiref) also exist");
    }
  }
 
  m_lines.clear();
  m_lines.reserve(kNumLayers);
  for (unsigned il = 0; il < kNumLayers; ++il) {
    int bin = il + 1;
    auto* line = new TLine(bin, 0, bin, 1.0);
    line->SetLineStyle(2);
    if (bin >= 8 && bin < 14)line->SetLineColor(kRed);  // U-type
    else if (bin >= 20 && bin < 26)line->SetLineColor(kGreen);  // V-type
    else if (bin >= 32 && bin < 38)line->SetLineColor(kRed);  // U-type
    else if (bin >= 44 && bin < 50)line->SetLineColor(kGreen);  // v-type
    else line->SetLineColor(kGray); // A-type
    m_lines.push_back(line);
  }
 
  c_hist_effphi = new TCanvas("CDC/c_hist_effphi", "c_hist_effphi", 500, 400);
  m_hist_effphi = new TH1D("CDC/hist_effphi", "m_hist_effphi", 360, -180.0, 180.0);
 
  c_hist_attach_eff[0] = new TCanvas("CDC/c_hist_attached_wires", "c_hist_attached_wires", 403, 400);
  c_hist_attach_eff[1] = new TCanvas("CDC/c_hist_expected_wires", "c_hist_expected_wires", 403, 400);
  c_hist_attach_eff[2] = new TCanvas("CDC/c_hist_attach_eff", "c_hist_attach_eff", 403, 400);
  c_hist_attach_eff[3] = new TCanvas("CDC/c_hist_attach_eff_1d", "c_hist_attach_eff_1d", 403, 400);
  if (m_doTH2PolyTrackingWireEff) {
    m_hist_attach_eff_Poly[0] = createEffiTH2Poly("CDC/hist_attachedWires",
                                                  "hist_attachedWires (backplate view);X [cm];Y [cm]; Track / bin");
    m_hist_attach_eff_Poly[0]->GetYaxis()->SetTitleOffset(1.4);
    m_hist_attach_eff_Poly[0]->SetDirectory(gDirectory);
    m_hist_attach_eff_Poly[1] = (TH2Poly*)m_hist_attach_eff_Poly[0]->Clone();
    m_hist_attach_eff_Poly[1]->SetNameTitle("CDC/hist_expectedWires", "hist_expectedWires (backplate view);X [cm];Y [cm]; Track / bin");
    m_hist_attach_eff_Poly[1]->SetDirectory(gDirectory);
    m_hist_attach_eff_Poly[2] = (TH2Poly*)m_hist_attach_eff_Poly[0]->Clone();
    m_hist_attach_eff_Poly[2]->SetNameTitle("CDC/hist_wireAttachEff", "hist_wireAttachEff (backplate view);X [cm];Y [cm]; Efficiency");
    m_hist_attach_eff_Poly[2]->SetDirectory(gDirectory);
  } else {
    static CDC::CDCGeometryPar& cdcgeo = CDC::CDCGeometryPar::Instance();
    int nSLayers = cdcgeo.getNumberOfSenseLayers();
    double maxLayerR = cdcgeo.senseWireR(nSLayers - 1);
    m_hist_attach_eff[0] = new TH2F("CDC/hist_attachedWires", "hist_attachedWires (backplate view);X [cm];Y [cm]; Track / bin",
                                    nSLayers * 6, -maxLayerR * 1.02, maxLayerR * 1.02,
                                    nSLayers * 6, -maxLayerR * 1.02, maxLayerR * 1.02);
    m_hist_attach_eff[0]->GetYaxis()->SetTitleOffset(1.4);
    m_hist_attach_eff[1] = (TH2F*)m_hist_attach_eff[0]->Clone();
    m_hist_attach_eff[1]->SetNameTitle("CDC/hist_expectedWires", "hist_expectedWires (backplate view);X [cm];Y [cm]; Track / bin");
    m_hist_attach_eff[2] = (TH2F*)m_hist_attach_eff[0]->Clone();
    m_hist_attach_eff[2]->SetNameTitle("CDC/hist_wireAttachEff", "hist_wireAttachEff (backplate view);X [cm];Y [cm]; Efficiency");
  }
  m_hist_wire_attach_eff_1d = new TH1F("CDC/hist_wire_attach_eff_1d", "hist_wire_attach_eff_1d;Wire Efficiency;Wire / bin",
                                       208, -0.02, 1.02);
  m_hist_wire_attach_eff_1d->GetYaxis()->SetTitleOffset(1.4);
 
  if (!hasDeltaPar(m_name_dir, m_hname_ladc))
    addDeltaPar(m_name_dir, m_hname_ladc, HistDelta::c_Entries, m_minevt, 1);
 
  if (!hasDeltaPar(m_name_dir, m_hname_badc))
    addDeltaPar(m_name_dir, m_hname_badc, HistDelta::c_Entries, m_minevt, 1);
 
  if (!hasDeltaPar(m_name_dir, m_hname_btdc))
    addDeltaPar(m_name_dir, m_hname_btdc, HistDelta::c_Entries, m_minevt, 1);
 
  if (!hasDeltaPar(m_name_dir, m_hname_idxphi))
    addDeltaPar(m_name_dir, m_hname_idxphi, HistDelta::c_Entries, m_minevt, 1);
 
  if (!hasDeltaPar(m_name_dir, m_hname_effphi))
    addDeltaPar(m_name_dir, m_hname_effphi, HistDelta::c_Entries, m_minevt, 1);
 
  if (!hasDeltaPar(m_name_dir, m_hname_hitsphi))
    addDeltaPar(m_name_dir, m_hname_hitsphi, HistDelta::c_Entries, m_minevt, 1);
 
  if (!hasDeltaPar(m_name_dir, m_histoTrackingWireEff))
    addDeltaPar(m_name_dir, m_histoTrackingWireEff, HistDelta::c_Events, m_minevt, 1);
 
  //creating box for normal adc and tdc windows, the real position is updated at begin run
  m_line_ladc = new TLine(0, m_minadc, 300, m_minadc);
  m_line_ladc->SetLineColor(kRed);
  m_line_ladc->SetLineWidth(2);
 
  m_line_hadc = new TLine(0, m_maxadc, 300, m_maxadc);
  m_line_hadc->SetLineColor(kRed);
  m_line_hadc->SetLineWidth(2);
 
  m_line_ltdc = new TLine(0, m_mintdc, 300, m_mintdc);
  m_line_ltdc->SetLineColor(kRed);
  m_line_ltdc->SetLineWidth(2);
 
  m_line_htdc = new TLine(0, m_maxtdc, 300, m_maxtdc);
  m_line_htdc->SetLineColor(kRed);
  m_line_htdc->SetLineWidth(2);
 
  registerEpicsPV(m_name_pvpfx + "cdcboards_wadc", "adcboards");
  registerEpicsPV(m_name_pvpfx + "cdcboards_wtdc", "tdcboards");
 
  registerEpicsPV(m_name_pvpfx + "adc_median_window", "adcmedianwindow");
  registerEpicsPV(m_name_pvpfx + "tdc_median_window", "tdcmedianwindow");
 
  registerEpicsPV(m_name_pvpfx + "phi_compare_window", "phicomparewindow");
 
  m_monObj = getMonitoringObject("cdc");
 
  B2DEBUG(20, "DQMHistAnalysisCDCEpics: initialized.");
}
 
void DQMHistAnalysisCDCEpicsModule::beginRun()
{
  double unused = 0;
  requestLimitsFromEpicsPVs("adcmedianwindow", unused, m_minadc, m_maxadc, unused);
  requestLimitsFromEpicsPVs("tdcmedianwindow", unused, m_mintdc, m_maxtdc, unused);
  requestLimitsFromEpicsPVs("phicomparewindow", m_phialarm, m_phiwarn, unused, unused);
 
  //in case if something is wrong in config file
  if (std::isnan(m_minadc)) m_minadc = 60.0;
  if (std::isnan(m_maxadc)) m_maxadc = 130.0;
  if (std::isnan(m_mintdc)) m_mintdc = 4600.0;
  if (std::isnan(m_maxtdc)) m_maxtdc = 5000.0;
 
  if (std::isnan(m_phiwarn)) m_phiwarn = 0.05; //>%5 is warning
  if (std::isnan(m_phialarm)) m_phialarm = 0.15; //>%15 is warning
 
  // Update Line position from Epics limits
  m_line_ladc->SetY1(m_minadc);
  m_line_ladc->SetY2(m_minadc);
  m_line_hadc->SetY1(m_maxadc);
  m_line_hadc->SetY2(m_maxadc);
  m_line_ltdc->SetY1(m_mintdc);
  m_line_ltdc->SetY2(m_mintdc);
  m_line_htdc->SetY1(m_maxtdc);
  m_line_htdc->SetY2(m_maxtdc);
 
  B2DEBUG(20, "DQMHistAnalysisCDCEpics: beginRun run called");
}
 
void DQMHistAnalysisCDCEpicsModule::event()
{
  //1. get adc median vs layer numbers
  auto m_delta_ladc = (TH2F*)getDelta(m_name_dir, m_hname_ladc, 0, true);
  if (m_delta_ladc) {
    m_histmd_ladc->Reset();
    for (unsigned il = 0; il < kNumLayers; ++il) {
      if (m_hists_lADC[il]) delete m_hists_lADC[il];
      m_hists_lADC[il] = m_delta_ladc->ProjectionY(Form("histmd_adc_layer%d", il + 1), il + 1, il + 1, "");
      m_hists_lADC[il]->SetTitle(Form("histmd_adc_layer%d", il));
      float md_ladc = getHistMedian(m_hists_lADC[il]);
      m_histmd_ladc->SetBinContent(il + 1, md_ladc);
    }
    // Draw canvas
    c_histmd_ladc->Clear();
    c_histmd_ladc->cd();
    getHistStyle(m_histmd_ladc, "layeradc", 0);
    double y_max = m_histmd_ladc->GetMaximum();
    m_histmd_ladc->SetFillColor(kYellow);
    m_histmd_ladc->SetMinimum(0);
    m_histmd_ladc->SetMaximum(y_max * 1.20);
    m_histmd_ladc->Draw("hist");
    for (auto* line : m_lines) {
      line->SetY2(y_max * 1.20);
      line->Draw("same");
    }
    c_histmd_ladc->Update();
    UpdateCanvas(c_histmd_ladc);
  }
 
  //2. get adc medians vs board ID
  auto m_delta_adc = (TH2F*)getDelta(m_name_dir, m_hname_badc, 0, true); //true=only if updated
  if (m_delta_adc) {
    m_hist_adc->Reset();
    int cadcgood = 0;
    int cadcbad = 0;
    double sumadcgood = 0;
    for (unsigned ic = 0; ic < kNumBoards; ++ic) {
      if (ic == 0) continue; //299 boards only
      if (m_hists_bADC[ic]) delete m_hists_bADC[ic];
      m_hists_bADC[ic] = m_delta_adc->ProjectionY(Form("histmd_tdc_board%d", ic + 1), ic + 1, ic + 1, "");
      m_hists_bADC[ic]->SetTitle(Form("histmd_adc_board%d", ic));
      float md_adc = getHistMedian(m_hists_bADC[ic]);
      m_hist_adc->SetBinContent(ic + 1, md_adc);
      if (md_adc >= m_minadc && md_adc <= m_maxadc) {
        sumadcgood = sumadcgood + md_adc;
        cadcgood++;
      } else cadcbad++;
    }
    double adcfrac = cadcgood / 2.99; // (100.0/299) in %
    setEpicsPV("adcboards", adcfrac);
    // Draw canvas
    c_hist_adc->Clear();
    c_hist_adc->cd();
    if (cadcgood > 0)sumadcgood = sumadcgood * 1.0 / cadcgood;
    getHistStyle(m_hist_adc, "adc", sumadcgood);
    m_hist_adc->SetTitle(Form("ADC Medians: Bad board count = %d (%0.01f%%)", cadcbad - 1, 100.0 - adcfrac));
    m_hist_adc->Draw("");
    m_line_ladc->Draw("same");
    m_line_hadc->Draw("same");
    c_hist_adc->Update();
    UpdateCanvas(c_hist_adc);
  }
 
  //3. get tdc medians vs board ID
  auto m_delta_tdc = (TH2F*)getDelta(m_name_dir, m_hname_btdc, 0, true);
  if (m_delta_tdc) {
    m_hist_tdc->Reset();
    int ctdcgood = 0;
    int ctdcbad = 0;
    double sumtdcgood = 0;
    for (unsigned ic = 0; ic < kNumBoards; ++ic) {
      if (ic == 0) continue; //299 boards only
      if (m_hists_bTDC[ic]) delete m_hists_bTDC[ic];
      m_hists_bTDC[ic] = m_delta_tdc->ProjectionY(Form("histmd_tdc_board%d", ic + 1), ic + 1, ic + 1, "");
      m_hists_bTDC[ic]->SetTitle(Form("histmd_tdc_board%d", ic));
      float md_tdc = getHistMedian(m_hists_bTDC[ic]);
      m_hist_tdc->SetBinContent(ic + 1, md_tdc);
      if (md_tdc >= m_mintdc && md_tdc <= m_maxtdc) {
        ctdcgood++;
        sumtdcgood = sumtdcgood + md_tdc;
      } else ctdcbad++;
    }
    double tdcfrac = ctdcgood / 2.99;
    setEpicsPV("tdcboards", tdcfrac);
    c_hist_tdc->Clear();
    c_hist_tdc->cd();
    if (ctdcgood > 0)sumtdcgood = sumtdcgood * 1.0 / ctdcgood;
    getHistStyle(m_hist_tdc, "tdc", sumtdcgood);
    m_hist_tdc->SetTitle(Form("TDC Medians: Bad board count = %d (%0.01f%%)", ctdcbad - 1, 100.0 - tdcfrac));
    m_hist_tdc->Draw("");
    m_line_ltdc->Draw("same");
    m_line_htdc->Draw("same");
    c_hist_tdc->Update();
    UpdateCanvas(c_hist_tdc);
  }
 
  //get phi plots for various options
  auto m_delta_skimphi = (TH2F*)getDelta(m_name_dir, m_hname_idxphi, 0, true); //true=only if updated
  if (m_delta_skimphi) {
    TString sip[2] = {"OffIP", "IP"};
    TString sname[4] = {"all", "bhabha", "hadron", "mumutrk"};
    for (int j = 0; j < 2; j++) { //ip selections
      for (int i = 0; i < 4; i++) { //skim selections
        int k = 4 * j + i; //0 to 7
        TString hname = TString::Format("histphi_%s_%sevt", sip[j].Data(), sname[i].Data());
        m_hist_skimphi[k] = m_delta_skimphi->ProjectionX(hname, k + 1, k + 1, "");
        m_hist_skimphi[k]->SetTitle(TString::Format("cdc-track #phi (%s, %s-events);#phi;entries", sip[j].Data(), sname[i].Data()));
        if (k < 4)m_hist_skimphi[k]->SetFillColor(kGray);
        else m_hist_skimphi[k]->SetFillColor(kCyan);
        c_hist_skimphi[k]->Clear();
        c_hist_skimphi[k]->cd();
        gPad->SetGridx(1);
        gPad->SetGridy(1);
        m_hist_skimphi[k]->Draw("hist");
      }
    }
  }
 
  //for CR shifter IP + all hadrons including alarm system
  if (m_delta_skimphi) {
    c_hist_crphi->Clear();
    bool isFew = false, isAlarm = false, isWarn = false;
    m_hist_crphi = m_delta_skimphi->ProjectionX("histphi_ip_hadrons", 7, 7, "");
    m_hist_crphi->SetTitle("cdc-track #phi (IP + hadrons);cdc-track #phi;norm entries");
    if (m_hist_crphi) {
      double maxnow = m_hist_crphi->Integral();
      if (maxnow > 0)m_hist_crphi->Scale(1.0 / maxnow);
      if (maxnow < 10000) {
        isFew = true;
      } else {
        if (m_histref_phiindex) {
          m_hist_refphi = m_histref_phiindex->ProjectionX("histphi_ip_hadronsref", 7, 7, "");
          double nbinref = m_hist_refphi->GetNbinsX();
          double nbinnow = m_hist_crphi->GetNbinsX();
          if (nbinref == nbinnow) { //same bins
            double maxref = m_hist_refphi->Integral();
            if (maxref > 0) {
              m_hist_refphi->Scale(1.0 / maxref);
              double maxphidiff = 0;
              double maxphidiff_angle = 0;
              for (int iphi = 0; iphi < nbinnow; iphi++) {
                double icnow = m_hist_crphi->GetBinContent(iphi + 1);
                double icref = m_hist_refphi->GetBinContent(iphi + 1);
                double phidiff = fabs(icnow - icref);
                if (phidiff > m_phiwarn)isWarn = true;
                if (phidiff > m_phialarm)isAlarm = true;
                if (phidiff > maxphidiff) {
                  maxphidiff = phidiff;
                  maxphidiff_angle = m_hist_crphi->GetBinLowEdge(iphi + 1) + m_hist_crphi->GetBinWidth(iphi + 1);
                }
              }
              m_hist_crphi->SetTitle(Form("%s (diff = %0.03f at %0.1f)", m_hist_crphi->GetTitle(), maxphidiff, maxphidiff_angle));
            }
          }
        }
      }
    }
    c_hist_crphi->cd();
    gPad->SetGridx(1);
    gPad->SetGridy(1);
    if (!m_histref_phiindex)m_hist_crphi->SetTitle(Form("%s (no ref file)", m_hist_crphi->GetTitle()));
    m_hist_crphi->Draw("hist");
    if (isFew) colorizeCanvas(c_hist_crphi, c_StatusTooFew);
    else if (isAlarm)colorizeCanvas(c_hist_crphi, c_StatusError);
    else if (isWarn)colorizeCanvas(c_hist_crphi, c_StatusWarning);
    else colorizeCanvas(c_hist_crphi, c_StatusGood);
    c_hist_crphi->Update();
    UpdateCanvas(c_hist_crphi);
  }
 
  //get tracking efficiency
  auto m_delta_effphi = (TH2F*)getDelta(m_name_dir, m_hname_effphi, 0, true); //true=only if updated
  if (m_delta_effphi) {
    c_hist_effphi->Clear();
    double eff = -1;
    const int all_phibins = m_delta_effphi->GetNbinsX();
    const int all_hitbins = m_delta_effphi->GetNbinsY();
    const int thr_hitbin = m_delta_effphi->GetYaxis()->FindBin(20);//min hits bin
    for (int iphi = 0; iphi < all_phibins; iphi++) {
      TH1D* temp = (TH1D*)m_delta_effphi->ProjectionY(Form("hhits_bin_%d", iphi + 1), iphi + 1, iphi + 1, "");
      Double_t num = temp->Integral(thr_hitbin, all_hitbins);
      Double_t den = temp->Integral();
      if (den > 0)eff = num * 100.0 / den;
      m_hist_effphi->SetBinContent(iphi + 1, eff);
      m_hist_effphi->SetBinError(iphi + 1, 0);
      delete temp;
    }
    m_hist_effphi->GetYaxis()->SetRangeUser(80.0, 110.0); //per efficiency
    m_hist_effphi->SetTitle("CDC track efficiency(cdchits>20/all); cdc-track #phi; tracking efficiency");
    c_hist_effphi->cd();
    gPad->SetGridx();
    gPad->SetGridy();
    m_hist_effphi->SetFillColor(kCyan);
    m_hist_effphi->Draw("hist");
    c_hist_effphi->Update();
    UpdateCanvas(c_hist_effphi);
  }
 
  //get cdc hits vs phi
  auto m_delta_hitphi = (TH2F*)getDelta(m_name_dir, m_hname_hitsphi, 0, true); //true=only if updated
  if (m_delta_hitphi) {
    c_hist_hitsphi->Clear();
    m_delta_hitphi->SetTitle("CDC track #phi vs cdchits; cdc-track #phi; nCDCHits");
    c_hist_hitsphi->cd();
    m_delta_hitphi->Draw("COLZ");
    c_hist_hitsphi->Update();
    UpdateCanvas(c_hist_hitsphi);
  }
 
  // get wire efficiency
  double meanWireAttachProb = 0;
  double fracWiresWithLowAttachProb = 0;
  double fracWiresWithHighAttachProb = 0;
  gStyle->SetNumberContours(100);
  auto m_delta_efflay = (TH2F*)getDelta(m_name_dir, m_histoTrackingWireEff, 0, true); //true=only if updated
  if (m_delta_efflay) {
    for (int ij = 0; ij < 4; ij++) c_hist_attach_eff[ij]->Clear();
    m_hist_wire_attach_eff_1d->Reset();
    int nEffiValues = 0;
    for (int ij = 1; ij <= m_delta_efflay->GetNbinsX(); ij++) {
      int halfYbin = m_delta_efflay->GetNbinsY() / 2;
      for (int jk = 0; jk < halfYbin; jk++) {
        if (m_delta_efflay->GetBinContent(ij, jk + 1) == 0) continue;
        double binEffi = m_delta_efflay->GetBinContent(ij, jk + halfYbin + 1) / m_delta_efflay->GetBinContent(ij, jk + 1);
        m_hist_wire_attach_eff_1d->Fill(binEffi);
        meanWireAttachProb += binEffi;
        nEffiValues++;
      }
    }
    if (nEffiValues) meanWireAttachProb /= nEffiValues;
    if (m_doTH2PolyTrackingWireEff)
      fillEffiTH2Poly(m_delta_efflay, m_hist_attach_eff_Poly[0], m_hist_attach_eff_Poly[1], m_hist_attach_eff_Poly[2]);
    else
      fillEffiTH2(m_delta_efflay, m_hist_attach_eff[0], m_hist_attach_eff[1], m_hist_attach_eff[2]);
    TLatex latex;
    latex.SetTextSize(0.025);
    for (int ij = 0; ij < 3; ij++) {
      c_hist_attach_eff[ij]->cd();
      if (m_doTH2PolyTrackingWireEff) {
        m_hist_attach_eff_Poly[ij]->SetStats(0);
        m_hist_attach_eff_Poly[ij]->Draw("COLZ");
      } else {
        m_hist_attach_eff[ij]->SetStats(0);
        m_hist_attach_eff[ij]->Draw("COLZ");
        int isl = 0;
        TEllipse* el[9];
        for (int ilay = 0; ilay < 56; ilay++) {
          int rmdr = int(abs(ilay - 2) % 6);
          if ((rmdr == 0 && ilay > 2) || ilay == 55) {
            isl++;
            el[isl] = new TEllipse(0, 0, m_lbinEdges[ilay], m_lbinEdges[ilay]);
            el[isl]->SetLineColor(kRed);
            el[isl]->SetLineWidth(2);
            el[isl]->SetFillStyle(0);
            el[isl]->Draw("same");
          }
        }
      }
      if (ij == 2)
        latex.DrawLatexNDC(0.12, 0.87, TString::Format("mean = %.3f%%", meanWireAttachProb * 100.0));
    }
    c_hist_attach_eff[3]->cd();
    if (nEffiValues) {
      int firstBoundaryBin = m_hist_wire_attach_eff_1d->GetXaxis()->FindBin(m_firstEffBoundary) - 1;
      fracWiresWithLowAttachProb = m_hist_wire_attach_eff_1d->Integral(1, firstBoundaryBin) / nEffiValues;
      int secondBoundaryBin = m_hist_wire_attach_eff_1d->GetXaxis()->FindBin(m_secondEffBoundary) - 1;
      fracWiresWithHighAttachProb =  m_hist_wire_attach_eff_1d->Integral(secondBoundaryBin + 1,
                                     m_hist_wire_attach_eff_1d->GetNbinsX()) / nEffiValues;
      m_hist_wire_attach_eff_1d->SetStats(0);
      m_hist_wire_attach_eff_1d->Draw();
      latex.DrawLatexNDC(0.15, 0.87, TString::Format("%06.3f%% wire : eff < %.2f",
                                                     fracWiresWithLowAttachProb * 100,
                                                     m_firstEffBoundary));
      latex.DrawLatexNDC(0.15, 0.84, TString::Format("%06.3f%% wire : %.2f < eff < %.2f",
                                                     (1. - fracWiresWithHighAttachProb - fracWiresWithLowAttachProb) * 100,
                                                     m_firstEffBoundary, m_secondEffBoundary));
      latex.DrawLatexNDC(0.15, 0.81, TString::Format("%06.3f%% wire : %.2f < eff",
                                                     fracWiresWithHighAttachProb * 100,
                                                     m_secondEffBoundary));
    }
    for (int ij = 0; ij < 4; ij++) {
      c_hist_attach_eff[ij]->Update();
      UpdateCanvas(c_hist_attach_eff[ij]);
    }
  }
 
  m_monObj->setVariable("meanWireAttachProb", meanWireAttachProb);
  m_monObj->setVariable("fracWiresWithLowAttachProb", fracWiresWithLowAttachProb);
  m_monObj->setVariable("fracWiresWithHighAttachProb", fracWiresWithHighAttachProb);
 
  B2DEBUG(20, "DQMHistAnalysisCDCEpics: end event");
}
 
//------------------------------------
void DQMHistAnalysisCDCEpicsModule::endRun()
{
  B2DEBUG(20, "DQMHistAnalysisCDCEpics: end run");
}
 
 
//------------------------------------
void DQMHistAnalysisCDCEpicsModule::terminate()
{
  for (auto* line : m_lines) delete line;
  m_lines.clear();
  B2DEBUG(20, "DQMHistAnalysisCDCEpics: terminate called");
}
 
 
//------------------------------------
float DQMHistAnalysisCDCEpicsModule::getHistMedian(TH1D* h) const
{
  TH1D* hist = (TH1D*)h->Clone();
  hist->SetBinContent(1, 0.0); // Exclude 0-th bin
  float median = 0.0;
  if (hist->GetMean() != 0) {
    // Avoid an error if only TCD/ADC=0 entries
    double quantiles[1] = {0.0}; // One element to store median
    double probSums[1] = {0.5}; // Median definition
    hist->GetQuantiles(1, quantiles, probSums);
    median = quantiles[0];
  }
  delete hist;
  return median;
}
 
 
void DQMHistAnalysisCDCEpicsModule::fillEffiTH2(TH2F* hist, TH2F* attached, TH2F* expected, TH2F* efficiency)
{
  static CDC::CDCGeometryPar& cdcgeo = CDC::CDCGeometryPar::Instance();
  int nSLayers = cdcgeo.getNumberOfSenseLayers();
 
  // Array to hold bin edges, with nSLayers + 1 edges needed for nSLayers bins
  std::vector<double> binEdges(nSLayers + 1);
  // Calculate r bin edges
  double firstR = cdcgeo.senseWireR(0);
  double secondR = cdcgeo.senseWireR(1);
  binEdges[0] = firstR - (secondR - firstR) / 2;
  m_lbinEdges[0] = firstR;
  for (int lay = 1; lay < nSLayers; lay++) {
    double prevR = cdcgeo.senseWireR(lay - 1);
    double currentR = cdcgeo.senseWireR(lay);
    binEdges[lay] = (prevR + currentR) / 2;
    m_lbinEdges[lay] = currentR;
  }
  double lastR = cdcgeo.senseWireR(nSLayers - 1);
  double secondLastR = cdcgeo.senseWireR(nSLayers - 2);
  binEdges[nSLayers] = lastR + (lastR - secondLastR) / 2;
  // convenient histogram to get layer number
  TH1F layerHist("layerHist", "Layer Histogram", nSLayers, binEdges.data());
 
  for (int binx = 1; binx <= efficiency->GetNbinsX(); binx++) {
    for (int biny = 1; biny <= efficiency->GetNbinsY(); biny++) {
      double bincenterx = efficiency->GetXaxis()->GetBinCenter(binx);
      double bincentery = efficiency->GetYaxis()->GetBinCenter(biny);
      double r = TMath::Sqrt(bincenterx * bincenterx + bincentery * bincentery);
      double phi = TMath::ATan2(bincentery, bincenterx);
      if (phi < 0) phi += 2 * TMath::Pi();
 
      int layerBin = layerHist.FindBin(r); // Get the bin corresponding to r
      if (layerBin < 1 || layerBin > nSLayers) continue;
      int layerExpected = layerBin - 1;
 
      int nWires = cdcgeo.nWiresInLayer(layerExpected);
      double offset = cdcgeo.offset(layerExpected);
      int wireExpected = phi * nWires / (2 * TMath::Pi()) - offset + 0.5;
      if (wireExpected < 0) wireExpected += nWires;
      if (wireExpected >= nWires) wireExpected -= nWires;
 
      int expBin = hist->GetYaxis()->FindBin(layerExpected);
      int obsBin = expBin + nSLayers;
      expected->SetBinContent(binx, biny, hist->GetBinContent(wireExpected + 1, expBin));
      attached->SetBinContent(binx, biny, hist->GetBinContent(wireExpected + 1, obsBin));
    }
  }
  efficiency->Divide(attached, expected);
}
 
TH2Poly* DQMHistAnalysisCDCEpicsModule::createEffiTH2Poly(const TString& name, const TString& title)
{
  static CDC::CDCGeometryPar& cdcgeo = CDC::CDCGeometryPar::Instance();
  int nSLayers = cdcgeo.getNumberOfSenseLayers();
  double maxLayerR = cdcgeo.senseWireR(nSLayers - 1);
  TH2Poly* hist = new TH2Poly(name, title, -maxLayerR * 1.02, maxLayerR * 1.02, -maxLayerR * 1.02, maxLayerR * 1.02);
  for (int lay = 0; lay < nSLayers; lay++) {
    int nWires = cdcgeo.nWiresInLayer(lay);
    double offset = cdcgeo.offset(lay);
    double layerR = cdcgeo.senseWireR(lay);
    double r_inner = 0;
    double r_outer = 0;
    if (lay == 0) {
      r_inner = layerR - (cdcgeo.senseWireR(1) - cdcgeo.senseWireR(0)) / 2;
      r_outer = layerR + (cdcgeo.senseWireR(1) - cdcgeo.senseWireR(0)) / 2;
    } else if (lay == nSLayers - 1) {
      r_inner = layerR - (cdcgeo.senseWireR(lay) - cdcgeo.senseWireR(lay - 1)) / 2;
      r_outer = layerR + (cdcgeo.senseWireR(lay) - cdcgeo.senseWireR(lay - 1)) / 2;
    } else {
      r_inner = layerR - (cdcgeo.senseWireR(lay) - cdcgeo.senseWireR(lay - 1)) / 2;
      r_outer = layerR + (cdcgeo.senseWireR(lay + 1) - cdcgeo.senseWireR(lay)) / 2;
    }
    for (int wire = 0; wire < nWires; wire++) {
      double phi_inner = (wire - 0.5 + offset) * 2 * TMath::Pi() / nWires;
      double phi_outer = (wire + 0.5 + offset) * 2 * TMath::Pi() / nWires;
      // Calculate the four corners of the bin
      double x0 = r_inner * TMath::Cos(phi_inner);
      double y0 = r_inner * TMath::Sin(phi_inner);
      double x1 = r_outer * TMath::Cos(phi_inner);
      double y1 = r_outer * TMath::Sin(phi_inner);
      double x2 = r_outer * TMath::Cos(phi_outer);
      double y2 = r_outer * TMath::Sin(phi_outer);
      double x3 = r_inner * TMath::Cos(phi_outer);
      double y3 = r_inner * TMath::Sin(phi_outer);
      double xx[] = {x0, x1, x2, x3};
      double yy[] = {y0, y1, y2, y3};
      hist->AddBin(4, xx, yy);
    }
  }
  return hist;
}
 
void DQMHistAnalysisCDCEpicsModule::fillEffiTH2Poly(TH2F* hist, TH2Poly* attached, TH2Poly* expected, TH2Poly* efficiency)
{
  attached->Reset("ICES");
  expected->Reset("ICES");
  static CDC::CDCGeometryPar& cdcgeo = CDC::CDCGeometryPar::Instance();
  int nSLayers = cdcgeo.getNumberOfSenseLayers();
  for (int lay = 0; lay < nSLayers; lay++) {
    int nWires = cdcgeo.nWiresInLayer(lay);
    double layerR = cdcgeo.senseWireR(lay);
    double offset = cdcgeo.offset(lay);
    int expBin = hist->GetYaxis()->FindBin(lay);
    int obsBin = expBin + nSLayers;
    // fill the bins for this layer
    for (int wire = 0; wire < nWires; wire++) {
      double phi = (wire + offset) * 2 * TMath::Pi() / nWires;
      double fillX = layerR * TMath::Cos(phi);
      double fillY = layerR * TMath::Sin(phi);
      attached->Fill(fillX, fillY, hist->GetBinContent(wire + 1, obsBin));
      expected->Fill(fillX, fillY, hist->GetBinContent(wire + 1, expBin));
    }
  }
  efficiency->Divide(attached, expected);
}