File indexing completed on 2025-12-17 09:21:17
0001 #include "InttMixupQA.h"
0002
0003 #include <ffarawobjects/InttRawHitContainer.h>
0004
0005 #include <qautils/QAHistManagerDef.h>
0006
0007 #include <ffarawobjects/InttRawHit.h>
0008
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>
0012
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h> // PHWHERE
0015
0016 #include <TCanvas.h>
0017 #include <TFile.h>
0018 #include <TGraph.h>
0019 #include <TH1.h>
0020 #include <TH2.h>
0021 #include <TLegend.h>
0022 #include <TLine.h>
0023 #include <TROOT.h>
0024 #include <TStyle.h>
0025 #include <TSystem.h>
0026
0027 #include <algorithm>
0028 #include <cstdint>
0029 #include <cstdlib>
0030 #include <format>
0031 #include <fstream>
0032 #include <iostream>
0033 #include <limits>
0034 #include <map>
0035 #include <sstream>
0036 #include <string>
0037 #include <utility>
0038
0039 InttMixupQA::InttMixupQA(const std::string &name, const int run_num, const int felix_num)
0040 : SubsysReco(name)
0041 , run_num_(run_num)
0042 , felix_num_(felix_num)
0043 {
0044 std::cout << "felix_num_=" << felix_num_ << "felix_num=" << felix_num << std::endl;
0045 }
0046
0047 int InttMixupQA::Init(PHCompositeNode * )
0048 {
0049 if (Verbosity() > 5)
0050 {
0051 std::cout << "Beginning Init in InttMixupQA" << std::endl;
0052 }
0053
0054 return 0;
0055 }
0056
0057 int InttMixupQA::InitRun(PHCompositeNode *topNode)
0058 {
0059 if (!topNode)
0060 {
0061 std::cout << "InttMixupQA::InitRun(PHCompositeNode* topNode)" << std::endl;
0062 std::cout << "\tCould not retrieve topNode; doing nothing" << std::endl;
0063
0064 return 1;
0065 }
0066
0067 std::cout << "felix_num_=" << felix_num_ << std::endl;
0068
0069
0070 for (int felix = 0; felix < kFelix_num_; felix++)
0071 {
0072 std::string name = std::format("{}_allmulti_intt{}", getHistoPrefix(), felix);
0073 std::string title = name + std::format("_Run{}", run_num_) + ": with clone cut";
0074 h_allmulti_[felix] = new TH1F(name.c_str(), title.c_str(), 200, 0, bin);
0075 h_allmulti_[felix]->SetXTitle("Multiplicity");
0076
0077
0078 name = std::format("{}_allclone_intt{}", getHistoPrefix(), felix);
0079 title = name + std::format("_Run{}", run_num_);
0080 h_allclone_[felix] = new TH1F(name.c_str(), title.c_str(), 200, 0, 200);
0081 h_allclone_[felix]->SetXTitle("Clone multiplicity");
0082 h_allclone_[felix]->SetLineColor(felix + 1);
0083
0084
0085 auto *hm = QAHistManagerDef::getHistoManager();
0086 if (!hm)
0087 {
0088 std::cerr << PHWHERE << "\n"
0089 << "\tQAHistManagerDef::getHistoManager() returned null\n"
0090 << "\texit(1)" << std::endl;
0091 gSystem->Exit(1);
0092 exit(1);
0093 }
0094
0095 name = std::format("{}_bco_full&0x7F_prev_vs_bco_intt{}", getHistoPrefix(), felix);
0096 title = name + std::format("_Run{}", run_num_);
0097 h_vsprefull_bco[felix] = new TH2F(name.c_str(), title.c_str(), 128, 0, 128, 128, 0, 128);
0098 h_vsprefull_bco[felix]->SetXTitle("BCO");
0099 h_vsprefull_bco[felix]->SetYTitle("BCO_FULL previous event &0x7F");
0100 hm->registerHisto(h_vsprefull_bco[felix]);
0101
0102 name = std::format("{}_bco_full&0x7F_vs_bco_intt{}", getHistoPrefix(), felix);
0103 title = name + std::format("_Run{}", run_num_);
0104 h_vsfull_bco[felix] = new TH2F(name.c_str(), title.c_str(), 128, 0, 128, 128, 0, 128);
0105 h_vsfull_bco[felix]->SetXTitle("BCO");
0106 h_vsfull_bco[felix]->SetYTitle("BCO_FULL &0x7F");
0107 hm->registerHisto(h_vsfull_bco[felix]);
0108
0109 name = std::format("{}_bco_full&0x7F_vs_bco_intt{}", getHistoPrefix(), felix);
0110 title = name + std::format("_Run{}", run_num_);
0111 h_prefull_bco[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0112 h_prefull_bco[felix]->SetXTitle("BCO_FULL previous event - BCO");
0113 h_prefull_bco[felix]->SetMinimum(0);
0114 hm->registerHisto(h_prefull_bco[felix]);
0115
0116 name = std::format("{}_bco_full&0x7F_bco_intt{}", getHistoPrefix(), felix);
0117 title = name + std::format("_Run{}", run_num_);
0118 h_full_bco[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0119 h_full_bco[felix]->SetXTitle("BCO_FULL - BCO");
0120 h_full_bco[felix]->SetMinimum(0);
0121 hm->registerHisto(h_full_bco[felix]);
0122
0123 name = std::format("{}_bco_full&0x7F_prev_bco_intt_all{}", getHistoPrefix(), felix);
0124 title = name + std::format("_Run{}", run_num_);
0125 h_prefull_bco_all[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0126 h_prefull_bco_all[felix]->SetXTitle("BCO_FULL previous event - BCO");
0127 h_prefull_bco_all[felix]->SetMinimum(0);
0128
0129 name = std::format("{}_bco_full&0x7F_prev_vs_bco_intt_all{}", getHistoPrefix(), felix);
0130 title = name + std::format("_Run{}", run_num_);
0131 h_vsprefull_bco_all[felix] = new TH2F(name.c_str(), title.c_str(), 128, 0, 128, 128, 0, 128);
0132 h_vsprefull_bco_all[felix]->SetXTitle("BCO");
0133 h_vsprefull_bco_all[felix]->SetYTitle("BCO_FULL previous event &0x7F");
0134
0135 std::string const name1 = "Mixup Multiplicity" + std::to_string(felix);
0136 std::string const name2 = "Mixup/all Multiplcity" + std::to_string(felix);
0137 std::string const title1 = name1 + std::format("_Run{}", run_num_);
0138 std::string const title2 = name2 + std::format("_Run{}", run_num_);
0139 for (int p = 0; p < divimul; p++)
0140 {
0141 h_mixupmulti[felix][p] = new TH1F((name1 + std::format("_{}", p)).c_str(), (title1 + std::format("_{}", p)).c_str(), 200, 0, bin);
0142 h_divmul[felix][p] = new TH1F((title1 + std::format("_{}", p)).c_str(), (title2 + std::format("_{}", p)).c_str(), 200, 0, bin);
0143 ;
0144 }
0145
0146 name = std::format("{}_Number of mixup hit{}", getHistoPrefix(), felix);
0147 title = name + std::format("_Run{}", run_num_);
0148 h_mixup[felix] = new TH1F(name.c_str(), title.c_str(), 400, 0, 400);
0149
0150 name = std::format("{}_prev_allhit vs Nmixup{}", getHistoPrefix(), felix);
0151 title = name + std::format("_Run{}", run_num_);
0152 h_prevsNmix[felix] = new TH2F(name.c_str(), title.c_str(), divimul + 10, 0, divimul + 10, bin, 0, bin);
0153 h_prevsNmix[felix]->GetXaxis()->SetNdivisions(405);
0154 h_prevsNmix[felix]->GetYaxis()->SetNdivisions(405);
0155 h_prevsNmix[felix]->SetXTitle("Number of Mixup hits");
0156 h_prevsNmix[felix]->SetYTitle("Number of previous event hits");
0157
0158 name = std::format("{}_bco_full_prev-bco No copyhit{}", getHistoPrefix(), felix);
0159 title = name + std::format("_Run{}", run_num_);
0160 h_nocopyhit[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0161 h_nocopyhit[felix]->SetXTitle("bco_full_prev-bco");
0162
0163 name = std::format("{}_bco_full_prev-bco copyhit{}", getHistoPrefix(), felix);
0164 title = name + std::format("_Run{}", run_num_);
0165 h_copyhit[felix] = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0166 h_copyhit[felix]->SetXTitle("bco_full_prev-bco");
0167
0168 name = std::format("{}_Mixup & copy hit {}", getHistoPrefix(), felix);
0169 title = name + std::format("_Run{}", run_num_);
0170 h_mixcopy[felix] = new TH1F(name.c_str(), title.c_str(), 50, 0, 50);
0171 h_mixcopy[felix]->SetXTitle("Number of Mixup copy hit");
0172
0173 name = std::format("{}_Nmixup vs Nclone{}", getHistoPrefix(), felix);
0174 title = name + std::format("_Run{}", run_num_);
0175 h_mixvscopy[felix] = new TH2F(name.c_str(), title.c_str(), 100, 0, 100, 50, 0, 50);
0176 h_mixvscopy[felix]->SetXTitle("Nmixup");
0177 h_mixvscopy[felix]->SetYTitle("Ncopy");
0178
0179 name = std::format("{}_Nmixup vs pre_Nhits{}", getHistoPrefix(), felix);
0180 title = name + std::format("_Run{}", run_num_);
0181 h_hitfra[felix] = new TH2F(name.c_str(), title.c_str(), bin, 0, bin, bin, 0, bin);
0182
0183 name = std::format("{}_Nmixup vs pre_Nhits others 4bin{}", getHistoPrefix(), felix);
0184 title = name + std::format("_Run{}", run_num_);
0185 h_bghit[felix] = new TH2F(name.c_str(), title.c_str(), bin, 0, bin, bin, 0, bin);
0186
0187 fNhit[felix].open(std::format("./txtfile/NHit_{}_intt{}.txt", run_num_, felix));
0188
0189 bcopeak_file[felix] = bcopeak_dir_ + std::format("bco_000{}_intt{}.root", run_num_, felix);
0190 }
0191
0192 std::string name = "bco_full - prev_bco_full";
0193 std::string title = name + std::format("_Run{}", run_num_);
0194 h_interval = new TH1F(name.c_str(), title.c_str(), 200, 0, bin2);
0195 h_interval->SetXTitle("bco_full - prev_bco_full");
0196
0197 name = "bco_full - prev_bco_full Mixup";
0198 title = name + std::format("_Run{}", run_num_);
0199 h_mixinterval = new TH1F(name.c_str(), title.c_str(), 200, 0, bin2);
0200
0201 name = "Interval normarize";
0202 title = name + std::format("_Run{}", run_num_);
0203 h_divinter = new TH1F(name.c_str(), title.c_str(), 200, 0, bin2);
0204 h_divinter->SetXTitle("bco_full - prev_bco_full");
0205
0206 name = "bco_full";
0207 title = name + std::format("_Run{}", run_num_);
0208 h_bcofull_7 = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0209 h_bcofull_7->SetXTitle("bco_full &0x7F");
0210
0211 name = "bco";
0212 title = name + std::format("_Run{}", run_num_);
0213 h_bco = new TH1F(name.c_str(), title.c_str(), 128, 0, 128);
0214 h_bco->SetXTitle("bco");
0215
0216 name = "felix vs NmixupEv";
0217 title = name + std::format("_Run{}", run_num_);
0218
0219 h_NmixEv = new TH1F(name.c_str(), title.c_str(), 8, 0, 8);
0220
0221 name = "All Event";
0222 title = name + std::format("_Run{}", run_num_);
0223 h_AllEv = new TH1F(name.c_str(), title.c_str(), 8, 0, 8);
0224
0225 g_evfraction = new TGraph(n);
0226
0227 g_cloevfraction = new TGraph(n);
0228
0229 g_hitfraction = new TGraph(n);
0230
0231 g_copyfraction = new TGraph(n);
0232
0233
0234
0235
0236
0237
0238
0239 tf_output_ = new TFile(output_root_.c_str(), "RECREATE");
0240
0241 GetBcopeak();
0242
0243 Readpeak();
0244
0245
0246
0247 return Fun4AllReturnCodes::EVENT_OK;
0248 }
0249
0250 int InttMixupQA::process_event(PHCompositeNode *topNode)
0251 {
0252 if (Verbosity() > 5)
0253 {
0254 std::cout << "prev_bcofull=" << prev_bcofull << " "
0255 << "pre_allhit=" << pre_allhit << std::endl;
0256 }
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277 std::string const m_InttRawNodeName = "INTTRAWHIT";
0278 InttRawHitContainer *inttcont = findNode::getClass<InttRawHitContainer>(topNode, m_InttRawNodeName);
0279 if (!inttcont)
0280 {
0281 std::cout << PHWHERE << std::endl;
0282 std::cout << "InttMixupQA::process_event(PHCompositeNode* topNode)" << std::endl;
0283 std::cout << "Could not get \"" << m_InttRawNodeName << "\" from Node Tree" << std::endl;
0284 std::cout << "Exiting" << std::endl;
0285 gSystem->Exit(1);
0286 exit(1);
0287 }
0288
0289 uint64_t const longbco_full = (0 < inttcont->get_nhits()) ? inttcont->get_hit(0)->get_bco() : std::numeric_limits<uint64_t>::max();
0290 uint64_t const difevent_bcofull = (longbco_full & bit) - (long_prev_bcofull & bit);
0291 h_interval->Fill(difevent_bcofull);
0292 uint64_t const bco_full = longbco_full & 0x7FU;
0293
0294
0295 ievent_++;
0296
0297 if (ievent_ == 1)
0298 {
0299 first_bcofull = longbco_full;
0300 }
0301
0302
0303
0304
0305
0306 if ((ievent_ % 100 == 0 && ievent_ < 1000) || ievent_ % 1000 == 0)
0307 {
0308 std::cout << "Process event #" << ievent_ << std::endl;
0309 }
0310
0311 int const nhits = inttcont->get_nhits();
0312 if (Verbosity() > 5)
0313 {
0314 std::cout << "Nhits = " << nhits << std::endl;
0315 }
0316
0317 std::map<int, int> map_hit;
0318 int nhit_fx[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0319 int ncln_fx[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0320
0321 int Nmixup[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0322 int Nother[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0323
0324 int Ncopy[8] = {0, 0, 0, 0, 0, 0, 0, 0};
0325
0326 for (int i = 0; i < nhits; i++)
0327 {
0328 InttRawHit *intthit = inttcont->get_hit(i);
0329
0330 int const fnum = intthit->get_packetid() - 3001;
0331 int const fchn = intthit->get_fee();
0332 int const adc = intthit->get_adc();
0333 int const chip = intthit->get_chip_id();
0334 int const chan = intthit->get_channel_id();
0335 int const bco = intthit->get_FPHX_BCO();
0336 int const hitID = 100000000 * (fnum + 1) + 1000000 * fchn + 10000 * chip + 10 * chan + adc;
0337
0338
0339
0340 auto itrHit = map_hit.find(hitID);
0341 if (itrHit == map_hit.end())
0342 {
0343 map_hit.insert(std::make_pair(hitID, 0));
0344 if (Verbosity() > 5)
0345 {
0346 std::cout << hitID << " " << fnum + 1 << " " << fchn << " " << chip << " " << chan << " " << adc << " " << bco << std::endl;
0347 }
0348
0349
0350
0351
0352 {
0353 nhit_fx[fnum]++;
0354
0355 h_bco->Fill(bco);
0356
0357 int bco_diff = bco_full - bco;
0358
0359 if (bco_diff < 0)
0360 {
0361 bco_diff = bco_diff + 128;
0362 }
0363
0364 int pre_bco_diff = prev_bcofull - bco;
0365
0366 if (pre_bco_diff < 0)
0367 {
0368 pre_bco_diff = pre_bco_diff + 128;
0369 }
0370
0371
0372
0373
0374
0375
0376 {
0377 if (!bcopar_[fnum].empty() && !bcopar_[fnum].contains(bco_diff))
0378 {
0379 if (bcopar_[fnum].contains(pre_bco_diff))
0380 {
0381 Nmixup[fnum]++;
0382
0383 auto pre_itrHit = premap_hit.find(hitID);
0384 if (pre_itrHit != premap_hit.end())
0385 {
0386 h_copyhit[fnum]->Fill(pre_bco_diff);
0387 Ncopy[fnum]++;
0388 }
0389 else
0390 {
0391
0392 h_nocopyhit[fnum]->Fill(pre_bco_diff);
0393 }
0394 }
0395 else if (otbcopar_[fnum].contains(pre_bco_diff))
0396 {
0397 Nother[fnum]++;
0398 }
0399
0400 h_prefull_bco[fnum]->Fill(pre_bco_diff);
0401 h_vsprefull_bco[fnum]->Fill(bco, prev_bcofull);
0402 }
0403 h_prefull_bco_all[fnum]->Fill(pre_bco_diff);
0404 h_vsprefull_bco_all[fnum]->Fill(bco, prev_bcofull);
0405
0406 h_full_bco[fnum]->Fill(bco_diff);
0407 h_vsfull_bco[fnum]->Fill(bco, bco_full);
0408 }
0409 }
0410
0411
0412
0413
0414 }
0415 else
0416 {
0417 ncln_fx[fnum]++;
0418 }
0419
0420 if (Verbosity() > 5)
0421 {
0422 std::cout << "bco_full=" << bco_full << " "
0423 << "bco=" << bco << " " << std::endl;
0424 }
0425 }
0426
0427 if (Verbosity() > 5)
0428 {
0429 for (int id = 0; id < 8; id++)
0430 {
0431 std::cout << "Felix#" << id << " "
0432 << "Nmixup=" << Nmixup[id] << " "
0433 << "Ncopy=" << Ncopy[id] << " "
0434 << "Nclone=" << ncln_fx[id] << std::endl;
0435 }
0436 }
0437
0438 for (int felix = 0; felix < kFelix_num_; felix++)
0439 {
0440
0441 h_allmulti_[felix]->Fill(pre_allhit[felix]);
0442 h_allclone_[felix]->Fill(ncln_fx[felix]);
0443 h_mixup[felix]->Fill(Nmixup[felix]);
0444 h_mixcopy[felix]->Fill(Ncopy[felix]);
0445 h_hitfra[felix]->Fill(Nmixup[felix], pre_allhit[felix]);
0446 h_bghit[felix]->Fill(Nother[felix], pre_allhit[felix]);
0447 for (int p = 0; p < divimul; p++)
0448 {
0449 if (Nmixup[felix] == p + 1)
0450
0451 {
0452 h_mixupmulti[felix][p]->Fill(pre_allhit[felix]);
0453
0454
0455
0456
0457
0458 }
0459 else
0460 {
0461 continue;
0462 }
0463
0464 h_divmul[felix][p]->Divide(h_mixupmulti[felix][p], h_allmulti_[felix]);
0465 }
0466
0467 mixupfraction[felix] = (double(Nmixup[felix]) / double(pre_allhit[felix] + Nmixup[felix]));
0468
0469
0470
0471 h_AllEv->Fill(felix);
0472 if (Nmixup[felix] > 0 && pre_allhit[felix] > 0)
0473
0474 {
0475
0476 h_mixinterval->Fill(difevent_bcofull);
0477 h_divinter->Divide(h_mixinterval, h_interval);
0478
0479 h_prevsNmix[felix]->Fill(Nmixup[felix], pre_allhit[felix]);
0480 NmixupEv[felix]++;
0481 h_NmixEv->Fill(felix);
0482
0483
0484
0485 mixupfraction_sum[felix] = mixupfraction_sum[felix] + mixupfraction[felix];
0486 Nmixup_sum[felix] = Nmixup_sum[felix] + Nmixup[felix];
0487 pre_allhit_sum[felix] = pre_allhit_sum[felix] + pre_allhit[felix];
0488 if (Ncopy[felix] > 0)
0489 {
0490 NmixcopyEv[felix]++;
0491 copyfraction_sum[felix] = copyfraction_sum[felix] + copyfraction[felix];
0492 }
0493 if (Verbosity() > 5)
0494 {
0495 std::cout << "Felix=" << felix << " "
0496 << "mixupfraction=" << mixupfraction[felix] << "Nmixup=" << Nmixup[felix] << "pre_allhit=" << pre_allhit[felix] << "Ncopy=" << Ncopy[felix] << "copyfraction=" << copyfraction[felix] << std::endl;
0497 }
0498 }
0499
0500 if (felix == felix_num_)
0501 {
0502
0503
0504 fNhit[felix] << nhit_fx[felix] << std::endl;
0505 }
0506 pre_allhit[felix] = nhit_fx[felix];
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516 }
0517
0518 long_prev_bcofull = longbco_full;
0519 prev_bcofull = bco_full;
0520 premap_hit = map_hit;
0521
0522 return Fun4AllReturnCodes::EVENT_OK;
0523 }
0524
0525 int InttMixupQA::End(PHCompositeNode * )
0526 {
0527
0528
0529 if (Verbosity() > 1)
0530 {
0531 std::cout << "Processing InttMixupQA done" << std::endl;
0532 }
0533
0534
0535
0536
0537 this->DrawHists();
0538
0539
0540
0541
0542
0543
0544
0545 return Fun4AllReturnCodes::EVENT_OK;
0546 }
0547
0548 int InttMixupQA::SetOutputDir(std::string const &dir)
0549 {
0550 output_dir_ = dir;
0551
0552 std::string const run_num_str = std::string(8 - std::to_string(run_num_).size(), '0') + std::to_string(run_num_);
0553 std::string const fnum_str = std::to_string(felix_num_);
0554
0555 output_root_ = output_dir_ + "root/" + output_basename_ + run_num_str + "intt" + fnum_str + ".root";
0556 output_pdf_ = output_dir_ + "plots/" + output_basename_ + run_num_str + "intt" + fnum_str + ".pdf";
0557
0558
0559
0560
0561 output_txt_ = output_dir_ + "txtfile/";
0562
0563 return Fun4AllReturnCodes::EVENT_OK;
0564 }
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590 void InttMixupQA::GetBcopeak()
0591 {
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609 for (int i = 0; i < 8; i++)
0610 {
0611 tf_bcopeak_[i] = new TFile(bcopeak_file[i].c_str(), "READ");
0612
0613 hbco[i] = dynamic_cast<TH1D *>(gROOT->FindObject(std::format("h2_bco_felix_{}", i).c_str()));
0614 if (!hbco[i])
0615 {
0616 if (Verbosity())
0617 {
0618 std::cerr << PHWHERE << std::endl;
0619 }
0620 continue;
0621 }
0622
0623
0624 int const maxbin_sub = hbco[i]->GetMaximumBin();
0625 DEFAULT_BCO_VALUE[i] = hbco[i]->GetBinLowEdge(maxbin_sub);
0626 }
0627
0628
0629 int peakbco[8];
0630 double bg_mean[8];
0631 double bg_rms[8];
0632
0633
0634 for (int id = 0; id < 8; id++)
0635 {
0636 double const max = hbco[id]->GetMaximum();
0637 int const maxbin = hbco[id]->GetMaximumBin();
0638 peakbco[id] = hbco[id]->GetBinLowEdge(maxbin);
0639 std::cout << id << " " << max << " " << maxbin << " " << peakbco[id] << std::endl;
0640
0641 hbcohist[id] = new TH1F(std::format("hbcohist_{}", id).c_str(), std::format("hbcohist_{}", id).c_str(), 200, 0, max * 1.1);
0642 int const N = hbco[id]->GetNbinsX();
0643 for (int ibin = 0; ibin < N; ibin++)
0644 {
0645 double const cont = hbco[id]->GetBinContent(ibin + 1);
0646 hbcohist[id]->Fill(cont);
0647 }
0648 }
0649
0650 for (int id = 0; id < 8; id++)
0651 {
0652 double const max = hbco[id]->GetMaximum();
0653 hbcohist2[id] = new TH1F(std::format("hbcohist2_{}", id).c_str(), std::format("hbcohist2_{}", id).c_str(), 200, 0, max * 1.1);
0654
0655 int const N = hbco[id]->GetNbinsX();
0656 double const mean = hbcohist[id]->GetMean();
0657 double const rms = hbcohist[id]->GetRMS();
0658
0659 for (int ibin = 0; ibin < N; ibin++)
0660 {
0661 double const cont = hbco[id]->GetBinContent(ibin + 1);
0662 if (mean + 3 * rms > cont)
0663 {
0664 hbcohist2[id]->Fill(cont);
0665 }
0666 }
0667 bg_mean[id] = hbcohist2[id]->GetMean();
0668 bg_rms[id] = hbcohist2[id]->GetRMS();
0669 }
0670
0671
0672
0673 {
0674 f_felixpeak.open(output_txt_ + std::format("bco_{}.txt", run_num_));
0675
0676 for (int id = 0; id < 8; id++)
0677 {
0678 f_felixpeak << id << " : ";
0679
0680 double const thre = 5 * bg_mean[id];
0681 ;
0682 int const maxbin = hbco[id]->GetMaximumBin();
0683
0684 int const N = hbco[id]->GetNbinsX();
0685 bool peak_found = false;
0686
0687 for (int ibin = 0; ibin < 10; ibin++)
0688 {
0689 int binid = ibin + maxbin - 10;
0690 if (binid < 1)
0691 {
0692 binid += N;
0693 }
0694 double const cont = hbco[id]->GetBinContent(binid + 1);
0695 std::cout << id << " : " << ibin << " " << binid << " " << maxbin;
0696 if (cont > thre)
0697 {
0698 int const peakbin = hbco[id]->GetBinLowEdge(binid + 1);
0699 std::cout << " exceed : " << peakbin;
0700 f_felixpeak << peakbin << " ";
0701 peak_found = true;
0702 }
0703 std::cout << std::endl;
0704 }
0705 for (int ibin = 0; ibin < 10; ibin++)
0706 {
0707 int binid = ibin + maxbin;
0708 if (binid >= N)
0709 {
0710 binid -= N;
0711 }
0712 double const cont = hbco[id]->GetBinContent(binid + 1);
0713 std::cout << id << " : " << ibin << " " << binid << " " << maxbin;
0714 if (cont > thre)
0715 {
0716 int const peakbin = hbco[id]->GetBinLowEdge(binid + 1);
0717 std::cout << " exceed : " << peakbin;
0718 f_felixpeak << peakbin << " ";
0719 peak_found = true;
0720 }
0721 std::cout << std::endl;
0722 }
0723 if (!peak_found)
0724 {
0725 f_felixpeak << DEFAULT_BCO_VALUE[id] << " ";
0726 }
0727
0728 f_felixpeak << std::endl;
0729 }
0730 f_felixpeak.close();
0731 }
0732
0733
0734 float minimum = 1000000;
0735 for (auto &id : hbco)
0736 {
0737 std::cout << id->GetMinimumBin() << std::endl;
0738 std::cout << id->GetBinContent(id->GetMinimumBin()) << std::endl;
0739 double const min = id->GetBinContent(id->GetMinimumBin());
0740 std::cout << min << std::endl;
0741 minimum = std::min<double>(min, minimum);
0742 }
0743
0744
0745
0746 TCanvas *c1 = new TCanvas("c1", "c1", 1200, 700);
0747 c1->Divide(4, 4);
0748 for (int id = 0; id < 8; id++)
0749 {
0750 c1->cd(id + 1);
0751 gPad->SetLogy();
0752 hbco[id]->SetMinimum(minimum * 0.2);
0753 hbco[id]->Draw();
0754
0755 TLine *la = new TLine(0, 5 * bg_mean[id], 128, 5 * bg_mean[id]);
0756 la->Draw("same");
0757 TLine *lb = new TLine(0, bg_mean[id] + 6 * bg_rms[id], 128, bg_mean[id] + 6 * bg_rms[id]);
0758 lb->SetLineColor(2);
0759 lb->Draw("same");
0760 }
0761
0762 for (int id = 0; id < 8; id++)
0763 {
0764 c1->cd(id + 9);
0765 gPad->SetLogy();
0766 double const max = hbcohist[id]->GetMaximum();
0767
0768 hbcohist[id]->Draw();
0769 hbcohist2[id]->SetLineColor(2);
0770 hbcohist2[id]->Draw("same");
0771
0772 TLine *lm = new TLine(bg_mean[id], 0.001, bg_mean[id], max * 1.1);
0773 lm->Draw();
0774 TLine *l5m = new TLine(5 * bg_mean[id], 0.001, 5 * bg_mean[id], max * 1.1);
0775 l5m->Draw();
0776 TLine *l = new TLine(bg_mean[id] + 6 * bg_rms[id], 0.001, bg_mean[id] + 6 * bg_rms[id], max * 1.1);
0777 l->SetLineColor(2);
0778 l->Draw();
0779 }
0780
0781 c1->Print((output_dir_ + "plots/" + std::format("check_bco_file_{}.pdf", run_num_)).c_str());
0782 }
0783
0784 void InttMixupQA::Readpeak()
0785 {
0786 std::string const bcofile_ = (output_txt_ + std::format("bco_{}.txt", run_num_));
0787
0788 if (!bcofile_.empty())
0789 {
0790 std::cout << "BCO file : " << bcofile_ << std::endl;
0791 std::ifstream file(bcofile_);
0792 if (!file.is_open())
0793 {
0794 std::cerr << "Failed to open file: " << bcofile_ << std::endl;
0795 return;
0796 }
0797
0798 std::string sline;
0799
0800 while (getline(file, sline))
0801 {
0802
0803 std::istringstream ss(sline);
0804 std::string sbuf;
0805
0806 int felix = -1;
0807 int idx = 0;
0808 while (getline(ss, sbuf, ':'))
0809 {
0810
0811 if (idx == 0)
0812 {
0813
0814 felix = std::stoi(sbuf);
0815 if (felix < 0 || felix > 7)
0816 {
0817 std::cout << "felixid out of range. id=" << felix << std::endl;
0818 break;
0819 }
0820 }
0821 else
0822 {
0823
0824 std::istringstream ss2(sbuf);
0825 std::string sbuf2;
0826 std::string sbuf3;
0827 std::string sbuf4;
0828 std::string sbuf5;
0829 std::string sbuf6;
0830 while (ss2 >> sbuf2)
0831 {
0832
0833 bcopar_[felix].insert(std::stoi(sbuf2));
0834 int const peak = std::stoi(sbuf2);
0835 int const peak1 = peak + 1;
0836 int const peak2 = peak + 2;
0837 int const peak_1 = peak - 1;
0838 int const peak_2 = peak - 2;
0839 sbuf3 = std::to_string(peak1);
0840 sbuf4 = std::to_string(peak2);
0841 sbuf5 = std::to_string(peak_1);
0842 sbuf6 = std::to_string(peak_2);
0843 otbcopar_[felix].insert(std::stoi(sbuf3));
0844 otbcopar_[felix].insert(std::stoi(sbuf4));
0845 otbcopar_[felix].insert(std::stoi(sbuf5));
0846 otbcopar_[felix].insert(std::stoi(sbuf6));
0847 }
0848 }
0849 idx++;
0850 }
0851 }
0852 file.close();
0853 }
0854
0855 std::cout << "BCO peaks " << std::endl;
0856 for (int i = 0; i < 8; i++)
0857 {
0858 std::cout << " felix " << i << " : ";
0859 for (int const itr : bcopar_[i])
0860 {
0861 std::cout << itr << " ";
0862 }
0863 std::cout << std::endl;
0864 }
0865
0866 std::cout << "Around BCO peaks " << std::endl;
0867 for (int i = 0; i < 8; i++)
0868 {
0869 std::cout << " felix " << i << " : ";
0870 for (int const itr : otbcopar_[i])
0871 {
0872 std::cout << itr << " ";
0873 }
0874 std::cout << std::endl;
0875 }
0876 }
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
0900
0901
0902
0903
0904
0905
0906
0907
0908
0909
0910 void InttMixupQA::DrawHists()
0911 {
0912 std::cout << "output pdf: " << output_pdf_ << std::endl;
0913
0914
0915 TCanvas *c0 = new TCanvas("c0", "c0");
0916 c0->Print((output_pdf_ + "[").c_str());
0917 TCanvas *c1[8];
0918 TCanvas *c2[8];
0919 TCanvas *c3[8];
0920
0921 TCanvas *c5[8];
0922 TCanvas *c6[8];
0923 TCanvas *c10[8];
0924 TCanvas *c7 = new TCanvas("canvas7", "interval", 1600, 1200);
0925 TCanvas *c8 = new TCanvas("canvas8", "fraction", 1600, 1200);
0926 TCanvas *c9 = new TCanvas("canvas9", "bco bco_full", 1600, 1200);
0927
0928 gStyle->SetOptStat(0);
0929
0930 for (int felix = 0; felix < kFelix_num_; felix++)
0931 {
0932 if (felix != felix_num_)
0933 {
0934 continue;
0935 }
0936 c1[felix] = new TCanvas(std::format("canvas1_{}", felix).c_str(), "bcofull&bco", 1600, 1200);
0937 c1[felix]->cd();
0938 gStyle->SetOptStat(0);
0939 c1[felix]->Divide(2, 2);
0940
0941 c1[felix]->cd(1);
0942 h_vsfull_bco[felix]->Draw("colz");
0943 c1[felix]->cd(2);
0944 h_vsprefull_bco[felix]->Draw("colz");
0945 c1[felix]->cd(3);
0946 h_full_bco[felix]->SetMinimum(0);
0947 h_full_bco[felix]->Draw();
0948 c1[felix]->cd(4);
0949 h_prefull_bco[felix]->SetMinimum(0);
0950 h_prefull_bco[felix]->Draw();
0951
0952 c2[felix] = new TCanvas(std::format("canvas2_{}", felix).c_str(), "Multiplicity", 1600, 1200);
0953 c2[felix]->cd();
0954 gStyle->SetOptStat(0);
0955 gPad->SetRightMargin(0.15);
0956 TLegend *legend = new TLegend(0.7, 0.45, 0.9, 0.9);
0957 gPad->SetLogy();
0958 h_allmulti_[felix]->GetXaxis()->SetNdivisions(405);
0959 h_allmulti_[felix]->GetYaxis()->SetNdivisions(405);
0960
0961 int const a = h_allmulti_[felix]->GetMaximum();
0962 h_allmulti_[felix]->GetYaxis()->SetRangeUser(1, a * 1.6);
0963 h_allmulti_[felix]->SetXTitle("Number of hits of previous event");
0964 h_allmulti_[felix]->SetYTitle("Entry");
0965 h_allmulti_[felix]->Draw();
0966 h_allmulti_[felix]->SetLineColor(1);
0967 legend->SetTextSize(0.04);
0968 legend->AddEntry(h_allmulti_[felix], "Number of", "");
0969 legend->AddEntry(h_allmulti_[felix], "Mixup hits", "");
0970 legend->AddEntry(h_allmulti_[felix], "All event", "l");
0971
0972 for (int p = 0; p < 10; p++)
0973 {
0974 h_mixupmulti[felix][p]->Draw("same");
0975 int const color = p + 2;
0976 if (color == 10)
0977 {
0978 h_mixupmulti[felix][p]->SetLineColor(46);
0979 }
0980 else
0981 {
0982 h_mixupmulti[felix][p]->SetLineColor(color);
0983 }
0984 legend->AddEntry(h_mixupmulti[felix][p], std::format("{}", p + 1).c_str(), "l");
0985 }
0986 legend->Draw();
0987
0988 c3[felix] = new TCanvas(std::format("canvas3_{}", felix).c_str(), "divide multiplicity", 1600, 1200);
0989 c3[felix]->cd();
0990 gStyle->SetOptStat(0);
0991 gPad->SetRightMargin(0.15);
0992 TLegend *leg = new TLegend(0.75, 0.5, 0.9, 0.9);
0993 gPad->SetLogy();
0994 for (int p = 0; p < 10; p++)
0995 {
0996 h_divmul[felix][p]->SetXTitle("Mixup/all hit Multiplicity");
0997 if (p == 0)
0998 {
0999 h_divmul[felix][p]->Draw();
1000 }
1001 else
1002 {
1003 h_divmul[felix][p]->Draw("same");
1004 }
1005 int const color2 = p + 2;
1006 if (color2 == 10)
1007 {
1008 h_divmul[felix][p]->SetLineColor(46);
1009 }
1010 else
1011 {
1012 h_divmul[felix][p]->SetLineColor(color2);
1013 }
1014 leg->AddEntry(h_divmul[felix][p], std::format("{}", p + 1).c_str(), "l");
1015 }
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030 c5[felix] = new TCanvas(std::format("canvas5_{}", felix).c_str(), "previous event hit vs Mixup hit", 1600, 1200);
1031 c5[felix]->cd();
1032 gStyle->SetOptStat(0);
1033 gPad->SetRightMargin(0.15);
1034 h_prevsNmix[felix]->Draw("colz");
1035
1036 c6[felix] = new TCanvas(std::format("canvas6_{}", felix).c_str(), "Mixup", 1600, 1200);
1037 c6[felix]->cd();
1038 c6[felix]->Divide(2, 2);
1039 c6[felix]->cd(1);
1040 gStyle->SetOptStat(0);
1041 gPad->SetLogy();
1042 h_mixup[felix]->Draw();
1043
1044 c6[felix]->cd(4);
1045 h_prefull_bco_all[felix]->Draw();
1046 c6[felix]->cd(2);
1047 h_vsprefull_bco_all[felix]->Draw();
1048
1049 c10[felix] = new TCanvas(std::format("canvas10_{}", felix).c_str(), "for fraction");
1050 c10[felix]->cd();
1051 c10[felix]->Divide(2, 1);
1052 c10[felix]->cd(1);
1053 h_hitfra[felix]->Draw("colz");
1054 c10[felix]->cd(2);
1055 h_bghit[felix]->Draw("colz");
1056
1057 c1[felix]->Print(output_pdf_.c_str());
1058 c2[felix]->Print(output_pdf_.c_str());
1059 c3[felix]->Print(output_pdf_.c_str());
1060
1061 c5[felix]->Print(output_pdf_.c_str());
1062 c6[felix]->Print(output_pdf_.c_str());
1063 c10[felix]->Print(output_pdf_.c_str());
1064 tf_output_->cd();
1065
1066
1067
1068
1069
1070
1071
1072 h_vsprefull_bco[felix]->Write();
1073 h_vsfull_bco[felix]->Write();
1074 h_prefull_bco[felix]->Write();
1075 h_full_bco[felix]->Write();
1076
1077
1078
1079
1080
1081 c2[felix]->Write();
1082 h_mixup[felix]->Write();
1083 h_prevsNmix[felix]->Write();
1084 h_hitfra[felix]->Write();
1085 h_bghit[felix]->Write();
1086 h_NmixEv->Write();
1087 h_AllEv->Write();
1088 h_divinter->Write();
1089
1090 }
1091
1092
1093
1094 c7->cd();
1095 gStyle->SetOptStat(0);
1096 TLegend *legsph3 = new TLegend(.55, .4, .65, .6);
1097 gPad->SetRightMargin(0.15);
1098 h_interval->Draw();
1099 h_mixinterval->Draw("same");
1100 h_mixinterval->SetLineColor(2);
1101 legsph3->AddEntry(h_interval, "All event", "l");
1102 legsph3->AddEntry(h_mixinterval, "Mixup event", "l");
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142 c9->cd();
1143 gStyle->SetOptStat(0);
1144 c9->Divide(2, 1);
1145 c9->cd(1);
1146 h_bcofull_7->Draw();
1147 c9->cd(2);
1148 h_bco->Draw();
1149
1150 tf_output_->cd();
1151 c7->Write();
1152
1153
1154
1155
1156
1157
1158
1159 tf_output_->Close();
1160
1161 c7->Print((output_pdf_).c_str());
1162 c8->Print((output_pdf_).c_str());
1163 c9->Print((output_pdf_).c_str());
1164
1165 c9->Print((output_pdf_ + "]").c_str());
1166 }