packages/CaloReco/CaloTowerTimeCalibration.cc

0001 #include "CaloTowerTimeCalibration.h"
0002
0003 #include <calobase/TowerInfo.h>
0004 #include <calobase/TowerInfoContainer.h>
0005
0006 #include <cdbobjects/CDBTTree.h>
0007
0008 #include <ffamodules/CDBInterface.h>
0009
0010 #include <globalvertex/GlobalVertex.h>
0011 #include <globalvertex/GlobalVertexMap.h>
0012
0013 #include <fun4all/Fun4AllHistoManager.h>
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015
0016 #include <qautils/QAHistManagerDef.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/PHObject.h>
0022 #include <phool/getClass.h>
0023 #include <phool/recoConsts.h>
0024
0025 #include <TAxis.h>
0026 #include <TH1.h>
0027 #include <TH2.h>
0028
0029 #include <cmath>
0030 #include <iostream>
0031 #include <limits>
0032 #include <stdexcept>
0033 #include <string>
0034 // Towers with invalid timing calibration are assigned NaN.
0035 // Downstream analyses should require std::isfinite(tower->get_time())
0036 // before using tower timing.
0037 //
0038 // ZS tower timing is not recalibrated; the standard calibrated
0039 // tower time is copied unchanged.
0040 namespace
0041 {
0042   constexpr float InvalidSentinelLimit = -998.5F;
0043
0044   bool IsUsableConstant(float value)
0045   {
0046     return std::isfinite(value) && value > InvalidSentinelLimit;
0047   }
0048
0049   std::string FoundText(const std::string &url)
0050   {
0051     return url.empty() ? " [missing]" : " [found]";
0052   }
0053
0054   template <class HistogramType>
0055   void StyleHistogram(HistogramType *histogram)
0056   {
0057     if (!histogram)
0058     {
0059       return;
0060     }
0061
0062     histogram->SetDirectory(nullptr);
0063
0064     // Store HIST as the default ROOT draw option for every QA histogram.
0065     // This keeps bin outlines/steps visible instead of point/error-bar drawing
0066     // when the histogram is opened without an explicit Draw() option.
0067     histogram->SetOption("HIST");
0068
0069     // ROOT font 62 is Helvetica bold. Use a slightly smaller stored axis
0070     // style so long labels like "(sample)" fit cleanly when users open the
0071     // histograms directly from the ROOT file.
0072     constexpr int boldFont = 62;
0073     constexpr float axisTitleSize = 0.045F;
0074     constexpr float axisLabelSize = 0.036F;
0075     constexpr float xTitleOffset = 0.95F;
0076     constexpr float yTitleOffset = 0.90F;
0077     constexpr float zTitleOffset = 0.95F;
0078
0079     TAxis *axes[] = {
0080         histogram->GetXaxis(),
0081         histogram->GetYaxis(),
0082         histogram->GetZaxis()};
0083
0084     for (TAxis *axis : axes)
0085     {
0086       if (!axis)
0087       {
0088         continue;
0089       }
0090
0091       axis->SetTitleFont(boldFont);
0092       axis->SetLabelFont(boldFont);
0093       axis->SetTitleSize(axisTitleSize);
0094       axis->SetLabelSize(axisLabelSize);
0095     }
0096
0097     if (histogram->GetXaxis())
0098     {
0099       histogram->GetXaxis()->SetTitleOffset(xTitleOffset);
0100     }
0101     if (histogram->GetYaxis())
0102     {
0103       histogram->GetYaxis()->SetTitleOffset(yTitleOffset);
0104     }
0105     if (histogram->GetZaxis())
0106     {
0107       histogram->GetZaxis()->SetTitleOffset(zTitleOffset);
0108     }
0109   }
0110 }  // namespace
0111
0112 CaloTowerTimeCalibration::CaloTowerTimeCalibration(const std::string &name)
0113   : SubsysReco(name)
0114   , m_qaEnergyThreshold(0.3F)
0115 {
0116 }
0117
0118 CaloTowerTimeCalibration::~CaloTowerTimeCalibration()
0119 {
0120   delete m_meanTimeCDB;
0121   delete m_timeCorrectionCDB;
0122 }
0123
0124 bool CaloTowerTimeCalibration::ResolveDetector()
0125 {
0126   if (m_detectorType == CaloTowerDefs::CEMC)
0127   {
0128     m_detector = "CEMC";
0129     return true;
0130   }
0131   if (m_detectorType == CaloTowerDefs::HCALIN)
0132   {
0133     m_detector = "HCALIN";
0134     return true;
0135   }
0136   if (m_detectorType == CaloTowerDefs::HCALOUT)
0137   {
0138     m_detector = "HCALOUT";
0139     return true;
0140   }
0141
0142   return false;
0143 }
0144
0145 void CaloTowerTimeCalibration::ResolveNames()
0146 {
0147  if (m_inputNodeName.empty())
0148  {
0149    m_inputNodeName = m_inputNodePrefix + m_detector;
0150  }
0151
0152  if (m_outputNodeName.empty())
0153  {
0154    m_outputNodeName = m_outputNodePrefix + m_detector;
0155  }
0156
0157  if (m_meanTimeCalibName.empty())
0158  {
0159    m_meanTimeCalibName = m_detector + "_meanTime";
0160  }
0161
0162  if (m_timeCorrectionCalibName.empty())
0163  {
0164    // Official custom tower timing calibration CDB domain.
0165    m_timeCorrectionCalibName =
0166        m_detector + "_towerTimeCalib";
0167  }
0168
0169 }
0170
0171
0172 int CaloTowerTimeCalibration::InitRun(PHCompositeNode *topNode)
0173 {
0174   if (!ResolveDetector())
0175   {
0176     std::cerr << Name() << ": unsupported detector type" << std::endl;
0177     return Fun4AllReturnCodes::ABORTRUN;
0178   }
0179
0180   ResolveNames();
0181
0182   delete m_meanTimeCDB;
0183   m_meanTimeCDB = nullptr;
0184   delete m_timeCorrectionCDB;
0185   m_timeCorrectionCDB = nullptr;
0186   m_timingInfo.clear();
0187   m_calibrationAvailable = false;
0188
0189   const std::string meanTimeURL =
0190     !m_directMeanTimeURL.empty()
0191           ? m_directMeanTimeURL
0192           : CDBInterface::instance()->getUrl(m_meanTimeCalibName);
0193
0194   const std::string timeCorrectionURL =
0195       !m_directTimeCorrectionURL.empty()
0196           ? m_directTimeCorrectionURL
0197           : CDBInterface::instance()->getUrl(
0198                 m_timeCorrectionCalibName);
0199
0200   if (!m_directMeanTimeURL.empty())
0201   {
0202     std::cout << Name() << "::" << m_detector
0203               << ": using direct mean-time payload "
0204               << meanTimeURL
0205               << std::endl;
0206   }
0207
0208   if (!m_directTimeCorrectionURL.empty())
0209   {
0210     std::cout << Name() << "::" << m_detector
0211               << ": using direct custom timing payload "
0212               << timeCorrectionURL
0213               << std::endl;
0214   }
0215
0216   if (meanTimeURL.empty() || timeCorrectionURL.empty())
0217   {
0218     std::cerr << Name() << "::" << m_detector
0219               << ": missing required timing calibration"
0220               << std::endl
0221               << "  mean-time CDB "
0222               << m_meanTimeCalibName
0223               << FoundText(meanTimeURL)
0224               << std::endl
0225               << "  custom timing CDB "
0226               << m_timeCorrectionCalibName
0227               << FoundText(timeCorrectionURL)
0228               << std::endl;
0229
0230     return Fun4AllReturnCodes::ABORTRUN;
0231   }
0232
0233   m_meanTimeCDB =
0234       new CDBTTree(meanTimeURL);
0235
0236   m_timeCorrectionCDB =
0237       new CDBTTree(timeCorrectionURL);
0238
0239   try
0240   {
0241     CreateNodeTree(topNode);
0242     CreateQAHistograms();
0243     LoadCalibration(topNode);
0244   }
0245   catch (const std::exception &error)
0246   {
0247     std::cerr << Name() << "::" << m_detector
0248               << ": " << error.what()
0249               << std::endl;
0250
0251     return Fun4AllReturnCodes::ABORTRUN;
0252   }
0253
0254   m_calibrationAvailable = true;
0255
0256   std::cout << Name() << "::" << m_detector
0257             << " input=" << m_inputNodeName
0258             << " output=" << m_outputNodeName
0259             << " mean_time_payload=" << meanTimeURL
0260             << " custom_payload=" << timeCorrectionURL
0261             << std::endl;
0262
0263   return Fun4AllReturnCodes::EVENT_OK;
0264 }
0265
0266 void CaloTowerTimeCalibration::CreateNodeTree(PHCompositeNode *topNode)
0267 {
0268   PHNodeIterator iterator(topNode);
0269
0270   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(
0271       iterator.findFirst("PHCompositeNode", "DST"));
0272   if (!dstNode)
0273   {
0274     throw std::runtime_error("DST node is missing");
0275   }
0276
0277   TowerInfoContainer *inputTowers =
0278       findNode::getClass<TowerInfoContainer>(topNode, m_inputNodeName);
0279   if (!inputTowers)
0280   {
0281     throw std::runtime_error(
0282         "missing input node " + m_inputNodeName
0283         + "; register this subsystem after standard CaloTowerCalib");
0284   }
0285
0286   PHNodeIterator dstIterator(dstNode);
0287
0288   PHCompositeNode *detectorNode = dynamic_cast<PHCompositeNode *>(
0289       dstIterator.findFirst("PHCompositeNode", m_detector));
0290   if (!detectorNode)
0291   {
0292     detectorNode = new PHCompositeNode(m_detector);
0293     dstNode->addNode(detectorNode);
0294   }
0295
0296   TowerInfoContainer *outputTowers =
0297       findNode::getClass<TowerInfoContainer>(topNode, m_outputNodeName);
0298
0299   if (!outputTowers)
0300   {
0301     outputTowers =
0302         dynamic_cast<TowerInfoContainer *>(inputTowers->CloneMe());
0303     if (!outputTowers)
0304     {
0305       throw std::runtime_error("failed to clone input tower container");
0306     }
0307
0308     auto *outputNode = new PHIODataNode<PHObject>(
0309         outputTowers,
0310         m_outputNodeName,
0311         "PHObject");
0312
0313     detectorNode->addNode(outputNode);
0314   }
0315 }
0316
0317 void CaloTowerTimeCalibration::LoadCalibration(PHCompositeNode *topNode)
0318 {
0319   TowerInfoContainer *inputTowers =
0320       findNode::getClass<TowerInfoContainer>(topNode, m_inputNodeName);
0321   if (!inputTowers)
0322   {
0323     throw std::runtime_error("missing input node " + m_inputNodeName);
0324   }
0325
0326   if (!m_meanTimeCDB || !m_timeCorrectionCDB)
0327   {
0328     throw std::runtime_error("timing CDB objects were not created");
0329   }
0330
0331   const unsigned int numberOfTowers = inputTowers->size();
0332   m_timingInfo.resize(numberOfTowers);
0333
0334   const int payloadTowerCount =
0335       m_timeCorrectionCDB->GetSingleIntValue("ntowers");
0336   const int formulaVersion =
0337       m_timeCorrectionCDB->GetSingleIntValue("formula_version");
0338   const int payloadRun =
0339       m_timeCorrectionCDB->GetSingleIntValue("runnumber");
0340
0341   recoConsts *recoConst = recoConsts::instance();
0342   const int requestedRun = recoConst->get_IntFlag("RUNNUMBER");
0343
0344   if (payloadTowerCount > 0
0345       && static_cast<unsigned int>(payloadTowerCount) != numberOfTowers)
0346   {
0347     throw std::runtime_error(
0348         "payload ntowers does not match " + m_inputNodeName);
0349   }
0350
0351   if (payloadRun > 0 && requestedRun > 0 && payloadRun != requestedRun)
0352   {
0353     throw std::runtime_error(
0354         "payload run " + std::to_string(payloadRun)
0355         + " does not match RUNNUMBER " + std::to_string(requestedRun));
0356   }
0357
0358   unsigned int invalidTowers = 0;
0359
0360   for (unsigned int channel = 0;
0361        channel < numberOfTowers;
0362        ++channel)
0363   {
0364     const unsigned int encodedKey = inputTowers->encode_key(channel);
0365     TimingInfo &timing = m_timingInfo[channel];
0366
0367     // Standard CaloTowerCalib used this value as:
0368     // standard_time = raw_time - mean_time.
0369     timing.meanTime =
0370         m_meanTimeCDB->GetFloatValue(encodedKey, "time");
0371
0372     // MakeAllTimingCDB_MV_test writes these values by flat channel.
0373     timing.phase1Shift =
0374         m_timeCorrectionCDB->GetFloatValue(channel, "phase1_shift");
0375     timing.sectorOffset =
0376         m_timeCorrectionCDB->GetFloatValue(channel, "sector_offset");
0377     timing.towerOffset =
0378         m_timeCorrectionCDB->GetFloatValue(channel, "tower_offset");
0379     timing.slewP0 =
0380         m_timeCorrectionCDB->GetFloatValue(channel, "slew_p0");
0381     timing.slewP1 =
0382         m_timeCorrectionCDB->GetFloatValue(channel, "slew_p1");
0383     timing.slewP2 =
0384         m_timeCorrectionCDB->GetFloatValue(channel, "slew_p2");
0385
0386     bool valid =
0387         IsUsableConstant(timing.meanTime)
0388         && IsUsableConstant(timing.phase1Shift)
0389         && IsUsableConstant(timing.sectorOffset)
0390         && IsUsableConstant(timing.towerOffset)
0391         && IsUsableConstant(timing.slewP0)
0392         && IsUsableConstant(timing.slewP1)
0393         && IsUsableConstant(timing.slewP2);
0394
0395     if (formulaVersion >= 3)
0396     {
0397       valid = valid
0398               && m_timeCorrectionCDB->GetIntValue(
0399                      channel,
0400                      "timing_valid") != 0;
0401     }
0402
0403     timing.valid = valid;
0404
0405     if (!valid)
0406     {
0407       ++invalidTowers;
0408     }
0409
0410   }
0411
0412   std::cout << Name() << "::" << m_detector
0413             << " loaded " << numberOfTowers
0414             << " channels; invalid=" << invalidTowers
0415             << ", formula_version=" << formulaVersion
0416             << std::endl;
0417 }
0418
0419 int CaloTowerTimeCalibration::process_event(PHCompositeNode *topNode)
0420 {
0421   if (!m_calibrationAvailable)
0422   {
0423     return Fun4AllReturnCodes::EVENT_OK;
0424   }
0425
0426   TowerInfoContainer *standardTowers =
0427       findNode::getClass<TowerInfoContainer>(topNode, m_inputNodeName);
0428   TowerInfoContainer *timingTowers =
0429       findNode::getClass<TowerInfoContainer>(topNode, m_outputNodeName);
0430
0431   if (!standardTowers || !timingTowers)
0432   {
0433     std::cerr << Name() << "::" << m_detector
0434               << ": required tower node is missing" << std::endl;
0435     return Fun4AllReturnCodes::ABORTEVENT;
0436   }
0437
0438   if (standardTowers->size() != timingTowers->size()
0439       || standardTowers->size() != m_timingInfo.size())
0440   {
0441     std::cerr << Name() << "::" << m_detector
0442               << ": tower-container size mismatch" << std::endl;
0443     return Fun4AllReturnCodes::ABORTEVENT;
0444   }
0445
0446   const unsigned int numberOfTowers = standardTowers->size();
0447
0448   // QA is filled only for events with a finite reconstructed global z vertex.
0449   // There is no |zvtx| magnitude cut and no zvtx-vs-time histogram.
0450   bool haveValidZVertex = false;
0451
0452   if (m_doQA)
0453   {
0454     GlobalVertexMap *vertexMap =
0455         findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0456
0457     if (vertexMap && !vertexMap->empty())
0458     {
0459       GlobalVertex *vertex = vertexMap->begin()->second;
0460       if (vertex)
0461       {
0462         haveValidZVertex = std::isfinite(vertex->get_z());
0463       }
0464     }
0465   }
0466
0467   for (unsigned int channel = 0;
0468        channel < numberOfTowers;
0469        ++channel)
0470   {
0471     TowerInfo *standardTower =
0472         standardTowers->get_tower_at_channel(channel);
0473     TowerInfo *timingTower =
0474         timingTowers->get_tower_at_channel(channel);
0475
0476     if (!standardTower || !timingTower)
0477     {
0478       continue;
0479     }
0480
0481     // Preserve calibrated energy and all metadata/status in the sidecar node.
0482     timingTower->copy_tower(standardTower);
0483
0484     // ZS timing is not custom-corrected. Fill dedicated ZS QA directly from
0485     // the standard calibrated tower, then leave the copied sidecar unchanged.
0486     // ZS QA is restricted to towers that pass get_isGood(). Do not apply
0487     // m_qaEnergyThreshold here: all energies of good ZS towers are retained.
0488     if (standardTower->get_isZS())
0489     {
0490       if (m_doQA && haveValidZVertex && standardTower->get_isGood())
0491       {
0492         const float zsTime = standardTower->get_time();
0493         const float zsEnergy = standardTower->get_energy();
0494
0495         if (std::isfinite(zsTime) && m_hZSTime)
0496         {
0497           m_hZSTime->Fill(zsTime);
0498         }
0499
0500         if (std::isfinite(zsEnergy) && m_hZSEnergy)
0501         {
0502           m_hZSEnergy->Fill(zsEnergy);
0503         }
0504
0505         if (std::isfinite(zsTime)
0506             && std::isfinite(zsEnergy)
0507             && m_hZSEnergyVsTime)
0508         {
0509           // Time is the x axis, matching the rest of the timing QA.
0510           m_hZSEnergyVsTime->Fill(zsTime, zsEnergy);
0511         }
0512       }
0513
0514       continue;
0515     }
0516
0517     const TimingInfo &timing = m_timingInfo[channel];
0518     if (!timing.valid)
0519     {
0520       timingTower->set_time(std::numeric_limits<float>::quiet_NaN());
0521       continue;
0522     }
0523
0524     const float standardTime = standardTower->get_time();
0525     const float calibratedEnergy = standardTower->get_energy();
0526
0527     if (!std::isfinite(standardTime)
0528         || !std::isfinite(calibratedEnergy))
0529     {
0530       timingTower->set_time(std::numeric_limits<float>::quiet_NaN());
0531       continue;
0532     }
0533
0534     // Standard CaloTowerCalib stores raw_time - official_mean_time.
0535     const float reconstructedRawTime =
0536         standardTime + timing.meanTime;
0537
0538     const float slewCorrection =
0539         timing.slewP0
0540         + timing.slewP1
0541               * std::exp(timing.slewP2 * calibratedEnergy);
0542
0543     const float correctedTime =
0544         reconstructedRawTime
0545         + timing.phase1Shift
0546         - timing.sectorOffset
0547         - timing.towerOffset
0548         - slewCorrection;
0549
0550     if (!std::isfinite(reconstructedRawTime)
0551         || !std::isfinite(slewCorrection)
0552         || !std::isfinite(correctedTime))
0553     {
0554       timingTower->set_time(std::numeric_limits<float>::quiet_NaN());
0555       continue;
0556     }
0557
0558     // Calibration is written before any QA quality definition is applied.
0559     timingTower->set_time(correctedTime);
0560
0561     if (!m_doQA
0562         || !haveValidZVertex
0563         || calibratedEnergy < m_qaEnergyThreshold)
0564     {
0565       continue;
0566     }
0567
0568     const auto fillTowerQA =
0569         [&](TowerQASet &qa)
0570         {
0571      if (qa.hStandardTime)
0572      {
0573      qa.hStandardTime->Fill(standardTime);
0574      }
0575
0576      if (qa.hCorrectedTime)
0577      {
0578      qa.hCorrectedTime->Fill(correctedTime);
0579      }
0580
0581      if (qa.hStandardEnergyVsTime)
0582      {
0583      qa.hStandardEnergyVsTime->Fill(standardTime, calibratedEnergy);
0584      }
0585
0586      if (qa.hCorrectedEnergyVsTime)
0587      {
0588      qa.hCorrectedEnergyVsTime->Fill(correctedTime, calibratedEnergy);
0589      }
0590
0591         };
0592
0593     // No get_isGood() requirement.
0594     fillTowerQA(m_qaAllTowers);
0595
0596     if (standardTower->get_isGood())
0597     {
0598       fillTowerQA(m_qaGoodTowers);
0599     }
0600   }
0601
0602   return Fun4AllReturnCodes::EVENT_OK;
0603 }
0604
0605 void CaloTowerTimeCalibration::CreateQAHistograms()
0606 {
0607   if (!m_doQA || m_qaHistogramsInitialized)
0608   {
0609     return;
0610   }
0611
0612   Fun4AllHistoManager *histogramManager =
0613       QAHistManagerDef::getHistoManager();
0614   if (!histogramManager)
0615   {
0616     throw std::runtime_error(
0617         Name() + "::" + m_detector
0618         + ": QAHistManagerDef returned a null histogram manager");
0619   }
0620
0621   const std::string prefix =
0622       "h_CaloTowerTimeCalibration_" + m_detector;
0623   const std::string title =
0624       "CaloTowerTimeCalibration " + m_detector;
0625
0626   constexpr int timeBins = 2000;
0627   constexpr float timeMin = -10.0F;
0628   constexpr float timeMax = 10.0F;
0629   constexpr int energyBins = 250;
0630   constexpr float energyMin = 0.0F;
0631   constexpr float energyMax = 25.0F;
0632
0633   // ZS towers are concentrated at low energy. Use 10 MeV (0.01 GeV) bins
0634   // over 0--1 GeV for the dedicated ZS energy QA.
0635   constexpr int zsEnergyBins = 100;
0636   constexpr float zsEnergyMin = 0.0F;
0637   constexpr float zsEnergyMax = 1.0F;
0638
0639   // Dedicated good-ZS QA. The timing sidecar intentionally does not modify ZS
0640   // tower times, so these are filled from standard towers passing get_isGood().
0641   m_hZSTime = new TH1F(
0642       (prefix + "_ZS_time").c_str(),
0643       (title + ", good ZS towers;ZS tower time (sample);towers").c_str(),
0644       timeBins,
0645       timeMin,
0646       timeMax);
0647
0648   m_hZSEnergy = new TH1F(
0649       (prefix + "_ZS_energy").c_str(),
0650       (title + ", good ZS towers;ZS tower calibrated energy (GeV);towers").c_str(),
0651       zsEnergyBins,
0652       zsEnergyMin,
0653       zsEnergyMax);
0654
0655   m_hZSEnergyVsTime = new TH2F(
0656       (prefix + "_ZS_energy_vs_time").c_str(),
0657       (title
0658        + ", good ZS towers;ZS tower time (sample);ZS tower calibrated energy (GeV)")
0659           .c_str(),
0660       timeBins,
0661       timeMin,
0662       timeMax,
0663       zsEnergyBins,
0664       zsEnergyMin,
0665       zsEnergyMax);
0666
0667   StyleHistogram(m_hZSTime);
0668   StyleHistogram(m_hZSEnergy);
0669   StyleHistogram(m_hZSEnergyVsTime);
0670
0671   histogramManager->registerHisto(m_hZSTime);
0672   histogramManager->registerHisto(m_hZSEnergy);
0673   histogramManager->registerHisto(m_hZSEnergyVsTime);
0674
0675   const auto bookTowerQASet =
0676       [&](TowerQASet &qa,
0677           const std::string &nameSuffix,
0678           const std::string &titleSuffix)
0679       {
0680         const std::string qaPrefix = prefix + nameSuffix;
0681         const std::string qaTitle = title + titleSuffix;
0682
0683         qa.hStandardTime = new TH1F(
0684             (qaPrefix + "_standard_time").c_str(),
0685             (qaTitle + ";standard time (sample);towers").c_str(),
0686             timeBins,
0687             timeMin,
0688             timeMax);
0689
0690         qa.hCorrectedTime = new TH1F(
0691             (qaPrefix + "_corrected_time").c_str(),
0692             (qaTitle + ";custom corrected time (sample);towers").c_str(),
0693             timeBins,
0694             timeMin,
0695             timeMax);
0696
0697         qa.hStandardEnergyVsTime = new TH2F(
0698             (qaPrefix + "_standard_energy_vs_time").c_str(),
0699             (qaTitle + ";standard time (sample);calibrated energy (GeV)").c_str(),
0700             timeBins,
0701             timeMin,
0702             timeMax,
0703             energyBins,
0704             energyMin,
0705             energyMax);
0706
0707         qa.hCorrectedEnergyVsTime = new TH2F(
0708             (qaPrefix + "_corrected_energy_vs_time").c_str(),
0709             (qaTitle + ";custom corrected time (sample);calibrated energy (GeV)").c_str(),
0710             timeBins,
0711             timeMin,
0712             timeMax,
0713             energyBins,
0714             energyMin,
0715             energyMax);
0716
0717         StyleHistogram(qa.hStandardTime);
0718         StyleHistogram(qa.hCorrectedTime);
0719         StyleHistogram(qa.hStandardEnergyVsTime);
0720         StyleHistogram(qa.hCorrectedEnergyVsTime);
0721
0722         histogramManager->registerHisto(qa.hStandardTime);
0723         histogramManager->registerHisto(qa.hCorrectedTime);
0724         histogramManager->registerHisto(qa.hStandardEnergyVsTime);
0725         histogramManager->registerHisto(qa.hCorrectedEnergyVsTime);
0726
0727         // No |zvtx| magnitude-cut histograms are booked. These inclusive
0728         // histograms are filled only when the event has a finite global z vertex.
0729       };
0730
0731   bookTowerQASet(
0732       m_qaAllTowers,
0733       "",
0734       ", all accepted towers (no get_isGood cut)");
0735
0736   bookTowerQASet(
0737       m_qaGoodTowers,
0738       "_isGood",
0739       ", strict TowerInfo::get_isGood() towers");
0740
0741   // Intentionally no Sumw2 calls; all histograms are filled with unit weight.
0742   m_qaHistogramsInitialized = true;
0743
0744   std::cout << Name() << "::" << m_detector
0745             << " booked all-tower and strict-get_isGood energy/time QA sets"
0746             << "; plus get_isGood ZS time, energy, and energy-vs-time QA"
0747             << "; QA requires a finite global z vertex but applies no |zvtx| cut"
0748             << std::endl;
0749 }