File indexing completed on 2026-05-23 08:08:41
0001 #include <algorithm>
0002 #include <cmath>
0003 #include <cstdlib>
0004 #include <glob.h>
0005 #include <iostream>
0006 #include <limits>
0007 #include <sstream>
0008 #include <set>
0009 #include <string>
0010 #include <utility>
0011 #include <vector>
0012
0013 #include "TAxis.h"
0014 #include "TCanvas.h"
0015 #include "TChain.h"
0016 #include "TColor.h"
0017 #include "TFile.h"
0018 #include "TH1D.h"
0019 #include "TLatex.h"
0020 #include "TLegend.h"
0021 #include "TLine.h"
0022 #include "TPad.h"
0023 #include "TROOT.h"
0024 #include "TStyle.h"
0025 #include "TSystem.h"
0026 #include "TTree.h"
0027
0028 bool HasWildcard(const std::string &path)
0029 {
0030 return path.find('*') != std::string::npos || path.find('?') != std::string::npos || path.find('[') != std::string::npos;
0031 }
0032
0033 std::vector<std::string> ResolveInputFiles(const std::string &inputPattern)
0034 {
0035 std::vector<std::string> files;
0036 if (inputPattern.empty())
0037 {
0038 return files;
0039 }
0040
0041 if (!HasWildcard(inputPattern))
0042 {
0043 files.push_back(inputPattern);
0044 return files;
0045 }
0046
0047 glob_t globResult;
0048 const int status = glob(inputPattern.c_str(), 0, nullptr, &globResult);
0049 if (status != 0)
0050 {
0051 globfree(&globResult);
0052 return files;
0053 }
0054
0055 files.reserve(globResult.gl_pathc);
0056 for (std::size_t i = 0; i < globResult.gl_pathc; ++i)
0057 {
0058 files.emplace_back(globResult.gl_pathv[i]);
0059 }
0060 globfree(&globResult);
0061 std::sort(files.begin(), files.end());
0062
0063 return files;
0064 }
0065
0066 std::vector<double> DerivePercentileBoundaries(const std::vector<int> &values, int nPercentileBins = 5)
0067 {
0068 std::vector<double> boundaries;
0069 if (values.empty() || nPercentileBins < 2)
0070 {
0071 return boundaries;
0072 }
0073
0074 std::vector<int> sorted_values = values;
0075 std::sort(sorted_values.begin(), sorted_values.end());
0076
0077 const double last_index = static_cast<double>(sorted_values.size() - 1);
0078 for (int i = 1; i < nPercentileBins; ++i)
0079 {
0080 const double fraction = static_cast<double>(i) / nPercentileBins;
0081 const double position = fraction * last_index;
0082 const auto low_index = static_cast<std::size_t>(position);
0083 const auto high_index = std::min(low_index + 1, sorted_values.size() - 1);
0084 const double weight = position - static_cast<double>(low_index);
0085 const double value = (1.0 - weight) * sorted_values[low_index] + weight * sorted_values[high_index];
0086 boundaries.push_back(value);
0087 }
0088
0089 return boundaries;
0090 }
0091
0092 std::vector<std::string> BuildPercentileBoundaryLabels(int nPercentileBins)
0093 {
0094 std::vector<std::string> labels;
0095 if (nPercentileBins < 2)
0096 {
0097 return labels;
0098 }
0099
0100 const double percentile_width = 100.0 / static_cast<double>(nPercentileBins);
0101 for (int i = 1; i < nPercentileBins; ++i)
0102 {
0103 const double low = 100.0 - (i + 1) * percentile_width;
0104 const double high = 100.0 - i * percentile_width;
0105
0106 auto format_percent = [](double value)
0107 {
0108 std::ostringstream out;
0109 if (std::fabs(value - std::round(value)) < 1e-6)
0110 {
0111 out << static_cast<int>(std::round(value));
0112 }
0113 else
0114 {
0115 out.setf(std::ios::fixed);
0116 out.precision(1);
0117 out << value;
0118 }
0119 return out.str();
0120 };
0121
0122 labels.push_back(format_percent(low) + "-" + format_percent(high) + "%");
0123 }
0124
0125 return labels;
0126 }
0127
0128 std::vector<std::pair<int, int>> BuildPercentileIntervals(const std::vector<double> &boundaries)
0129 {
0130 std::vector<std::pair<int, int>> intervals;
0131 const int nPercentileBins = static_cast<int>(boundaries.size()) + 1;
0132 if (nPercentileBins < 2)
0133 {
0134 return intervals;
0135 }
0136
0137 std::vector<int> integer_boundaries;
0138 integer_boundaries.reserve(boundaries.size());
0139 for (double boundary : boundaries)
0140 {
0141 integer_boundaries.push_back(static_cast<int>(std::ceil(boundary)));
0142 }
0143
0144 intervals.reserve(nPercentileBins);
0145 intervals.emplace_back(integer_boundaries.empty() ? 0 : integer_boundaries.back(), std::numeric_limits<int>::max());
0146 for (int i = static_cast<int>(integer_boundaries.size()) - 1; i > 0; --i)
0147 {
0148 intervals.emplace_back(integer_boundaries[i - 1], integer_boundaries[i]);
0149 }
0150 intervals.emplace_back(0, integer_boundaries.empty() ? std::numeric_limits<int>::max() : integer_boundaries.front());
0151
0152 return intervals;
0153 }
0154
0155 std::vector<std::string> BuildPercentileIntervalLabels(const std::vector<std::pair<int, int>> &intervals)
0156 {
0157 std::vector<std::string> labels;
0158 labels.reserve(intervals.size());
0159 for (const auto &[low, high] : intervals)
0160 {
0161 std::ostringstream out;
0162 out << "[" << low << ", ";
0163 if (high == std::numeric_limits<int>::max())
0164 {
0165 out << "INT_MAX";
0166 }
0167 else
0168 {
0169 out << high;
0170 }
0171 out << ")";
0172 labels.push_back(out.str());
0173 }
0174 return labels;
0175 }
0176
0177 std::vector<double> BuildIntervalLabelPositions(const std::vector<double> &boundaries, double xMin, double xMax)
0178 {
0179 std::vector<double> positions;
0180 if (xMax <= xMin)
0181 {
0182 return positions;
0183 }
0184
0185 const double left_anchor = std::max(xMin, 0.0);
0186 if (boundaries.empty())
0187 {
0188 positions.push_back(0.5 * (left_anchor + xMax));
0189 return positions;
0190 }
0191
0192 positions.reserve(boundaries.size() + 1);
0193 positions.push_back(0.5 * (boundaries.back() + xMax));
0194 for (int i = static_cast<int>(boundaries.size()) - 1; i > 0; --i)
0195 {
0196 positions.push_back(0.5 * (boundaries[i - 1] + boundaries[i]));
0197 }
0198 positions.push_back(0.5 * (left_anchor + boundaries.front()));
0199
0200 return positions;
0201 }
0202
0203 void SavePercentileIntervals(TFile *outfile, const char *tree_name, const std::vector<std::pair<int, int>> &intervals, int nPercentileBins)
0204 {
0205 if (!outfile)
0206 {
0207 return;
0208 }
0209
0210 TTree *tree = new TTree(tree_name, tree_name);
0211 int percentile_bin = -1;
0212 double percentile_low = 0.0;
0213 double percentile_high = 0.0;
0214 int cluster_low = 0;
0215 int cluster_high = 0;
0216
0217 tree->Branch("percentile_bin", &percentile_bin, "percentile_bin/I");
0218 tree->Branch("percentile_low", &percentile_low, "percentile_low/D");
0219 tree->Branch("percentile_high", &percentile_high, "percentile_high/D");
0220 tree->Branch("cluster_low", &cluster_low, "cluster_low/I");
0221 tree->Branch("cluster_high", &cluster_high, "cluster_high/I");
0222
0223 const double percentile_width = 100.0 / static_cast<double>(nPercentileBins);
0224 for (std::size_t i = 0; i < intervals.size(); ++i)
0225 {
0226 percentile_bin = static_cast<int>(i);
0227 percentile_low = static_cast<double>(i) * percentile_width;
0228 percentile_high = static_cast<double>(i + 1) * percentile_width;
0229 cluster_low = intervals[i].first;
0230 cluster_high = intervals[i].second;
0231 tree->Fill();
0232 }
0233
0234 outfile->cd();
0235 tree->Write();
0236 }
0237
0238 void Draw1DHistogramsWithPercentiles(const std::vector<TH1 *> &histograms,
0239 const std::vector<std::vector<double>> &percentile_boundaries,
0240 const std::string plotinfo,
0241 const std::vector<std::string> &hist_labels,
0242 const std::string &x_axis_title,
0243 const std::string &outname,
0244 const std::vector<std::string> &colors = {"#3f90da", "#bd1f01", "#ffa90e", "#6CA651", "#832db6"},
0245 const std::vector<double> &label_boundaries = {},
0246 const std::vector<std::string> &custom_boundary_labels = {},
0247 bool normalize = true,
0248 bool logy = true
0249 )
0250 {
0251 if (histograms.empty())
0252 {
0253 std::cout << "[Draw1DHistogramsWithPercentiles] No histograms provided.\n";
0254 return;
0255 }
0256 if (percentile_boundaries.size() != histograms.size())
0257 {
0258 std::cout << "[Draw1DHistogramsWithPercentiles] Boundary set count does not match histogram count.\n";
0259 return;
0260 }
0261
0262 gStyle->SetOptStat(0);
0263
0264 TCanvas *c = new TCanvas("c_NInttClusterCrossing", "c_NInttClusterCrossing", 900, 700);
0265 c->cd();
0266 gPad->SetLeftMargin(0.14);
0267 gPad->SetRightMargin(0.08);
0268 gPad->SetTopMargin(0.08);
0269 if (logy)
0270 {
0271 c->SetLogy();
0272 }
0273
0274 std::vector<TH1 *> draw_hists;
0275 draw_hists.reserve(histograms.size());
0276
0277 double y_max = 0.0;
0278 double min_positive = 0.0;
0279 for (std::size_t i = 0; i < histograms.size(); ++i)
0280 {
0281 if (!histograms[i])
0282 {
0283 continue;
0284 }
0285
0286 TH1 *hist = static_cast<TH1 *>(histograms[i]->Clone(Form("%s_draw_%zu", histograms[i]->GetName(), i)));
0287 hist->SetDirectory(nullptr);
0288 if (normalize && hist->Integral(0, hist->GetNbinsX() + 1) > 0.0)
0289 {
0290 hist->Scale(1.0 / hist->Integral(0, hist->GetNbinsX() + 1));
0291 }
0292
0293 hist->SetLineColor(TColor::GetColor(colors[i % colors.size()].c_str()));
0294 hist->SetMarkerColor(TColor::GetColor(colors[i % colors.size()].c_str()));
0295 hist->SetLineWidth(2);
0296 hist->GetXaxis()->SetTitle(x_axis_title.c_str());
0297 hist->GetYaxis()->SetTitle(normalize ? "Normalized entries" : "Entries");
0298
0299 y_max = std::max(y_max, hist->GetMaximum());
0300 const double hist_min_positive = hist->GetMinimum(0.0);
0301 if (hist_min_positive > 0.0)
0302 {
0303 min_positive = (min_positive > 0.0) ? std::min(min_positive, hist_min_positive) : hist_min_positive;
0304 }
0305
0306 draw_hists.push_back(hist);
0307 }
0308
0309 if (draw_hists.empty())
0310 {
0311 std::cout << "[Draw1DHistogramsWithPercentiles] All histograms were null.\n";
0312 delete c;
0313 return;
0314 }
0315
0316 draw_hists.front()->GetYaxis()->SetRangeUser(logy ? std::max(min_positive * 0.5, 1e-6) : 0.0, y_max * (logy ? 15.0 : 1.3));
0317 draw_hists.front()->Draw("hist");
0318 for (std::size_t i = 1; i < draw_hists.size(); ++i)
0319 {
0320 draw_hists[i]->Draw("hist same");
0321 }
0322
0323 TLegend *leg = new TLegend((1 - gPad->GetRightMargin()) - 0.3,
0324 1 - gPad->GetTopMargin() - 0.03 - 0.045 * (static_cast<double>(draw_hists.size()) + 1),
0325 (1 - gPad->GetRightMargin()) - 0.15,
0326 1 - gPad->GetTopMargin() - 0.03
0327 );
0328 leg->SetHeader(plotinfo.c_str());
0329 leg->SetFillStyle(0);
0330 leg->SetBorderSize(0);
0331 leg->SetTextSize(0.035);
0332 for (std::size_t i = 0; i < draw_hists.size(); ++i)
0333 {
0334 const std::string legend_label = (i < hist_labels.size() && !hist_labels[i].empty()) ? hist_labels[i] : draw_hists[i]->GetName();
0335 leg->AddEntry(draw_hists[i], legend_label.c_str(), "l");
0336 }
0337 leg->Draw();
0338
0339
0340 for (std::size_t i = 0; i < 1; ++i)
0341 {
0342 const auto &boundaries = percentile_boundaries[i];
0343 const double label_y = y_max * (logy ? (0.60 - 0.12 * static_cast<double>(i)) : (0.85 - 0.08 * static_cast<double>(i)));
0344
0345 for (std::size_t j = 0; j < boundaries.size(); ++j)
0346 {
0347 TLine *line = new TLine(boundaries[j], draw_hists.front()->GetYaxis()->GetXmin(), boundaries[j], label_y);
0348 line->SetLineColor(TColor::GetColor(colors[i % colors.size()].c_str()));
0349 line->SetLineStyle(kDashed);
0350 line->SetLineWidth(2);
0351 line->Draw("same");
0352 }
0353 }
0354
0355 const std::vector<double> &boundaries_for_labels = label_boundaries.empty() ? percentile_boundaries.front() : label_boundaries;
0356 const std::vector<std::string> labels_for_boundaries = custom_boundary_labels.empty() ? BuildPercentileBoundaryLabels(static_cast<int>(boundaries_for_labels.size()) + 1) : custom_boundary_labels;
0357 const std::vector<double> label_x_positions = (labels_for_boundaries.size() == boundaries_for_labels.size() + 1) ? BuildIntervalLabelPositions(boundaries_for_labels, draw_hists.front()->GetXaxis()->GetXmin(), draw_hists.front()->GetXaxis()->GetXmax()) : boundaries_for_labels;
0358 const double label_y = y_max * (logy ? 0.24 : 0.72);
0359
0360 for (std::size_t i = 0; i < label_x_positions.size() && i < labels_for_boundaries.size(); ++i)
0361 {
0362 TLatex *label = new TLatex(label_x_positions[i], label_y, labels_for_boundaries[i].c_str());
0363 label->SetTextAngle(90);
0364 label->SetTextAlign(13);
0365 label->SetTextSize(0.025);
0366 label->SetTextColor(TColor::GetColor(colors[0].c_str()));
0367 label->Draw("same");
0368 }
0369
0370 c->SaveAs(Form("%s.pdf", outname.c_str()));
0371 c->SaveAs(Form("%s.png", outname.c_str()));
0372
0373 for (TH1 *hist : draw_hists)
0374 {
0375 delete hist;
0376 }
0377 delete leg;
0378 delete c;
0379 }
0380
0381
0382
0383
0384
0385 void NInttClusterCrossing(const std::string inputfile = "/sphenix/tg/tg01/hf/hjheng/ppg-dNdEta-OOpp/simulation-nopileup-ntuple/ACTS/ntuple_*.root",
0386 const bool isMC = true,
0387 int nPercentileBins = 5
0388 )
0389 {
0390 if (inputfile.empty())
0391 {
0392 std::cout << "[NInttClusterCrossing] Please provide an input ROOT file.\n";
0393 return;
0394 }
0395
0396 std::string plotdir = "./figure/figure-NInttClusterCrossing";
0397 system(("mkdir -p " + plotdir).c_str());
0398 const std::string outputPrefix = isMC ? "MC_" : "Data_";
0399
0400 std::vector<int> NInttCluster_TrgCrossing;
0401 std::vector<int> NInttCluster_NonTrgCrossing;
0402
0403 TH1 *hM_NInttCluster_TrgCrossing = new TH1D("hM_NInttCluster_TrgCrossing", "N_{INTT} clusters in trigger crossing;N_{INTT} clusters;Counts", 200, -0.5, 999.5);
0404 TH1 *hM_NInttCluster_NonTrgCrossing = new TH1D("hM_NInttCluster_NonTrgCrossing", "N_{INTT} clusters in non-trigger crossing;N_{INTT} clusters;Counts", 200, -0.5, 999.5);
0405
0406 TChain *t = new TChain("VTX");
0407 const auto inputFiles = ResolveInputFiles(inputfile);
0408 if (inputFiles.empty())
0409 {
0410 std::cout << "[NInttClusterCrossing] No files matched: " << inputfile << "\n";
0411 delete t;
0412 return;
0413 }
0414
0415 const std::size_t maxFiles = isMC ? 100 : inputFiles.size();
0416 int nAddedFiles = 0;
0417 for (std::size_t i = 0; i < inputFiles.size() && i < maxFiles; ++i)
0418 {
0419 if (t->Add(inputFiles[i].c_str()) > 0)
0420 {
0421 ++nAddedFiles;
0422 }
0423 }
0424
0425 if (isMC && inputFiles.size() > maxFiles)
0426 {
0427 std::cout << "[NInttClusterCrossing] MC mode: using first " << maxFiles << " files out of " << inputFiles.size() << " matched files.\n";
0428 }
0429
0430 if (nAddedFiles == 0)
0431 {
0432 std::cout << "[NInttClusterCrossing] Failed to add input files to TChain.\n";
0433 delete t;
0434 return;
0435 }
0436
0437 std::cout << "[NInttClusterCrossing] Added " << nAddedFiles << " files to TChain.\n";
0438
0439 int counter = 0;
0440 std::vector<int> *firedTriggers = nullptr;
0441 std::vector<unsigned int> *cluster_layer = nullptr;
0442 std::vector<int> *cluster_timeBucketID = nullptr;
0443 t->SetBranchAddress("counter", &counter);
0444 t->SetBranchAddress("firedTriggers", &firedTriggers);
0445 t->SetBranchAddress("cluster_layer", &cluster_layer);
0446 t->SetBranchAddress("cluster_timeBucketID", &cluster_timeBucketID);
0447
0448 const Long64_t nEntries = t->GetEntries();
0449 for (Long64_t i = 0; i < nEntries; ++i)
0450 {
0451 t->GetEntry(i);
0452
0453 std::set<int> unique_time_buckets(cluster_timeBucketID->begin(), cluster_timeBucketID->end());
0454 for (int time_bucket : unique_time_buckets)
0455 {
0456 int n_clusters_in_bucket = 0;
0457 for (std::size_t j = 0; j < cluster_timeBucketID->size(); ++j)
0458 {
0459 if (cluster_timeBucketID->at(j) == time_bucket)
0460 {
0461 ++n_clusters_in_bucket;
0462 }
0463 }
0464
0465 if (time_bucket == 0)
0466 {
0467 hM_NInttCluster_TrgCrossing->Fill(n_clusters_in_bucket);
0468 NInttCluster_TrgCrossing.push_back(n_clusters_in_bucket);
0469 }
0470 else
0471 {
0472 hM_NInttCluster_NonTrgCrossing->Fill(n_clusters_in_bucket);
0473 NInttCluster_NonTrgCrossing.push_back(n_clusters_in_bucket);
0474 }
0475 }
0476 }
0477
0478 const auto trg_percentile_boundaries = DerivePercentileBoundaries(NInttCluster_TrgCrossing, nPercentileBins);
0479 const auto nontrg_percentile_boundaries = DerivePercentileBoundaries(NInttCluster_NonTrgCrossing, nPercentileBins);
0480 const auto trg_percentile_intervals = BuildPercentileIntervals(trg_percentile_boundaries);
0481 const auto trg_percentile_interval_labels = BuildPercentileIntervalLabels(trg_percentile_intervals);
0482
0483 std::cout << "[NInttClusterCrossing] Trigger crossing percentile boundaries:";
0484 for (double boundary : trg_percentile_boundaries)
0485 {
0486 std::cout << " " << boundary;
0487 }
0488 std::cout << "\n";
0489
0490 std::cout << "[NInttClusterCrossing] Non-trigger crossing percentile boundaries:";
0491 for (double boundary : nontrg_percentile_boundaries)
0492 {
0493 std::cout << " " << boundary;
0494 }
0495 std::cout << "\n";
0496
0497 std::cout << "[NInttClusterCrossing] Trigger crossing percentile intervals:\n";
0498 for (std::size_t i = 0; i < trg_percentile_interval_labels.size(); ++i)
0499 {
0500 std::cout << " bin " << i << ": " << trg_percentile_interval_labels[i] << "\n";
0501 }
0502
0503 Draw1DHistogramsWithPercentiles({hM_NInttCluster_TrgCrossing, hM_NInttCluster_NonTrgCrossing},
0504 {trg_percentile_boundaries, nontrg_percentile_boundaries},
0505 (isMC ? "Simulation" : "Data, O+O 82405"),
0506 {"Trigger crossing", "Non-trigger crossing"},
0507 "N_{INTT clusters}",
0508 Form("%s/%shM_NInttCluster_CrossingComparison_wPercentile", plotdir.c_str(), outputPrefix.c_str()),
0509 {"#3f90da", "#bd1f01"},
0510 trg_percentile_boundaries, BuildPercentileBoundaryLabels(nPercentileBins)
0511 );
0512
0513 TFile *fout = TFile::Open(Form("%s/%sNInttClusterPercentileBoundaries.root", plotdir.c_str(), outputPrefix.c_str()), "RECREATE");
0514 SavePercentileIntervals(fout, "trigger_percentile_intervals", trg_percentile_intervals, nPercentileBins);
0515 fout->Close();
0516 delete fout;
0517
0518 delete t;
0519 delete hM_NInttCluster_TrgCrossing;
0520 delete hM_NInttCluster_NonTrgCrossing;
0521 }
