File indexing completed on 2025-08-06 08:18:41
0001 #include "GlobalQA.h"
0002
0003
0004 #include <calobase/TowerInfo.h>
0005 #include <calobase/TowerInfoContainer.h>
0006 #include <calobase/TowerInfoDefs.h>
0007
0008 #include <mbd/MbdPmtContainer.h>
0009 #include <mbd/MbdPmtHit.h>
0010
0011 #include <zdcinfo/Zdcinfo.h>
0012
0013 #include <globalvertex/MbdVertex.h>
0014 #include <globalvertex/MbdVertexMap.h>
0015
0016 #include <qautils/QAHistManagerDef.h>
0017
0018 #include <ffarawobjects/Gl1Packet.h>
0019
0020 #include <ffamodules/CDBInterface.h>
0021
0022 #include <cdbobjects/CDBTTree.h> // for CDBTTree
0023
0024 #include <fun4all/Fun4AllHistoManager.h>
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h> // for PHWHERE
0029
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <TProfile2D.h>
0033 #include <TSystem.h>
0034
0035 #include <boost/format.hpp>
0036
0037 #include <algorithm> // for sort
0038 #include <cassert>
0039 #include <cmath> // for log10, pow, sqrt, abs, M_PI
0040 #include <cstdint>
0041 #include <iostream> // for operator<<, endl, basic_...
0042 #include <limits>
0043 #include <map> // for operator!=, _Rb_tree_con...
0044 #include <string>
0045 #include <utility> // for pair
0046
0047 GlobalQA::GlobalQA(const std::string &name)
0048 : SubsysReco(name)
0049 , detector("HCALIN")
0050 {
0051 }
0052
0053 GlobalQA::~GlobalQA()
0054 {
0055 }
0056
0057 int GlobalQA::Init(PHCompositeNode * )
0058 {
0059 if (m_debug)
0060 {
0061 std::cout << "In GlobalQA::Init" << std::endl;
0062 }
0063
0064 createHistos();
0065
0066 if (m_debug)
0067 {
0068 std::cout << "Leaving GlobalQA::Init" << std::endl;
0069 }
0070
0071 return Fun4AllReturnCodes::EVENT_OK;
0072 }
0073
0074 int GlobalQA::process_event(PHCompositeNode *topNode)
0075 {
0076 _eventcounter++;
0077 process_towers(topNode);
0078
0079 return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 int GlobalQA::process_towers(PHCompositeNode *topNode)
0083 {
0084 if (m_debug)
0085 {
0086 std::cout << _eventcounter << std::endl;
0087 }
0088
0089
0090 Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, 14001);
0091 if (!gl1PacketInfo)
0092 {
0093 gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0094 if (!gl1PacketInfo)
0095 {
0096 std::cout << PHWHERE << "GlobalQA::process_event: GL1Packet node is missing" << std::endl;
0097 }
0098 }
0099
0100 uint64_t triggervec = 0;
0101 if (gl1PacketInfo)
0102 {
0103 triggervec = gl1PacketInfo->getScaledVector();
0104 for (int i = 0; i < 64; i++)
0105 {
0106 bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0107
0108 if (trig_decision)
0109 {
0110 h_GlobalQA_triggerVec->Fill(i);
0111 }
0112 triggervec = (triggervec >> 1U) & 0xffffffffU;
0113 }
0114 triggervec = gl1PacketInfo->getScaledVector();
0115 }
0116
0117 if ((triggervec >> 0xAU) & 0x1U)
0118 {
0119
0120 MbdVertexMap *mbdmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0121 MbdVertex *bvertex = nullptr;
0122 float mbd_zvtx = std::numeric_limits<float>::quiet_NaN();
0123 if (mbdmap)
0124 {
0125 for (MbdVertexMap::ConstIter mbditer = mbdmap->begin(); mbditer != mbdmap->end(); ++mbditer)
0126 {
0127 bvertex = mbditer->second;
0128 }
0129 if (bvertex)
0130 {
0131 mbd_zvtx = bvertex->get_z();
0132 }
0133 }
0134 h_GlobalQA_mbd_zvtx->Fill(mbd_zvtx);
0135 h_GlobalQA_mbd_zvtx_wide->Fill(mbd_zvtx);
0136 if (!std::isfinite(mbd_zvtx))
0137 {
0138 h_GlobalQA_mbd_zvtxq->SetBinContent(1, h_GlobalQA_mbd_zvtxq->GetBinContent(1) + 1);
0139 }
0140 else
0141 {
0142 h_GlobalQA_mbd_zvtxq->SetBinContent(2, h_GlobalQA_mbd_zvtxq->GetBinContent(2) + 1);
0143 }
0144 }
0145
0146
0147
0148 if (triggervec & mbdtrig)
0149 {
0150
0151 TowerInfoContainer *_sepd_towerinfo = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_SEPD");
0152 unsigned int ntowers = 0;
0153 if (_sepd_towerinfo)
0154 {
0155 ntowers = _sepd_towerinfo->size();
0156 }
0157 if (ntowers != 744)
0158 {
0159 std::cout << "sEPD container has unexpected size - exiting now!" << std::endl;
0160 gSystem->Exit(1);
0161 }
0162
0163 float sepdsouthadcsum = 0.;
0164 float sepdnorthadcsum = 0.;
0165 if (_sepd_towerinfo)
0166 {
0167 for (unsigned int i = 0; i < ntowers; i++)
0168 {
0169 TowerInfo *_tower = _sepd_towerinfo->get_tower_at_channel(i);
0170 float _e = _tower->get_energy();
0171 bool isZS = _tower->get_isZS();
0172 unsigned int key = TowerInfoDefs::encode_epd(i);
0173 int arm = TowerInfoDefs::get_epd_arm(key);
0174 float rbin = (float) (TowerInfoDefs::get_epd_rbin(key));
0175 float phibin = (float) (TowerInfoDefs::get_epd_phibin(key));
0176 if (!isZS)
0177 {
0178 if (arm == 0)
0179 {
0180 sepdsouthadcsum += _e;
0181 h2Profile_GlobalQA_sEPD_tiles_south->Fill(rbin, phibin, _e);
0182 }
0183 else if (arm == 1)
0184 {
0185 sepdnorthadcsum += _e;
0186 h2Profile_GlobalQA_sEPD_tiles_north->Fill(rbin, phibin, _e);
0187 }
0188 }
0189 }
0190
0191 h_GlobalQA_sEPD_adcsum_s->Fill(sepdsouthadcsum);
0192 h_GlobalQA_sEPD_adcsum_n->Fill(sepdnorthadcsum);
0193 h2_GlobalQA_sEPD_adcsum_ns->Fill(sepdsouthadcsum, sepdnorthadcsum);
0194 }
0195 }
0196
0197 if ((triggervec >> 0x3U) & 0x1U)
0198 {
0199
0200
0201
0202 Zdcinfo *zdcinfo = findNode::getClass<Zdcinfo>(topNode, "Zdcinfo");
0203 float totalzdcsouthcalib = 0.;
0204 float totalzdcnorthcalib = 0.;
0205 if (zdcinfo)
0206 {
0207 totalzdcsouthcalib = zdcinfo->get_zdc_energy(0);
0208 totalzdcnorthcalib = zdcinfo->get_zdc_energy(1);
0209 zdc_zvtx = zdcinfo->get_zvertex();
0210 h_GlobalQA_zdc_zvtx->Fill(zdc_zvtx);
0211 h_GlobalQA_zdc_zvtx_wide->Fill(zdc_zvtx);
0212 h_GlobalQA_zdc_energy_s->Fill(totalzdcsouthcalib);
0213 h_GlobalQA_zdc_energy_n->Fill(totalzdcnorthcalib);
0214 }
0215 }
0216
0217 if ((triggervec >> 0xAU) & 0x1U)
0218 {
0219
0220 MbdPmtContainer *bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0221 if (!bbcpmts)
0222 {
0223 std::cout << "GlobalQA::process_event: Could not find MbdPmtContainer,"
0224 << std::endl;
0225
0226 }
0227
0228
0229 int hits_n = 0;
0230 int hits_s = 0;
0231 int hits_n_t = 0;
0232 int hits_s_t = 0;
0233 std::vector<float> time_sum_s;
0234 std::vector<float> time_sum_n;
0235 float sum_s = 0.;
0236 float sum_n = 0.;
0237 float sum_s2 = 0.;
0238 float sum_n2 = 0.;
0239 float tot_charge_s = 0.;
0240 float tot_charge_n = 0.;
0241
0242 float charge_thresh = 0.4;
0243 if (bbcpmts)
0244 {
0245 int nPMTs = bbcpmts->get_npmt();
0246 for (int i = 0; i < nPMTs; i++)
0247 {
0248 MbdPmtHit *mbdpmt = bbcpmts->get_pmt(i);
0249 float q = mbdpmt->get_q();
0250 float t = mbdpmt->get_time();
0251 if (i < 64)
0252 {
0253 tot_charge_s += q;
0254 if (q > charge_thresh)
0255 {
0256 hits_s++;
0257 }
0258 if (i == 56 || std::isnan(t))
0259 {
0260 continue;
0261 }
0262 hits_s_t++;
0263 time_sum_s.push_back(t);
0264 sum_s += t;
0265 sum_s2 += t * t;
0266 }
0267 else if (i >= 64)
0268 {
0269 tot_charge_n += q;
0270 if (q > charge_thresh)
0271 {
0272 hits_n++;
0273 }
0274 if (i == 120 || std::isnan(t))
0275 {
0276 continue;
0277 }
0278 hits_n_t++;
0279 time_sum_n.push_back(t);
0280 sum_n += t;
0281 sum_n2 += t * t;
0282 }
0283
0284
0285
0286 }
0287 }
0288
0289
0290 std::sort(time_sum_n.begin(), time_sum_n.end());
0291 std::sort(time_sum_s.begin(), time_sum_s.end());
0292 unsigned length_s = time_sum_s.size();
0293 unsigned length_n = time_sum_n.size();
0294 float mean_north = 999;
0295 float mean_south = 999;
0296 int central_cut = 4;
0297 float sigma_cut = 1.5;
0298
0299 if (hits_s_t >= central_cut)
0300 {
0301 mean_south = sum_s / static_cast<float>(hits_s_t);
0302 float rms_s = std::sqrt((sum_s2 / static_cast<float>(hits_s_t)) - (mean_south * mean_south));
0303 int nhit_s_center = 0;
0304 float sum_s_center = 0.;
0305
0306 for (unsigned int is = 0; is < length_s; is++)
0307 {
0308 if (std::fabs(time_sum_s.at(is) - mean_south) < sigma_cut * rms_s)
0309 {
0310 sum_s_center += time_sum_s.at(is);
0311 nhit_s_center++;
0312 }
0313 }
0314
0315 if (nhit_s_center > 0)
0316 {
0317 float mean_south_center =
0318 sum_s_center / static_cast<float>(nhit_s_center);
0319 mean_south = mean_south_center;
0320 }
0321 }
0322 else if (hits_s >= 2 && (hits_s_t >= 1))
0323 {
0324 mean_south = sum_s / static_cast<float>(hits_s_t);
0325 }
0326
0327 if (hits_n_t >= central_cut)
0328 {
0329 mean_north = sum_n / static_cast<float>(hits_n_t);
0330 float rms_n = std::sqrt((sum_n2 / static_cast<float>(hits_n_t)) - (mean_north * mean_north));
0331 int nhit_n_center = 0;
0332 float sum_n_center = 0.;
0333
0334 for (unsigned int ino = 0; ino < length_n; ino++)
0335 {
0336 if (std::abs(time_sum_n.at(ino) - mean_north) < sigma_cut * rms_n)
0337 {
0338 sum_n_center += time_sum_n.at(ino);
0339 nhit_n_center++;
0340 }
0341 }
0342
0343 if (nhit_n_center > 0)
0344 {
0345 float mean_north_center = sum_n_center / static_cast<float>(nhit_n_center);
0346 mean_north = mean_north_center;
0347 }
0348 }
0349 else if (hits_n >= 2 && hits_n_t >= 1)
0350 {
0351 mean_north = sum_n / static_cast<float>(hits_n_t);
0352 }
0353 float calc_zvtx;
0354 if (mean_north != 999 && mean_south != 999)
0355 {
0356 calc_zvtx = 15 * (mean_south - mean_north);
0357 }
0358 else
0359 {
0360 calc_zvtx = 999;
0361 }
0362
0363 h_GlobalQA_calc_zvtx->Fill(calc_zvtx);
0364 h_GlobalQA_calc_zvtx_wide->Fill(calc_zvtx);
0365 h_GlobalQA_mbd_charge_s->Fill(tot_charge_s);
0366 h_GlobalQA_mbd_charge_n->Fill(tot_charge_n);
0367 h_GlobalQA_mbd_nhit_s->Fill(hits_s);
0368 h_GlobalQA_mbd_nhit_n->Fill(hits_n);
0369 h_GlobalQA_mbd_charge_sum->Fill(tot_charge_s + tot_charge_n);
0370 h2_GlobalQA_mbd_charge_NS_correlation->Fill(tot_charge_s, tot_charge_n);
0371 h2_GlobalQA_mbd_nhits_NS_correlation->Fill(hits_s, hits_n);
0372 }
0373
0374 return Fun4AllReturnCodes::EVENT_OK;
0375 }
0376
0377 void GlobalQA::createHistos()
0378 {
0379 auto *hm = QAHistManagerDef::getHistoManager();
0380 assert(hm);
0381
0382
0383 h_GlobalQA_mbd_zvtxq =
0384 new TH1D("h_GlobalQA_mbd_zvtxq",
0385 ";Scaled Trigger 10: MBD Coincidence Has zvtx?;percentage", 2,
0386 -0.5, 1.5);
0387 h_GlobalQA_mbd_zvtx = new TH1D(
0388 "h_GlobalQA_mbd_zvtx", ";Scaled Trigger 10: MBD Coincidence zvtx [cm]",
0389 100, -50, 50);
0390 h_GlobalQA_mbd_zvtx_wide = new TH1D(
0391 "h_GlobalQA_mbd_zvtx_wide",
0392 ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -300, 300);
0393 h_GlobalQA_calc_zvtx = new TH1D(
0394 "h_GlobalQA_calc_zvtx",
0395 ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -50, 50);
0396 h_GlobalQA_calc_zvtx_wide = new TH1D(
0397 "h_GlobalQA_calc_zvtx_wide",
0398 ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -300, 300);
0399 h_GlobalQA_mbd_charge_s =
0400 new TH1D("h_GlobalQA_mbd_charge_s",
0401 ";Scaled Trigger 10: MBD Coincidence charge", 100, 0, 10);
0402 h_GlobalQA_mbd_charge_n =
0403 new TH1D("h_GlobalQA_mbd_charge_n",
0404 ";Scaled Trigger 10: MBD Coincidence charge", 100, 0, 10);
0405 h_GlobalQA_mbd_nhit_s =
0406 new TH1D("h_GlobalQA_mbd_nhit_s",
0407 ";Scaled Trigger 10: MBD Coincidence nhit", 30, -0.5, 29.5);
0408 h_GlobalQA_mbd_nhit_n =
0409 new TH1D("h_GlobalQA_mbd_nhit_n",
0410 ";Scaled Trigger 10: MBD Coincidence nhit", 30, -0.5, 29.5);
0411
0412 h_GlobalQA_mbd_charge_sum =
0413 new TH1F("h_GlobalQA_mbd_charge_sum", " ; MBD Total Charge ; Counts",
0414 100, 0., 20);
0415
0416 h2_GlobalQA_mbd_charge_NS_correlation = new TH2F(
0417 "h2_GlobalQA_mbd_charge_NS_correlation",
0418 "MBD Charge Correlation ; Total Charge (South); Total Charge (North)",
0419 150, 0, 1500, 150, 0, 1500);
0420
0421 h2_GlobalQA_mbd_nhits_NS_correlation =
0422 new TH2F("h2_GlobalQA_mbd_nhits_NS_correlation",
0423 "MBD Number Of Hits Correlation ; Number Of Hits (South); "
0424 "Number Of Hits (North)",
0425 70, 0., 70, 70, 0., 70);
0426
0427 hm->registerHisto(h_GlobalQA_mbd_zvtx);
0428 hm->registerHisto(h_GlobalQA_mbd_zvtxq);
0429 hm->registerHisto(h_GlobalQA_mbd_zvtx_wide);
0430 hm->registerHisto(h_GlobalQA_calc_zvtx);
0431 hm->registerHisto(h_GlobalQA_calc_zvtx_wide);
0432 hm->registerHisto(h_GlobalQA_mbd_charge_s);
0433 hm->registerHisto(h_GlobalQA_mbd_charge_n);
0434 hm->registerHisto(h_GlobalQA_mbd_nhit_s);
0435 hm->registerHisto(h_GlobalQA_mbd_nhit_n);
0436 hm->registerHisto(h_GlobalQA_mbd_charge_sum);
0437 hm->registerHisto(h2_GlobalQA_mbd_charge_NS_correlation);
0438 hm->registerHisto(h2_GlobalQA_mbd_nhits_NS_correlation);
0439
0440
0441 h_GlobalQA_sEPD_adcsum_s = new TH1D(
0442 "h_GlobalQA_sEPD_adcsum_s", " ; sEPD ADC sum ; Counts", 200, 0, 6e6);
0443
0444 h_GlobalQA_sEPD_adcsum_n = new TH1D(
0445 "h_GlobalQA_sEPD_adcsum_n", " ; sEPD ADC sum ; Counts", 200, 0, 6e6);
0446
0447 h2_GlobalQA_sEPD_adcsum_ns =
0448 new TH2D("h2_GlobalQA_sEPD_adcsum_ns",
0449 "sEPD NS ADC sum correlation ; ADC sum (south); ADC sum (north)",
0450 500, 0, 6e6, 500, 0, 6e6);
0451
0452 h2Profile_GlobalQA_sEPD_tiles_north =
0453 new TProfile2D("h2Profile_GlobalQA_sEPD_tiles_north",
0454 "sEPD Tile Mean Energy (north); #eta; #phi", 16, -0.5,
0455 15.5, 24, -0.5, 23.5);
0456
0457 h2Profile_GlobalQA_sEPD_tiles_south =
0458 new TProfile2D("h2Profile_GlobalQA_sEPD_tiles_south",
0459 "sEPD Tile Mean Energy (south); #eta; #phi", 16, -0.5,
0460 15.5, 24, -0.5, 23.5);
0461
0462 h2_GlobalQA_sEPD_ADC_channel_south =
0463 new TH2D("h2_GlobalQA_sEPD_ADC_channel_south",
0464 "h2_GlobalQA_sEPD_ADC_channel_south ; #eta; #phi", 16, -0.5,
0465 15.5, 24, -0.5, 23.5);
0466 h2_GlobalQA_sEPD_ADC_channel_north =
0467 new TH2D("h2_GlobalQA_sEPD_ADC_channel_north",
0468 "h2_GlobalQA_sEPD_ADC_channel_north ; #eta; #phi", 16, -0.5,
0469 15.5, 24, -0.5, 23.5);
0470
0471
0472 hm->registerHisto(h_GlobalQA_sEPD_adcsum_s);
0473 hm->registerHisto(h_GlobalQA_sEPD_adcsum_n);
0474 hm->registerHisto(h2_GlobalQA_sEPD_adcsum_ns);
0475 hm->registerHisto(h2Profile_GlobalQA_sEPD_tiles_south);
0476 hm->registerHisto(h2Profile_GlobalQA_sEPD_tiles_north);
0477 hm->registerHisto(h2_GlobalQA_sEPD_ADC_channel_south);
0478 hm->registerHisto(h2_GlobalQA_sEPD_ADC_channel_north);
0479
0480
0481 h_GlobalQA_zdc_zvtx = new TH1D(
0482 "h_GlobalQA_zdc_zvtx", ";Scaled Trigger 3: ZDC Coincidence zvtx [cm]",
0483 100, -1000, 1000);
0484 h_GlobalQA_zdc_zvtx_wide = new TH1D(
0485 "h_GlobalQA_zdc_zvtx_wide",
0486 ";Scaled Trigger 3: ZDC Coincidence zvtx [cm]", 100, -2000, 2000);
0487 h_GlobalQA_zdc_energy_s = new TH1D(
0488 "h_GlobalQA_zdc_energy_s",
0489 ";Scaled Trigger 3: ZDC Coincidence Energy [Gev]", 100, 10, 510);
0490 h_GlobalQA_zdc_energy_n = new TH1D(
0491 "h_GlobalQA_zdc_energy_n",
0492 ";Scaled Trigger 3: ZDC Coincidence Energy [Gev]", 100, 10, 510);
0493 hm->registerHisto(h_GlobalQA_zdc_zvtx);
0494 hm->registerHisto(h_GlobalQA_zdc_zvtx_wide);
0495 hm->registerHisto(h_GlobalQA_zdc_energy_s);
0496 hm->registerHisto(h_GlobalQA_zdc_energy_n);
0497
0498 h_GlobalQA_triggerVec = new TH1F("h_GlobalQA_triggerVec", "", 65, -0.5, 64.5);
0499 h_GlobalQA_triggerVec->SetDirectory(nullptr);
0500 hm->registerHisto(h_GlobalQA_triggerVec);
0501 }
