KLMDQMModule.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.                  *
 **************************************************************************/
 
/* Own header. */
#include <klm/modules/KLMDQM/KLMDQMModule.h>
 
/* KLM headers. */
#include <klm/dataobjects/KLMChannelIndex.h>
#include <klm/dataobjects/KLMDigitRaw.h>
 
/* ROOT headers. */
#include <TDirectory.h>
 
#include <set>
 
using namespace Belle2;
 
REG_MODULE(KLMDQM);
 
KLMDQMModule::KLMDQMModule() :
  HistoModule(),
  m_DAQInclusion{nullptr},
  m_TimeRPC{nullptr},
  m_TimeScintillatorBKLM{nullptr},
  m_TimeScintillatorEKLM{nullptr},
  m_ChargeScintillatorBKLM{nullptr},
  m_ChargeScintillatorEKLM{nullptr},
  m_ChargeScintillatorBKLM_SingleStrip{nullptr},
  m_ChargeScintillatorEKLM_SingleStrip{nullptr},
  m_ChargeScintillatorBKLM_MultiStrip{nullptr},
  m_ChargeScintillatorEKLM_MultiStrip{nullptr},
  m_ChargeClusterBKLM{nullptr},
  m_ChargeClusterEKLM{nullptr},
  m_AverageChargeClusterBKLM{nullptr},
  m_AverageChargeClusterEKLM{nullptr},
  m_PlaneBKLMPhi{nullptr},
  m_PlaneBKLMZ{nullptr},
  m_PlaneEKLM{nullptr},
  m_DigitsKLM{nullptr},
  m_DigitsRPC{nullptr},
  m_DigitsScintillatorBKLM{nullptr},
  m_DigitsScintillatorEKLM{nullptr},
  m_DigitsMultiStripBKLM{nullptr},
  m_DigitsMultiStripEKLM{nullptr},
  m_TriggerBitsBKLM{nullptr},
  m_TriggerBitsEKLM{nullptr},
  m_EventBackgroundTriggerSummary{nullptr},
  m_DigitsAfterLERInj{nullptr},
  m_TriggersLERInj{nullptr},
  m_DigitsAfterHERInj{nullptr},
  m_TriggersHERInj{nullptr},
  m_FE_BKLM_Layer_0(nullptr),
  m_FE_BKLM_Layer_1(nullptr),
  m_FE_EKLM_Plane_0(nullptr),
  m_FE_EKLM_Plane_1(nullptr),
  m_ChannelArrayIndex{&(KLMChannelArrayIndex::Instance())},
  m_SectorArrayIndex{&(KLMSectorArrayIndex::Instance())},
  m_ElementNumbers{&(KLMElementNumbers::Instance())},
  m_eklmElementNumbers{&(EKLMElementNumbers::Instance())},
  m_klmTime{&(KLMTime::Instance())}
{
  setDescription("KLM data quality monitor.");
  setPropertyFlags(c_ParallelProcessingCertified);
  addParam("histogramDirectoryName", m_HistogramDirectoryName,
           "Directory for KLM DQM histograms in ROOT file.",
           std::string("KLM"));
  addParam("RPCTimeMin", m_RPCTimeMin,
           "Min time for RPC time histogram.", double(-1223.5));
  addParam("RPCTimeMax", m_RPCTimeMax,
           "Max time for RPC time histogram.", double(-199.5));
  addParam("BKLMScintTimeMin", m_BKLMScintTimeMin,
           "Min time for BKLM scintillator time histogram.", double(-5300));
  addParam("BKLMScintTimeMax", m_BKLMScintTimeMax,
           "Max time for BKLM scintillator time histogram.", double(-4300));
  addParam("EKLMScintTimeMin", m_EKLMScintTimeMin,
           "Min time for EKLM scintillator time histogram.", double(-5300));
  addParam("EKLMScintTimeMax", m_EKLMScintTimeMax,
           "Max time for EKLM scintillator time histogram.", double(-4300));
  addParam("Revo9DCArrivalTimeMin", m_Revo9DCArrivalTimeMin,
           "Min time for RPC hit revo9DCArrivalTime histogram.", double(-10000));
  addParam("Revo9DCArrivalTimeMax", m_Revo9DCArrivalTimeMax,
           "Max time for RPC hit revo9DCArrivalTime histogram.", double(3000));
}
 
KLMDQMModule::~KLMDQMModule()
{
  KLMChannelIndex klmIndex(KLMChannelIndex::c_IndexLevelSector);
  for (KLMChannelIndex& klmSector : klmIndex) {
    KLMSectorNumber sector = klmSector.getKLMSectorNumber();
    KLMSectorNumber sectorIndex = m_SectorArrayIndex->getIndex(sector);
    if (m_ChannelHits[sectorIndex] != nullptr)
      delete[] m_ChannelHits[sectorIndex];
  }
  klmIndex.setIndexLevel(KLMChannelIndex::c_IndexLevelSection);
  for (KLMChannelIndex& klmSection : klmIndex) {
    KLMSubdetectorNumber subdetector = klmSection.getSubdetector();
    if (subdetector == KLMElementNumbers::c_EKLM) {
      KLMSectionNumber section = klmSection.getSection();
      if (m_Spatial2DHitsEKLM[section - 1] != nullptr)
        delete[] m_Spatial2DHitsEKLM[section - 1];
    }
  }
}
 
void KLMDQMModule::defineHisto()
{
  TDirectory* oldDirectory, *newDirectory;
  oldDirectory = gDirectory;
  newDirectory = oldDirectory->mkdir(m_HistogramDirectoryName.c_str());
  newDirectory->cd();
  /* DAQ inclusion. */
  m_DAQInclusion = new TH1F("daq_inclusion", "Is KLM included in DAQ?", 2, 0.0, 2.0);
  m_DAQInclusion->GetXaxis()->SetBinLabel(1, "No");
  m_DAQInclusion->GetXaxis()->SetBinLabel(2, "Yes");
  /* Time histograms. */
  m_TimeRPC = new TH1F("time_rpc", "RPC hit time", 128, m_RPCTimeMin, m_RPCTimeMax);
  m_TimeRPC->GetXaxis()->SetTitle("Time, ns");
  m_TimeScintillatorBKLM =
    new TH1F("time_scintillator_bklm", "Scintillator hit time (BKLM)",
             100, m_BKLMScintTimeMin, m_BKLMScintTimeMax);
  m_TimeScintillatorBKLM->GetXaxis()->SetTitle("Time, ns");
  m_TimeScintillatorEKLM =
    new TH1F("time_scintillator_eklm", "Scintillator hit time (EKLM)",
             100, m_EKLMScintTimeMin, m_EKLMScintTimeMax);
  m_TimeScintillatorEKLM->GetXaxis()->SetTitle("Time, ns");
  m_ChargeScintillatorBKLM =
    new TH1F("charge_scintillator_bklm", "Scintillator charge (BKLM)",
             256, 0., 512.);
  m_ChargeScintillatorBKLM->GetXaxis()->SetTitle("Charge");
  m_ChargeScintillatorEKLM =
    new TH1F("charge_scintillator_eklm", "Scintillator charge (EKLM)",
             256, 0., 512.);
  m_ChargeScintillatorEKLM->GetXaxis()->SetTitle("Charge");
  m_ChargeScintillatorBKLM_SingleStrip =
    new TH1F("charge_scintillator_bklm_single_strip", "Scintillator charge (BKLM) single strip",
             256, 0., 512.);
  m_ChargeScintillatorBKLM_SingleStrip->GetXaxis()->SetTitle("Charge");
  m_ChargeScintillatorEKLM_SingleStrip =
    new TH1F("charge_scintillator_eklm_single_strip", "Scintillator charge (EKLM) single strip",
             256, 0., 512.);
  m_ChargeScintillatorEKLM_SingleStrip->GetXaxis()->SetTitle("Charge");
  m_ChargeScintillatorBKLM_MultiStrip =
    new TH1F("charge_scintillator_bklm_multi_strip", "Scintillator charge (BKLM) multi strip",
             256, 0., 512.);
  m_ChargeScintillatorBKLM_MultiStrip->GetXaxis()->SetTitle("Charge");
  m_ChargeScintillatorEKLM_MultiStrip =
    new TH1F("charge_scintillator_eklm_multi_strip", "Scintillator charge (EKLM) multi strip",
             256, 0., 512.);
  m_ChargeScintillatorEKLM_MultiStrip->GetXaxis()->SetTitle("Charge");
  m_ChargeClusterBKLM =
    new TH1F("charge_cluster_bklm", "BKLM Scintillator Clusters",
             500, 0., 10000.);
  m_ChargeClusterBKLM->GetXaxis()->SetTitle("Cluster Charge");
  m_ChargeClusterEKLM =
    new TH1F("charge_cluster_eklm", "EKLM Clusters",
             500, 0., 10000.);
  m_ChargeClusterEKLM->GetXaxis()->SetTitle("Cluster Charge");
  m_AverageChargeClusterBKLM =
    new TH1F("avg_charge_cluster_bklm", "Mean Cluster Charge (BKLM)",
             256, 0., 512.);
  m_AverageChargeClusterBKLM->GetXaxis()->SetTitle("Cluster Charge / Number of Digits in Cluster");
  m_AverageChargeClusterEKLM =
    new TH1F("avg_charge_cluster_eklm", "Mean Cluster Charge (EKLM)",
             256, 0., 512.);
  m_AverageChargeClusterEKLM->GetXaxis()->SetTitle("Cluster Charge / Number of Digits in Cluster");
  m_TimeRevo9DCArrivalTime = new TH1F("time_revo9dc_arrival_time", "DC arrival hit time (RPC)",
                                      500, m_Revo9DCArrivalTimeMin, m_Revo9DCArrivalTimeMax);
  m_TimeRevo9DCArrivalTime->GetXaxis()->SetTitle("Time, ns");
  /* Number of hits per plane. */
  m_PlaneBKLMPhi = new TH1F("plane_bklm_phi",
                            "BKLM plane occupancy (#phi readout)",
                            240, 0.5, 240.5);
  m_PlaneBKLMPhi->GetXaxis()->SetTitle("Layer number");
  m_PlaneBKLMZ = new TH1F("plane_bklm_z",
                          "BKLM plane occupancy (z readout)",
                          240, 0.5, 240.5);
  m_PlaneBKLMZ->GetXaxis()->SetTitle("Layer number");
  m_PlaneEKLM = new TH1F("plane_eklm", "EKLM plane occupancy (both readouts)", 208, 0.5, 208.5);
  m_PlaneEKLM->GetXaxis()->SetTitle("Plane number");
  /* Number of hits per channel. */
  int nChannelHistograms =
    BKLMElementNumbers::getMaximalSectorGlobalNumber() *
    m_ChannelHitHistogramsBKLM +
    EKLMElementNumbers::getMaximalSectorGlobalNumberKLMOrder() *
    m_ChannelHitHistogramsEKLM;
  KLMChannelNumber* firstChannelNumbers =
    new KLMChannelNumber[nChannelHistograms + 1];
  int i = 0;
  KLMChannelIndex klmIndex(KLMChannelIndex::c_IndexLevelSector);
  for (KLMChannelIndex& klmSector : klmIndex) {
    KLMChannelIndex klmChannel(klmSector);
    klmChannel.setIndexLevel(KLMChannelIndex::c_IndexLevelStrip);
    KLMChannelNumber channel = klmChannel.getKLMChannelNumber();
    firstChannelNumbers[i] = m_ChannelArrayIndex->getIndex(channel);
    if (klmSector.getSubdetector() == KLMElementNumbers::c_BKLM) {
      channel = m_ElementNumbers->channelNumberBKLM(
                  klmChannel.getSection(), klmChannel.getSector(), 8, 0, 1);
      firstChannelNumbers[i + 1] = m_ChannelArrayIndex->getIndex(channel);
      i += 2;
    } else {
      int layerIncrease = (klmSector.getSection() == 1) ? 4 : 5;
      channel = m_ElementNumbers->channelNumberEKLM(
                  klmChannel.getSection(), klmChannel.getSector(),
                  1 + layerIncrease, 1, 1);
      firstChannelNumbers[i + 1] = m_ChannelArrayIndex->getIndex(channel);
      channel = m_ElementNumbers->channelNumberEKLM(
                  klmChannel.getSection(), klmChannel.getSector(),
                  1 + layerIncrease * 2, 1, 1);
      firstChannelNumbers[i + 2] = m_ChannelArrayIndex->getIndex(channel);
      i += 3;
    }
  }
  firstChannelNumbers[nChannelHistograms] = m_ChannelArrayIndex->getNElements();
  i = 0;
  klmIndex.setIndexLevel(KLMChannelIndex::c_IndexLevelSector);
  for (KLMChannelIndex& klmSector : klmIndex) {
    int nHistograms;
    if (klmSector.getSubdetector() == KLMElementNumbers::c_BKLM)
      nHistograms = m_ChannelHitHistogramsBKLM;
    else
      nHistograms = m_ChannelHitHistogramsEKLM;
    KLMSectorNumber sector = klmSector.getKLMSectorNumber();
    KLMSectorNumber sectorIndex = m_SectorArrayIndex->getIndex(sector);
    m_ChannelHits[sectorIndex] = new TH1F*[nHistograms];
    for (int j = 0; j < nHistograms; j++) {
      std::string name =
        "strip_hits_subdetector_" + std::to_string(klmSector.getSubdetector()) +
        "_section_" + std::to_string(klmSector.getSection()) +
        "_sector_" + std::to_string(klmSector.getSector()) +
        "_" + std::to_string(j);
      std::string title = "Sector " + std::to_string(klmSector.getSector()) + " -- " +
                          m_ElementNumbers->getSectorDAQName(klmSector.getSubdetector(), klmSector.getSection(), klmSector.getSector());
      m_ChannelHits[sectorIndex][j] = new TH1F(
        name.c_str(), title.c_str(),
        firstChannelNumbers[i + 1] - firstChannelNumbers[i],
        firstChannelNumbers[i] - 0.5, firstChannelNumbers[i + 1] - 0.5);
      i++;
    }
  }
  delete[] firstChannelNumbers;
  /* Number of digits. */
  m_DigitsKLM = new TH1F("digits_klm", "Number of KLM digits",
                         250.0, 0.0, 250.0);
  m_DigitsKLM->GetXaxis()->SetTitle("Number of digits");
  m_DigitsRPC = new TH1F("digits_rpc", "Number of RPC digits",
                         250.0, 0.0, 250.0);
  m_DigitsRPC->GetXaxis()->SetTitle("Number of digits");
  m_DigitsScintillatorBKLM = new TH1F("digits_scintillator_bklm", "Number of scintillator digits (BKLM)",
                                      250.0, 0.0, 250.0);
  m_DigitsScintillatorBKLM->GetXaxis()->SetTitle("Number of digits");
  m_DigitsScintillatorEKLM = new TH1F("digits_scintillator_eklm", "Number of scintillator digits (EKLM)",
                                      250.0, 0.0, 250.0);
  m_DigitsScintillatorEKLM->GetXaxis()->SetTitle("Number of digits");
  m_DigitsMultiStripBKLM = new TH1F("digits_multi_bklm", "Number of multi-strip digits (BKLM)",
                                    50.0, 0.0, 50.0);
  m_DigitsMultiStripBKLM->GetXaxis()->SetTitle("Number of multi-strip digits");
  m_DigitsMultiStripEKLM = new TH1F("digits_multi_eklm", "Number of multi-strip digits (EKLM)",
                                    50.0, 0.0, 50.0);
  m_DigitsMultiStripEKLM->GetXaxis()->SetTitle("Number of multi-strip digits");
  /* Trigger bits. */
  m_TriggerBitsBKLM = new TH1F("trigger_bits_bklm", "Trigger bits of multi-strip digits (BKLM)",
                               (double)c_0x1, (double)c_0x8, (double)c_0x1 + 1.0);
  m_TriggerBitsBKLM->GetXaxis()->SetBinLabel(c_0x8, "0x8");
  m_TriggerBitsBKLM->GetXaxis()->SetBinLabel(c_0x4, "0x4");
  m_TriggerBitsBKLM->GetXaxis()->SetBinLabel(c_0x2, "0x2");
  m_TriggerBitsBKLM->GetXaxis()->SetBinLabel(c_0x1, "0x1");
  m_TriggerBitsEKLM = new TH1F("trigger_bits_eklm", "Trigger bits of multi-strip digits (EKLM)",
                               (double)c_0x1, (double)c_0x8, (double)c_0x1 + 1.0);
  m_TriggerBitsEKLM->GetXaxis()->SetBinLabel(c_0x8, "0x8");
  m_TriggerBitsEKLM->GetXaxis()->SetBinLabel(c_0x4, "0x4");
  m_TriggerBitsEKLM->GetXaxis()->SetBinLabel(c_0x2, "0x2");
  m_TriggerBitsEKLM->GetXaxis()->SetBinLabel(c_0x1, "0x1");
 
  /* Event-level L1 trigger bits (TRGSummary). */
  int nTimingBits = (int)c_KlmL1Triggers.size();          // TTYP_DPHY, TTYP_RAND, TTYP_POIS
  m_EventBackgroundTriggerSummary = new TH1F("event_background_trigger_summary",
                                             "Event background trigger summary;Trigger Decision;Events",
                                             nTimingBits, 0.5, 0.5 + nTimingBits);
  /* Keep labels explicit to match the actual trigger bits. */
  m_EventBackgroundTriggerSummary->GetXaxis()->SetBinLabel(1, "TTYP_DPHY");
  m_EventBackgroundTriggerSummary->GetXaxis()->SetBinLabel(2, "TTYP_RAND");
  m_EventBackgroundTriggerSummary->GetXaxis()->SetBinLabel(3, "TTYP_POIS");
 
  /* Number of digits after injection */
  /* For the histograms below, we use the same style as for other subdetectors. */
  m_DigitsAfterLERInj = new TH1F("KLMOccInjLER", "KLM digits after LER injection / Time;Time [#mus];Number of KLM digits / (5 #mus)",
                                 4000, 0, 20000);
  m_TriggersLERInj = new TH1F("KLMTrigInjLER", "Triggers after KER injection / Time;Time [#mus];Number of triggers / (5 #mus)",
                              4000, 0, 20000);
  m_DigitsAfterHERInj = new TH1F("KLMOccInjHER", "KLM digits after HER injection / Time;Time [#mus];Number of KLM digits / (5 #mus)",
                                 4000, 0, 20000);
  m_TriggersHERInj = new TH1F("KLMTrigInjHER", "Triggers after HER injection / Time;Time [#mus];Number of triggers / (5 #mus)",
                              4000, 0, 20000);
  /* Spatial distribution of EKLM 2d hits per layer. */
  klmIndex.setIndexLevel(KLMChannelIndex::c_IndexLevelSection);
  for (KLMChannelIndex& klmSection : klmIndex) {
    KLMSubdetectorNumber subdetector = klmSection.getSubdetector();
    if (subdetector == KLMElementNumbers::c_EKLM) {
      KLMSectionNumber section = klmSection.getSection();
      int maximalLayerNumber = m_eklmElementNumbers->getMaximalDetectorLayerNumber(section);
      m_Spatial2DHitsEKLM[section - 1] = new TH2F*[maximalLayerNumber];
      std::string sectionName = (section == EKLMElementNumbers::c_ForwardSection) ? "Forward" : "Backward";
      for (int j = 1; j <= maximalLayerNumber; ++j) {
        std::string name = "spatial_2d_hits_subdetector_" + std::to_string(subdetector) +
                           "_section_" + std::to_string(section) +
                           "_layer_" + std::to_string(j);
        std::string title = "Endcap " + sectionName + " , Layer " + std::to_string(j);
        /* Use bins with a size of 10 cm per side. */
        m_Spatial2DHitsEKLM[section - 1][j - 1] = new TH2F(name.c_str(), title.c_str(),
                                                           340 * 2 / 10, -340, 340,
                                                           340 * 2 / 10, -340, 340);
        m_Spatial2DHitsEKLM[section - 1][j - 1]->GetXaxis()->SetTitle("X coordinate [cm]");
        m_Spatial2DHitsEKLM[section - 1][j - 1]->GetYaxis()->SetTitle("Y coordinate [cm]");
      }
    }
  }
  // Feature extraction status histogram for BKLM
  int bklmSectors = BKLMElementNumbers::getMaximalSectorGlobalNumber(); // 16
  int eklmPlanes = EKLMElementNumbers::getMaximalPlaneGlobalNumber(); // 208
 
  m_FE_BKLM_Layer_0 = new TH1F("feStatus_bklm_scintillator_layers_0",
                               "BKLM Scintillator Standard Readout;FEE Card", bklmSectors * 2, 0.5, 0.5 + bklmSectors * 2);
  m_FE_BKLM_Layer_1 = new TH1F("feStatus_bklm_scintillator_layers_1",
                               "BKLM Scintillator Feature Extraction;FEE Card", bklmSectors * 2, 0.5, 0.5 + bklmSectors * 2);
  m_FE_EKLM_Plane_0 = new TH1F("feStatus_eklm_plane_0",
                               "EKLM Standard Readout;Plane number", eklmPlanes, 0.5, 0.5 + eklmPlanes);
  m_FE_EKLM_Plane_1 = new TH1F("feStatus_eklm_plane_1",
                               "EKLM Feature Extraction;Plane number", eklmPlanes, 0.5, 0.5 + eklmPlanes);
  oldDirectory->cd();
}
 
void KLMDQMModule::initialize()
{
  REG_HISTOGRAM;
  m_RawFtsws.isOptional();
  m_RawKlms.isOptional();
  m_trgSummary.isOptional();
  m_Digits.isOptional();
  m_BklmHit1ds.isOptional();
  m_Hit2ds.isOptional();
  m_KLMClusters.isOptional();
}
 
void KLMDQMModule::beginRun()
{
  /* DAQ inclusion. */
  m_DAQInclusion->Reset();
  if (m_RawKlms.isValid())
    m_DAQInclusion->Fill("Yes", 1);
  else
    m_DAQInclusion->Fill("No", 1);
  /* Time. */
  m_TimeRPC->Reset();
  m_TimeScintillatorBKLM->Reset();
  m_TimeScintillatorEKLM->Reset();
  m_ChargeScintillatorBKLM->Reset();
  m_ChargeScintillatorEKLM->Reset();
  m_ChargeScintillatorBKLM_SingleStrip->Reset();
  m_ChargeScintillatorEKLM_SingleStrip->Reset();
  m_ChargeScintillatorBKLM_MultiStrip->Reset();
  m_ChargeScintillatorEKLM_MultiStrip->Reset();
  m_ChargeClusterBKLM->Reset();
  m_ChargeClusterEKLM->Reset();
  m_AverageChargeClusterBKLM->Reset();
  m_AverageChargeClusterEKLM->Reset();
  m_TimeRevo9DCArrivalTime->Reset();
  m_klmTime->updateConstants(); //to get correct CTime
  /* Plane hits. */
  m_PlaneEKLM->Reset();
  m_PlaneBKLMPhi->Reset();
  m_PlaneBKLMZ->Reset();
  /* Channel hits. */
  KLMChannelIndex klmIndex(KLMChannelIndex::c_IndexLevelSector);
  for (KLMChannelIndex& klmSector : klmIndex) {
    int nHistograms;
    if (klmSector.getSubdetector() == KLMElementNumbers::c_BKLM)
      nHistograms = m_ChannelHitHistogramsBKLM;
    else
      nHistograms = m_ChannelHitHistogramsEKLM;
    KLMSectorNumber sector = klmSector.getKLMSectorNumber();
    KLMSectorNumber sectorIndex = m_SectorArrayIndex->getIndex(sector);
    for (int j = 0; j < nHistograms; j++)
      m_ChannelHits[sectorIndex][j]->Reset();
  }
  /* Digits. */
  m_DigitsKLM->Reset();
  m_DigitsRPC->Reset();
  m_DigitsScintillatorBKLM->Reset();
  m_DigitsScintillatorEKLM->Reset();
  m_DigitsMultiStripBKLM->Reset();
  m_DigitsMultiStripEKLM->Reset();
  /* Trigger bits. */
  m_TriggerBitsBKLM->Reset();
  m_TriggerBitsEKLM->Reset();
  /* Event-level background trigger summary. */
  if (m_EventBackgroundTriggerSummary != nullptr)
    m_EventBackgroundTriggerSummary->Reset();
  /* Injection information. */
  m_DigitsAfterLERInj->Reset();
  m_TriggersLERInj->Reset();
  m_DigitsAfterHERInj->Reset();
  m_TriggersHERInj->Reset();
  /* Spatial 2D hits distributions. */
  klmIndex.setIndexLevel(KLMChannelIndex::c_IndexLevelSection);
  for (KLMChannelIndex& klmSection : klmIndex) {
    KLMSubdetectorNumber subdetector = klmSection.getSubdetector();
    if (subdetector == KLMElementNumbers::c_EKLM) {
      KLMSectionNumber section = klmSection.getSection();
      int maximalLayerNumber =
        m_eklmElementNumbers->getMaximalDetectorLayerNumber(section);
      for (int j = 1; j <= maximalLayerNumber; ++j)
        m_Spatial2DHitsEKLM[section - 1][j - 1]->Reset();
    }
  }
  /* Feature extraction. */
  m_FE_BKLM_Layer_0->Reset();
  m_FE_BKLM_Layer_1->Reset();
  m_FE_EKLM_Plane_0->Reset();
  m_FE_EKLM_Plane_1->Reset();
}
 
void KLMDQMModule::event()
{
  /* Event-level background trigger summary, filled once per event. */
  if (m_trgSummary.isValid() && m_EventBackgroundTriggerSummary != nullptr) {
    const int nTimingBits = (int)c_KlmL1Triggers.size();
    // GDL Background triggers: TTYP_DPHY, TTYP_RAND, TTYP_POIS (event-level timing source)
    for (int i = 0; i < nTimingBits; ++i) {
      if (m_trgSummary->getTimType() == c_KlmL1Triggers[i])
        m_EventBackgroundTriggerSummary->Fill((double)i + 1.0);
    }
  }
 
  int nDigits = m_Digits.getEntries();
  int nDigitsRPC = 0, nDigitsScintillatorBKLM = 0, nDigitsScintillatorEKLM = 0;
  int nDigitsMultiStripBKLM = 0, nDigitsMultiStripEKLM = 0;
  for (const KLMDigit& digit : m_Digits) {
    /*
     * Reject digits that are below the threshold (such digits may appear
     * for simulated events).
     */
 
    if (!digit.isGood())
      continue;
    KLMDigitRaw* digitRaw = digit.getRelated<KLMDigitRaw>();
    if (digit.getSubdetector() == KLMElementNumbers::c_EKLM) {
      nDigitsScintillatorEKLM++;
      int section = digit.getSection();
      int sector = digit.getSector();
      int layer = digit.getLayer();
      int plane = digit.getPlane();
      int strip = digit.getStrip();
      if (not digit.isMultiStrip()) {
        KLMSectorNumber klmSector = m_ElementNumbers->sectorNumberEKLM(section, sector);
        KLMSectorNumber klmSectorIndex = m_SectorArrayIndex->getIndex(klmSector);
        KLMChannelNumber channel = m_ElementNumbers->channelNumberEKLM(section, sector, layer, plane, strip);
        KLMChannelNumber channelIndex = m_ChannelArrayIndex->getIndex(channel);
        for (int j = 0; j < m_ChannelHitHistogramsEKLM; j++) {
          double xMin = m_ChannelHits[klmSectorIndex][j]->GetXaxis()->GetXmin();
          double xMax = m_ChannelHits[klmSectorIndex][j]->GetXaxis()->GetXmax();
          if ((xMin > channelIndex) || (xMax < channelIndex))
            continue;
          m_ChannelHits[klmSectorIndex][j]->Fill(channelIndex);
        }
      } else
        nDigitsMultiStripEKLM++;
      int planeGlobal = m_eklmElementNumbers->planeNumber(section, layer, sector, plane);
      m_PlaneEKLM->Fill(planeGlobal);
      m_TimeScintillatorEKLM->Fill(digit.getTime());
      if (digit.isMultiStrip()) {
        if (digitRaw) {
          uint16_t triggerBits = digitRaw->getTriggerBits();
          if ((triggerBits & 0x1) != 0)
            m_TriggerBitsEKLM->Fill(c_0x1);
          if ((triggerBits & 0x2) != 0)
            m_TriggerBitsEKLM->Fill(c_0x2);
          if ((triggerBits & 0x4) != 0)
            m_TriggerBitsEKLM->Fill(c_0x4);
          if ((triggerBits & 0x8) != 0)
            m_TriggerBitsEKLM->Fill(c_0x8);
        }
      }
      if (digitRaw) {
        // Extract m_FE from m_word4
        uint16_t feStatus = digitRaw->getFEStatus(); // Extract the most significant bit
        if (feStatus != 0) {
          m_FE_EKLM_Plane_1->Fill(planeGlobal);
          m_ChargeScintillatorEKLM->Fill(digit.getCharge());
          uint16_t triggerBits = digitRaw->getTriggerBits();
          if ((triggerBits & 0x10) != 0) {
            m_ChargeScintillatorEKLM_MultiStrip->Fill(digit.getCharge());
          } else if ((triggerBits & 0x10) == 0) {
            m_ChargeScintillatorEKLM_SingleStrip->Fill(digit.getCharge());
          }
        } else {
          m_FE_EKLM_Plane_0->Fill(planeGlobal);
        }
      }
    } else if (digit.getSubdetector() == KLMElementNumbers::c_BKLM) {
      int section = digit.getSection();
      int sector = digit.getSector();
      int layer = digit.getLayer();
      int plane = digit.getPlane();
      int strip = digit.getStrip();
 
      KLMSectorNumber klmSector = m_ElementNumbers->sectorNumberBKLM(section, sector);
      KLMSectorNumber klmSectorIndex = m_SectorArrayIndex->getIndex(klmSector);
 
      if (not digit.isMultiStrip()) {
        KLMChannelNumber channel = m_ElementNumbers->channelNumberBKLM(section, sector, layer, plane, strip);
        KLMChannelNumber channelIndex = m_ChannelArrayIndex->getIndex(channel);
        for (int j = 0; j < m_ChannelHitHistogramsBKLM; j++) {
          double xMin = m_ChannelHits[klmSectorIndex][j]->GetXaxis()->GetXmin();
          double xMax = m_ChannelHits[klmSectorIndex][j]->GetXaxis()->GetXmax();
          if ((xMin > channelIndex) || (xMax < channelIndex))
            continue;
          m_ChannelHits[klmSectorIndex][j]->Fill(channelIndex);
        }
      } else
        nDigitsMultiStripBKLM++;
      if (digit.inRPC()) {
        nDigitsRPC++;
        m_TimeRPC->Fill(digit.getTime());
        m_TimeRevo9DCArrivalTime->Fill(digit.getRevo9DCArrivalTime() * m_klmTime->getCTimePeriod());
      } else {
        nDigitsScintillatorBKLM++;
        m_TimeScintillatorBKLM->Fill(digit.getTime());
        if (digitRaw) {
          uint16_t feStatus = digitRaw->getFEStatus(); // Extract the most significant bit
          if (feStatus != 0) {
            m_FE_BKLM_Layer_1->Fill((klmSectorIndex) * 2 + layer);
            m_ChargeScintillatorBKLM->Fill(digit.getCharge());
            uint16_t triggerBits = digitRaw->getTriggerBits();
            if ((triggerBits & 0x10) != 0) {
              m_ChargeScintillatorBKLM_MultiStrip->Fill(digit.getCharge());
            } else if ((triggerBits & 0x10) == 0) {
              m_ChargeScintillatorBKLM_SingleStrip->Fill(digit.getCharge());
            }
          } else {
            m_FE_BKLM_Layer_0->Fill((klmSectorIndex) * 2 + layer);
          }
        }
      }
      if (digit.isMultiStrip()) {
        if (digitRaw) {
          uint16_t triggerBits = digitRaw->getTriggerBits();
          if ((triggerBits & 0x1) != 0)
            m_TriggerBitsBKLM->Fill(c_0x1);
          if ((triggerBits & 0x2) != 0)
            m_TriggerBitsBKLM->Fill(c_0x2);
          if ((triggerBits & 0x4) != 0)
            m_TriggerBitsBKLM->Fill(c_0x4);
          if ((triggerBits & 0x8) != 0)
            m_TriggerBitsBKLM->Fill(c_0x8);
        }
      }
    } else
      B2FATAL("Not a BKLM or a EKLM digit, something went really wrong.");
  }
  for (const BKLMHit1d& hit1d : m_BklmHit1ds) {
    int section = hit1d.getSection();
    int sector = hit1d.getSector();
    int layer = hit1d.getLayer();
    int layerGlobal = BKLMElementNumbers::layerGlobalNumber(
                        section, sector, layer);
    if (hit1d.isPhiReadout())
      m_PlaneBKLMPhi->Fill(layerGlobal);
    else
      m_PlaneBKLMZ->Fill(layerGlobal);
  }
  /* Digits. */
  m_DigitsKLM->Fill((double)nDigits);
  m_DigitsRPC->Fill((double)nDigitsRPC);
  m_DigitsScintillatorBKLM->Fill((double)nDigitsScintillatorBKLM);
  m_DigitsScintillatorEKLM->Fill((double)nDigitsScintillatorEKLM);
  if (nDigitsMultiStripBKLM > 0)
    m_DigitsMultiStripBKLM->Fill((double)nDigitsMultiStripBKLM);
  if (nDigitsMultiStripEKLM > 0)
    m_DigitsMultiStripEKLM->Fill((double)nDigitsMultiStripEKLM);
  /* Injection information. */
  for (RawFTSW& rawFtsw : m_RawFtsws) {
    unsigned int difference = rawFtsw.GetTimeSinceLastInjection(0);
    if (difference != 0x7FFFFFFF) {
      /* 127 MHz clock ticks to us, inexact rounding. */
      float differenceInUs = difference / 127.;
      if (rawFtsw.GetIsHER(0)) {
        m_DigitsAfterHERInj->Fill(differenceInUs, nDigits);
        m_TriggersHERInj->Fill(differenceInUs);
      } else {
        m_DigitsAfterLERInj->Fill(differenceInUs, nDigits);
        m_TriggersLERInj->Fill(differenceInUs);
      }
    }
    /*
     * Usually, only one RawFTSW object is stored per event.
     * If there are more, ignore the others.
     */
    break;
  }
  /* Spatial 2D hits distributions. */
  for (const KLMHit2d& hit2d : m_Hit2ds) {
    if (hit2d.getSubdetector() != KLMElementNumbers::c_EKLM)
      continue;
    int section = hit2d.getSection();
    int layer = hit2d.getLayer();
    m_Spatial2DHitsEKLM[section - 1][layer - 1]->Fill(hit2d.getPositionX(), hit2d.getPositionY());
  }
  /* Cluster charge from related scintillator digits (unique per cluster). */
  if (m_KLMClusters.isValid()) {
    for (KLMCluster& cluster : m_KLMClusters) {
      std::set<const KLMDigit*> digitsSeen;
      double sumChargeBKLM = 0.;
      double sumChargeEKLM = 0.;
      int nDigitsClusterBKLM = 0;
      int nDigitsClusterEKLM = 0;
      auto addDigits = [&](const RelationVector<KLMDigit>& digits) {
        for (unsigned id = 0; id < digits.size(); ++id) {
          const KLMDigit* digit = digits[id];
          if (!digit->isGood())
            continue;
          if (!digitsSeen.insert(digit).second)
            continue;
          if (digit->getSubdetector() == KLMElementNumbers::c_EKLM && digit->getCharge() > 35.) {
            sumChargeEKLM += digit->getCharge();
            ++nDigitsClusterEKLM;
          } else if (digit->getSubdetector() == KLMElementNumbers::c_BKLM && !digit->inRPC() && digit->getCharge() > 35.) {
            sumChargeBKLM += digit->getCharge();
            ++nDigitsClusterBKLM;
          }
        }
      };
      RelationVector<KLMHit2d> clusterHits = cluster.getRelationsTo<KLMHit2d>();
      for (unsigned ih = 0; ih < clusterHits.size(); ++ih) {
        KLMHit2d* hit2d = clusterHits[ih];
        addDigits(hit2d->getRelationsTo<KLMDigit>());
        RelationVector<BKLMHit1d> hit1ds = hit2d->getRelationsTo<BKLMHit1d>();
        for (unsigned i1 = 0; i1 < hit1ds.size(); ++i1)
          addDigits(hit1ds[i1]->getRelationsTo<KLMDigit>());
      }
      if (clusterHits.size() > 1 && sumChargeBKLM > 0.) {
        m_ChargeClusterBKLM->Fill(sumChargeBKLM);
        m_AverageChargeClusterBKLM->Fill(sumChargeBKLM /
                                         static_cast<double>(nDigitsClusterBKLM));
      }
      if (clusterHits.size() > 1 && sumChargeEKLM > 0.) {
        m_ChargeClusterEKLM->Fill(sumChargeEKLM);
        m_AverageChargeClusterEKLM->Fill(sumChargeEKLM /
                                         static_cast<double>(nDigitsClusterEKLM));
      }
    }
  }
}
 
void KLMDQMModule::endRun()
{
}
 
void KLMDQMModule::terminate()
{
}