File indexing completed on 2026-07-16 08:12:43
0001 #include "sEPDTowerPairQA.h"
0002
0003 #include <calobase/TowerInfo.h>
0004 #include <calobase/TowerInfoContainer.h>
0005 #include <calobase/TowerInfoDefs.h>
0006
0007 #include <cdbobjects/CDBTTree.h>
0008
0009 #include <ffarawobjects/CaloPacket.h>
0010 #include <ffarawobjects/CaloPacketContainer.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phool/getClass.h>
0015
0016 #include <TFile.h>
0017 #include <TH1F.h>
0018 #include <TH2F.h>
0019 #include <TProfile.h>
0020 #include <TProfile2D.h>
0021 #include <TString.h>
0022 #include <TTree.h>
0023
0024 #include <algorithm>
0025 #include <cmath>
0026 #include <iostream>
0027
0028 sEPDTowerPairQA::sEPDTowerPairQA(const std::string &name)
0029 : SubsysReco(name)
0030 {
0031 }
0032
0033 bool sEPDTowerPairQA::sepd_skip_channel(int ich, int pid)
0034 {
0035 int sector = (ich + 1) / 32;
0036 int emptych = -999;
0037 if (sector == 0 && pid == 9001)
0038 {
0039 emptych = 1;
0040 }
0041 else
0042 {
0043 emptych = 14 + 32 * sector;
0044 }
0045 return (ich == emptych);
0046 }
0047
0048 void sEPDTowerPairQA::build_index_map()
0049 {
0050 m_idxpid.clear();
0051 m_idxch.clear();
0052 for (int pid = 9001; pid <= 9006; pid++)
0053 {
0054 for (int ch = 0; ch < 128; ch++)
0055 {
0056 if (sepd_skip_channel(ch, pid))
0057 {
0058 continue;
0059 }
0060 m_idxpid.push_back(pid);
0061 m_idxch.push_back(ch);
0062 }
0063 }
0064 }
0065
0066 TProfile *sEPDTowerPairQA::clone_template(const std::string &file, const char *newname)
0067 {
0068 if (file.empty())
0069 {
0070 return nullptr;
0071 }
0072 TFile *ft = TFile::Open(file.c_str(), "READ");
0073 if (!ft || ft->IsZombie())
0074 {
0075 std::cout << "sEPDTowerPairQA: WARNING could not open template file for cloning: " << file << std::endl;
0076 delete ft;
0077 m_file->cd();
0078 return nullptr;
0079 }
0080 TProfile *src = dynamic_cast<TProfile *>(ft->Get("waveform_template"));
0081 TProfile *clone = nullptr;
0082 if (src)
0083 {
0084 m_file->cd();
0085 clone = dynamic_cast<TProfile *>(src->Clone(newname));
0086 if (clone)
0087 {
0088 clone->SetDirectory(m_file);
0089 }
0090 }
0091 else
0092 {
0093 std::cout << "sEPDTowerPairQA: WARNING no waveform_template in " << file << std::endl;
0094 }
0095 ft->Close();
0096 delete ft;
0097 m_file->cd();
0098 return clone;
0099 }
0100
0101 int sEPDTowerPairQA::Init(PHCompositeNode * )
0102 {
0103 m_file = new TFile(m_outputFile.c_str(), "RECREATE");
0104
0105 const std::string tA = m_labelA;
0106 const std::string tB = m_labelB;
0107
0108 h_dchi2 = new TH1F("h_dchi2",
0109 Form("Paired #Delta#chi^{2} = #chi^{2}_{%s} - #chi^{2}_{%s};#Delta#chi^{2} [ADC^{2}];towers",
0110 tA.c_str(), tB.c_str()),
0111 1000, -5000, 5000);
0112 h_dchi2_frac = new TH1F("h_dchi2_frac",
0113 Form("Fractional improvement (#chi^{2}_{%s}-#chi^{2}_{%s})/#chi^{2}_{%s};fraction;towers",
0114 tA.c_str(), tB.c_str(), tA.c_str()),
0115 300, -2.0, 1.0);
0116 h_chi2_scatter = new TH2F("h_chi2_scatter",
0117 Form("#chi^{2} per tower;#chi^{2}_{%s} [ADC^{2}];#chi^{2}_{%s} [ADC^{2}]",
0118 tA.c_str(), tB.c_str()),
0119 300, 0, 3000, 300, 0, 3000);
0120 p_dchi2_vs_amp = new TProfile("p_dchi2_vs_amp",
0121 Form("<#Delta#chi^{2}> vs amplitude;amplitude [ADC];<#chi^{2}_{%s}-#chi^{2}_{%s}>",
0122 tA.c_str(), tB.c_str()),
0123 100, 0, 3000);
0124 p_dchi2frac_vs_amp = new TProfile("p_dchi2frac_vs_amp",
0125 "<fractional improvement> vs amplitude;amplitude [ADC];<fraction>",
0126 100, 0, 3000);
0127 p_dchi2frac_map = new TProfile2D("p_dchi2frac_map",
0128 "<fractional improvement> vs arm/ring;arm;ring;<fraction>",
0129 2, -0.5, 1.5, NRING, -0.5, NRING - 0.5);
0130
0131 h_dtime = new TH1F("h_dtime",
0132 Form("Paired #Delta t = t_{%s} - t_{%s};#Delta t [samples];towers", tA.c_str(), tB.c_str()),
0133 300, -3, 3);
0134 h_amp_scatter = new TH2F("h_amp_scatter",
0135 Form("amplitude per tower;amp_{%s} [ADC];amp_{%s} [ADC]", tA.c_str(), tB.c_str()),
0136 300, 0, 3000, 300, 0, 3000);
0137
0138 m_tree = new TTree("pair_tree", "sEPD paired template-fit comparison");
0139 m_tree->Branch("ich", &t_ich);
0140 m_tree->Branch("arm", &t_arm);
0141 m_tree->Branch("ring", &t_ring);
0142 m_tree->Branch("ampA", &t_ampA);
0143 m_tree->Branch("ampB", &t_ampB);
0144 m_tree->Branch("chi2A", &t_chi2A);
0145 m_tree->Branch("chi2B", &t_chi2B);
0146 m_tree->Branch("timeA", &t_timeA);
0147 m_tree->Branch("timeB", &t_timeB);
0148
0149 if (m_captureEvents > 0)
0150 {
0151 build_index_map();
0152 if (!m_channelMapFile.empty())
0153 {
0154 m_sepdmap = new CDBTTree(m_channelMapFile);
0155 }
0156 else
0157 {
0158 std::cout << "sEPDTowerPairQA: WARNING capture requested but no channel-map file set "
0159 << "- raw-sample capture will be skipped" << std::endl;
0160 }
0161 m_tmplA = clone_template(m_templateAFile, "template_A");
0162 m_tmplB = clone_template(m_templateBFile, "template_B");
0163
0164 m_file->cd();
0165 m_dtree = new TTree("display_tree", "sEPD raw waveform + both production template fits");
0166 m_dtree->Branch("evt", &d_evt);
0167 m_dtree->Branch("ich", &d_ich);
0168 m_dtree->Branch("arm", &d_arm);
0169 m_dtree->Branch("ring", &d_ring);
0170 m_dtree->Branch("ns", &d_ns);
0171 m_dtree->Branch("raw", d_raw, "raw[16]/F");
0172 m_dtree->Branch("ampA", &d_ampA);
0173 m_dtree->Branch("timeA", &d_timeA);
0174 m_dtree->Branch("chi2A", &d_chi2A);
0175 m_dtree->Branch("ampB", &d_ampB);
0176 m_dtree->Branch("timeB", &d_timeB);
0177 m_dtree->Branch("chi2B", &d_chi2B);
0178 std::cout << "sEPDTowerPairQA::Init - example-fit capture ENABLED for first " << m_captureEvents
0179 << " events (index map " << m_idxpid.size() << " channels)" << std::endl;
0180 }
0181
0182 std::cout << "sEPDTowerPairQA::Init - comparing node A=" << m_nodeA << " (" << m_labelA
0183 << ") vs node B=" << m_nodeB << " (" << m_labelB << "), amp cut " << m_amp_cut
0184 << ", output " << m_outputFile << std::endl;
0185
0186 return Fun4AllReturnCodes::EVENT_OK;
0187 }
0188
0189 int sEPDTowerPairQA::process_event(PHCompositeNode *topNode)
0190 {
0191 m_event++;
0192
0193 TowerInfoContainer *towersA = findNode::getClass<TowerInfoContainer>(topNode, m_nodeA);
0194 TowerInfoContainer *towersB = findNode::getClass<TowerInfoContainer>(topNode, m_nodeB);
0195 if (!towersA || !towersB)
0196 {
0197 if (m_event <= 3)
0198 {
0199 std::cout << "sEPDTowerPairQA::process_event - missing node "
0200 << (towersA ? m_nodeB : m_nodeA) << "; skipping event" << std::endl;
0201 }
0202 return Fun4AllReturnCodes::EVENT_OK;
0203 }
0204
0205 unsigned int nA = towersA->size();
0206 unsigned int nB = towersB->size();
0207 if (nA != nB && m_event <= 3)
0208 {
0209 std::cout << "sEPDTowerPairQA::process_event - WARNING size mismatch " << m_nodeA << "=" << nA
0210 << " " << m_nodeB << "=" << nB << std::endl;
0211 }
0212 unsigned int n = std::min(nA, nB);
0213
0214 const bool capturing = (m_dtree != nullptr && m_sepdmap != nullptr && m_event <= m_captureEvents);
0215 CaloPacket *pkts[6] = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
0216 if (capturing)
0217 {
0218 if (!m_mapValidated)
0219 {
0220 m_mapValidated = true;
0221 m_mapOK = (m_idxpid.size() == n);
0222 if (!m_mapOK)
0223 {
0224 std::cout << "sEPDTowerPairQA: WARNING index-map size " << m_idxpid.size()
0225 << " != container size " << n << " - example-fit capture disabled" << std::endl;
0226 }
0227 }
0228 if (m_mapOK)
0229 {
0230 CaloPacketContainer *cont = findNode::getClass<CaloPacketContainer>(topNode, "SEPDPackets");
0231 for (int pid = 9001; pid <= 9006; pid++)
0232 {
0233 pkts[pid - 9001] = cont ? cont->getPacketbyId(pid) : findNode::getClass<CaloPacket>(topNode, pid);
0234 }
0235 }
0236 }
0237
0238 for (unsigned int i = 0; i < n; i++)
0239 {
0240 TowerInfo *tA = towersA->get_tower_at_channel(i);
0241 TowerInfo *tB = towersB->get_tower_at_channel(i);
0242 if (!tA || !tB)
0243 {
0244 continue;
0245 }
0246
0247 float ampA = tA->get_energy();
0248 float ampB = tB->get_energy();
0249 float chi2A = tA->get_chi2();
0250 float chi2B = tB->get_chi2();
0251 float timeA = tA->get_time();
0252 float timeB = tB->get_time();
0253
0254 if (!std::isfinite(chi2A) || !std::isfinite(chi2B) || chi2A <= 0 || ampA < m_amp_cut)
0255 {
0256 continue;
0257 }
0258
0259 unsigned int key = towersA->encode_key(i);
0260 int arm = TowerInfoDefs::get_epd_arm(key);
0261 int ring = TowerInfoDefs::get_epd_rbin(key);
0262
0263 double dchi2 = chi2A - chi2B;
0264 double frac = dchi2 / chi2A;
0265
0266 h_dchi2->Fill(dchi2);
0267 h_dchi2_frac->Fill(frac);
0268 h_chi2_scatter->Fill(chi2A, chi2B);
0269 p_dchi2_vs_amp->Fill(ampA, dchi2);
0270 p_dchi2frac_vs_amp->Fill(ampA, frac);
0271 p_dchi2frac_map->Fill(arm, ring, frac);
0272 h_dtime->Fill(timeA - timeB);
0273 h_amp_scatter->Fill(ampA, ampB);
0274
0275 t_ich = (int) i;
0276 t_arm = arm;
0277 t_ring = ring;
0278 t_ampA = ampA;
0279 t_ampB = ampB;
0280 t_chi2A = chi2A;
0281 t_chi2B = chi2B;
0282 t_timeA = timeA;
0283 t_timeB = timeB;
0284 m_tree->Fill();
0285
0286 if (capturing && m_mapOK)
0287 {
0288 int idx = m_sepdmap->GetIntValue((int) i, m_mapFieldName);
0289 if (idx >= 0 && idx < (int) m_idxpid.size())
0290 {
0291 CaloPacket *cp = pkts[m_idxpid[idx] - 9001];
0292 int ch = m_idxch[idx];
0293 if (cp && !cp->iValue(ch, "SUPPRESSED"))
0294 {
0295 int ns = cp->iValue(0, "SAMPLES");
0296 if (ns > MAXSAMP)
0297 {
0298 ns = MAXSAMP;
0299 }
0300 for (int s = 0; s < MAXSAMP; s++)
0301 {
0302 d_raw[s] = 0.f;
0303 }
0304 for (int s = 0; s < ns; s++)
0305 {
0306 d_raw[s] = (float) cp->iValue(s, ch);
0307 }
0308 d_evt = m_event;
0309 d_ich = (int) i;
0310 d_arm = arm;
0311 d_ring = ring;
0312 d_ns = ns;
0313 d_ampA = ampA;
0314 d_timeA = timeA;
0315 d_chi2A = chi2A;
0316 d_ampB = ampB;
0317 d_timeB = timeB;
0318 d_chi2B = chi2B;
0319 m_dtree->Fill();
0320 }
0321 }
0322 }
0323
0324 m_sum_dchi2 += dchi2;
0325 m_sum_dchi2_frac += frac;
0326 m_npaired++;
0327 if (chi2B < chi2A)
0328 {
0329 m_nimproved++;
0330 }
0331 }
0332
0333 return Fun4AllReturnCodes::EVENT_OK;
0334 }
0335
0336 int sEPDTowerPairQA::End(PHCompositeNode * )
0337 {
0338 if (m_npaired > 0)
0339 {
0340 double mean_dchi2 = m_sum_dchi2 / (double) m_npaired;
0341 double mean_frac = m_sum_dchi2_frac / (double) m_npaired;
0342 double frac_improved = (double) m_nimproved / (double) m_npaired;
0343 std::cout << "sEPDTowerPairQA::End - paired towers (" << m_labelA << " vs " << m_labelB << ", amp>"
0344 << m_amp_cut << "): " << m_npaired << std::endl;
0345 std::cout << " mean Delta(chi2) = " << mean_dchi2 << " ADC^2" << std::endl;
0346 std::cout << " mean fractional improvement = " << mean_frac << std::endl;
0347 std::cout << " fraction of towers improved = " << frac_improved << std::endl;
0348 }
0349 else
0350 {
0351 std::cout << "sEPDTowerPairQA::End - WARNING no paired towers above amp cut " << m_amp_cut << std::endl;
0352 }
0353
0354 m_file->cd();
0355 h_dchi2->Write();
0356 h_dchi2_frac->Write();
0357 h_chi2_scatter->Write();
0358 p_dchi2_vs_amp->Write();
0359 p_dchi2frac_vs_amp->Write();
0360 p_dchi2frac_map->Write();
0361 h_dtime->Write();
0362 h_amp_scatter->Write();
0363 m_tree->Write();
0364
0365 if (m_dtree)
0366 {
0367 if (m_tmplA)
0368 {
0369 m_tmplA->Write();
0370 }
0371 if (m_tmplB)
0372 {
0373 m_tmplB->Write();
0374 }
0375 m_dtree->Write();
0376 std::cout << "sEPDTowerPairQA::End - captured " << m_dtree->GetEntries()
0377 << " example waveforms (display_tree) + template_A/template_B" << std::endl;
0378 }
0379
0380 m_file->Close();
0381 delete m_file;
0382 m_file = nullptr;
0383
0384 delete m_sepdmap;
0385 m_sepdmap = nullptr;
0386
0387 return Fun4AllReturnCodes::EVENT_OK;
0388 }
