File indexing completed on 2025-08-05 08:20:23
0001
0002
0003
0004
0005
0006 #include "CemcMon.h"
0007
0008 #include <onlmon/OnlMon.h> // for OnlMon
0009 #include <onlmon/OnlMonServer.h>
0010 #include <onlmon/pseudoRunningMean.h>
0011 #include <onlmon/triggerEnum.h>
0012
0013 #include <calobase/TowerInfoDefs.h>
0014 #include <caloreco/CaloWaveformFitting.h>
0015
0016 #include <cdbobjects/CDBTTree.h>
0017
0018 #include <Event/Event.h>
0019 #include <Event/EventTypes.h>
0020 #include <Event/eventReceiverClient.h>
0021 #include <Event/msg_profile.h>
0022
0023 #include <TH1.h>
0024 #include <TH2.h>
0025 #include <TProfile.h>
0026 #include <TProfile2D.h>
0027
0028 #include <cmath>
0029 #include <cstdio> // for printf
0030 #include <fstream>
0031 #include <iostream>
0032 #include <sstream>
0033 #include <string> // for allocator, string, char_traits
0034
0035 enum
0036 {
0037 TRGMESSAGE = 1,
0038 FILLMESSAGE = 2
0039 };
0040
0041 const int depth = 50000;
0042
0043
0044 const float hit_threshold = 100;
0045 const float waveform_hit_threshold = 100;
0046 const float chi2_check_threshold = 3000;
0047
0048 using namespace std;
0049
0050 CemcMon::CemcMon(const std::string &name)
0051 : OnlMon(name)
0052 , eventCounter(0)
0053 {
0054
0055
0056
0057 return;
0058 }
0059
0060 CemcMon::~CemcMon()
0061 {
0062 for (auto iter : rm_vector_twr)
0063 {
0064 delete iter;
0065 }
0066 for (auto iter : rm_vector_twrhits)
0067 {
0068 delete iter;
0069 }
0070 for (auto iter : rm_vector_twrhits_alltrig)
0071 {
0072 delete iter;
0073 }
0074
0075 delete WaveformProcessingFast;
0076 delete WaveformProcessingTemp;
0077 delete erc;
0078 return;
0079 }
0080
0081 int CemcMon::Init()
0082 {
0083
0084 const char *cemccalib = getenv("CEMCCALIB");
0085 if (!cemccalib)
0086 {
0087 std::cout << "CEMCCALIB environment variable not set" << std::endl;
0088 exit(1);
0089 }
0090 std::string fullfile = std::string(cemccalib) + "/" + "CemcMonData.dat";
0091 std::ifstream calib(fullfile);
0092 calib.close();
0093
0094
0095 printf("CemcMon::Init()\n");
0096 std::string cdbfile = std::string(cemccalib) + "/cemc_adcskipmask_40430.root";
0097 cdbttree = new CDBTTree(cdbfile.c_str());
0098
0099
0100 for (int itrig = 0; itrig < 64; itrig++)
0101 {
0102 h2_cemc_hits_trig[itrig] = new TH2F(Form("h2_cemc_hits_trig_bit_%d", itrig), "", 96, 0, 96, 256, 0, 256);
0103 }
0104 p2_zsFrac_etaphi = new TProfile2D("p2_zsFrac_etaphi", "", 96, 0, 96, 256, 0, 256);
0105 p2_zsFrac_etaphi_all = new TProfile2D("p2_zsFrac_etaphi_all", "", 96, 0, 96, 256, 0, 256);
0106 h1_cemc_trig = new TH1F("h1_cemc_trig", "", 64, -0.5, 63.5);
0107 h1_packet_event = new TH1F("h1_packet_event", "", 8, packetlow - 0.5, packethigh + 0.5);
0108 h2_caloPack_gl1_clock_diff = new TH2F("h2_caloPack_gl1_clock_diff", "", 8, packetlow - 0.5, packethigh + 0.5, 65536, 0, 65536);
0109 h_evtRec = new TProfile("h_evtRec", "", 1, 0, 1);
0110
0111
0112 h2_cemc_hits = new TH2F("h2_cemc_hits", "", 96, 0, 96, 256, 0, 256);
0113 h2_cemc_rm = new TH2F("h2_cemc_rm", "", 96, 0, 96, 256, 0, 256);
0114 h2_cemc_rmhits = new TH2F("h2_cemc_rmhits", "", 96, 0, 96, 256, 0, 256);
0115 h2_cemc_rmhits_alltrig = new TH2F("h2_cemc_rmhits_alltrig", "", 96, 0, 96, 256, 0, 256);
0116 h2_cemc_mean = new TH2F("h2_cemc_mean", "", 96, 0, 96, 256, 0, 256);
0117 h1_cemc_adc = new TH1F("h1_cemc_adc", "", 1000, 0.5, pow(2, 14));
0118
0119 h1_event = new TH1F("h1_event", "", 1, 0, 1);
0120
0121
0122
0123 h2_waveform_twrAvg = new TH2F("h2_waveform_twrAvg", "", 12, -0.5, 11.5, 1000, 0, 15000);
0124 h1_waveform_time = new TH1F("h1_waveform_time", "", 12, -0.5, 11.5);
0125 h1_waveform_pedestal = new TH1F("h1_waveform_pedestal", "", 5000, 1, 5000);
0126
0127
0128 h1_cemc_fitting_sigDiff = new TH1F("h1_fitting_sigDiff", "", 50, 0, 2);
0129 h1_cemc_fitting_pedDiff = new TH1F("h1_fitting_pedDiff", "", 50, 0, 2);
0130 h1_cemc_fitting_timeDiff = new TH1F("h1_fitting_timeDiff", "", 50, -10, 10);
0131
0132
0133 h1_packet_number = new TH1F("h1_packet_number", "", 128, 6000.5, 6128.5);
0134 h1_packet_length = new TH1F("h1_packet_length", "", 128, 6000.5, 6128.5);
0135 h1_packet_chans = new TH1F("h1_packet_chans", "", 128, 6000.5, 6128.5);
0136
0137 p2_bad_chi2 = new TProfile2D("p2_bad_chi2", "", 96, 0, 96, 256, 0, 256);
0138
0139 p2_pre_post = new TProfile2D("p2_pre_post", "", 96, 0, 96, 256, 0, 256, "S");
0140
0141
0142
0143 for (int i = 0; i < Ntower; i++)
0144 {
0145 rm_vector_twr.push_back(new pseudoRunningMean(1, depth));
0146 rm_vector_twrhits.push_back(new pseudoRunningMean(1, depth));
0147 rm_vector_twrhits_alltrig.push_back(new pseudoRunningMean(1, depth));
0148 }
0149
0150
0151
0152 OnlMonServer *se = OnlMonServer::instance();
0153
0154
0155
0156 for (int itrig = 0; itrig < 64; itrig++)
0157 {
0158 se->registerHisto(this, h2_cemc_hits_trig[itrig]);
0159 }
0160 se->registerHisto(this, p2_zsFrac_etaphi);
0161 se->registerHisto(this, p2_zsFrac_etaphi_all);
0162 se->registerHisto(this, p2_bad_chi2);
0163 se->registerHisto(this, h1_cemc_trig);
0164 se->registerHisto(this, h1_packet_event);
0165 se->registerHisto(this, h2_caloPack_gl1_clock_diff);
0166 se->registerHisto(this, h_evtRec);
0167
0168 se->registerHisto(this, h2_cemc_hits);
0169 se->registerHisto(this, h2_cemc_rm);
0170 se->registerHisto(this, h2_cemc_rmhits);
0171 se->registerHisto(this, h2_cemc_rmhits_alltrig);
0172 se->registerHisto(this, p2_pre_post);
0173 se->registerHisto(this, h2_cemc_mean);
0174 se->registerHisto(this, h1_event);
0175
0176 se->registerHisto(this, h2_waveform_twrAvg);
0177 se->registerHisto(this, h1_waveform_time);
0178 se->registerHisto(this, h1_waveform_pedestal);
0179 se->registerHisto(this, h1_cemc_fitting_sigDiff);
0180 se->registerHisto(this, h1_cemc_fitting_pedDiff);
0181 se->registerHisto(this, h1_cemc_fitting_timeDiff);
0182 se->registerHisto(this, h1_packet_number);
0183 se->registerHisto(this, h1_packet_length);
0184 se->registerHisto(this, h1_packet_chans);
0185 se->registerHisto(this, h1_cemc_adc);
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201 WaveformProcessingFast = new CaloWaveformFitting();
0202
0203 WaveformProcessingTemp = new CaloWaveformFitting();
0204
0205 std::string cemctemplate;
0206 if (getenv("CEMCCALIB"))
0207 {
0208 cemctemplate = getenv("CEMCCALIB");
0209 }
0210 else
0211 {
0212 cemctemplate = ".";
0213 }
0214 cemctemplate += std::string("/testbeam_cemc_template.root");
0215 WaveformProcessingTemp->initialize_processing(cemctemplate);
0216
0217 if (anaGL1)
0218 {
0219 erc = new eventReceiverClient("gl1daq");
0220 }
0221
0222 return 0;
0223 }
0224
0225 int CemcMon::BeginRun(const int )
0226 {
0227
0228
0229
0230
0231 std::vector<runningMean *>::iterator rm_it;
0232 for (rm_it = rm_vector_twr.begin(); rm_it != rm_vector_twr.end(); ++rm_it)
0233 {
0234 (*rm_it)->Reset();
0235 }
0236 for (rm_it = rm_vector_twrhits.begin(); rm_it != rm_vector_twrhits.end(); ++rm_it)
0237 {
0238 (*rm_it)->Reset();
0239 }
0240 for (rm_it = rm_vector_twrhits_alltrig.begin(); rm_it != rm_vector_twrhits_alltrig.end(); ++rm_it)
0241 {
0242 (*rm_it)->Reset();
0243 }
0244 if (anaGL1)
0245 {
0246 OnlMonServer *se = OnlMonServer::instance();
0247 se->UseGl1();
0248 }
0249 return 0;
0250 }
0251
0252
0253 std::vector<float> CemcMon::getSignal(Packet *p, const int channel)
0254 {
0255 double baseline = 0;
0256 for (int s = 0; s < 3; s++)
0257 {
0258 baseline += p->iValue(s, channel);
0259 }
0260 baseline /= 3.;
0261
0262 double signal = 0;
0263 float x = 0;
0264 for (int s = 3; s < p->iValue(0, "SAMPLES"); s++)
0265 {
0266 x++;
0267 signal += p->iValue(s, channel) - baseline;
0268 }
0269
0270 signal /= x;
0271
0272 std::vector<float> result;
0273 result.push_back(signal);
0274 result.push_back(2);
0275 result.push_back(1);
0276 return result;
0277 }
0278
0279 std::vector<float> CemcMon::anaWaveformFast(Packet *p, const int channel)
0280 {
0281 std::vector<float> waveform;
0282
0283
0284 if (p->iValue(channel, "SUPPRESSED"))
0285 {
0286 waveform.push_back(p->iValue(channel, "PRE"));
0287 waveform.push_back(p->iValue(channel, "POST"));
0288 }
0289 else
0290 {
0291 int nSamples = p->iValue(0, "SAMPLES");
0292 waveform.reserve(nSamples);
0293 for (int s = 0; s < nSamples; s++)
0294 {
0295 waveform.push_back(p->iValue(s, channel));
0296 }
0297 }
0298
0299 std::vector<std::vector<float>> multiple_wfs;
0300 multiple_wfs.push_back(waveform);
0301
0302 std::vector<std::vector<float>> fitresults_cemc;
0303 fitresults_cemc = WaveformProcessingFast->calo_processing_fast(multiple_wfs);
0304
0305 std::vector<float> result;
0306 result = fitresults_cemc.at(0);
0307 return result;
0308 }
0309
0310 std::vector<float> CemcMon::anaWaveformTemp(Packet *p, const int channel)
0311 {
0312 std::vector<float> waveform;
0313
0314 if (p->iValue(channel, "SUPPRESSED"))
0315 {
0316 waveform.push_back(p->iValue(channel, "PRE"));
0317 waveform.push_back(p->iValue(channel, "POST"));
0318 }
0319 else
0320 {
0321 waveform.reserve(p->iValue(0, "SAMPLES"));
0322 for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0323 {
0324 waveform.push_back(p->iValue(s, channel));
0325 }
0326 }
0327 std::vector<std::vector<float>> multiple_wfs;
0328 multiple_wfs.push_back(waveform);
0329
0330 std::vector<std::vector<float>> fitresults_cemc;
0331 fitresults_cemc = WaveformProcessingTemp->process_waveform(multiple_wfs);
0332
0333 std::vector<float> result;
0334 result = fitresults_cemc.at(0);
0335 return result;
0336 }
0337
0338 int CemcMon::process_event(Event *e )
0339 {
0340 float sectorAvg[Nsector] = {0};
0341 unsigned int towerNumber = 0;
0342 bool fillhist = true;
0343 std::vector<bool> trig_bools;
0344 trig_bools.resize(64);
0345 long long int gl1_clock = 0;
0346 bool have_gl1 = false;
0347 if (anaGL1)
0348 {
0349 int evtnr = e->getEvtSequence();
0350 Event *gl1Event = erc->getEvent(evtnr);
0351 if (gl1Event)
0352 {
0353 OnlMonServer *se = OnlMonServer::instance();
0354 se->IncrementGl1FoundCounter();
0355 have_gl1 = true;
0356 Packet *p = gl1Event->getPacket(14001);
0357 h_evtRec->Fill(0.0, 1.0);
0358 if (p)
0359 {
0360 gl1_clock = p->lValue(0, "BCO");
0361 uint64_t triggervec = p->lValue(0, "ScaledVector");
0362 for (int i = 0; i < 64; i++)
0363 {
0364 bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0365 trig_bools[i] = trig_decision;
0366 if (trig_decision)
0367 {
0368 h1_cemc_trig->Fill(i);
0369 }
0370 triggervec = (triggervec >> 1U);
0371 }
0372 delete p;
0373 }
0374 delete gl1Event;
0375 }
0376 else
0377 {
0378 std::cout << "GL1 event is null" << std::endl;
0379 h_evtRec->Fill(0.0, 0.0);
0380 }
0381
0382
0383
0384 if(usembdtrig){
0385 if(trig_bools.at(TriggerEnum::BitCodes::MBD_NS2_ZVRTX10) == 0 && trig_bools.at(TriggerEnum::BitCodes::MBD_NS2_ZVRTX30) == 0 && trig_bools.at(TriggerEnum::BitCodes::MBD_NS2_ZVRTX150) == 0){
0386 fillhist = false;
0387 }
0388 }
0389
0390 }
0391
0392
0393 eventCounter++;
0394 int one = 1;
0395 int zero = 0;
0396 for (int packet = packetlow; packet <= packethigh; packet++)
0397 {
0398 Packet *p = e->getPacket(packet);
0399
0400 if (p)
0401 {
0402 unsigned int adc_skip_mask = cdbttree->GetIntValue(packet, "adcskipmask");
0403 h1_packet_number->Fill(packet);
0404 h1_packet_length->SetBinContent(packet - 6000, h1_packet_length->GetBinContent(packet - 6000) + p->getLength());
0405
0406 h1_packet_event->SetBinContent(packet - 6000, p->lValue(0, "CLOCK"));
0407
0408 if (have_gl1)
0409 {
0410 long long int p_clock = p->lValue(0, "CLOCK");
0411 long long int diff = (p_clock - gl1_clock) % 65536;
0412 h2_caloPack_gl1_clock_diff->Fill(packet, diff);
0413 }
0414 int nChannels = p->iValue(0, "CHANNELS");
0415 int skiped_channel = 0;
0416 if (nChannels > m_nChannels)
0417 {
0418 return -1;
0419 }
0420
0421 for (int c = 0; c < nChannels; c++)
0422 {
0423 h1_packet_chans->Fill(packet);
0424
0425 if(c % 64 ==0){
0426 unsigned int adcboard = (unsigned int) c / 64;
0427 if ((adc_skip_mask >> adcboard) & 0x1U){
0428
0429 towerNumber += 64;
0430 skiped_channel += 64;
0431 }
0432 }
0433
0434
0435 towerNumber++;
0436
0437 unsigned int key = TowerInfoDefs::encode_emcal(towerNumber - 1);
0438 unsigned int phi_bin = TowerInfoDefs::getCaloTowerPhiBin(key);
0439 unsigned int eta_bin = TowerInfoDefs::getCaloTowerEtaBin(key);
0440
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455 std::vector<float> resultFast = anaWaveformFast(p, c);
0456 float signalFast = resultFast.at(0);
0457 float timeFast = resultFast.at(1);
0458 float pedestalFast = resultFast.at(2);
0459 int bin = h2_cemc_mean->FindBin(eta_bin + 0.5, phi_bin + 0.5);
0460 float prepost = p->iValue(c, "PRE") - p->iValue(c, "POST");
0461
0462 if (fillhist)
0463 {
0464 if (p->iValue(c, "SUPPRESSED"))
0465 {
0466 p2_zsFrac_etaphi->Fill(eta_bin, phi_bin, 0);
0467 }
0468 else
0469 {
0470 p2_zsFrac_etaphi->Fill(eta_bin, phi_bin, 1);
0471 }
0472
0473 h1_waveform_pedestal->Fill(pedestalFast);
0474
0475 rm_vector_twr[towerNumber - 1]->Add(&signalFast);
0476
0477 if (signalFast > hit_threshold)
0478 {
0479 rm_vector_twrhits[towerNumber - 1]->Add(&one);
0480 h2_cemc_hits->SetBinContent(bin, h2_cemc_hits->GetBinContent(bin) + 1);
0481 }
0482 else
0483 {
0484 rm_vector_twrhits[towerNumber - 1]->Add(&zero);
0485 }
0486 h2_cemc_mean->SetBinContent(bin, h2_cemc_mean->GetBinContent(bin) + signalFast);
0487 h2_cemc_rm->SetBinContent(bin, rm_vector_twr[towerNumber - 1]->getMean(0));
0488 h2_cemc_rmhits->SetBinContent(bin, rm_vector_twrhits[towerNumber - 1]->getMean(0));
0489 h1_cemc_adc->Fill(signalFast);
0490 }
0491
0492
0493 if (p->iValue(c, "SUPPRESSED"))
0494 {
0495 p2_zsFrac_etaphi_all->Fill(eta_bin, phi_bin, 0);
0496 }
0497 else
0498 {
0499 p2_zsFrac_etaphi_all->Fill(eta_bin, phi_bin, 1);
0500 }
0501
0502
0503
0504 if (signalFast > waveform_hit_threshold)
0505 {
0506
0507 if (!isHottower(packet, c))
0508 {
0509 for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0510 {
0511 h2_waveform_twrAvg->Fill(s, p->iValue(s, c) - pedestalFast);
0512 }
0513 h1_waveform_time->Fill(timeFast);
0514 }
0515 }
0516 if (signalFast > hit_threshold)
0517 {
0518
0519 if (have_gl1)
0520 {
0521 for (int itrig = 0; itrig < 64; itrig++)
0522 {
0523 if (trig_bools[itrig])
0524 {
0525 h2_cemc_hits_trig[itrig]->Fill(eta_bin + 0.5, phi_bin + 0.5);
0526 }
0527 }
0528 }
0529 rm_vector_twrhits_alltrig[towerNumber - 1]->Add(&one);
0530 }
0531 else
0532 {
0533 rm_vector_twrhits_alltrig[towerNumber - 1]->Add(&zero);
0534 }
0535 if(prepost > 0)
0536 {
0537 p2_pre_post->Fill(eta_bin, phi_bin, prepost);
0538 p2_pre_post->Fill(eta_bin, phi_bin, -prepost);
0539 }
0540
0541 h2_cemc_rmhits_alltrig->SetBinContent(bin, rm_vector_twrhits_alltrig[towerNumber - 1]->getMean(0));
0542
0543 if (signalFast > chi2_check_threshold)
0544 {
0545 std::vector<float> resultTemp = anaWaveformTemp(p, c);
0546 float chi2 = resultTemp.at(3);
0547 float signalTemp = resultTemp.at(0);
0548 if(chi2 > signalTemp*signalTemp / 50. )
0549 {
0550 p2_bad_chi2->Fill(eta_bin, phi_bin, 1);
0551 }
0552 else
0553 {
0554 p2_bad_chi2->Fill(eta_bin, phi_bin, 0);
0555 }
0556 }
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573 }
0574 if ((nChannels + skiped_channel) < m_nChannels)
0575 {
0576
0577
0578 for (int channel = 0; channel < m_nChannels - (nChannels + skiped_channel); channel++)
0579 {
0580 towerNumber++;
0581
0582 unsigned int key = TowerInfoDefs::encode_emcal(towerNumber - 1);
0583 unsigned int phi_bin = TowerInfoDefs::getCaloTowerPhiBin(key);
0584 unsigned int eta_bin = TowerInfoDefs::getCaloTowerEtaBin(key);
0585
0586 int sectorNumber = phi_bin / 8 + 1;
0587
0588 int bin = h2_cemc_mean->FindBin(eta_bin + 0.5, phi_bin + 0.5);
0589
0590 sectorAvg[sectorNumber - 1] += 0.;
0591
0592 float signalFast = 0.0;
0593
0594 rm_vector_twr[towerNumber - 1]->Add(&signalFast);
0595
0596 h2_cemc_rm->SetBinContent(bin, rm_vector_twr[towerNumber - 1]->getMean(0));
0597 h2_cemc_rmhits->SetBinContent(bin, rm_vector_twrhits[towerNumber - 1]->getMean(0));
0598
0599 h2_cemc_mean->SetBinContent(bin, h2_cemc_mean->GetBinContent(bin));
0600 }
0601 }
0602 delete p;
0603 }
0604 else
0605 {
0606 for (int channel = 0; channel < m_nChannels; channel++)
0607 {
0608 towerNumber++;
0609 unsigned int key = TowerInfoDefs::encode_emcal(towerNumber - 1);
0610 unsigned int phi_bin = TowerInfoDefs::getCaloTowerPhiBin(key);
0611 unsigned int eta_bin = TowerInfoDefs::getCaloTowerEtaBin(key);
0612
0613 int sectorNumber = phi_bin / 8 + 1;
0614
0615 int bin = h2_cemc_mean->FindBin(eta_bin + 0.5, phi_bin + 0.5);
0616
0617 sectorAvg[sectorNumber - 1] += 0;
0618
0619 float signalFast = 0;
0620
0621 rm_vector_twr[towerNumber - 1]->Add(&signalFast);
0622
0623 h2_cemc_rm->SetBinContent(bin, rm_vector_twr[towerNumber - 1]->getMean(0));
0624 h2_cemc_rmhits->SetBinContent(bin, rm_vector_twrhits[towerNumber - 1]->getMean(0));
0625
0626 h2_cemc_mean->SetBinContent(bin, h2_cemc_mean->GetBinContent(bin));
0627 }
0628 }
0629 }
0630
0631 h1_event->Fill(0);
0632
0633 eventCounter++;
0634 return 0;
0635 }
0636
0637 int CemcMon::Reset()
0638 {
0639
0640 eventCounter = 0;
0641 idummy = 0;
0642 return 0;
0643 }