File indexing completed on 2025-12-17 09:21:16
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 int GlobalQA::Init(PHCompositeNode * )
0054 {
0055 if (m_debug)
0056 {
0057 std::cout << "In GlobalQA::Init" << std::endl;
0058 }
0059
0060 createHistos();
0061
0062 if (m_debug)
0063 {
0064 std::cout << "Leaving GlobalQA::Init" << std::endl;
0065 }
0066
0067 return Fun4AllReturnCodes::EVENT_OK;
0068 }
0069
0070 int GlobalQA::process_event(PHCompositeNode *topNode)
0071 {
0072 _eventcounter++;
0073 if (m_debug)
0074 {
0075 std::cout << _eventcounter << std::endl;
0076 }
0077
0078
0079 Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, 14001);
0080 if (!gl1PacketInfo)
0081 {
0082 gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0083 if (!gl1PacketInfo)
0084 {
0085 std::cout << PHWHERE << "GlobalQA::process_event: GL1Packet node is missing" << std::endl;
0086 }
0087 }
0088
0089 triggervec = 0;
0090 if (gl1PacketInfo)
0091 {
0092 triggervec = gl1PacketInfo->getScaledVector();
0093 for (int i = 0; i < 64; i++)
0094 {
0095 bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0096
0097 if (trig_decision)
0098 {
0099 h_GlobalQA_triggerVec->Fill(i);
0100 }
0101 triggervec = (triggervec >> 1U) & 0xffffffffU;
0102 }
0103 triggervec = gl1PacketInfo->getScaledVector();
0104 }
0105
0106 process_towers(topNode);
0107
0108 if ( triggervec & mbdtrig )
0109 {
0110 process_mbd(topNode);
0111 }
0112
0113 return Fun4AllReturnCodes::EVENT_OK;
0114 }
0115
0116 int GlobalQA::process_towers(PHCompositeNode *topNode)
0117 {
0118
0119
0120
0121 TowerInfoContainer *_sepd_towerinfo = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_SEPD");
0122 unsigned int ntowers = 0;
0123 if (_sepd_towerinfo)
0124 {
0125 ntowers = _sepd_towerinfo->size();
0126 if (ntowers != 744)
0127 {
0128 std::cout << "sEPD container has unexpected size - skipping sEPD!" << std::endl;
0129
0130 }
0131 }
0132 if ( (triggervec & mbdtrig) && (ntowers == 744) )
0133 {
0134
0135
0136 float sepdsouthadcsum = 0.;
0137 float sepdnorthadcsum = 0.;
0138 if (_sepd_towerinfo)
0139 {
0140 for (unsigned int i = 0; i < ntowers; i++)
0141 {
0142 TowerInfo *_tower = _sepd_towerinfo->get_tower_at_channel(i);
0143 float _e = _tower->get_energy();
0144 bool isZS = _tower->get_isZS();
0145 unsigned int key = TowerInfoDefs::encode_epd(i);
0146 int arm = TowerInfoDefs::get_epd_arm(key);
0147 float rbin = (float) (TowerInfoDefs::get_epd_rbin(key));
0148 float phibin = (float) (TowerInfoDefs::get_epd_phibin(key));
0149 if (!isZS)
0150 {
0151 if (arm == 0)
0152 {
0153 sepdsouthadcsum += _e;
0154 h2Profile_GlobalQA_sEPD_tiles_south->Fill(rbin, phibin, _e);
0155 }
0156 else if (arm == 1)
0157 {
0158 sepdnorthadcsum += _e;
0159 h2Profile_GlobalQA_sEPD_tiles_north->Fill(rbin, phibin, _e);
0160 }
0161 }
0162 }
0163
0164 h_GlobalQA_sEPD_adcsum_s->Fill(sepdsouthadcsum);
0165 h_GlobalQA_sEPD_adcsum_n->Fill(sepdnorthadcsum);
0166 h2_GlobalQA_sEPD_adcsum_ns->Fill(sepdsouthadcsum, sepdnorthadcsum);
0167 }
0168 }
0169
0170
0171
0172
0173 Zdcinfo *zdcinfo = findNode::getClass<Zdcinfo>(topNode, "Zdcinfo");
0174 float totalzdcsouthcalib = 0.;
0175 float totalzdcnorthcalib = 0.;
0176 if (zdcinfo)
0177 {
0178 totalzdcsouthcalib = zdcinfo->get_zdc_energy(0);
0179 totalzdcnorthcalib = zdcinfo->get_zdc_energy(1);
0180 zdc_zvtx = zdcinfo->get_zvertex();
0181 h_GlobalQA_zdc_zvtx->Fill(zdc_zvtx);
0182 h_GlobalQA_zdc_zvtx_wide->Fill(zdc_zvtx);
0183 h_GlobalQA_zdc_energy_s->Fill(totalzdcsouthcalib);
0184 h_GlobalQA_zdc_energy_n->Fill(totalzdcnorthcalib);
0185 }
0186
0187
0188 return Fun4AllReturnCodes::EVENT_OK;
0189 }
0190
0191
0192 int GlobalQA::process_mbd(PHCompositeNode *topNode)
0193 {
0194
0195 MbdVertexMap *mbdmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0196 MbdVertex *bvertex = nullptr;
0197 float mbd_zvtx = std::numeric_limits<float>::quiet_NaN();
0198 if (mbdmap)
0199 {
0200 for (MbdVertexMap::ConstIter mbditer = mbdmap->begin(); mbditer != mbdmap->end(); ++mbditer)
0201 {
0202 bvertex = mbditer->second;
0203 }
0204 if (bvertex)
0205 {
0206 mbd_zvtx = bvertex->get_z();
0207 }
0208 }
0209 h_GlobalQA_mbd_zvtx->Fill(mbd_zvtx);
0210 h_GlobalQA_mbd_zvtx_wide->Fill(mbd_zvtx);
0211 if (!std::isfinite(mbd_zvtx))
0212 {
0213 h_GlobalQA_mbd_zvtxq->SetBinContent(1, h_GlobalQA_mbd_zvtxq->GetBinContent(1) + 1);
0214 }
0215 else
0216 {
0217 h_GlobalQA_mbd_zvtxq->SetBinContent(2, h_GlobalQA_mbd_zvtxq->GetBinContent(2) + 1);
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 return Fun4AllReturnCodes::EVENT_OK;
0374 }
0375
0376
0377 void GlobalQA::createHistos()
0378 {
0379 auto *hm = QAHistManagerDef::getHistoManager();
0380 assert(hm);
0381
0382
0383 h_GlobalQA_mbd_zvtxq = new TH1D("h_GlobalQA_mbd_zvtxq", ";Scaled Trigger 10: MBD Coincidence Has zvtx?;percentage", 2, -0.5, 1.5);
0384
0385 h_GlobalQA_mbd_zvtx = new TH1D( "h_GlobalQA_mbd_zvtx", ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -50, 50);
0386
0387 h_GlobalQA_mbd_zvtx_wide = new TH1D( "h_GlobalQA_mbd_zvtx_wide", ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -300, 300);
0388
0389 h_GlobalQA_calc_zvtx = new TH1D("h_GlobalQA_calc_zvtx", ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -50, 50);
0390
0391 h_GlobalQA_calc_zvtx_wide = new TH1D("h_GlobalQA_calc_zvtx_wide", ";Scaled Trigger 10: MBD Coincidence zvtx [cm]", 100, -300, 300);
0392
0393 h_GlobalQA_mbd_charge_s = new TH1D("h_GlobalQA_mbd_charge_s", ";Scaled Trigger 10: MBD Coincidence charge", 1500, 0, 1500);
0394
0395 h_GlobalQA_mbd_charge_n = new TH1D("h_GlobalQA_mbd_charge_n", ";Scaled Trigger 10: MBD Coincidence charge", 1500, 0, 1500);
0396
0397 h_GlobalQA_mbd_nhit_s = new TH1D("h_GlobalQA_mbd_nhit_s", ";Scaled Trigger 10: MBD Coincidence nhit", 65, -0.5, 64.5);
0398
0399 h_GlobalQA_mbd_nhit_n = new TH1D("h_GlobalQA_mbd_nhit_n", ";Scaled Trigger 10: MBD Coincidence nhit", 65, -0.5, 64.5);
0400
0401 h_GlobalQA_mbd_charge_sum = new TH1F("h_GlobalQA_mbd_charge_sum", " ; MBD Total Charge ; Counts", 2500, 0., 2500);
0402
0403 h2_GlobalQA_mbd_charge_NS_correlation = new TH2F("h2_GlobalQA_mbd_charge_NS_correlation",
0404 "MBD Charge Correlation ; Total Charge (South); Total Charge (North)",
0405 500, 0, 1500, 500, 0, 1500);
0406
0407 h2_GlobalQA_mbd_nhits_NS_correlation = new TH2F("h2_GlobalQA_mbd_nhits_NS_correlation",
0408 "MBD Number Of Hits Correlation ; Number Of Hits (South); " "Number Of Hits (North)",
0409 65, -0.5, 64.5, 65, -0.5, 64.5);
0410
0411 hm->registerHisto(h_GlobalQA_mbd_zvtx);
0412 hm->registerHisto(h_GlobalQA_mbd_zvtxq);
0413 hm->registerHisto(h_GlobalQA_mbd_zvtx_wide);
0414 hm->registerHisto(h_GlobalQA_calc_zvtx);
0415 hm->registerHisto(h_GlobalQA_calc_zvtx_wide);
0416 hm->registerHisto(h_GlobalQA_mbd_charge_s);
0417 hm->registerHisto(h_GlobalQA_mbd_charge_n);
0418 hm->registerHisto(h_GlobalQA_mbd_nhit_s);
0419 hm->registerHisto(h_GlobalQA_mbd_nhit_n);
0420 hm->registerHisto(h_GlobalQA_mbd_charge_sum);
0421 hm->registerHisto(h2_GlobalQA_mbd_charge_NS_correlation);
0422 hm->registerHisto(h2_GlobalQA_mbd_nhits_NS_correlation);
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437 h_GlobalQA_sEPD_adcsum_s = new TH1D(
0438 "h_GlobalQA_sEPD_adcsum_s", " ; sEPD ADC sum ; Counts", 200, 0, 2e4);
0439
0440 h_GlobalQA_sEPD_adcsum_n = new TH1D(
0441 "h_GlobalQA_sEPD_adcsum_n", " ; sEPD ADC sum ; Counts", 200, 0, 2e4);
0442
0443 h2_GlobalQA_sEPD_adcsum_ns =
0444 new TH2D("h2_GlobalQA_sEPD_adcsum_ns",
0445 "sEPD NS ADC sum correlation ; ADC sum (south); ADC sum (north)",
0446 500, 0, 2e4, 500, 0, 2e4);
0447
0448 h2Profile_GlobalQA_sEPD_tiles_north =
0449 new TProfile2D("h2Profile_GlobalQA_sEPD_tiles_north",
0450 "sEPD Tile Mean Energy (north); #eta; #phi", 16, -0.5,
0451 15.5, 24, -0.5, 23.5);
0452
0453 h2Profile_GlobalQA_sEPD_tiles_south =
0454 new TProfile2D("h2Profile_GlobalQA_sEPD_tiles_south",
0455 "sEPD Tile Mean Energy (south); #eta; #phi", 16, -0.5,
0456 15.5, 24, -0.5, 23.5);
0457
0458 h2_GlobalQA_sEPD_ADC_channel_south =
0459 new TH2D("h2_GlobalQA_sEPD_ADC_channel_south",
0460 "h2_GlobalQA_sEPD_ADC_channel_south ; #eta; #phi", 16, -0.5,
0461 15.5, 24, -0.5, 23.5);
0462 h2_GlobalQA_sEPD_ADC_channel_north =
0463 new TH2D("h2_GlobalQA_sEPD_ADC_channel_north",
0464 "h2_GlobalQA_sEPD_ADC_channel_north ; #eta; #phi", 16, -0.5,
0465 15.5, 24, -0.5, 23.5);
0466
0467 hm->registerHisto(h_GlobalQA_sEPD_adcsum_s);
0468 hm->registerHisto(h_GlobalQA_sEPD_adcsum_n);
0469 hm->registerHisto(h2_GlobalQA_sEPD_adcsum_ns);
0470 hm->registerHisto(h2Profile_GlobalQA_sEPD_tiles_south);
0471 hm->registerHisto(h2Profile_GlobalQA_sEPD_tiles_north);
0472 hm->registerHisto(h2_GlobalQA_sEPD_ADC_channel_south);
0473 hm->registerHisto(h2_GlobalQA_sEPD_ADC_channel_north);
0474
0475
0476 h_GlobalQA_zdc_zvtx = new TH1D(
0477 "h_GlobalQA_zdc_zvtx", ";Scaled Trigger 3: ZDC Coincidence zvtx [cm]",
0478 100, -1000, 1000);
0479 h_GlobalQA_zdc_zvtx_wide = new TH1D(
0480 "h_GlobalQA_zdc_zvtx_wide",
0481 ";Scaled Trigger 3: ZDC Coincidence zvtx [cm]", 100, -2000, 2000);
0482 h_GlobalQA_zdc_energy_s = new TH1D(
0483 "h_GlobalQA_zdc_energy_s",
0484 ";Scaled Trigger 3: ZDC Coincidence Energy [Gev]", 100, 10, 510);
0485 h_GlobalQA_zdc_energy_n = new TH1D(
0486 "h_GlobalQA_zdc_energy_n",
0487 ";Scaled Trigger 3: ZDC Coincidence Energy [Gev]", 100, 10, 510);
0488 hm->registerHisto(h_GlobalQA_zdc_zvtx);
0489 hm->registerHisto(h_GlobalQA_zdc_zvtx_wide);
0490 hm->registerHisto(h_GlobalQA_zdc_energy_s);
0491 hm->registerHisto(h_GlobalQA_zdc_energy_n);
0492
0493 h_GlobalQA_triggerVec = new TH1F("h_GlobalQA_triggerVec", "", 65, -0.5, 64.5);
0494 h_GlobalQA_triggerVec->SetDirectory(nullptr);
0495 hm->registerHisto(h_GlobalQA_triggerVec);
0496 }
