File indexing completed on 2026-07-16 08:12:43
0001 #include "sEPDWaveform.h"
0002
0003 #include <calobase/TowerInfoContainer.h>
0004 #include <calobase/TowerInfo.h>
0005 #include <calobase/TowerInfoDefs.h>
0006
0007 #include <Event/Event.h>
0008 #include <Event/packet.h>
0009
0010 #include <ffarawobjects/CaloPacket.h>
0011 #include <ffarawobjects/CaloPacketContainer.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>
0015 #include <fun4all/Fun4AllServer.h>
0016 #include <fun4all/Fun4AllHistoManager.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/PHObject.h>
0023 #include <phool/getClass.h>
0024
0025 #include <vector>
0026 #include <map>
0027 #include <iomanip>
0028
0029 #include <TFile.h>
0030 #include <TGraph.h>
0031 #include <TProfile2D.h>
0032
0033 using namespace std;
0034
0035 sEPDWaveform::sEPDWaveform(const std::string &name) : SubsysReco(name)
0036 {
0037 }
0038
0039 sEPDWaveform::~sEPDWaveform()
0040 {
0041 }
0042
0043 int sEPDWaveform::Init(PHCompositeNode *topNode)
0044 {
0045 WaveformProcessing = nullptr;
0046
0047 WaveformProcessing2 = new CaloWaveformProcessing();
0048 WaveformProcessing2->set_processing_type(CaloWaveformProcessing::FAST);
0049 WaveformProcessing2->initialize_processing();
0050
0051 m_file = new TFile(m_outputFileName.c_str(), "RECREATE");
0052
0053 if (m_fill_geometry_histos)
0054 {
0055
0056 for (int arm = 0; arm < NARM; arm++)
0057 {
0058 for (int ring = 0; ring < NRING; ring++)
0059 {
0060 hpwaveform_1[arm][ring] = new TProfile(Form("hpwaveform_1_arm%d_ring%d", arm, ring),
0061 Form("Low amplitude waveform arm%d ring%d", arm, ring),
0062 300, 0, 30);
0063 hpwaveform_2[arm][ring] = new TProfile(Form("hpwaveform_2_arm%d_ring%d", arm, ring),
0064 Form("Medium amplitude waveform arm%d ring%d", arm, ring),
0065 300, 0, 30);
0066 hpwaveform_3[arm][ring] = new TProfile(Form("hpwaveform_3_arm%d_ring%d", arm, ring),
0067 Form("High amplitude waveform arm%d ring%d", arm, ring),
0068 300, 0, 30);
0069
0070 h2waveform_1[arm][ring] = new TH2F(Form("h2waveform_1_arm%d_ring%d", arm, ring),
0071 Form("Low amplitude waveform 2D arm%d ring%d;Sample;ADC", arm, ring),
0072 m_n_samples, 0, m_n_samples, 200, -100, 1500);
0073 h2waveform_2[arm][ring] = new TH2F(Form("h2waveform_2_arm%d_ring%d", arm, ring),
0074 Form("Medium amplitude waveform 2D arm%d ring%d;Sample;ADC", arm, ring),
0075 m_n_samples, 0, m_n_samples, 200, -100, 1500);
0076 h2waveform_3[arm][ring] = new TH2F(Form("h2waveform_3_arm%d_ring%d", arm, ring),
0077 Form("High amplitude waveform 2D arm%d ring%d;Sample;ADC", arm, ring),
0078 m_n_samples, 0, m_n_samples, 200, -100, 1500);
0079 }
0080 }
0081
0082 hpwaveform_1_sum = new TProfile("hpwaveform_1_sum", "Low amplitude sum", 300, 0, 30);
0083 hpwaveform_2_sum = new TProfile("hpwaveform_2_sum", "Medium amplitude sum", 300, 0, 30);
0084 hpwaveform_3_sum = new TProfile("hpwaveform_3_sum", "High amplitude sum", 300, 0, 30);
0085
0086 h2waveform_1_sum = new TH2F("h2waveform_1_sum", "Low amplitude sum 2D;Sample;ADC", m_n_samples, 0, m_n_samples, 200, -100, 1500);
0087 h2waveform_2_sum = new TH2F("h2waveform_2_sum", "Medium amplitude sum 2D;Sample;ADC", m_n_samples, 0, m_n_samples, 200, -100, 1500);
0088 h2waveform_3_sum = new TH2F("h2waveform_3_sum", "High amplitude sum 2D;Sample;ADC", m_n_samples, 0, m_n_samples, 200, -100, 1500);
0089 }
0090
0091 if (m_fill_global_histos)
0092 {
0093 hpwaveform_global_1 = new TProfile("hpwaveform_global_1", "Global Low Amplitude Template", 300, 0, 30);
0094 hpwaveform_global_2 = new TProfile("hpwaveform_global_2", "Global Medium Amplitude Template", 300, 0, 30);
0095 hpwaveform_global_3 = new TProfile("hpwaveform_global_3", "Global High Amplitude Template", 300, 0, 30);
0096
0097 h2waveform_global_1 = new TH2F("h2waveform_global_1", "Global Low Amplitude 2D;Sample;ADC", m_n_samples, 0, m_n_samples, 200, -100, 1500);
0098 h2waveform_global_2 = new TH2F("h2waveform_global_2", "Global Medium Amplitude 2D;Sample;ADC", m_n_samples, 0, m_n_samples, 200, -100, 1500);
0099 h2waveform_global_3 = new TH2F("h2waveform_global_3", "Global High Amplitude 2D;Sample;ADC", m_n_samples, 0, m_n_samples, 200, -100, 1500);
0100 }
0101
0102 waveform_template = new TProfile("waveform_template", "sEPD Waveform Template", m_template_nbins, 0, 30);
0103
0104 h_sig_chi2 = new TH2F("h_sig_chi2","Signal vs Chi2",5500,-100,1e4,5000,0,5e6);
0105 h_sig_chi2_nt = new TH2F("h_sig_chi2_nt","Signal vs Chi2 (no template)",1100,-100,1e4,1000,0,1e5);
0106 h_sig1_sig2 = new TH2F("h_sig1_sig2","Signal1 vs Signal2",1100,-100,1e3,1100,-100,1e3);
0107 h_sigP1_sigP2 = new TH2F("h_sigP1_sigP2","Peak1 vs Peak2",320,-1,31,320,-1,31);
0108 h_chi1_chi2 = new TH2F("h_chi1_chi2","Chi2_1 vs Chi2_2",1e3,0,1e4,1e3,0,1e4);
0109 h_ped1_ped2 = new TH2F("h_ped1_ped2","Pedestal1 vs Pedestal2",1e3,0,1e4,1e3,0,1e4);
0110 h_dif_peak_time = new TH2F("h_dif_peak_time","Peak time difference",100,-5,5,100,-5,5);
0111 h_chi2 = new TH1F("h_chi2","Chi2 distribution",1000,0,1e5);
0112
0113 p_chi2 = new TProfile2D("p_chi2","Chi2 vs arm/ring",NARM,0,NARM,NRING,0,NRING);
0114 h_etaphi_ped = new TProfile2D("h_etaphi_ped","Pedestal vs arm/ring",NARM,0,NARM,NRING,0,NRING);
0115 h_etaphi_ped2 = new TProfile2D("h_etaphi_ped2","Pedestal2 vs arm/ring",NARM,0,NARM,NRING,0,NRING);
0116
0117 h_time1 = new TH1F("h_time1","Peak time distribution 1",100,0,30);
0118 h_time2 = new TH1F("h_time2","Peak time distribution 2",100,0,30);
0119
0120 m_tree = new TTree("waveform_tree", "sEPD Waveform Tree");
0121 m_tree->Branch("arm", &m_arm);
0122 m_tree->Branch("ring", &m_ring);
0123 m_tree->Branch("phibin", &m_phibin);
0124 m_tree->Branch("peak", &m_peak);
0125 m_tree->Branch("peaktime", &m_peaktime);
0126 m_tree->Branch("pedestal", &m_pedestal);
0127 m_tree->Branch("ADC_peak", &m_ADC_peak);
0128 m_tree->Branch("ADC_peaktime", &m_ADC_peaktime);
0129 m_tree->Branch("waveform", &m_waveform);
0130
0131 return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133
0134 int sEPDWaveform::process_event(PHCompositeNode *topNode)
0135 {
0136
0137 m_arm.clear();
0138 m_ring.clear();
0139 m_phibin.clear();
0140 m_peak.clear();
0141 m_peaktime.clear();
0142 m_pedestal.clear();
0143 m_ADC_peak.clear();
0144 m_ADC_peaktime.clear();
0145 m_waveform.clear();
0146
0147 auto processPacket = [&](auto *packet, int pid)
0148 {
0149 if (!packet)
0150 {
0151 return;
0152 }
0153
0154 int nchannels = packet->iValue(0, "CHANNELS");
0155 if (nchannels > m_nchannels)
0156 {
0157 std::cout << "Skipping packet " << pid << " - nchannels " << nchannels
0158 << " exceeds expected " << m_nchannels << std::endl;
0159 return;
0160 }
0161
0162 for (int channel = 0; channel < nchannels; channel++)
0163 {
0164 if (packet->iValue(channel, "SUPPRESSED"))
0165 {
0166 continue;
0167 }
0168
0169 std::vector<double> waveform;
0170 for (int samp = 0; samp < m_n_samples; samp++)
0171 {
0172 waveform.push_back(packet->iValue(samp, channel));
0173 }
0174
0175 if (waveform.size() != (unsigned int)m_n_samples)
0176 {
0177 continue;
0178 }
0179
0180 int global_channel = (pid - m_packet_low) * m_nchannels + channel;
0181
0182 unsigned int tower_key = TowerInfoDefs::encode_epd(global_channel);
0183
0184 int arm = TowerInfoDefs::get_epd_arm(tower_key);
0185 int ring = TowerInfoDefs::get_epd_rbin(tower_key);
0186 int phi = TowerInfoDefs::get_epd_phibin(tower_key);
0187
0188 m_arm.push_back(arm);
0189 m_ring.push_back(ring);
0190 m_phibin.push_back(phi);
0191 m_waveform.push_back(waveform);
0192
0193 double peak = 0;
0194 double peaktime = 0;
0195 double pedestal = 0;
0196
0197 int peakbin = 0;
0198 for (int i = 0; i < m_n_samples; i++)
0199 {
0200 if (waveform[i] > peak)
0201 {
0202 peak = waveform[i];
0203 peakbin = i;
0204 }
0205 }
0206 peaktime = peakbin;
0207
0208 if (peakbin > 0 && peakbin < m_n_samples - 1)
0209 {
0210 double ym1 = waveform[peakbin - 1];
0211 double y0 = waveform[peakbin];
0212 double yp1 = waveform[peakbin + 1];
0213 double denom = ym1 - (2.0 * y0) + yp1;
0214 if (denom < 0)
0215 {
0216 double delta = 0.5 * (ym1 - yp1) / denom;
0217 if (delta > -0.5 && delta < 0.5)
0218 {
0219 peaktime = peakbin + delta;
0220 }
0221 }
0222 }
0223
0224 for (int i = 0; i < 3 && i < m_n_samples; i++)
0225 {
0226 pedestal += waveform[i];
0227 }
0228 pedestal /= 3.0;
0229
0230 m_peak.push_back(peak);
0231 m_peaktime.push_back(peaktime);
0232 m_pedestal.push_back(pedestal);
0233 m_ADC_peak.push_back(peak - pedestal);
0234 m_ADC_peaktime.push_back(peaktime);
0235
0236 if (arm < NARM && ring < NRING)
0237 {
0238
0239 int max_phi = (ring == 0) ? NPHI_RING0 : NPHI;
0240 if (phi >= max_phi)
0241 {
0242 std::cout << "Warning: Invalid phi bin " << phi << " for ring " << ring << std::endl;
0243 continue;
0244 }
0245
0246 double amplitude = peak - pedestal;
0247
0248 if (amplitude > 50)
0249 {
0250
0251 double thr = pedestal + 0.5 * amplitude;
0252 double t_hm = -1.0;
0253 for (int i = 1; i <= peakbin; i++)
0254 {
0255 if (waveform[i - 1] < thr && waveform[i] >= thr && waveform[i] != waveform[i - 1])
0256 {
0257 t_hm = (i - 1) + (thr - waveform[i - 1]) / (waveform[i] - waveform[i - 1]);
0258 break;
0259 }
0260 }
0261
0262 double area = 0.0;
0263 for (int i = 0; i < m_n_samples; i++)
0264 {
0265 double y = waveform[i] - pedestal;
0266 if (y > 0) { area += y; }
0267 }
0268 if (t_hm >= 0 && area > 0)
0269 {
0270 double shift = m_target_hm - t_hm;
0271 for (int i = 0; i < m_n_samples; i++)
0272 {
0273 double ttmpl = i + shift;
0274 if (ttmpl > 0 && ttmpl < 30)
0275 {
0276 waveform_template->Fill(ttmpl, (waveform[i] - pedestal) / area);
0277 }
0278 }
0279 }
0280 }
0281
0282 for (int i = 0; i < m_n_samples; i++)
0283 {
0284 double time = i - peaktime + 6.0;
0285
0286 if (amplitude > 50 && amplitude < 500)
0287 {
0288
0289 if (m_fill_geometry_histos)
0290 {
0291 h2waveform_1[arm][ring]->Fill(i, waveform[i] - pedestal);
0292 h2waveform_1_sum->Fill(i, waveform[i] - pedestal);
0293
0294 if (time > 0 && time < 30)
0295 {
0296 hpwaveform_1[arm][ring]->Fill(time, waveform[i] - pedestal);
0297 hpwaveform_1_sum->Fill(time, waveform[i] - pedestal);
0298 }
0299 }
0300
0301 if (m_fill_global_histos)
0302 {
0303 h2waveform_global_1->Fill(i, waveform[i] - pedestal);
0304 if (time > 0 && time < 30)
0305 {
0306 hpwaveform_global_1->Fill(time, waveform[i] - pedestal);
0307 }
0308 }
0309 }
0310 else if (amplitude >= 500 && amplitude < 1000)
0311 {
0312
0313 if (m_fill_geometry_histos)
0314 {
0315 h2waveform_2[arm][ring]->Fill(i, waveform[i] - pedestal);
0316 h2waveform_2_sum->Fill(i, waveform[i] - pedestal);
0317
0318 if (time > 0 && time < 30)
0319 {
0320 hpwaveform_2[arm][ring]->Fill(time, waveform[i] - pedestal);
0321 hpwaveform_2_sum->Fill(time, waveform[i] - pedestal);
0322 }
0323 }
0324
0325 if (m_fill_global_histos)
0326 {
0327 h2waveform_global_2->Fill(i, waveform[i] - pedestal);
0328 if (time > 0 && time < 30)
0329 {
0330 hpwaveform_global_2->Fill(time, waveform[i] - pedestal);
0331 }
0332 }
0333 }
0334 else if (amplitude >= 1000)
0335 {
0336
0337 if (m_fill_geometry_histos)
0338 {
0339 h2waveform_3[arm][ring]->Fill(i, waveform[i] - pedestal);
0340 h2waveform_3_sum->Fill(i, waveform[i] - pedestal);
0341
0342 if (time > 0 && time < 30)
0343 {
0344 hpwaveform_3[arm][ring]->Fill(time, waveform[i] - pedestal);
0345 hpwaveform_3_sum->Fill(time, waveform[i] - pedestal);
0346 }
0347 }
0348
0349 if (m_fill_global_histos)
0350 {
0351 h2waveform_global_3->Fill(i, waveform[i] - pedestal);
0352 if (time > 0 && time < 30)
0353 {
0354 hpwaveform_global_3->Fill(time, waveform[i] - pedestal);
0355 }
0356 }
0357 }
0358 }
0359 }
0360 }
0361
0362 };
0363
0364 if (m_offline)
0365 {
0366
0367 CaloPacketContainer *cont = findNode::getClass<CaloPacketContainer>(topNode, "SEPDPackets");
0368 bool found = false;
0369 for (int pid = m_packet_low; pid <= m_packet_high; pid++)
0370 {
0371 CaloPacket *cp = cont ? cont->getPacketbyId(pid) : findNode::getClass<CaloPacket>(topNode, pid);
0372 if (cp)
0373 {
0374 found = true;
0375
0376 if (eventCount == 0)
0377 {
0378 int nch = cp->iValue(0, "CHANNELS");
0379 int nunsupp = 0;
0380 for (int c = 0; c < nch; c++)
0381 {
0382 if (!cp->iValue(c, "SUPPRESSED")) { nunsupp++; }
0383 }
0384 std::cout << "[sEPDWaveform] packet " << pid
0385 << ": SAMPLES=" << cp->iValue(0, "SAMPLES")
0386 << " CHANNELS=" << nch
0387 << " unsuppressed=" << nunsupp << std::endl;
0388 }
0389 }
0390 processPacket(cp, pid);
0391 }
0392 if (!found)
0393 {
0394 if (eventCount < 3)
0395 {
0396 std::cout << "sEPDWaveform::process_event - no sEPD CaloPacket nodes (9001-9006) found" << std::endl;
0397 }
0398 return Fun4AllReturnCodes::DISCARDEVENT;
0399 }
0400 }
0401 else
0402 {
0403
0404 Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0405 if (!_event)
0406 {
0407 std::cout << "sEPDWaveform::process_event - Event not found" << std::endl;
0408 return Fun4AllReturnCodes::DISCARDEVENT;
0409 }
0410 for (int pid = m_packet_low; pid <= m_packet_high; pid++)
0411 {
0412 Packet *packet = _event->getPacket(pid);
0413 processPacket(packet, pid);
0414 delete packet;
0415 }
0416 }
0417
0418 if (m_arm.size() > 0)
0419 {
0420 m_tree->Fill();
0421 }
0422
0423 eventCount++;
0424
0425 if (m_show_progress && eventCount % m_progress_interval == 0)
0426 {
0427
0428 std::cout << "\rProcessed " << eventCount << " events...";
0429 std::cout.flush();
0430 }
0431
0432 return Fun4AllReturnCodes::EVENT_OK;
0433 }
0434
0435 int sEPDWaveform::ResetEvent(PHCompositeNode *topNode)
0436 {
0437 return Fun4AllReturnCodes::EVENT_OK;
0438 }
0439
0440 int sEPDWaveform::End(PHCompositeNode *topNode)
0441 {
0442 if (m_show_progress)
0443 {
0444 std::cout << std::endl;
0445 }
0446
0447 std::cout << "sEPDWaveform::End - Writing output to " << m_outputFileName << std::endl;
0448
0449 m_file->cd();
0450
0451 for (int arm = 0; arm < NARM; arm++)
0452 {
0453 for (int ring = 0; ring < NRING; ring++)
0454 {
0455 if (hpwaveform_1[arm][ring]) hpwaveform_1[arm][ring]->Write();
0456 if (hpwaveform_2[arm][ring]) hpwaveform_2[arm][ring]->Write();
0457 if (hpwaveform_3[arm][ring]) hpwaveform_3[arm][ring]->Write();
0458 if (h2waveform_1[arm][ring]) h2waveform_1[arm][ring]->Write();
0459 if (h2waveform_2[arm][ring]) h2waveform_2[arm][ring]->Write();
0460 if (h2waveform_3[arm][ring]) h2waveform_3[arm][ring]->Write();
0461 }
0462 }
0463
0464 if (hpwaveform_1_sum) hpwaveform_1_sum->Write();
0465 if (hpwaveform_2_sum) hpwaveform_2_sum->Write();
0466 if (hpwaveform_3_sum) hpwaveform_3_sum->Write();
0467 if (h2waveform_1_sum) h2waveform_1_sum->Write();
0468 if (h2waveform_2_sum) h2waveform_2_sum->Write();
0469 if (h2waveform_3_sum) h2waveform_3_sum->Write();
0470
0471 if (hpwaveform_global_1) hpwaveform_global_1->Write();
0472 if (hpwaveform_global_2) hpwaveform_global_2->Write();
0473 if (hpwaveform_global_3) hpwaveform_global_3->Write();
0474 if (h2waveform_global_1) h2waveform_global_1->Write();
0475 if (h2waveform_global_2) h2waveform_global_2->Write();
0476 if (h2waveform_global_3) h2waveform_global_3->Write();
0477
0478 if (waveform_template && waveform_template->GetMaximum() > 0)
0479 {
0480 std::cout << "sEPDWaveform::End - Normalizing template from max=" << waveform_template->GetMaximum() << " to 1.0" << std::endl;
0481 waveform_template->Scale(1.0 / waveform_template->GetMaximum());
0482 waveform_template->Write();
0483 }
0484 else if (waveform_template)
0485 {
0486 std::cout << "Warning: waveform_template has zero maximum, writing unnormalized" << std::endl;
0487 waveform_template->Write();
0488 }
0489
0490 h_sig_chi2->Write();
0491 h_sig_chi2_nt->Write();
0492 h_sig1_sig2->Write();
0493 h_sigP1_sigP2->Write();
0494 h_chi1_chi2->Write();
0495 h_ped1_ped2->Write();
0496 h_dif_peak_time->Write();
0497 h_chi2->Write();
0498 p_chi2->Write();
0499 h_etaphi_ped->Write();
0500 h_etaphi_ped2->Write();
0501 h_time1->Write();
0502 h_time2->Write();
0503
0504 m_tree->Write();
0505
0506 m_file->Close();
0507 delete m_file;
0508
0509 if (WaveformProcessing) delete WaveformProcessing;
0510 if (WaveformProcessing2) delete WaveformProcessing2;
0511
0512 std::cout << "sEPDWaveform::End - Processed " << eventCount << " events" << std::endl;
0513
0514 return Fun4AllReturnCodes::EVENT_OK;
0515 }
0516
0517 void sEPDWaveform::Print(const std::string &what) const
0518 {
0519 }
0520
0521 void sEPDWaveform::ShowProgress(int current, int total)
0522 {
0523 if (total <= 0) return;
0524
0525 float percent = (float)current / (float)total * 100.0;
0526 int barWidth = 50;
0527 int pos = barWidth * current / total;
0528
0529 std::cout << "\r[";
0530 for (int i = 0; i < barWidth; ++i)
0531 {
0532 if (i < pos) std::cout << "=";
0533 else if (i == pos) std::cout << ">";
0534 else std::cout << " ";
0535 }
0536 std::cout << "] " << std::fixed << std::setprecision(1) << percent << "% ("
0537 << current << "/" << total << " events)";
0538 std::cout.flush();
0539
0540 if (current == total)
0541 {
0542 std::cout << std::endl;
0543 }
0544 }
