File indexing completed on 2026-05-23 08:08:42
0001 #include <algorithm>
0002 #include <cctype>
0003 #include <cmath>
0004 #include <cstdlib>
0005 #include <iostream>
0006 #include <iomanip>
0007 #include <limits>
0008 #include <map>
0009 #include <set>
0010 #include <sstream>
0011 #include <string>
0012 #include <vector>
0013
0014 #include "TCanvas.h"
0015 #include "TColor.h"
0016 #include "TFile.h"
0017 #include "TGaxis.h"
0018 #include "TH1D.h"
0019 #include "TH2D.h"
0020 #include "THStack.h"
0021 #include "TKey.h"
0022 #include "TLatex.h"
0023 #include "TLegend.h"
0024 #include "TLine.h"
0025 #include "TPad.h"
0026 #include "TString.h"
0027 #include "TStyle.h"
0028 #include "TSystem.h"
0029 #include "TTree.h"
0030
0031
0032
0033
0034 void CanvasPartition(TCanvas *c, int Nx, int Ny, float lMargin, float rMargin, float bMargin, float tMargin)
0035 {
0036 if (!c)
0037 return;
0038
0039 const float hStep = (1.f - lMargin - rMargin) / Nx;
0040 const float vStep = (1.f - bMargin - tMargin) / Ny;
0041
0042 for (int i = 0; i < Nx; ++i)
0043 {
0044 float hposl = (i == 0) ? 0.f : lMargin + i * hStep;
0045 float hposr = (i == Nx - 1) ? 1.f : lMargin + (i + 1) * hStep;
0046 float hmarl = (i == 0) ? lMargin / (hposr - hposl) : 0.f;
0047 float hmarr = (i == Nx - 1) ? rMargin / (hposr - hposl) : 0.f;
0048
0049 for (int j = 0; j < Ny; ++j)
0050 {
0051 float vposd = (j == 0) ? 0.f : bMargin + j * vStep;
0052 float vposu = (j == Ny - 1) ? 1.f : bMargin + (j + 1) * vStep;
0053 float vmard = (j == 0) ? bMargin / (vposu - vposd) : 0.f;
0054 float vmaru = (j == Ny - 1) ? tMargin / (vposu - vposd) : 0.f;
0055
0056 c->cd(0);
0057 TString name = TString::Format("pad_%d_%d", i, j);
0058 delete c->FindObject(name.Data());
0059
0060 auto *pad = new TPad(name, "", hposl, vposd, hposr, vposu);
0061 pad->SetLeftMargin(hmarl);
0062 pad->SetRightMargin(hmarr);
0063 pad->SetBottomMargin(vmard);
0064 pad->SetTopMargin(vmaru);
0065 pad->SetBorderMode(0);
0066 pad->SetBorderSize(0);
0067 pad->SetFrameBorderMode(0);
0068 pad->Draw();
0069 }
0070 }
0071 }
0072
0073
0074
0075
0076 static double XtoPad(double x)
0077 {
0078 double lm = gPad->GetLeftMargin(), rm = gPad->GetRightMargin();
0079 return lm + x * (1. - lm - rm);
0080 }
0081 static double YtoPad(double y)
0082 {
0083 double bm = gPad->GetBottomMargin(), tm = gPad->GetTopMargin();
0084 return bm + y * (1. - bm - tm);
0085 }
0086
0087 template <typename T> static std::string to_string_with_precision(const T a_value, const int n = 2)
0088 {
0089 std::ostringstream out;
0090 out.precision(n);
0091 out << std::fixed << a_value;
0092 return out.str();
0093 }
0094
0095 static std::string EtaPhiRangeNameSuffix(const TH2D *h2, int ieta, int iphi)
0096 {
0097 const auto *etaAx = h2->GetXaxis();
0098 const auto *phiAx = h2->GetYaxis();
0099 const double etaLow = etaAx->GetBinLowEdge(ieta + 1);
0100 const double etaHigh = etaAx->GetBinUpEdge(ieta + 1);
0101 const double phiLow = phiAx->GetBinLowEdge(iphi + 1);
0102 const double phiHigh = phiAx->GetBinUpEdge(iphi + 1);
0103
0104 std::ostringstream nm;
0105 nm << "eta" << to_string_with_precision(etaLow) << "-" << to_string_with_precision(etaHigh) << "_phi" << to_string_with_precision(phiLow) << "-" << to_string_with_precision(phiHigh);
0106 return nm.str();
0107 }
0108
0109
0110
0111
0112 static TH1D *GetCellHist(TFile *fin, const char *dirname, const std::string &prefix, const std::string &mode, const std::string &tag, int ieta, int iphi, const TH2D *hGrid)
0113 {
0114 if (!fin || !hGrid)
0115 return nullptr;
0116
0117 std::vector<std::string> dir_candidates;
0118 dir_candidates.push_back(dirname);
0119 dir_candidates.push_back(prefix + "_" + dirname);
0120
0121 for (const auto &dir_name : dir_candidates)
0122 {
0123 auto *dir = fin->GetDirectory(dir_name.c_str());
0124 if (!dir)
0125 continue;
0126
0127 auto tryGet = [&](const std::string &m) -> TH1D *
0128 {
0129 std::string n = prefix + "_" + m + "_" + tag + "_" + EtaPhiRangeNameSuffix(hGrid, ieta, iphi);
0130 return dynamic_cast<TH1D *>(dir->Get(n.c_str()));
0131 };
0132
0133 if (auto *h = tryGet(mode))
0134 return h;
0135 if (mode != "dR")
0136 if (auto *h = tryGet("dR"))
0137 return h;
0138 }
0139 return nullptr;
0140 }
0141
0142 static TH2D *GetPrimaryGrid(TFile *fin, const std::string &prefix)
0143 {
0144 if (!fin)
0145 return nullptr;
0146 for (const char *suffix : {"_Sig", "_Bkg", "_Sub"})
0147 {
0148 const std::string hname = prefix + suffix;
0149 if (auto *h = dynamic_cast<TH2D *>(fin->Get(hname.c_str())))
0150 {
0151 return h;
0152 }
0153 }
0154 return nullptr;
0155 }
0156 static TH2D *GetRawCountGrid(TFile *fin, const std::string &prefix)
0157 {
0158 if (!fin)
0159 return nullptr;
0160 return dynamic_cast<TH2D *>(fin->Get((prefix + "_Sub_AbsDphiCount").c_str()));
0161 }
0162
0163 static constexpr double kDirectRawCountCut = 0.015;
0164
0165
0166
0167
0168
0169
0170
0171 static double IntegralInRange(TH1D *h, double xmin, double xmax)
0172 {
0173 if (!h || xmax <= xmin)
0174 return 0.;
0175 int b1 = h->GetXaxis()->FindBin(xmin + 1e-9);
0176 int b2 = h->GetXaxis()->FindBin(xmax - 1e-9);
0177 return (b2 >= b1) ? h->Integral(b1, b2) : 0.;
0178 }
0179
0180 static TH1D *MakeFractionalHist(TH1D *h, TH1D *hRef, const std::string &newname)
0181 {
0182 if (!h || !hRef)
0183 return nullptr;
0184
0185 TH1D *hFrac = dynamic_cast<TH1D *>(h->Clone(newname.c_str()));
0186 if (!hFrac)
0187 return nullptr;
0188 hFrac->SetDirectory(nullptr);
0189
0190 const int nbins = hFrac->GetNbinsX();
0191 for (int i = 1; i <= nbins; ++i)
0192 {
0193 const double denom = hRef->GetBinContent(i);
0194 if (denom > 0.0)
0195 {
0196 hFrac->SetBinContent(i, h->GetBinContent(i) / denom);
0197 hFrac->SetBinError(i, h->GetBinError(i) / denom);
0198 }
0199 else
0200 {
0201 hFrac->SetBinContent(i, 0.0);
0202 hFrac->SetBinError(i, 0.0);
0203 }
0204 }
0205 return hFrac;
0206 }
0207
0208 static TH1D *MakeInverseHist(TH1D *h, const std::string &newname)
0209 {
0210 if (!h)
0211 return nullptr;
0212
0213 TH1D *hInv = dynamic_cast<TH1D *>(h->Clone(newname.c_str()));
0214 if (!hInv)
0215 return nullptr;
0216 hInv->SetDirectory(nullptr);
0217
0218 const int nbins = hInv->GetNbinsX();
0219 for (int i = 1; i <= nbins; ++i)
0220 {
0221 const double value = h->GetBinContent(i);
0222 const double error = h->GetBinError(i);
0223 if (value != 0.0)
0224 {
0225 hInv->SetBinContent(i, 1.0 / value);
0226 hInv->SetBinError(i, error / (value * value));
0227 }
0228 else
0229 {
0230 hInv->SetBinContent(i, 0.0);
0231 hInv->SetBinError(i, 0.0);
0232 }
0233 }
0234 return hInv;
0235 }
0236
0237 static void draw2Dhistogram_adjZaxis(TH2 *hist,
0238 bool logz,
0239 const std::string &xtitle,
0240 const std::string &ytitle,
0241 const std::string &ztitle,
0242 std::pair<double, double> zrange,
0243 const std::vector<std::string> &addinfo,
0244 const std::string &drawoption = "colz text",
0245 const std::string &paintTextFmt = ".3g",
0246 bool show_bin_value_with_error = true,
0247 bool show_bin_value_only = false,
0248 const std::string &filename = "plot2d")
0249 {
0250 if (!hist)
0251 return;
0252
0253 gStyle->SetPalette(kLightTemperature);
0254
0255
0256
0257
0258 TCanvas *c = new TCanvas("c_2d", "c_2d", 900, 800);
0259 c->cd();
0260 gPad->SetRightMargin(0.21);
0261 c->SetLogz(logz);
0262
0263 gStyle->SetPaintTextFormat(paintTextFmt.c_str());
0264 hist->GetXaxis()->SetTitle(xtitle.c_str());
0265 hist->GetYaxis()->SetTitle(ytitle.c_str());
0266 hist->GetZaxis()->SetTitle(ztitle.c_str());
0267 hist->GetZaxis()->SetTitleOffset(1.5);
0268 if (zrange.first < zrange.second)
0269 hist->GetZaxis()->SetRangeUser(zrange.first, zrange.second);
0270
0271 std::string drawopt = drawoption;
0272 if (show_bin_value_with_error)
0273 {
0274 const std::vector<std::string> tokens_to_strip = {"TEXT", "text"};
0275 for (const auto &tok : tokens_to_strip)
0276 {
0277 size_t pos = std::string::npos;
0278 while ((pos = drawopt.find(tok)) != std::string::npos)
0279 drawopt.erase(pos, tok.size());
0280 }
0281 if (drawopt.find("colz") == std::string::npos && drawopt.find("COLZ") == std::string::npos)
0282 drawopt += " colz";
0283 }
0284 else if (show_bin_value_only)
0285 {
0286 const std::vector<std::string> tokens_to_strip = {"TEXT", "text"};
0287 for (const auto &tok : tokens_to_strip)
0288 {
0289 size_t pos = std::string::npos;
0290 while ((pos = drawopt.find(tok)) != std::string::npos)
0291 drawopt.erase(pos, tok.size());
0292 }
0293 if (drawopt.find("colz") == std::string::npos && drawopt.find("COLZ") == std::string::npos)
0294 drawopt += " colz";
0295 }
0296 hist->Draw(drawopt.c_str());
0297
0298 if (show_bin_value_with_error)
0299 {
0300 TLatex tbin_val;
0301 tbin_val.SetTextAlign(22);
0302 tbin_val.SetTextFont(42);
0303 tbin_val.SetTextSize(0.018);
0304
0305 TLatex tbin_err;
0306 tbin_err.SetTextAlign(22);
0307 tbin_err.SetTextFont(42);
0308 tbin_err.SetTextSize(0.015);
0309
0310
0311 for (int ibx = 1; ibx <= hist->GetNbinsX(); ++ibx)
0312 {
0313 for (int iby = 1; iby <= hist->GetNbinsY(); ++iby)
0314 {
0315 const double val = hist->GetBinContent(ibx, iby);
0316 const double err = hist->GetBinError(ibx, iby);
0317 const double x = hist->GetXaxis()->GetBinCenter(ibx);
0318 const double y = hist->GetYaxis()->GetBinCenter(iby);
0319 const double y_offset = 0.18 * hist->GetYaxis()->GetBinWidth(iby);
0320
0321 tbin_val.DrawLatex(x, y + y_offset, Form("%.2f", val));
0322 tbin_err.DrawLatex(x, y - y_offset, Form("#pm%.3g", err));
0323 }
0324 }
0325 }
0326 else if (show_bin_value_only)
0327 {
0328 TLatex tbin_val;
0329 tbin_val.SetTextAlign(22);
0330 tbin_val.SetTextFont(42);
0331 tbin_val.SetTextSize(0.014);
0332
0333 for (int ibx = 1; ibx <= hist->GetNbinsX(); ++ibx)
0334 {
0335 for (int iby = 1; iby <= hist->GetNbinsY(); ++iby)
0336 {
0337 const double val = hist->GetBinContent(ibx, iby);
0338 const double x = hist->GetXaxis()->GetBinCenter(ibx);
0339 const double y = hist->GetYaxis()->GetBinCenter(iby);
0340 tbin_val.DrawLatex(x, y, Form("%.3g", val));
0341 }
0342 }
0343 }
0344
0345
0346 TLatex latex;
0347 latex.SetNDC();
0348 latex.SetTextFont(42);
0349 latex.SetTextSize(0.035);
0350 for (size_t i = 0; i < addinfo.size(); ++i)
0351 {
0352 latex.DrawLatex(gPad->GetLeftMargin(), (1 - gPad->GetTopMargin() + 0.01) - 0.045 * i, addinfo[i].c_str());
0353 }
0354
0355 c->SaveAs((filename + ".png").c_str());
0356 c->SaveAs((filename + ".pdf").c_str());
0357 delete c;
0358 }
0359
0360 static void draw1Dhistogram(TH1 *hist,
0361 const bool logy,
0362 const std::string &xtitle,
0363 const std::string &ytitle,
0364 const std::vector<std::string> &addinfo,
0365 const std::string &drawoption = "hist e1",
0366 const std::string &filename = "plot1d")
0367 {
0368 if (!hist)
0369 return;
0370
0371 TCanvas *c = new TCanvas("c_1d", "c_1d", 900, 800);
0372 c->SetLogy(logy);
0373 c->cd();
0374 gPad->SetLeftMargin(0.14);
0375 gPad->SetRightMargin(0.05);
0376
0377 gPad->SetTopMargin(0.08);
0378
0379 hist->SetTitle("");
0380 hist->GetXaxis()->SetTitle(xtitle.c_str());
0381 hist->GetYaxis()->SetTitle(ytitle.c_str());
0382 hist->GetYaxis()->SetTitleOffset(1.4);
0383 hist->GetYaxis()->SetRangeUser((logy ? hist->GetMinimum(0) * 0.5 : 0.0), hist->GetMaximum() * (logy ? 10.0 : 1.4));
0384 hist->SetLineWidth(2);
0385 hist->SetLineColor(kBlack);
0386 hist->SetMarkerStyle(kFullCircle);
0387 hist->SetMarkerSize(1.0);
0388 hist->SetMarkerColor(kBlack);
0389 hist->Draw(drawoption.c_str());
0390
0391 TLatex latex;
0392 latex.SetNDC();
0393 latex.SetTextFont(42);
0394 latex.SetTextSize(0.035);
0395
0396 latex.SetTextAlign(kHAlignRight + kVAlignBottom);
0397 for (size_t i = 0; i < addinfo.size(); ++i)
0398 {
0399 latex.DrawLatex(1 - gPad->GetRightMargin(), (1 - gPad->GetTopMargin() + 0.01) - 0.045 * i, addinfo[i].c_str());
0400 }
0401
0402 c->SaveAs((filename + ".png").c_str());
0403 c->SaveAs((filename + ".pdf").c_str());
0404 delete c;
0405 }
0406
0407 static void drawComparisonWithRatio(std::pair<TH1D *, std::string> hDenInput,
0408 std::pair<TH1D *, std::string> hNumInput,
0409 const bool logy,
0410 const std::string &xtitle,
0411 const std::string &ytitle,
0412 const std::string &ratio_ytitle,
0413 const std::vector<std::string> &addinfo,
0414 const std::string &filename)
0415 {
0416 if (!hDenInput.first || !hNumInput.first)
0417 return;
0418
0419 TH1D *hDen = dynamic_cast<TH1D *>(hDenInput.first->Clone((filename + "_den").c_str()));
0420 TH1D *hNum = dynamic_cast<TH1D *>(hNumInput.first->Clone((filename + "_num").c_str()));
0421 TH1D *hRatio = dynamic_cast<TH1D *>(hNumInput.first->Clone((filename + "_ratio").c_str()));
0422 if (!hDen || !hNum || !hRatio)
0423 {
0424 delete hDen;
0425 delete hNum;
0426 delete hRatio;
0427 return;
0428 }
0429
0430 hDen->SetDirectory(nullptr);
0431 hNum->SetDirectory(nullptr);
0432 hRatio->SetDirectory(nullptr);
0433 hRatio->Reset("ICES");
0434 hRatio->Divide(hNum, hDen, 1.0, 1.0, "B");
0435
0436 hDen->SetLineColor(kBlack);
0437 hDen->SetMarkerColor(kBlack);
0438 hDen->SetMarkerStyle(kFullCircle);
0439 hDen->SetMarkerSize(0.8);
0440 hDen->SetLineWidth(2);
0441
0442 hNum->SetLineColor(kRed + 1);
0443 hNum->SetMarkerColor(kRed + 1);
0444 hNum->SetMarkerStyle(kOpenCircle);
0445 hNum->SetMarkerSize(0.8);
0446 hNum->SetLineWidth(2);
0447
0448 hRatio->SetLineColor(kRed + 1);
0449 hRatio->SetMarkerColor(kRed + 1);
0450 hRatio->SetMarkerStyle(kOpenCircle);
0451 hRatio->SetMarkerSize(0.8);
0452 hRatio->SetLineWidth(2);
0453
0454 double ymax = std::max(hDen->GetMaximum(), hNum->GetMaximum());
0455 if (ymax <= 0.0)
0456 ymax = 1.0;
0457
0458 double ymin = std::min(hDen->GetMinimum(0), hNum->GetMinimum(0));
0459
0460 TCanvas *c = new TCanvas((filename + "_canvas").c_str(), (filename + "_canvas").c_str(), 900, 800);
0461 c->cd();
0462 TPad *pad_top = new TPad("pad_top", "pad_top", 0.0, 0.32, 1.0, 1.0);
0463 TPad *pad_bottom = new TPad("pad_bottom", "pad_bottom", 0.0, 0.0, 1.0, 0.32);
0464 pad_top->SetLeftMargin(0.14);
0465 pad_top->SetRightMargin(0.05);
0466 pad_top->SetTopMargin(0.05);
0467 pad_top->SetBottomMargin(0.025);
0468 pad_bottom->SetLeftMargin(0.14);
0469 pad_bottom->SetRightMargin(0.05);
0470 pad_bottom->SetTopMargin(0.025);
0471 pad_bottom->SetBottomMargin(0.35);
0472 pad_top->Draw();
0473 pad_bottom->Draw();
0474
0475 float toptitlesize = 0.05;
0476 pad_top->cd();
0477 pad_top->SetLogy(logy);
0478 hDen->SetTitle("");
0479 hDen->GetYaxis()->SetTitle(ytitle.c_str());
0480 hDen->GetYaxis()->SetTitleOffset(1.35);
0481 hDen->GetYaxis()->SetTitleSize(toptitlesize);
0482 hDen->GetYaxis()->SetLabelSize(toptitlesize);
0483 hDen->GetXaxis()->SetLabelSize(0.0);
0484 hDen->GetXaxis()->SetTitleSize(0.0);
0485
0486 hDen->SetMaximum((logy ? 100.0 : 1.5) * ymax);
0487 hDen->Draw("e1");
0488 hNum->Draw("e1 same");
0489
0490 TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.03,
0491 1 - gPad->GetTopMargin() - 0.2,
0492 gPad->GetLeftMargin() + 0.15,
0493 1 - gPad->GetTopMargin() - 0.05
0494 );
0495 leg->SetHeader(addinfo[0].c_str(), "L");
0496 leg->SetBorderSize(0);
0497 leg->SetFillStyle(0);
0498 leg->SetTextSize(0.04);
0499 leg->AddEntry(hDen, hDenInput.second.c_str(), "pel");
0500 leg->AddEntry(hNum, hNumInput.second.c_str(), "pel");
0501 leg->Draw();
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513 float textsize_scalefactor = 2.1;
0514 pad_bottom->cd();
0515 hRatio->SetTitle("");
0516 hRatio->GetXaxis()->SetTitle(xtitle.c_str());
0517 hRatio->GetYaxis()->SetTitle(ratio_ytitle.c_str());
0518 hRatio->GetXaxis()->SetTitleSize(toptitlesize * textsize_scalefactor);
0519 hRatio->GetXaxis()->SetLabelSize(toptitlesize * textsize_scalefactor);
0520
0521 hRatio->GetYaxis()->SetTitleSize(toptitlesize * textsize_scalefactor);
0522 hRatio->GetYaxis()->SetLabelSize(toptitlesize * textsize_scalefactor);
0523 hRatio->GetYaxis()->SetTitleOffset(0.6);
0524 hRatio->GetYaxis()->SetNdivisions(505);
0525
0526 hRatio->SetMaximum(1);
0527 hRatio->Draw("e1");
0528
0529 TLine *line = new TLine(hRatio->GetXaxis()->GetXmin(), 1.0, hRatio->GetXaxis()->GetXmax(), 1.0);
0530 line->SetLineColor(kBlack);
0531 line->SetLineStyle(2);
0532 line->SetLineWidth(2);
0533 line->Draw("same");
0534
0535 c->SaveAs((filename + ".png").c_str());
0536 c->SaveAs((filename + ".pdf").c_str());
0537
0538 delete line;
0539 delete leg;
0540 delete hRatio;
0541 delete hNum;
0542 delete hDen;
0543 delete c;
0544 }
0545
0546
0547 static void draw_comparison(std::vector<TH1 *> histograms,
0548 const std::vector<std::string> &colors,
0549 std::vector<std::string> labels,
0550 const std::string &xAxisTitle,
0551 const std::string &yAxisTitle,
0552 const std::pair<float, float> &yRange,
0553 const bool logy,
0554 const bool normalize,
0555 const int xAxisNdivisions,
0556 const std::pair<bool, float> &referenceline,
0557 const std::string &outputFileName
0558 )
0559 {
0560 if (histograms.empty())
0561 return;
0562
0563 if (normalize)
0564 {
0565 for (auto *hist : histograms)
0566 {
0567 if (hist && hist->Integral(-1, -1) > 0.0)
0568 hist->Scale(1. / hist->Integral(-1, -1));
0569 }
0570 }
0571
0572 float yMin = yRange.first, yMax = yRange.second;
0573 if (yRange.first < 0)
0574 {
0575 float histMin = 1E6;
0576 for (auto *hist : histograms)
0577 {
0578 if (hist && hist->GetMinimum(0) < histMin)
0579 histMin = hist->GetMinimum(0);
0580 }
0581 yMin = logy ? histMin * 0.5f : 0.f;
0582 }
0583
0584 if (yRange.second < 0)
0585 {
0586 float histMax = -1E6;
0587 for (auto *hist : histograms)
0588 {
0589 if (hist && hist->GetMaximum() > histMax)
0590 histMax = hist->GetMaximum();
0591 }
0592 yMax = logy ? histMax * 250.f : histMax * 1.5f;
0593 }
0594
0595 TCanvas *c = new TCanvas("c_comparison", "Comparison", 800, 700);
0596 c->SetLogy(logy);
0597 c->cd();
0598 gPad->SetLeftMargin(0.16);
0599 gPad->SetRightMargin(0.05);
0600 gPad->SetTopMargin(0.07);
0601 for (size_t i = 0; i < histograms.size(); ++i)
0602 {
0603 if (!histograms[i])
0604 continue;
0605 histograms[i]->SetLineColor(TColor::GetColor(colors[i].c_str()));
0606 histograms[i]->SetMarkerColor(TColor::GetColor(colors[i].c_str()));
0607 histograms[i]->SetMarkerStyle(20 + static_cast<int>(i));
0608 histograms[i]->SetMarkerSize(1.0);
0609 histograms[i]->SetLineWidth(2);
0610 if (i == 0)
0611 {
0612 histograms[i]->SetTitle("");
0613 histograms[i]->GetXaxis()->SetTitle(xAxisTitle.c_str());
0614 histograms[i]->GetXaxis()->SetTitleOffset(1.4);
0615 histograms[i]->GetXaxis()->SetNdivisions(xAxisNdivisions);
0616 histograms[i]->GetYaxis()->SetTitle(yAxisTitle.c_str());
0617 histograms[i]->GetYaxis()->SetRangeUser(yMin, yMax);
0618 histograms[i]->GetYaxis()->SetTitleOffset(1.5);
0619 histograms[i]->Draw("HIST E1");
0620 }
0621 else
0622 {
0623 histograms[i]->Draw("HIST E1 SAME");
0624 }
0625 }
0626
0627
0628 if (referenceline.first)
0629 {
0630 TLine *line = new TLine(histograms[0]->GetXaxis()->GetXmin(), referenceline.second, histograms[0]->GetXaxis()->GetXmax(), referenceline.second);
0631 line->SetLineColor(kGray + 1);
0632 line->SetLineStyle(2);
0633 line->SetLineWidth(2);
0634 line->Draw("same");
0635 }
0636
0637 TLegend *legend = new TLegend(gPad->GetLeftMargin() + 0.04,
0638 (1 - gPad->GetTopMargin()) - 0.035 - 0.04 * (histograms.size() + 1),
0639 gPad->GetLeftMargin() + 0.14,
0640 (1 - gPad->GetTopMargin()) - 0.035
0641 );
0642 legend->SetTextAlign(kHAlignLeft + kVAlignCenter);
0643 legend->SetTextSize(0.035);
0644 legend->SetFillStyle(0);
0645 legend->SetBorderSize(0);
0646 for (size_t i = 0; i < labels.size() && i < histograms.size(); ++i)
0647 {
0648 if (histograms[i])
0649 legend->AddEntry(histograms[i], labels[i].c_str(), "pel");
0650 }
0651 legend->Draw();
0652
0653 c->SaveAs(Form("%s.png", outputFileName.c_str()));
0654 c->SaveAs(Form("%s.pdf", outputFileName.c_str()));
0655 delete legend;
0656 delete c;
0657 }
0658
0659
0660
0661 static void draw_comparison_alt(std::vector<TH1 *> histograms,
0662 const std::vector<std::string> &colors,
0663 const std::vector<int> &markers,
0664 const std::vector<float> &markerSizes,
0665 std::vector<std::string> labels,
0666 const std::string &xAxisTitle,
0667 const std::string &yAxisTitle,
0668 const std::pair<float, float> &yRange,
0669 const bool logy,
0670 const bool normalize,
0671 const int xAxisNdivisions,
0672 const std::pair<bool, float> &referenceline,
0673 const std::string &outputFileName
0674 )
0675 {
0676
0677 TH1 *hDummy = new TH1D("hDummy", "", 1, histograms[0]->GetXaxis()->GetXmin(), histograms[0]->GetXaxis()->GetXmax());
0678
0679 if (histograms.empty())
0680 return;
0681
0682 if (normalize)
0683 {
0684 for (auto *hist : histograms)
0685 {
0686 if (hist && hist->Integral(-1, -1) > 0.0)
0687 hist->Scale(1. / hist->Integral(-1, -1));
0688 }
0689 }
0690
0691 float yMin = yRange.first, yMax = yRange.second;
0692 if (yRange.first < 0)
0693 {
0694 float histMin = std::numeric_limits<float>::max();
0695 for (auto *hist : histograms)
0696 {
0697 std::cout << "hist min: " << (hist ? hist->GetMinimum(0) : -1) << std::endl;
0698 if (hist && hist->GetMinimum(0) < histMin)
0699 histMin = hist->GetMinimum(0);
0700 }
0701 yMin = logy ? histMin * 0.5f : 0.f;
0702 }
0703
0704 if (yRange.second < 0)
0705 {
0706 float histMax = -std::numeric_limits<float>::max();
0707 for (auto *hist : histograms)
0708 {
0709 if (hist && hist->GetMaximum() > histMax)
0710 histMax = hist->GetMaximum();
0711 }
0712 yMax = logy ? histMax * 2.f : histMax * 1.5f;
0713 }
0714
0715 std::cout << "yMin: " << yMin << ", yMax: " << yMax << std::endl;
0716
0717 TCanvas *c = new TCanvas("c_comparison", "Comparison", 800, 700);
0718 c->SetLogy(logy);
0719 c->cd();
0720 gPad->SetLeftMargin(0.16);
0721 gPad->SetRightMargin(0.05);
0722 gPad->SetTopMargin(0.07);
0723
0724 hDummy->SetTitle("");
0725 hDummy->GetXaxis()->SetTitle(xAxisTitle.c_str());
0726 hDummy->GetXaxis()->SetTitleOffset(1.4);
0727 hDummy->GetXaxis()->SetNdivisions(xAxisNdivisions);
0728 hDummy->GetYaxis()->SetTitle(yAxisTitle.c_str());
0729 hDummy->GetYaxis()->SetRangeUser(yMin, yMax);
0730 hDummy->GetYaxis()->SetTitleOffset(1.5);
0731 hDummy->Draw("AXIS");
0732
0733 for (size_t i = 0; i < histograms.size(); ++i)
0734 {
0735 if (!histograms[i])
0736 continue;
0737 histograms[i]->SetLineColor(TColor::GetColor(colors[i].c_str()));
0738 histograms[i]->SetMarkerColor(TColor::GetColor(colors[i].c_str()));
0739 histograms[i]->SetMarkerStyle(markers[i]);
0740 histograms[i]->SetMarkerSize(markerSizes[i]);
0741 histograms[i]->SetLineWidth(2);
0742 histograms[i]->Draw("E1 SAME");
0743 }
0744
0745
0746 if (referenceline.first)
0747 {
0748 TLine *line = new TLine(histograms[0]->GetXaxis()->GetXmin(), referenceline.second, histograms[0]->GetXaxis()->GetXmax(), referenceline.second);
0749 line->SetLineColor(kGray + 1);
0750 line->SetLineStyle(2);
0751 line->SetLineWidth(2);
0752 line->Draw("same");
0753 }
0754
0755 TLegend *legend = new TLegend(gPad->GetLeftMargin() + 0.25,
0756 (1 - gPad->GetTopMargin()) - 0.035 - 0.04 * (histograms.size() + 1),
0757 gPad->GetLeftMargin() + 0.35,
0758 (1 - gPad->GetTopMargin()) - 0.035
0759 );
0760 legend->SetTextAlign(kHAlignLeft + kVAlignCenter);
0761 legend->SetTextSize(0.035);
0762 legend->SetFillStyle(0);
0763 legend->SetBorderSize(0);
0764 for (size_t i = 0; i < labels.size() && i < histograms.size(); ++i)
0765 {
0766 if (histograms[i])
0767 legend->AddEntry(histograms[i], labels[i].c_str(), "pel");
0768 }
0769 legend->Draw();
0770
0771 c->SaveAs(Form("%s.png", outputFileName.c_str()));
0772 c->SaveAs(Form("%s.pdf", outputFileName.c_str()));
0773 delete legend;
0774 delete c;
0775 }
0776
0777 struct ProjectionRatioOutputs
0778 {
0779 std::string prefix;
0780 std::string label;
0781 TH1D *eta_ratio = nullptr;
0782 TH1D *phi_ratio = nullptr;
0783 TH1D *eta_ratio_good = nullptr;
0784 TH1D *phi_ratio_good = nullptr;
0785 };
0786
0787 static bool StartsWith(const std::string &value, const std::string &prefix) { return value.rfind(prefix, 0) == 0; }
0788
0789 static bool EndsWith(const std::string &value, const std::string &suffix) { return value.size() >= suffix.size() && value.compare(value.size() - suffix.size(), suffix.size(), suffix) == 0; }
0790
0791 static int ExtractMultiplicityBinNumber(const std::string &base_prefix, const std::string &object_name)
0792 {
0793 const std::string probe_prefix = base_prefix + "_multPercentileBin";
0794 const std::string suffix = "_Sub";
0795 if (!StartsWith(object_name, probe_prefix) || !EndsWith(object_name, suffix))
0796 return -1;
0797
0798 const std::string bin_text = object_name.substr(probe_prefix.size(), object_name.size() - probe_prefix.size() - suffix.size());
0799 if (bin_text.empty())
0800 return -1;
0801 for (char ch : bin_text)
0802 {
0803 if (!std::isdigit(static_cast<unsigned char>(ch)))
0804 return -1;
0805 }
0806 return std::stoi(bin_text);
0807 }
0808
0809 static std::vector<std::string> CollectAvailablePrefixes(TFile *fin, const std::string &base_prefix)
0810 {
0811 std::vector<std::pair<int, std::string>> multiplicity_prefixes;
0812 if (!fin)
0813 return {base_prefix};
0814
0815 TIter next(fin->GetListOfKeys());
0816 while (auto *key = dynamic_cast<TKey *>(next()))
0817 {
0818 const std::string name = key->GetName();
0819 const int bin_number = ExtractMultiplicityBinNumber(base_prefix, name);
0820 if (bin_number >= 0)
0821 multiplicity_prefixes.push_back({bin_number, base_prefix + "_multPercentileBin" + std::to_string(bin_number)});
0822 }
0823
0824 std::sort(multiplicity_prefixes.begin(), multiplicity_prefixes.end(), [](const auto &lhs, const auto &rhs) { return lhs.first < rhs.first; });
0825 multiplicity_prefixes.erase(std::unique(multiplicity_prefixes.begin(), multiplicity_prefixes.end(), [](const auto &lhs, const auto &rhs) { return lhs.first == rhs.first; }), multiplicity_prefixes.end());
0826
0827 std::vector<std::string> prefixes = {base_prefix};
0828 for (const auto &entry : multiplicity_prefixes)
0829 prefixes.push_back(entry.second);
0830 return prefixes;
0831 }
0832
0833 static std::string PrefixSuffixForOutput(const std::string &base_prefix, const std::string ¤t_prefix)
0834 {
0835 if (current_prefix == base_prefix)
0836 return "";
0837 if (StartsWith(current_prefix, base_prefix))
0838 return current_prefix.substr(base_prefix.size());
0839 return "_" + current_prefix;
0840 }
0841
0842 static std::string PrefixDisplayLabel(const std::string &base_prefix, const std::string ¤t_prefix)
0843 {
0844 if (current_prefix == base_prefix)
0845 return "Inclusive";
0846 if (StartsWith(current_prefix, base_prefix + "_"))
0847 return current_prefix.substr(base_prefix.size() + 1);
0848 return current_prefix;
0849 }
0850
0851 static int ExtractMultiplicityBinFromPrefix(const std::string &base_prefix, const std::string ¤t_prefix)
0852 {
0853 const std::string probe_prefix = base_prefix + "_multPercentileBin";
0854 if (!StartsWith(current_prefix, probe_prefix))
0855 return -1;
0856
0857 const std::string bin_text = current_prefix.substr(probe_prefix.size());
0858 if (bin_text.empty())
0859 return -1;
0860 for (char ch : bin_text)
0861 {
0862 if (!std::isdigit(static_cast<unsigned char>(ch)))
0863 return -1;
0864 }
0865 return std::stoi(bin_text);
0866 }
0867
0868 static std::string FormatPercentileValue(double value)
0869 {
0870 std::ostringstream out;
0871 if (std::fabs(value - std::round(value)) < 1e-6)
0872 out << static_cast<int>(std::round(value));
0873 else
0874 out << std::fixed << std::setprecision(1) << value;
0875 return out.str();
0876 }
0877
0878 static std::map<int, std::pair<double, double>> LoadPercentileBoundariesByBin(const bool is_simulation)
0879 {
0880 std::map<int, std::pair<double, double>> boundaries;
0881 const std::string multiplicity_boundary_prefix = is_simulation ? "MC_" : "Data_";
0882 const std::string input_file = "figure/figure-NInttClusterCrossing/" + multiplicity_boundary_prefix + "NInttClusterPercentileBoundaries.root";
0883
0884 TFile fin(input_file.c_str(), "READ");
0885 if (fin.IsZombie())
0886 return boundaries;
0887
0888 auto *tree = dynamic_cast<TTree *>(fin.Get("trigger_percentile_intervals"));
0889 if (!tree)
0890 return boundaries;
0891
0892 int percentile_bin = -1;
0893 double percentile_low = 0.0;
0894 double percentile_high = 0.0;
0895 tree->SetBranchAddress("percentile_bin", &percentile_bin);
0896 tree->SetBranchAddress("percentile_low", &percentile_low);
0897 tree->SetBranchAddress("percentile_high", &percentile_high);
0898
0899 const Long64_t nentries = tree->GetEntries();
0900 for (Long64_t i = 0; i < nentries; ++i)
0901 {
0902 tree->GetEntry(i);
0903 boundaries[percentile_bin] = {percentile_low, percentile_high};
0904 }
0905
0906 return boundaries;
0907 }
0908
0909 static std::string PrefixDisplayLabel(const std::string &base_prefix, const std::string ¤t_prefix, const std::map<int, std::pair<double, double>> &percentile_boundaries)
0910 {
0911 if (current_prefix == base_prefix)
0912 return "Inclusive";
0913
0914 const int bin_number = ExtractMultiplicityBinFromPrefix(base_prefix, current_prefix);
0915 if (bin_number >= 0)
0916 {
0917 auto it = percentile_boundaries.find(bin_number);
0918 if (it != percentile_boundaries.end())
0919 {
0920 return "N_{INTT clusters} " + FormatPercentileValue(it->second.first) + "-" + FormatPercentileValue(it->second.second) + "%";
0921 }
0922 }
0923
0924 if (StartsWith(current_prefix, base_prefix + "_"))
0925 return current_prefix.substr(base_prefix.size() + 1);
0926 return current_prefix;
0927 }
0928
0929 static double DecodeTruthVtxZHistToken(const std::string &token, bool &is_infinite, bool &ok)
0930 {
0931 is_infinite = false;
0932 ok = false;
0933 if (token == "mInf")
0934 {
0935 is_infinite = true;
0936 ok = true;
0937 return -std::numeric_limits<double>::infinity();
0938 }
0939 if (token == "pInf")
0940 {
0941 is_infinite = true;
0942 ok = true;
0943 return std::numeric_limits<double>::infinity();
0944 }
0945
0946 std::string value = token;
0947 double sign = 1.0;
0948 if (!value.empty() && value[0] == 'm')
0949 {
0950 sign = -1.0;
0951 value.erase(0, 1);
0952 }
0953 else if (!value.empty() && value[0] == 'p')
0954 {
0955 value.erase(0, 1);
0956 }
0957
0958 std::replace(value.begin(), value.end(), 'p', '.');
0959 char *endptr = nullptr;
0960 const double parsed = std::strtod(value.c_str(), &endptr);
0961 if (endptr && *endptr == 0)
0962 {
0963 ok = true;
0964 return sign * parsed;
0965 }
0966 return 0.0;
0967 }
0968
0969 struct TruthVtxZDCAHistBin
0970 {
0971 int bin_index = -1;
0972 double xlow = 0.0;
0973 double xhigh = 0.0;
0974 TH1D *hist = nullptr;
0975 };
0976
0977 static bool ParseTruthVtxZDCAHistName(const std::string &hist_name, const std::string &prefix, TruthVtxZDCAHistBin &info)
0978 {
0979 const std::string probe_prefix = prefix + "_SilSeedDCA3D_truthVtxZ_bin";
0980 if (!StartsWith(hist_name, probe_prefix))
0981 return false;
0982
0983 const std::string remainder = hist_name.substr(probe_prefix.size());
0984 if (remainder.size() < 4 || !std::isdigit(static_cast<unsigned char>(remainder[0])) || !std::isdigit(static_cast<unsigned char>(remainder[1])) || remainder[2] != '_')
0985 return false;
0986
0987 const int bin_index = std::stoi(remainder.substr(0, 2));
0988 const std::string encoded_bounds = remainder.substr(3);
0989 const size_t sep = encoded_bounds.find('_');
0990 if (sep == std::string::npos)
0991 return false;
0992
0993 bool low_inf = false;
0994 bool high_inf = false;
0995 bool low_ok = false;
0996 bool high_ok = false;
0997 const double low_raw = DecodeTruthVtxZHistToken(encoded_bounds.substr(0, sep), low_inf, low_ok);
0998 const double high_raw = DecodeTruthVtxZHistToken(encoded_bounds.substr(sep + 1), high_inf, high_ok);
0999 if (!low_ok || !high_ok)
1000 return false;
1001
1002 const double low = low_inf ? (high_raw - 10.0) : low_raw;
1003 const double high = high_inf ? (low_raw + 10.0) : high_raw;
1004 if (!(high > low))
1005 return false;
1006
1007 info.bin_index = bin_index;
1008 info.xlow = low;
1009 info.xhigh = high;
1010 return true;
1011 }
1012
1013 static std::vector<TruthVtxZDCAHistBin> CollectTruthVtxZDCAHistBins(TFile *fin, const std::string &prefix)
1014 {
1015 std::vector<TruthVtxZDCAHistBin> bins;
1016 if (!fin)
1017 return bins;
1018
1019 TIter next(fin->GetListOfKeys());
1020 while (auto *key = dynamic_cast<TKey *>(next()))
1021 {
1022 const std::string hist_name = key->GetName();
1023 TruthVtxZDCAHistBin info;
1024 if (!ParseTruthVtxZDCAHistName(hist_name, prefix, info))
1025 continue;
1026
1027 info.hist = dynamic_cast<TH1D *>(fin->Get(hist_name.c_str()));
1028 if (!info.hist)
1029 continue;
1030 bins.push_back(info);
1031 }
1032
1033 std::sort(bins.begin(), bins.end(), [](const TruthVtxZDCAHistBin &lhs, const TruthVtxZDCAHistBin &rhs) { return lhs.bin_index < rhs.bin_index; });
1034 return bins;
1035 }
1036
1037 static TH2D *BuildTruthVtxZDCAHeatmap(TFile *fin, const std::string &prefix)
1038 {
1039 std::vector<TruthVtxZDCAHistBin> bins = CollectTruthVtxZDCAHistBins(fin, prefix);
1040 if (bins.empty())
1041 return nullptr;
1042
1043 TH1D *href = bins.front().hist;
1044 if (!href)
1045 return nullptr;
1046
1047 std::vector<double> xedges;
1048 xedges.reserve(bins.size() + 1);
1049 xedges.push_back(bins.front().xlow);
1050 for (const auto &entry : bins)
1051 xedges.push_back(entry.xhigh);
1052
1053 std::vector<double> yedges;
1054 yedges.reserve(href->GetNbinsX() + 1);
1055 for (int ibin = 1; ibin <= href->GetNbinsX(); ++ibin)
1056 yedges.push_back(href->GetXaxis()->GetBinLowEdge(ibin));
1057 yedges.push_back(href->GetXaxis()->GetBinUpEdge(href->GetNbinsX()));
1058
1059 TH2D *h2 = new TH2D((prefix + "_SilSeedDCA3D_truthVtxZ_2D").c_str(), (prefix + "_SilSeedDCA3D_truthVtxZ_2D").c_str(), static_cast<int>(bins.size()), xedges.data(), href->GetNbinsX(), yedges.data());
1060 h2->SetDirectory(nullptr);
1061 h2->Sumw2();
1062
1063 for (size_t ix = 0; ix < bins.size(); ++ix)
1064 {
1065 TH1D *h1 = bins[ix].hist;
1066 if (!h1)
1067 continue;
1068
1069 const double norm = h1->Integral(1, h1->GetNbinsX());
1070 for (int iy = 1; iy <= h1->GetNbinsX(); ++iy)
1071 {
1072 const double value = (norm > 0.0) ? (h1->GetBinContent(iy) / norm) : 0.0;
1073 const double error = (norm > 0.0) ? (h1->GetBinError(iy) / norm) : 0.0;
1074 h2->SetBinContent(static_cast<int>(ix) + 1, iy, value);
1075 h2->SetBinError(static_cast<int>(ix) + 1, iy, error);
1076 }
1077 }
1078
1079 return h2;
1080 }
1081
1082 static std::vector<std::pair<TH1D *, std::string>> BuildMergedAbsTruthVtxZDCAHists(TFile *fin, const std::string &prefix)
1083 {
1084 std::vector<std::pair<TH1D *, std::string>> merged_hists;
1085 const std::vector<TruthVtxZDCAHistBin> bins = CollectTruthVtxZDCAHistBins(fin, prefix);
1086 if (bins.empty())
1087 return merged_hists;
1088
1089 auto format_abs_bound = [](double value) -> std::string
1090 {
1091 std::ostringstream out;
1092 out << std::fixed << std::setprecision((std::fabs(value - std::round(value)) < 1e-6) ? 0 : 1) << value;
1093 return out.str();
1094 };
1095
1096 const size_t npairs = (bins.size() + 1) / 2;
1097 merged_hists.reserve(npairs);
1098 for (size_t i = 0; i < npairs; ++i)
1099 {
1100 const size_t j = bins.size() - 1 - i;
1101 TH1D *href = bins[i].hist;
1102 if (!href)
1103 continue;
1104
1105 TH1D *hmerge = dynamic_cast<TH1D *>(href->Clone((prefix + "_SilSeedDCA3D_absTruthVtxZ_bin" + std::to_string(i)).c_str()));
1106 if (!hmerge)
1107 continue;
1108 hmerge->SetDirectory(nullptr);
1109 if (j != i && bins[j].hist)
1110 hmerge->Add(bins[j].hist);
1111
1112 const double abs_low = (j != i) ? std::min(std::fabs(bins[i].xhigh), std::fabs(bins[j].xlow)) : std::min(std::fabs(bins[i].xlow), std::fabs(bins[i].xhigh));
1113 const double abs_high = (j != i) ? std::max(std::fabs(bins[i].xlow), std::fabs(bins[j].xhigh)) : std::max(std::fabs(bins[i].xlow), std::fabs(bins[i].xhigh));
1114 merged_hists.push_back({hmerge, "|z_{truth}|=" + format_abs_bound(abs_low) + "-" + format_abs_bound(abs_high) + " cm"});
1115 }
1116
1117 return merged_hists;
1118 }
1119
1120 static ProjectionRatioOutputs plot_tklCombinatoric_single(const std::string &inputfile = "test_OO_combinatoric.root",
1121 const std::string &prefix = "tkl_Combinatoric",
1122 const std::string &mode = "absdPhi",
1123 bool is_simulation = false,
1124 double zvtx_cut_min = -10.0,
1125 double zvtx_cut_max = 10.0,
1126 double raw_count_xmin = 0.0,
1127 double raw_count_xmax = kDirectRawCountCut,
1128 double plot_xmin = 0.0,
1129 double plot_xmax = 0.2,
1130 const std::string &outstem = "plot_tklCombinatoric",
1131 const std::string &plotdir = "./figure/figure-tklCombinatoric",
1132 bool zero_x_error = true,
1133 const std::string &ratio_draw_option = "colz text",
1134 const std::string &ratio_text_format = ".3f",
1135 bool use_rounded_zmax = true,
1136 const std::string &display_label = "")
1137 {
1138 ProjectionRatioOutputs ratio_outputs;
1139 ratio_outputs.prefix = prefix;
1140 ratio_outputs.label = display_label.empty() ? prefix : display_label;
1141
1142 gStyle->SetOptStat(0);
1143 TGaxis::SetMaxDigits(3);
1144 if (zero_x_error)
1145 gStyle->SetErrorX(0.f);
1146
1147 gSystem->mkdir(plotdir.c_str(), true);
1148 const std::string single_canvas_dir = plotdir + "/single-canvas";
1149 gSystem->mkdir(single_canvas_dir.c_str(), true);
1150
1151 auto *fin = TFile::Open(inputfile.c_str(), "READ");
1152 if (!fin || fin->IsZombie())
1153 {
1154 std::cerr << "[plot_tklCombinatoric] Cannot open: " << inputfile << "\n";
1155 return ratio_outputs;
1156 }
1157
1158 auto *hGrid = GetPrimaryGrid(fin, prefix);
1159 if (!hGrid)
1160 {
1161 std::cerr << "[plot_tklCombinatoric] Cannot load " << prefix << "_Sig/_Bkg/_Sub for grid axes.\n";
1162 fin->Close();
1163 return ratio_outputs;
1164 }
1165 const int nEta = hGrid->GetXaxis()->GetNbins();
1166 const int nPhi = hGrid->GetYaxis()->GetNbins();
1167
1168 auto *hSub = dynamic_cast<TH2D *>(fin->Get((prefix + "_Sub").c_str()));
1169 if (!hSub)
1170 {
1171 std::cerr << "[plot_tklCombinatoric] Missing " << prefix << "_Sub.\n";
1172 fin->Close();
1173 return ratio_outputs;
1174 }
1175
1176 auto *hSilSeed = dynamic_cast<TH2D *>(fin->Get((prefix + "_SilSeedEtaPhi").c_str()));
1177 if (!hSilSeed)
1178 {
1179 std::cerr << "[plot_tklCombinatoric] Missing " << prefix << "_SilSeedEtaPhi.\n";
1180 fin->Close();
1181 return ratio_outputs;
1182 }
1183 auto *hSilSeed_nMVTX3nINTT2 = dynamic_cast<TH2D *>(fin->Get((prefix + "_SilSeedEtaPhi_nMVTX3nINTT2").c_str()));
1184 auto *hSilSeed_nMVTX3nINTT1 = dynamic_cast<TH2D *>(fin->Get((prefix + "_SilSeedEtaPhi_nMVTX3nINTT1").c_str()));
1185 auto *hSilSeed_nMVTX2nINTT2 = dynamic_cast<TH2D *>(fin->Get((prefix + "_SilSeedEtaPhi_nMVTX2nINTT2").c_str()));
1186 auto *hSilSeed_nMVTX2nINTT1 = dynamic_cast<TH2D *>(fin->Get((prefix + "_SilSeedEtaPhi_nMVTX2nINTT1").c_str()));
1187 auto *hSilSeed_nMVTXnINTTOther = dynamic_cast<TH2D *>(fin->Get((prefix + "_SilSeedEtaPhi_nMVTXnINTTOther").c_str()));
1188 auto *hINTTClusterVsAssociatedSilSeed = dynamic_cast<TH2D *>(fin->Get((prefix + "_INTTClusterVsAssociatedSilSeed").c_str()));
1189 auto *hINTTClusterVsAllSilSeed = dynamic_cast<TH2D *>(fin->Get((prefix + "_INTTClusterVsAllSilSeed").c_str()));
1190 auto *hSeedCrossing_selectedCrossings = dynamic_cast<TH1D *>(fin->Get((prefix + "_SeedCrossing_selectedCrossings").c_str()));
1191 auto *hSeedCrossing_selectedCrossings_largeDiffNClusterSeeds = dynamic_cast<TH1D *>(fin->Get((prefix + "_SeedCrossing_selectedCrossings_largeDiffNClusterSeeds").c_str()));
1192 auto *hDCA3D_Sig = dynamic_cast<TH1D *>(fin->Get((prefix + "_DCA3D_Sig").c_str()));
1193 auto *hDCA3D_Bkg = dynamic_cast<TH1D *>(fin->Get((prefix + "_DCA3D_Bkg").c_str()));
1194 auto *hSilSeedDCA3D = dynamic_cast<TH1D *>(fin->Get((prefix + "_SilSeedDCA3D").c_str()));
1195 TH2D *hSilSeedDCA3D_truthVtxZ_2D = is_simulation ? BuildTruthVtxZDCAHeatmap(fin, prefix) : nullptr;
1196 std::vector<std::pair<TH1D *, std::string>> hSilSeedDCA3D_truthVtxZ_absMerged = is_simulation ? BuildMergedAbsTruthVtxZDCAHists(fin, prefix) : std::vector<std::pair<TH1D *, std::string>>{};
1197
1198 if (!hSilSeed_nMVTX3nINTT2 || !hSilSeed_nMVTX3nINTT1 || !hSilSeed_nMVTX2nINTT2 || !hSilSeed_nMVTX2nINTT1 || !hSilSeed_nMVTXnINTTOther)
1199 {
1200 std::cerr << "[plot_tklCombinatoric] Missing one or more silicon-seed category histograms." << std::endl;
1201 fin->Close();
1202 return ratio_outputs;
1203 }
1204
1205
1206 auto *hSilSeed_good = dynamic_cast<TH2D *>(hSilSeed_nMVTX3nINTT2->Clone((prefix + "_SilSeedEtaPhi_good").c_str()));
1207 if (!hSilSeed_good)
1208 {
1209 std::cerr << "[plot_tklCombinatoric] Failed to build " << prefix << "_SilSeedEtaPhi_good." << std::endl;
1210 fin->Close();
1211 return ratio_outputs;
1212 }
1213 hSilSeed_good->SetDirectory(nullptr);
1214 hSilSeed_good->Add(hSilSeed_nMVTX3nINTT1);
1215 hSilSeed_good->Add(hSilSeed_nMVTX2nINTT2);
1216
1217
1218 auto *hSeed = dynamic_cast<TH2D *>(fin->Get((prefix + "_SeedEtaPhi").c_str()));
1219 if (!hSeed)
1220 {
1221 std::cerr << "[plot_tklCombinatoric] Missing " << prefix << "_SeedEtaPhi.\n";
1222 fin->Close();
1223 return ratio_outputs;
1224 }
1225
1226 TH2D *hPrimaryCharged = nullptr;
1227 if (is_simulation)
1228 {
1229 hPrimaryCharged = dynamic_cast<TH2D *>(fin->Get((prefix + "_PrimaryChargedEtaPhi").c_str()));
1230 if (!hPrimaryCharged)
1231 {
1232 std::cerr << "[plot_tklCombinatoric] Missing " << prefix << "_PrimaryChargedEtaPhi for simulation mode.\n";
1233 fin->Close();
1234 return ratio_outputs;
1235 }
1236 }
1237
1238 bool use_direct_raw_count_grid = false;
1239 TH2D *h2_nraw = nullptr;
1240 if (auto *hRawDirect = GetRawCountGrid(fin, prefix))
1241 {
1242 h2_nraw = dynamic_cast<TH2D *>(hRawDirect->Clone((prefix + "_NRawEtaPhi").c_str()));
1243 if (h2_nraw)
1244 {
1245 h2_nraw->SetDirectory(nullptr);
1246 use_direct_raw_count_grid = true;
1247 }
1248 }
1249 if (!h2_nraw)
1250 {
1251 h2_nraw = dynamic_cast<TH2D *>(hSub->Clone((prefix + "_NRawEtaPhi").c_str()));
1252 h2_nraw->SetDirectory(nullptr);
1253 h2_nraw->Reset("ICES");
1254 std::cerr << "[plot_tklCombinatoric] Missing " << prefix << "_Sub_AbsDphiCount. Falling back to per-cell integration for h2_nraw.\n";
1255 }
1256 if (h2_nraw->GetNbinsX() != hSilSeed->GetNbinsX() || h2_nraw->GetNbinsY() != hSilSeed->GetNbinsY())
1257 {
1258 std::cerr << "[plot_tklCombinatoric] Incompatible binning between " << prefix << "_Sub and " << prefix << "_SilSeedEtaPhi.\n";
1259 delete h2_nraw;
1260 fin->Close();
1261 return ratio_outputs;
1262 }
1263 if (is_simulation && (h2_nraw->GetNbinsX() != hPrimaryCharged->GetNbinsX() || h2_nraw->GetNbinsY() != hPrimaryCharged->GetNbinsY()))
1264 {
1265 std::cerr << "[plot_tklCombinatoric] Incompatible binning between " << prefix << "_NRawEtaPhi and " << prefix << "_PrimaryChargedEtaPhi.\n";
1266 delete h2_nraw;
1267 fin->Close();
1268 return ratio_outputs;
1269 }
1270 auto check_same_binning = [](const TH2D *ha, const TH2D *hb) -> bool { return ha && hb && ha->GetNbinsX() == hb->GetNbinsX() && ha->GetNbinsY() == hb->GetNbinsY(); };
1271 if (!check_same_binning(hSilSeed_nMVTX3nINTT2, hSilSeed) || !check_same_binning(hSilSeed_nMVTX3nINTT1, hSilSeed) || !check_same_binning(hSilSeed_nMVTX2nINTT2, hSilSeed) || !check_same_binning(hSilSeed_nMVTX2nINTT1, hSilSeed) || !check_same_binning(hSilSeed_nMVTXnINTTOther, hSilSeed))
1272 {
1273 std::cerr << "[plot_tklCombinatoric] Incompatible binning between silicon-seed categories and inclusive silicon-seed histogram.\n";
1274 delete h2_nraw;
1275 fin->Close();
1276 return ratio_outputs;
1277 }
1278
1279 const bool etaAscending = hGrid->GetXaxis()->GetBinCenter(1) <= hGrid->GetXaxis()->GetBinCenter(nEta);
1280 const bool phiAscending = hGrid->GetYaxis()->GetBinCenter(1) <= hGrid->GetYaxis()->GetBinCenter(nPhi);
1281 const bool use_custom_xrange = std::isfinite(plot_xmin) && std::isfinite(plot_xmax) && (plot_xmax > plot_xmin);
1282
1283
1284 double global_ymax = 0.0;
1285 double ref_bin_size = -1.0;
1286 for (int ieta = 0; ieta < nEta; ++ieta)
1287 {
1288 for (int iphi = 0; iphi < nPhi; ++iphi)
1289 {
1290 auto *hSig = GetCellHist(fin, "dR_Sig", prefix, mode, "sig", ieta, iphi, hGrid);
1291 auto *hBkg = GetCellHist(fin, "dR_Bkg", prefix, mode, "bkg", ieta, iphi, hGrid);
1292
1293 if (hSig)
1294 {
1295 global_ymax = std::max(global_ymax, hSig->GetMaximum());
1296 if (ref_bin_size < 0.0)
1297 ref_bin_size = hSig->GetXaxis()->GetBinWidth(1);
1298 }
1299 if (hBkg)
1300 {
1301 global_ymax = std::max(global_ymax, hBkg->GetMaximum());
1302 if (ref_bin_size < 0.0)
1303 ref_bin_size = hBkg->GetXaxis()->GetBinWidth(1);
1304 }
1305 }
1306 }
1307 const double shared_ymax = (global_ymax > 0.0) ? 1.3 * global_ymax : 1.0;
1308 const std::string yAxisTitle = (ref_bin_size > 0.0) ? Form("Counts / (%.2g radian)", ref_bin_size) : "Counts / (bin size) [radian]";
1309
1310 auto *c = new TCanvas("c_tklCombinatoric", "c_tklCombinatoric", 5500, 5000);
1311 const float lM = 0.03f, rM = 0.02f, bM = 0.03f, tM = 0.02f;
1312 CanvasPartition(c, nEta, nPhi, lM, rM, bM, tM);
1313
1314 auto *pad00 = dynamic_cast<TPad *>(c->FindObject("pad_0_0"));
1315 if (!pad00)
1316 {
1317 std::cerr << "[plot_tklCombinatoric] pad_0_0 not found.\n";
1318 delete h2_nraw;
1319 fin->Close();
1320 delete c;
1321 return ratio_outputs;
1322 }
1323
1324 std::vector<TH1D *> frames;
1325 std::vector<TLine *> cutLines;
1326
1327 int nDrawn = 0;
1328
1329 auto styleHist = [](TH1D *h, Style_t ms, Color_t col, float markerSize)
1330 {
1331 h->SetMarkerStyle(ms);
1332 h->SetMarkerSize(markerSize);
1333 h->SetLineColor(col);
1334 h->SetMarkerColor(col);
1335 h->Draw("E1 SAME");
1336 };
1337
1338 auto createAndDrawCutLine = [&](double x, int linewidth) -> TLine *
1339 {
1340 auto *line = new TLine(x, 0.0, x, shared_ymax);
1341 line->SetLineStyle(3);
1342 line->SetLineWidth(linewidth);
1343 line->SetLineColor(kBlue + 2);
1344 line->Draw("SAME");
1345 cutLines.push_back(line);
1346 return line;
1347 };
1348
1349 auto SaveSingleCellPlot = [&](bool isdata, TH1D *hSig, TH1D *hBkg, int ieta, int iphi, double etaLow, double etaHigh, double phiLow, double phiHigh, double nraw, std::string str_mode, std::pair<double, double> cutrange)
1350 {
1351 if (!hSig && !hBkg)
1352 return;
1353
1354 const std::string range_suffix = EtaPhiRangeNameSuffix(hGrid, ieta, iphi);
1355 const std::string cname = "c_single_" + range_suffix;
1356 TCanvas csingle(cname.c_str(), cname.c_str(), 800, 700);
1357 csingle.cd();
1358 gPad->SetLeftMargin(0.16);
1359 gPad->SetRightMargin(0.05);
1360
1361 gPad->SetTopMargin(0.06);
1362 gPad->SetTicks(1, 1);
1363
1364 TH1D *frame_single = (TH1D *)(hSig ? hSig : hBkg)->Clone((cname + "_frame").c_str());
1365 frame_single->SetDirectory(nullptr);
1366 frame_single->Reset("ICES");
1367
1368 frame_single->SetMaximum(shared_ymax);
1369 if (use_custom_xrange)
1370 frame_single->GetXaxis()->SetRangeUser(plot_xmin, plot_xmax);
1371 frame_single->GetXaxis()->SetTitle("INTT doublet |#Delta#phi| [radian]");
1372 frame_single->GetYaxis()->SetTitle(yAxisTitle.c_str());
1373 frame_single->GetYaxis()->SetTitleOffset(1.5);
1374 frame_single->Draw("AXIS");
1375
1376 if (hSig)
1377 styleHist(hSig, kFullCircle, kBlack, 1.0f);
1378 if (hBkg)
1379 styleHist(hBkg, kOpenCircle, kRed + 1, 1.0f);
1380
1381 TLegend leg(gPad->GetLeftMargin() + 0.05,
1382 1. - gPad->GetTopMargin() - 0.22,
1383 gPad->GetLeftMargin() + 0.2,
1384 1. - gPad->GetTopMargin() - 0.04
1385 );
1386 leg.SetHeader(Form("%s. #eta=[%s,%s], #phi=[%s,%s]", isdata ? "Data" : "Simulation", to_string_with_precision(etaLow).c_str(), to_string_with_precision(etaHigh).c_str(), to_string_with_precision(phiLow).c_str(), to_string_with_precision(phiHigh).c_str()));
1387 leg.SetBorderSize(0);
1388 leg.SetFillStyle(0);
1389 leg.SetTextSize(0.035);
1390 if (hSig)
1391 leg.AddEntry(hSig, "Unrotated", "pel");
1392 if (hBkg)
1393 leg.AddEntry(hBkg, "Rotated", "pel");
1394 leg.AddEntry((TObject *)0, Form("N^{Uncorr.}_{Doublets}=%.1f (%s=[%s,%s])", nraw, str_mode.c_str(), to_string_with_precision(cutrange.first).c_str(), to_string_with_precision(cutrange.second).c_str()), "");
1395 leg.Draw();
1396
1397 createAndDrawCutLine(cutrange.first, 2);
1398 createAndDrawCutLine(cutrange.second, 2);
1399
1400 const std::string outbase = single_canvas_dir + "/" + outstem + "_" + range_suffix;
1401 csingle.SaveAs((outbase + ".pdf").c_str());
1402 csingle.SaveAs((outbase + ".png").c_str());
1403 delete frame_single;
1404 };
1405
1406
1407 auto modeStringLambda = [&]() -> std::string
1408 {
1409 if (mode == "absdPhi")
1410 return "|#Delta#phi|";
1411 else if (mode == "dR")
1412 return "#DeltaR";
1413 else
1414 return mode;
1415 };
1416 std::string modeString = modeStringLambda();
1417 std::string multiplicity_text = "inclusive";
1418 const std::string multiplicity_prefix_1 = "Cluster multiplicity percentile ";
1419 const std::string multiplicity_prefix_2 = "N_{INTT clusters} ";
1420 if (StartsWith(ratio_outputs.label, multiplicity_prefix_1))
1421 {
1422 multiplicity_text = ratio_outputs.label.substr(multiplicity_prefix_1.size());
1423 }
1424 else if (StartsWith(ratio_outputs.label, multiplicity_prefix_2))
1425 {
1426 multiplicity_text = ratio_outputs.label.substr(multiplicity_prefix_2.size());
1427 }
1428 else if (!ratio_outputs.label.empty() && ratio_outputs.label != "Inclusive")
1429 {
1430 multiplicity_text = ratio_outputs.label;
1431 }
1432 const std::string sample_label = Form("%s, %s#leqZ_{vtx}^{silicon}#leq%s [cm], %s", is_simulation ? "Simulation" : "Data", to_string_with_precision(zvtx_cut_min, 1).c_str(), to_string_with_precision(zvtx_cut_max, 1).c_str(), multiplicity_text.c_str());
1433
1434 for (int ieta = 0; ieta < nEta; ++ieta)
1435 {
1436 for (int iphi = 0; iphi < nPhi; ++iphi)
1437 {
1438
1439 c->cd(0);
1440 const int padEta = etaAscending ? ieta : (nEta - 1 - ieta);
1441 const int padPhi = phiAscending ? iphi : (nPhi - 1 - iphi);
1442
1443 auto *pad = dynamic_cast<TPad *>(c->FindObject(TString::Format("pad_%d_%d", padEta, padPhi).Data()));
1444 if (!pad)
1445 continue;
1446
1447 pad->cd();
1448 pad->SetTicks(1, 1);
1449
1450 float xFactor = pad00->GetAbsWNDC() / pad->GetAbsWNDC();
1451 float yFactor = pad00->GetAbsHNDC() / pad->GetAbsHNDC();
1452
1453 auto *hSig = GetCellHist(fin, "dR_Sig", prefix, mode, "sig", ieta, iphi, hGrid);
1454 auto *hBkg = GetCellHist(fin, "dR_Bkg", prefix, mode, "bkg", ieta, iphi, hGrid);
1455 if (!hSig && !hBkg)
1456 continue;
1457
1458
1459 TH1D *baseHist = hSig ? hSig : hBkg;
1460 auto *frame = (TH1D *)baseHist->Clone(Form("hframe_ieta%d_iphi%d", ieta, iphi));
1461 frame->SetDirectory(nullptr);
1462 frame->Reset("ICES");
1463 frames.push_back(frame);
1464
1465
1466 frame->SetMaximum(shared_ymax);
1467 frame->SetTitle("");
1468 frame->GetXaxis()->SetNdivisions(505);
1469 frame->GetYaxis()->SetNdivisions(505);
1470 if (use_custom_xrange)
1471 {
1472 frame->GetXaxis()->SetRangeUser(plot_xmin, plot_xmax);
1473 }
1474 frame->GetXaxis()->SetLabelFont(43);
1475 frame->GetXaxis()->SetLabelSize(34);
1476 frame->GetYaxis()->SetLabelFont(43);
1477 frame->GetYaxis()->SetLabelSize(34);
1478 frame->GetXaxis()->SetTickLength(yFactor * 0.06f / xFactor);
1479 frame->GetYaxis()->SetTickLength(xFactor * 0.04f / yFactor);
1480 frame->GetXaxis()->SetTitle("");
1481 frame->GetYaxis()->SetTitle("");
1482 frame->Draw("AXIS");
1483
1484 createAndDrawCutLine(raw_count_xmin, 1);
1485 createAndDrawCutLine(raw_count_xmax, 1);
1486
1487 if (hSig)
1488 styleHist(hSig, kFullCircle, kBlack, 0.0f);
1489 if (hBkg)
1490 styleHist(hBkg, kOpenCircle, kRed + 1, 0.0f);
1491
1492 double nraw = h2_nraw->GetBinContent(ieta + 1, iphi + 1);
1493 double nraw_err = h2_nraw->GetBinError(ieta + 1, iphi + 1);
1494 if (!use_direct_raw_count_grid)
1495 {
1496 nraw = IntegralInRange(hSig, raw_count_xmin, raw_count_xmax) - IntegralInRange(hBkg, raw_count_xmin, raw_count_xmax);
1497 nraw_err = (nraw > 0.0) ? std::sqrt(nraw) : 0.0;
1498 h2_nraw->SetBinContent(ieta + 1, iphi + 1, nraw);
1499 h2_nraw->SetBinError(ieta + 1, iphi + 1, nraw_err);
1500 }
1501
1502 const double etaLow = hGrid->GetXaxis()->GetBinLowEdge(ieta + 1);
1503 const double etaHigh = hGrid->GetXaxis()->GetBinUpEdge(ieta + 1);
1504 const double phiLow = hGrid->GetYaxis()->GetBinLowEdge(iphi + 1);
1505 const double phiHigh = hGrid->GetYaxis()->GetBinUpEdge(iphi + 1);
1506
1507 TLatex lbl;
1508 lbl.SetNDC();
1509 lbl.SetTextFont(43);
1510 lbl.SetTextSize(28);
1511 lbl.DrawLatex(XtoPad(0.05), YtoPad(0.90), Form("#eta=[%s,%s]", to_string_with_precision(etaLow).c_str(), to_string_with_precision(etaHigh).c_str()));
1512 lbl.DrawLatex(XtoPad(0.05), YtoPad(0.81), Form("#phi=[%s,%s]", to_string_with_precision(phiLow).c_str(), to_string_with_precision(phiHigh).c_str()));
1513 lbl.DrawLatex(XtoPad(0.05), YtoPad(0.70), use_direct_raw_count_grid ? Form("N_{raw}^{stored}=%.1f", nraw) : Form("N_{raw}(%.3f<|#Delta#phi|<%.3f)=%.1f", raw_count_xmin, raw_count_xmax, nraw));
1514
1515
1516
1517 pad->Modified();
1518 ++nDrawn;
1519 }
1520 }
1521
1522
1523 c->cd();
1524 TLatex tx;
1525 tx.SetNDC();
1526 tx.SetTextFont(43);
1527
1528 tx.SetTextSize(50);
1529 tx.SetTextAlign(kHAlignCenter + kVAlignTop);
1530 tx.DrawLatex(0.5, bM / 2., "|#Delta#phi| [radian]");
1531
1532 tx.SetTextAngle(90);
1533 tx.SetTextAlign(kHAlignCenter + kVAlignBottom);
1534 tx.SetTextSize(50);
1535 tx.DrawLatex(lM / 2., 0.5, yAxisTitle.c_str());
1536
1537 std::cout << "[plot_tklCombinatoric] Pads with histograms: " << nDrawn << " / " << nEta * nPhi << "\n";
1538
1539 auto build_and_draw_ratio = [&](const TH2D *hNum,
1540 const TH2D *hDen,
1541 const std::string &ratio_hist_name,
1542 const std::string &ztitle,
1543 const std::string &out_suffix,
1544 int round_zmax_mode = -1,
1545 const std::string &text_format = "") -> TH2D *
1546 {
1547 double ratio_zmax = -1;
1548 auto *h2_ratio = dynamic_cast<TH2D *>(h2_nraw->Clone(ratio_hist_name.c_str()));
1549 h2_ratio->SetDirectory(nullptr);
1550 h2_ratio->Reset("ICES");
1551
1552 for (int ibx = 1; ibx <= h2_ratio->GetNbinsX(); ++ibx)
1553 {
1554 for (int iby = 1; iby <= h2_ratio->GetNbinsY(); ++iby)
1555 {
1556 const double num = hNum->GetBinContent(ibx, iby);
1557 const double den = hDen->GetBinContent(ibx, iby);
1558 const double ratio = (den > 0.0) ? (num / den) : 0.0;
1559 if (ratio > ratio_zmax)
1560 ratio_zmax = ratio;
1561
1562
1563 double ratio_err = 0.0;
1564 if (den > 0.0 && num > 0.0)
1565 {
1566 const double sigma_num = std::sqrt(num);
1567 const double sigma_den = std::sqrt(den);
1568 const double rel_num = sigma_num / num;
1569 const double rel_den = sigma_den / den;
1570 ratio_err = ratio * std::sqrt(rel_num * rel_num + rel_den * rel_den);
1571 }
1572 h2_ratio->SetBinContent(ibx, iby, ratio);
1573 h2_ratio->SetBinError(ibx, iby, ratio_err);
1574 }
1575 }
1576
1577
1578 const bool round_zmax = (round_zmax_mode < 0) ? use_rounded_zmax : (round_zmax_mode != 0);
1579 double ratio_zmax_plot = ratio_zmax;
1580 if (round_zmax)
1581 {
1582 ratio_zmax_plot = std::ceil(ratio_zmax / 5.0) * 5.0;
1583 if (ratio_zmax_plot < ratio_zmax)
1584 ratio_zmax_plot += 5.0;
1585 }
1586 const std::string ratio_paint_fmt = text_format.empty() ? ratio_text_format : text_format;
1587
1588 draw2Dhistogram_adjZaxis(h2_ratio,
1589 false,
1590 "#eta",
1591 "#phi [radian]",
1592 ztitle,
1593 std::pair<double, double>(0, ratio_zmax_plot), {sample_label},
1594 ratio_draw_option,
1595 ratio_paint_fmt,
1596 true,
1597 false,
1598 plotdir + "/" + outstem + out_suffix);
1599 return h2_ratio;
1600 };
1601
1602 auto draw_input_2d = [&](TH2D *h2,
1603 const std::string &ztitle,
1604 const std::string &out_suffix,
1605 const std::string &drawoption = "colz",
1606 bool show_bin_value_with_error = false,
1607 bool show_bin_value_only = false,
1608 const std::string &xtitle = "#eta",
1609 const std::string &ytitle = "#phi [radian]",
1610 bool logz = false)
1611 {
1612 if (!h2)
1613 return;
1614
1615
1616 double zmax = h2->GetMaximum();
1617 double zmax_plot = zmax;
1618 if (use_rounded_zmax)
1619 {
1620 zmax_plot = std::ceil(zmax / 5.0) * 5.0;
1621 if (zmax_plot < zmax)
1622 zmax_plot += 5.0;
1623 }
1624
1625 draw2Dhistogram_adjZaxis(h2,
1626 logz,
1627 xtitle,
1628 ytitle,
1629 ztitle,
1630 std::pair<double, double>(0, zmax_plot),
1631 {sample_label},
1632 drawoption,
1633 ratio_text_format,
1634 show_bin_value_with_error,
1635 show_bin_value_only,
1636 plotdir + "/" + outstem + out_suffix);
1637 };
1638
1639 auto make_projection = [&](TH2D *h2, bool project_x, const std::string &hist_name) -> TH1D *
1640 {
1641 if (!h2)
1642 return nullptr;
1643
1644 TH1D *hproj = project_x ? h2->ProjectionX(hist_name.c_str(), 1, h2->GetNbinsY()) : h2->ProjectionY(hist_name.c_str(), 1, h2->GetNbinsX());
1645 if (!hproj)
1646 return nullptr;
1647 hproj->SetDirectory(nullptr);
1648 return hproj;
1649 };
1650
1651 auto draw_projection = [&](TH2D *h2, bool project_x, const std::string &axis_title, const std::string &ytitle, const std::string &out_suffix)
1652 {
1653 if (!h2)
1654 return;
1655
1656 TH1D *hproj = make_projection(h2, project_x, std::string(h2->GetName()) + (project_x ? "_projX" : "_projY"));
1657 if (!hproj)
1658 return;
1659 hproj->GetXaxis()->SetTitle(axis_title.c_str());
1660 hproj->GetYaxis()->SetTitle(ytitle.c_str());
1661 draw1Dhistogram(hproj, false, axis_title, ytitle, {sample_label}, "hist e1", plotdir + "/" + outstem + out_suffix);
1662 delete hproj;
1663 };
1664
1665 auto build_projection_ratio = [&](std::pair<TH2D *, std::string> hNum2D, std::pair<TH2D *, std::string> hDen2D, bool project_x, const std::string &axis_title, const std::string &ratio_ytitle, const std::string &out_suffix, const std::string &ratio_hist_name) -> TH1D *
1666 {
1667 if (!hNum2D.first || !hDen2D.first)
1668 return nullptr;
1669
1670
1671 TH1D *hNumProj = project_x ? hNum2D.first->ProjectionX((ratio_hist_name + "_num").c_str(), 1, hNum2D.first->GetNbinsY()) : hNum2D.first->ProjectionY((ratio_hist_name + "_num").c_str(), 1, hNum2D.first->GetNbinsX());
1672 TH1D *hDenProj = project_x ? hDen2D.first->ProjectionX((ratio_hist_name + "_den").c_str(), 1, hDen2D.first->GetNbinsY()) : hDen2D.first->ProjectionY((ratio_hist_name + "_den").c_str(), 1, hDen2D.first->GetNbinsX());
1673 if (!hNumProj || !hDenProj)
1674 {
1675 delete hNumProj;
1676 delete hDenProj;
1677 return nullptr;
1678 }
1679
1680 hNumProj->SetDirectory(nullptr);
1681 hDenProj->SetDirectory(nullptr);
1682 drawComparisonWithRatio(std::make_pair(hDenProj, hDen2D.second), std::make_pair(hNumProj, hNum2D.second), false, axis_title, "Counts", ratio_ytitle, {sample_label}, plotdir + "/" + outstem + out_suffix + "_components");
1683
1684 TH1D *hRatio = dynamic_cast<TH1D *>(hNumProj->Clone(ratio_hist_name.c_str()));
1685 if (hRatio)
1686 {
1687 hRatio->SetDirectory(nullptr);
1688 hRatio->Reset("ICES");
1689 hRatio->Divide(hNumProj, hDenProj, 1.0, 1.0);
1690 draw1Dhistogram(hRatio, false, axis_title, ratio_ytitle, {sample_label}, "hist e1", plotdir + "/" + outstem + out_suffix);
1691 }
1692
1693 delete hNumProj;
1694 delete hDenProj;
1695 return hRatio;
1696 };
1697
1698 auto draw_projection_stack = [&](TH2D *hInclusive, const std::vector<std::pair<TH2D *, std::string>> &components, bool project_x, const std::string &axis_title, const std::string &out_suffix)
1699 {
1700 if (!hInclusive)
1701 return;
1702
1703 TH1D *hInclusiveProj = make_projection(hInclusive, project_x, std::string(hInclusive->GetName()) + (project_x ? "_projX_stack_inclusive" : "_projY_stack_inclusive"));
1704 if (!hInclusiveProj)
1705 return;
1706
1707 const float alpha = 0.5f;
1708 const std::vector<std::string> fill_colors = {"#ffa90e", "#6CA651", "#832db6", "#bd1f01", "#3f90da"};
1709
1710 TCanvas *cstack = new TCanvas((std::string("c") + out_suffix).c_str(), (std::string("c") + out_suffix).c_str(), 800, 700);
1711 cstack->cd();
1712 gPad->SetLeftMargin(0.14);
1713 gPad->SetRightMargin(0.05);
1714 gPad->SetTopMargin(0.08);
1715
1716 THStack *hs = new THStack((std::string(hInclusive->GetName()) + (project_x ? "_projX_stack" : "_projY_stack")).c_str(), (std::string(";") + axis_title + ";Counts").c_str());
1717 std::vector<TH1D *> component_proj;
1718 std::vector<std::string> component_labels;
1719 component_proj.reserve(components.size());
1720 component_labels.reserve(components.size());
1721
1722 for (size_t i = 0; i < components.size(); ++i)
1723 {
1724 TH2D *h2 = components[i].first;
1725 if (!h2)
1726 continue;
1727
1728 TH1D *hproj = make_projection(h2, project_x, std::string(h2->GetName()) + (project_x ? "_projX_stack_component" : "_projY_stack_component"));
1729 if (!hproj)
1730 continue;
1731
1732 const Color_t color = TColor::GetColor(fill_colors[i % fill_colors.size()].c_str());
1733 hproj->SetFillColorAlpha(color, alpha);
1734 hproj->SetLineColor(color);
1735 hproj->SetLineWidth(1);
1736 hs->Add(hproj);
1737 component_proj.push_back(hproj);
1738 component_labels.push_back(components[i].second);
1739 }
1740
1741 hs->Draw("hist");
1742 hs->GetXaxis()->SetTitle(axis_title.c_str());
1743 hs->GetYaxis()->SetTitle("Counts");
1744 hs->GetYaxis()->SetTitleOffset(1.4);
1745 hs->SetMaximum(hInclusiveProj->GetMaximum() * 1.5);
1746 gPad->Modified();
1747 gPad->Update();
1748
1749 hInclusiveProj->SetLineColor(kBlack);
1750 hInclusiveProj->SetLineWidth(2);
1751 hInclusiveProj->Draw("hist SAME");
1752
1753 TLegend *leg = new TLegend(gPad->GetLeftMargin() + 0.02, 1 - gPad->GetTopMargin() - 0.24, gPad->GetLeftMargin() + 0.34, 1 - gPad->GetTopMargin() - 0.04);
1754 leg->SetBorderSize(0);
1755 leg->SetFillStyle(0);
1756 leg->SetTextSize(0.035);
1757 leg->AddEntry(hInclusiveProj, "Inclusive", "l");
1758 for (size_t i = 0; i < component_proj.size(); ++i)
1759 leg->AddEntry(component_proj[i], component_labels[i].c_str(), "f");
1760 leg->Draw();
1761
1762 TLatex latex;
1763 latex.SetNDC();
1764 latex.SetTextFont(42);
1765 latex.SetTextSize(0.035);
1766 latex.SetTextAlign(kHAlignRight + kVAlignBottom);
1767 latex.DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, sample_label.c_str());
1768
1769 cstack->SaveAs((plotdir + "/" + outstem + out_suffix + ".png").c_str());
1770 cstack->SaveAs((plotdir + "/" + outstem + out_suffix + ".pdf").c_str());
1771
1772 delete leg;
1773 for (auto *h : component_proj)
1774 delete h;
1775 delete hs;
1776 delete hInclusiveProj;
1777 delete cstack;
1778 };
1779
1780 auto draw_projection_stack_fractional = [&](TH2D *hInclusive, const std::vector<std::pair<TH2D *, std::string>> &components, bool project_x, const std::string &axis_title, const std::string &out_suffix)
1781 {
1782 if (!hInclusive)
1783 return;
1784
1785 TH1D *hInclusiveProj = make_projection(hInclusive, project_x, std::string(hInclusive->GetName()) + (project_x ? "_projX_frac_inclusive" : "_projY_frac_inclusive"));
1786 if (!hInclusiveProj)
1787 return;
1788
1789 const float alpha = 0.5f;
1790 const std::vector<std::string> fill_colors = {"#ffa90e", "#6CA651", "#832db6", "#bd1f01", "#3f90da"};
1791
1792 TCanvas *cFrac = new TCanvas((std::string("c") + out_suffix).c_str(), (std::string("c") + out_suffix).c_str(), 800, 700);
1793 cFrac->cd();
1794 gPad->SetLeftMargin(0.14);
1795 gPad->SetRightMargin(0.05);
1796 gPad->SetTopMargin(0.23);
1797
1798 THStack *hsFrac = new THStack((std::string(hInclusive->GetName()) + (project_x ? "_projX_stack_fractional" : "_projY_stack_fractional")).c_str(), (std::string(";") + axis_title + ";Fractional contribution").c_str());
1799 std::vector<TH1D *> component_frac;
1800 std::vector<std::string> component_labels;
1801 component_frac.reserve(components.size());
1802 component_labels.reserve(components.size());
1803
1804 for (size_t i = 0; i < components.size(); ++i)
1805 {
1806 TH2D *h2 = components[i].first;
1807 if (!h2)
1808 continue;
1809
1810 TH1D *hproj = make_projection(h2, project_x, std::string(h2->GetName()) + (project_x ? "_projX_frac_component" : "_projY_frac_component"));
1811 if (!hproj)
1812 continue;
1813
1814 TH1D *hFrac = MakeFractionalHist(hproj, hInclusiveProj, std::string(hproj->GetName()) + "_fractional");
1815 delete hproj;
1816 if (!hFrac)
1817 continue;
1818
1819 const Color_t color = TColor::GetColor(fill_colors[i % fill_colors.size()].c_str());
1820 hFrac->SetFillColorAlpha(color, alpha);
1821 hFrac->SetLineColor(color);
1822 hFrac->SetLineWidth(2);
1823 hsFrac->Add(hFrac);
1824 component_frac.push_back(hFrac);
1825 component_labels.push_back(components[i].second);
1826 }
1827
1828 hsFrac->Draw("hist");
1829 hsFrac->GetXaxis()->SetTitle(axis_title.c_str());
1830 hsFrac->GetYaxis()->SetTitle("Fractional contribution");
1831 hsFrac->GetYaxis()->SetTitleOffset(1.5);
1832 hsFrac->SetMaximum(1.0);
1833 gPad->Modified();
1834 gPad->Update();
1835
1836 TLine *line1 = new TLine(hInclusiveProj->GetXaxis()->GetXmin(), 1.0, hInclusiveProj->GetXaxis()->GetXmax(), 1.0);
1837 line1->SetLineColor(kBlack);
1838 line1->SetLineWidth(2);
1839
1840 TLegend *legFrac = new TLegend(gPad->GetLeftMargin() + 0.02, 1 - gPad->GetTopMargin() + 0.01, gPad->GetLeftMargin() + 0.34, 0.99);
1841 legFrac->SetBorderSize(0);
1842 legFrac->SetFillStyle(0);
1843 legFrac->SetTextSize(0.035);
1844 legFrac->AddEntry(line1, "Inclusive", "l");
1845 for (size_t i = 0; i < component_frac.size(); ++i)
1846 legFrac->AddEntry(component_frac[i], component_labels[i].c_str(), "f");
1847 legFrac->Draw();
1848
1849 TLatex latex;
1850 latex.SetNDC();
1851 latex.SetTextFont(42);
1852 latex.SetTextSize(0.03);
1853 latex.SetTextAlign(kHAlignRight + kVAlignBottom);
1854 latex.DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, sample_label.c_str());
1855
1856 cFrac->SaveAs((plotdir + "/" + outstem + out_suffix + ".png").c_str());
1857 cFrac->SaveAs((plotdir + "/" + outstem + out_suffix + ".pdf").c_str());
1858
1859 delete legFrac;
1860 delete line1;
1861 for (auto *h : component_frac)
1862 delete h;
1863 delete hsFrac;
1864 delete hInclusiveProj;
1865 delete cFrac;
1866 };
1867
1868 draw_input_2d(h2_nraw, "N_{Uncorrected}^{INTT Doublets}", "_h2_nraw");
1869 draw_input_2d(hSilSeed, "N_{Silicon seeds}", "_hSilSeed", "colz text", false, true);
1870 draw_input_2d(hSilSeed_nMVTX3nINTT2, "N_{Silicon seeds}^{nMVTX3nINTT2}", "_hSilSeed_nMVTX3nINTT2", "colz text", false, true);
1871 draw_input_2d(hSilSeed_nMVTX3nINTT1, "N_{Silicon seeds}^{nMVTX3nINTT1}", "_hSilSeed_nMVTX3nINTT1", "colz text", false, true);
1872 draw_input_2d(hSilSeed_nMVTX2nINTT2, "N_{Silicon seeds}^{nMVTX2nINTT2}", "_hSilSeed_nMVTX2nINTT2", "colz text", false, true);
1873 draw_input_2d(hSilSeed_nMVTX2nINTT1, "N_{Silicon seeds}^{nMVTX2nINTT1}", "_hSilSeed_nMVTX2nINTT1", "colz text", false, true);
1874 draw_input_2d(hSilSeed_nMVTXnINTTOther, "N_{Silicon seeds}^{nMVTXnINTTOther}", "_hSilSeed_nMVTXnINTTOther", "colz text", false, true);
1875 draw_input_2d(hSilSeed_good, "N_{Silicon seeds}^{good}", "_hSilSeed_good", "colz text", false, true);
1876 draw_input_2d(hINTTClusterVsAssociatedSilSeed, "Entries", "_hINTTClusterVsAssociatedSilSeed", "colz", false, false, "N_{INTT clusters} (selected crossing)", "N_{Silicon seeds}^{assoc. w. vertex} (selected crossing)", true);
1877 draw_input_2d(hINTTClusterVsAllSilSeed, "Entries", "_hINTTClusterVsAllSilSeed", "colz", false, false, "N_{INTT clusters} (selected crossing)", "N_{Silicon seeds}^{all in crossing} (selected crossing)", true);
1878 drawComparisonWithRatio(std::make_pair(hSeedCrossing_selectedCrossings, "Selected crossings"),
1879 std::make_pair(hSeedCrossing_selectedCrossings_largeDiffNClusterSeeds, "Selected crossings with N_{INTT clusters}>4#timesN_{Silicon seeds}"),
1880 true,
1881 "Crossing",
1882 "Entries",
1883 "Ratio",
1884 {sample_label},
1885 plotdir + "/" + outstem + "_hSeedCrossing_selectedCrossings_comparison"
1886 );
1887 draw_input_2d(hSeed, "N_{Cluster}^{Layer 3+4}", "_hSeed");
1888 if (is_simulation)
1889 draw_input_2d(hPrimaryCharged, "N_{sPHENIX primary charged}", "_hPrimaryCharged");
1890
1891 draw_projection(hSilSeed, false, "#phi [radian]", "Counts", "_hSilSeed_phiProjection");
1892 draw_projection(hSilSeed, true, "#eta", "Counts", "_hSilSeed_etaProjection");
1893 draw_projection(hSilSeed_nMVTX3nINTT2, false, "#phi [radian]", "Counts", "_hSilSeed_nMVTX3nINTT2_phiProjection");
1894 draw_projection(hSilSeed_nMVTX3nINTT2, true, "#eta", "Counts", "_hSilSeed_nMVTX3nINTT2_etaProjection");
1895 draw_projection(hSilSeed_nMVTX3nINTT1, false, "#phi [radian]", "Counts", "_hSilSeed_nMVTX3nINTT1_phiProjection");
1896 draw_projection(hSilSeed_nMVTX3nINTT1, true, "#eta", "Counts", "_hSilSeed_nMVTX3nINTT1_etaProjection");
1897 draw_projection(hSilSeed_nMVTX2nINTT2, false, "#phi [radian]", "Counts", "_hSilSeed_nMVTX2nINTT2_phiProjection");
1898 draw_projection(hSilSeed_nMVTX2nINTT2, true, "#eta", "Counts", "_hSilSeed_nMVTX2nINTT2_etaProjection");
1899 draw_projection(hSilSeed_nMVTX2nINTT1, false, "#phi [radian]", "Counts", "_hSilSeed_nMVTX2nINTT1_phiProjection");
1900 draw_projection(hSilSeed_nMVTX2nINTT1, true, "#eta", "Counts", "_hSilSeed_nMVTX2nINTT1_etaProjection");
1901 draw_projection(hSilSeed_nMVTXnINTTOther, false, "#phi [radian]", "Counts", "_hSilSeed_nMVTXnINTTOther_phiProjection");
1902 draw_projection(hSilSeed_nMVTXnINTTOther, true, "#eta", "Counts", "_hSilSeed_nMVTXnINTTOther_etaProjection");
1903 draw_projection_stack(hSilSeed,
1904 {
1905
1906 {hSilSeed_nMVTX3nINTT2, "MVTX+INTT: 3+2"},
1907 {hSilSeed_nMVTX3nINTT1, "MVTX+INTT: 3+1"},
1908 {hSilSeed_nMVTX2nINTT2, "MVTX+INTT: 2+2"},
1909 {hSilSeed_nMVTX2nINTT1, "MVTX+INTT: 2+1"},
1910 {hSilSeed_nMVTXnINTTOther, "MVTX+INTT: Other"},
1911 },
1912 false,
1913 "#phi [radian]",
1914 "_hSilSeed_phiProjection_stack"
1915 );
1916 draw_projection_stack(hSilSeed,
1917 {
1918
1919 {hSilSeed_nMVTX3nINTT2, "MVTX+INTT: 3+2"},
1920 {hSilSeed_nMVTX3nINTT1, "MVTX+INTT: 3+1"},
1921 {hSilSeed_nMVTX2nINTT2, "MVTX+INTT: 2+2"},
1922 {hSilSeed_nMVTX2nINTT1, "MVTX+INTT: 2+1"},
1923 {hSilSeed_nMVTXnINTTOther, "MVTX+INTT: Other"},
1924 },
1925 true,
1926 "#eta",
1927 "_hSilSeed_etaProjection_stack"
1928 );
1929 draw_projection_stack_fractional(hSilSeed,
1930 {
1931
1932 {hSilSeed_nMVTX3nINTT2, "MVTX+INTT: 3+2"},
1933 {hSilSeed_nMVTX3nINTT1, "MVTX+INTT: 3+1"},
1934 {hSilSeed_nMVTX2nINTT2, "MVTX+INTT: 2+2"},
1935 {hSilSeed_nMVTX2nINTT1, "MVTX+INTT: 2+1"},
1936 {hSilSeed_nMVTXnINTTOther, "MVTX+INTT: Other"},
1937 },
1938 false,
1939 "#phi [radian]",
1940 "_hSilSeed_phiProjection_stack_fractional"
1941 );
1942 draw_projection_stack_fractional(hSilSeed,
1943 {
1944
1945 {hSilSeed_nMVTX3nINTT2, "MVTX+INTT: 3+2"},
1946 {hSilSeed_nMVTX3nINTT1, "MVTX+INTT: 3+1"},
1947 {hSilSeed_nMVTX2nINTT2, "MVTX+INTT: 2+2"},
1948 {hSilSeed_nMVTX2nINTT1, "MVTX+INTT: 2+1"},
1949 {hSilSeed_nMVTXnINTTOther, "MVTX+INTT: Other"},
1950 },
1951 true,
1952 "#eta",
1953 "_hSilSeed_etaProjection_stack_fractional"
1954 );
1955
1956 if (hDCA3D_Sig && hDCA3D_Bkg)
1957 {
1958 TH1D *hDCA3D_Sub = dynamic_cast<TH1D *>(hDCA3D_Sig->Clone((prefix + "_DCA3D_Sub").c_str()));
1959 if (hDCA3D_Sub && hSilSeedDCA3D)
1960 {
1961 hDCA3D_Sub->SetDirectory(nullptr);
1962 hDCA3D_Sub->Add(hDCA3D_Bkg, -1.0);
1963
1964 double dca_ymin = std::min({hDCA3D_Sig->GetMinimum(0), hDCA3D_Bkg->GetMinimum(0), hDCA3D_Sub->GetMinimum(0), hSilSeedDCA3D->GetMinimum(0)});
1965 double dca_ymax = std::max({hDCA3D_Sig->GetMaximum(), hDCA3D_Bkg->GetMaximum(), hDCA3D_Sub->GetMaximum(), hSilSeedDCA3D->GetMaximum()});
1966 if (dca_ymax <= dca_ymin)
1967 dca_ymax = dca_ymin + 1.0;
1968 const double dca_margin = 0.15 * (dca_ymax - dca_ymin);
1969
1970 draw_comparison_alt({hDCA3D_Sig, hDCA3D_Bkg, hDCA3D_Sub},
1971 {"#1f77b4", "#67B2D8", "#000000"},
1972 {24, 25, 20},
1973 {0.5, 0.5, 1},
1974 {"INTT doublets - Unrotated", "INTT doublets - Rotated", "INTT doublets - Background subtracted"},
1975 "DCA_{w.r.t vertex} [cm]",
1976 "Counts",
1977 {0, std::max({hDCA3D_Sig->GetMaximum(), hDCA3D_Bkg->GetMaximum(), hDCA3D_Sub->GetMaximum()})*1.3},
1978 false,
1979 false,
1980 510,
1981 {false, 0.f},
1982 plotdir + "/" + outstem + "_hDCA3D_comparison"
1983 );
1984
1985 draw_comparison_alt({hDCA3D_Sig, hDCA3D_Bkg, hDCA3D_Sub, hSilSeedDCA3D},
1986 {"#1f77b4", "#67B2D8", "#000000", "#bd1f01"},
1987 {24, 25, 20, 21},
1988 {0.5, 0.5, 1, 1},
1989 {"INTT doublets - Unrotated", "INTT doublets - Rotated", "INTT doublets - Background subtracted", "Silicon seeds"},
1990 "DCA_{w.r.t vertex} [cm]",
1991 "Counts",
1992 {-1, -1},
1993 true,
1994 false,
1995 510,
1996 {false, 0.f},
1997 plotdir + "/" + outstem + "_hDCA3D_SilSeed_comparison"
1998 );
1999
2000 }
2001 }
2002
2003 if (hSilSeedDCA3D)
2004 draw1Dhistogram(hSilSeedDCA3D, true, "DCA_{Silicon seeds}^{w.r.t vertex} [cm]", "Counts", {sample_label}, "hist e1", plotdir + "/" + outstem + "_hSilSeedDCA3D");
2005
2006 if (hSilSeedDCA3D_truthVtxZ_2D)
2007 {
2008 const double heatmap_zmax = std::max(0.05, hSilSeedDCA3D_truthVtxZ_2D->GetMaximum() * 1.1);
2009 hSilSeedDCA3D_truthVtxZ_2D->GetYaxis()->SetRangeUser(0, 3);
2010 draw2Dhistogram_adjZaxis(hSilSeedDCA3D_truthVtxZ_2D,
2011 true,
2012 "Truth vertex z [cm]",
2013 "DCA_{Silicon seeds}^{w.r.t vertex} [cm]",
2014 "Normalized entries",
2015 std::pair<double, double>(0.0, heatmap_zmax),
2016 {sample_label},
2017 "colz",
2018 ".3g",
2019 false,
2020 false,
2021 plotdir + "/" + outstem + "_hSilSeedDCA3D_truthVtxZ_2D");
2022 }
2023
2024 if (!hSilSeedDCA3D_truthVtxZ_absMerged.empty())
2025 {
2026 std::vector<TH1 *> hists;
2027 std::vector<std::string> labels;
2028 hists.reserve(hSilSeedDCA3D_truthVtxZ_absMerged.size());
2029 labels.reserve(hSilSeedDCA3D_truthVtxZ_absMerged.size());
2030 for (const auto &entry : hSilSeedDCA3D_truthVtxZ_absMerged)
2031 {
2032 hists.push_back(entry.first);
2033 labels.push_back(entry.second);
2034 }
2035
2036 draw_comparison(hists,
2037 {"#4477aa", "#228833", "#ccbb44", "#ee8866", "#cc3311"},
2038 labels,
2039 "DCA_{Silicon seeds}^{w.r.t vertex} [cm]",
2040 "Normalized counts",
2041 {-1, -1},
2042 true,
2043 true,
2044 510,
2045 std::make_pair(false, 0.0f),
2046 plotdir + "/" + outstem + "_hSilSeedDCA3D_truthVtxZAbs_comparison");
2047 }
2048
2049
2050 TH2D *h2_ratio_nraw_over_silseed = build_and_draw_ratio(h2_nraw, hSilSeed,
2051 prefix + "_NRawOverSilSeedEtaPhi",
2052 "N_{Uncorrected}^{INTT Doublets} / N_{Silicon seeds}",
2053 "_ratio_nraw_over_silseed");
2054 TH2D *h2_ratio_nraw_over_silseed_good = build_and_draw_ratio(h2_nraw, hSilSeed_good,
2055 prefix + "_NRawOverSilSeedGoodEtaPhi",
2056 "N_{Uncorrected}^{INTT Doublets} / N_{Silicon seeds}^{good}",
2057 "_ratio_nraw_over_silseed_good");
2058
2059 ratio_outputs.eta_ratio = build_projection_ratio(std::make_pair(h2_nraw, "N_{Uncorrected}^{INTT Doublets}"),
2060 std::make_pair(hSilSeed, "N_{Silicon seeds}"),
2061 true,
2062 "#eta",
2063 "N_{Uncorrected}^{INTT Doublets} / N_{Silicon seeds}",
2064 "_ratio1D_nraw_over_silseed_eta",
2065 prefix + "_NRawOverSilSeedEtaProjection"
2066 );
2067 ratio_outputs.phi_ratio = build_projection_ratio(std::make_pair(h2_nraw, "N_{Uncorrected}^{INTT Doublets}"),
2068 std::make_pair(hSilSeed, "N_{Silicon seeds}"),
2069 false,
2070 "#phi [radian]",
2071 "N_{Uncorrected}^{INTT Doublets} / N_{Silicon seeds}",
2072 "_ratio1D_nraw_over_silseed_phi",
2073 prefix + "_NRawOverSilSeedPhiProjection"
2074 );
2075 ratio_outputs.eta_ratio_good = build_projection_ratio(std::make_pair(h2_nraw, "N_{Uncorrected}^{INTT Doublets}"),
2076 std::make_pair(hSilSeed_good, "N_{Silicon seeds}^{good}"),
2077 true,
2078 "#eta",
2079 "N_{Uncorrected}^{INTT Doublets} / N_{Silicon seeds}^{good}",
2080 "_ratio1D_nraw_over_silseed_good_eta",
2081 prefix + "_NRawOverSilSeedGoodEtaProjection"
2082 );
2083 ratio_outputs.phi_ratio_good = build_projection_ratio(std::make_pair(h2_nraw, "N_{Uncorrected}^{INTT Doublets}"),
2084 std::make_pair(hSilSeed_good, "N_{Silicon seeds}^{good}"),
2085 false,
2086 "#phi [radian]",
2087 "N_{Uncorrected}^{INTT Doublets} / N_{Silicon seeds}^{good}",
2088 "_ratio1D_nraw_over_silseed_good_phi",
2089 prefix + "_NRawOverSilSeedGoodPhiProjection"
2090 );
2091
2092
2093
2094 TH2D *h2_ratio_primary_over_nraw = nullptr;
2095 TH2D *h2_ratio_primary_over_silseed = nullptr;
2096 if (is_simulation)
2097 {
2098 h2_ratio_primary_over_nraw = build_and_draw_ratio(hPrimaryCharged, h2_nraw,
2099 prefix + "_PrimaryChargedOverNRawEtaPhi",
2100 "N_{sPHENIX primary charged} / N_{raw}^{INTT Doublets}",
2101 "_ratio_primarycharged_over_nraw");
2102 h2_ratio_primary_over_silseed = build_and_draw_ratio(hPrimaryCharged, hSilSeed,
2103 prefix + "_PrimaryChargedOverSilSeedEtaPhi",
2104 "N_{sPHENIX primary charged} / N_{Silicon seeds}",
2105 "_ratio_primarycharged_over_silseed");
2106 }
2107
2108 TH2D *h2_ratio_silseed_nMVTX3nINTT2_over_inclusive = build_and_draw_ratio(hSilSeed_nMVTX3nINTT2, hSilSeed,
2109 prefix + "_SilSeed_nMVTX3nINTT2_OverInclusive",
2110 "N_{SilSeed}^{nMVTX3nINTT2} / N_{SilSeed}^{Inclusive}",
2111 "_ratio_silseed_nMVTX3nINTT2_over_inclusive",
2112 false,
2113 ".2g");
2114 TH2D *h2_ratio_silseed_nMVTX3nINTT1_over_inclusive = build_and_draw_ratio(hSilSeed_nMVTX3nINTT1, hSilSeed,
2115 prefix + "_SilSeed_nMVTX3nINTT1_OverInclusive",
2116 "N_{SilSeed}^{nMVTX3nINTT1} / N_{SilSeed}^{Inclusive}",
2117 "_ratio_silseed_nMVTX3nINTT1_over_inclusive",
2118 false,
2119 ".2g");
2120 TH2D *h2_ratio_silseed_nMVTX2nINTT2_over_inclusive = build_and_draw_ratio(hSilSeed_nMVTX2nINTT2, hSilSeed,
2121 prefix + "_SilSeed_nMVTX2nINTT2_OverInclusive",
2122 "N_{SilSeed}^{nMVTX2nINTT2} / N_{SilSeed}^{Inclusive}",
2123 "_ratio_silseed_nMVTX2nINTT2_over_inclusive",
2124 false,
2125 ".2g");
2126 TH2D *h2_ratio_silseed_nMVTX2nINTT1_over_inclusive = build_and_draw_ratio(hSilSeed_nMVTX2nINTT1, hSilSeed,
2127 prefix + "_SilSeed_nMVTX2nINTT1_OverInclusive",
2128 "N_{SilSeed}^{nMVTX2nINTT1} / N_{SilSeed}^{Inclusive}",
2129 "_ratio_silseed_nMVTX2nINTT1_over_inclusive",
2130 false,
2131 ".2g");
2132 TH2D *h2_ratio_silseed_nMVTXnINTTOther_over_inclusive = build_and_draw_ratio(hSilSeed_nMVTXnINTTOther, hSilSeed,
2133 prefix + "_SilSeed_nMVTXnINTTOther_OverInclusive",
2134 "N_{SilSeed}^{nMVTXnINTTOther} / N_{SilSeed}^{Inclusive}",
2135 "_ratio_silseed_nMVTXnINTTOther_over_inclusive",
2136 false,
2137 ".2g");
2138
2139 c->SaveAs((plotdir + "/" + outstem + ".pdf").c_str());
2140
2141
2142
2143 for (auto *f : frames)
2144 delete f;
2145 for (auto *ln : cutLines)
2146 delete ln;
2147 delete h2_ratio_nraw_over_silseed;
2148 delete h2_ratio_nraw_over_silseed_good;
2149 delete h2_ratio_silseed_nMVTX3nINTT2_over_inclusive;
2150 delete h2_ratio_silseed_nMVTX3nINTT1_over_inclusive;
2151 delete h2_ratio_silseed_nMVTX2nINTT2_over_inclusive;
2152 delete h2_ratio_silseed_nMVTX2nINTT1_over_inclusive;
2153 delete h2_ratio_silseed_nMVTXnINTTOther_over_inclusive;
2154 if (h2_ratio_primary_over_nraw)
2155 delete h2_ratio_primary_over_nraw;
2156 if (h2_ratio_primary_over_silseed)
2157 delete h2_ratio_primary_over_silseed;
2158 delete hSilSeedDCA3D_truthVtxZ_2D;
2159 for (auto &entry : hSilSeedDCA3D_truthVtxZ_absMerged)
2160 delete entry.first;
2161 delete hSilSeed_good;
2162 delete h2_nraw;
2163 fin->Close();
2164 delete c;
2165 return ratio_outputs;
2166 }
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185 void plot_tklCombinatoric(const std::string &inputfile = "/sphenix/tg/tg01/hf/hjheng/ppg-dNdEta-OOpp/TrackletAna-HistOutput/simulation/ACTS/cotThetaMax2p9/merged_hist.root",
2186 const std::string &prefix = "tkl_Combinatoric",
2187 const std::string &mode = "absdPhi",
2188 bool is_simulation = true,
2189 double zvtx_cut_min = -10.0,
2190 double zvtx_cut_max = 10.0,
2191 double raw_count_xmin = 0.0,
2192 double raw_count_xmax = kDirectRawCountCut,
2193 double plot_xmin = 0.0,
2194 double plot_xmax = 0.2,
2195 const std::string &outstem = "plot_tklCombinatoric",
2196 const std::string &plotdir = "./figure/figure-tklCombinatoric-simulation-cotThetaMax2p9",
2197 bool zero_x_error = true,
2198 const std::string &ratio_draw_option = "colz text",
2199 const std::string &ratio_text_format = ".3f",
2200 bool use_rounded_zmax = true
2201 )
2202 {
2203 TFile *fin = TFile::Open(inputfile.c_str(), "READ");
2204 if (!fin || fin->IsZombie())
2205 {
2206 std::cerr << "[plot_tklCombinatoric] Cannot open: " << inputfile << "\n";
2207 delete fin;
2208 return;
2209 }
2210
2211 const std::vector<std::string> prefixes = CollectAvailablePrefixes(fin, prefix);
2212 const auto percentile_boundaries = LoadPercentileBoundariesByBin(is_simulation);
2213 fin->Close();
2214 delete fin;
2215
2216 const std::vector<std::string> colors = {"#4477aa", "#228833", "#ccbb44", "#ee8866", "#cc3311", "#aa3377"};
2217 std::vector<ProjectionRatioOutputs> ratio_sets;
2218 ratio_sets.reserve(prefixes.size());
2219
2220 for (const auto ¤t_prefix : prefixes)
2221 {
2222 const std::string suffix = PrefixSuffixForOutput(prefix, current_prefix);
2223 ProjectionRatioOutputs outputs = plot_tklCombinatoric_single(inputfile, current_prefix, mode, is_simulation, zvtx_cut_min, zvtx_cut_max, raw_count_xmin, raw_count_xmax, plot_xmin, plot_xmax, outstem + suffix, plotdir, zero_x_error, ratio_draw_option, ratio_text_format, use_rounded_zmax, PrefixDisplayLabel(prefix, current_prefix, percentile_boundaries));
2224 if (outputs.eta_ratio || outputs.phi_ratio || outputs.eta_ratio_good || outputs.phi_ratio_good)
2225 ratio_sets.push_back(outputs);
2226 }
2227
2228 std::vector<TH1 *> eta_hists;
2229 std::vector<TH1 *> phi_hists;
2230 std::vector<TH1 *> eta_hists_good;
2231 std::vector<TH1 *> phi_hists_good;
2232 std::vector<TH1 *> eta_hists_inverse;
2233 std::vector<TH1 *> phi_hists_inverse;
2234 std::vector<TH1 *> eta_hists_good_inverse;
2235 std::vector<TH1 *> phi_hists_good_inverse;
2236 std::vector<std::string> labels;
2237 std::vector<std::string> labels_good;
2238 std::vector<std::string> comparison_colors;
2239 std::vector<std::string> comparison_colors_good;
2240 for (size_t i = 0; i < ratio_sets.size(); ++i)
2241 {
2242 if (!ratio_sets[i].eta_ratio || !ratio_sets[i].phi_ratio)
2243 continue;
2244 eta_hists.push_back(ratio_sets[i].eta_ratio);
2245 phi_hists.push_back(ratio_sets[i].phi_ratio);
2246 if (TH1D *hInv = MakeInverseHist(ratio_sets[i].eta_ratio, ratio_sets[i].prefix + "_NRawOverSilSeedEtaProjection_inverse"))
2247 eta_hists_inverse.push_back(hInv);
2248 if (TH1D *hInv = MakeInverseHist(ratio_sets[i].phi_ratio, ratio_sets[i].prefix + "_NRawOverSilSeedPhiProjection_inverse"))
2249 phi_hists_inverse.push_back(hInv);
2250 if (ratio_sets[i].eta_ratio_good && ratio_sets[i].phi_ratio_good)
2251 {
2252 eta_hists_good.push_back(ratio_sets[i].eta_ratio_good);
2253 phi_hists_good.push_back(ratio_sets[i].phi_ratio_good);
2254 if (TH1D *hInv = MakeInverseHist(ratio_sets[i].eta_ratio_good, ratio_sets[i].prefix + "_NRawOverSilSeedGoodEtaProjection_inverse"))
2255 eta_hists_good_inverse.push_back(hInv);
2256 if (TH1D *hInv = MakeInverseHist(ratio_sets[i].phi_ratio_good, ratio_sets[i].prefix + "_NRawOverSilSeedGoodPhiProjection_inverse"))
2257 phi_hists_good_inverse.push_back(hInv);
2258 labels_good.push_back(ratio_sets[i].label);
2259 comparison_colors_good.push_back((i == 0) ? "#000000" : colors[(i - 1) % colors.size()]);
2260 }
2261 labels.push_back(ratio_sets[i].label);
2262 comparison_colors.push_back((i == 0) ? "#000000" : colors[(i - 1) % colors.size()]);
2263 }
2264
2265 if (!eta_hists.empty())
2266 {
2267 draw_comparison(eta_hists,
2268 comparison_colors,
2269 labels,
2270 "#eta",
2271 "N_{Uncorrected}^{INTT Doublets}/N_{Silicon seeds}",
2272 {0, (!is_simulation) ? 2 : 3},
2273 false,
2274 false,
2275 510,
2276 std::make_pair(true, 1.0),
2277 plotdir + "/" + outstem + "_ratio1D_nraw_over_seed_eta_comparison"
2278 );
2279
2280 draw_comparison(eta_hists_inverse,
2281 comparison_colors,
2282 labels,
2283 "#eta",
2284 "N_{Silicon seeds}/N_{Uncorrected}^{INTT Doublets}",
2285 {0, (!is_simulation) ? 5 : 2},
2286 false,
2287 false,
2288 510,
2289 std::make_pair(true, 1.0),
2290 plotdir + "/" + outstem + "_ratio1D_seed_over_nraw_eta_comparison"
2291 );
2292 }
2293
2294 if (!phi_hists.empty())
2295 {
2296 draw_comparison(phi_hists,
2297 comparison_colors,
2298 labels,
2299 "#phi [radian]",
2300 "N_{Uncorrected}^{INTT Doublets}/N_{Silicon seeds}",
2301 {0, (!is_simulation) ? 5 : 1.7},
2302 false,
2303 false,
2304 510,
2305 std::make_pair(true, 1.0),
2306 plotdir + "/" + outstem + "_ratio1D_nraw_over_seed_phi_comparison"
2307 );
2308
2309 draw_comparison(phi_hists_inverse,
2310 comparison_colors,
2311 labels,
2312 "#phi [radian]",
2313 "N_{Silicon seeds}/N_{Uncorrected}^{INTT Doublets}",
2314 {0, (!is_simulation) ? 5 : 2},
2315 false,
2316 false,
2317 510,
2318 std::make_pair(true, 1.0),
2319 plotdir + "/" + outstem + "_ratio1D_seed_over_nraw_phi_comparison"
2320 );
2321 }
2322
2323 if (!eta_hists_good.empty())
2324 {
2325 draw_comparison(eta_hists_good,
2326 comparison_colors_good,
2327 labels_good,
2328 "#eta",
2329 "N_{Uncorrected}^{INTT Doublets}/N_{Silicon seeds}^{good}",
2330 {0, (!is_simulation) ? 2 : 3},
2331 false,
2332 false,
2333 510,
2334 std::make_pair(true, 1.0),
2335 plotdir + "/" + outstem + "_ratio1D_nraw_over_seed_good_eta_comparison"
2336 );
2337
2338 draw_comparison(eta_hists_good_inverse,
2339 comparison_colors_good,
2340 labels_good,
2341 "#eta",
2342 "N_{Silicon seeds}^{good}/N_{Uncorrected}^{INTT Doublets}",
2343 {0, (!is_simulation) ? 2 : 2},
2344 false,
2345 false,
2346 510,
2347 std::make_pair(true, 1.0),
2348 plotdir + "/" + outstem + "_ratio1D_seed_good_over_nraw_eta_comparison"
2349 );
2350 }
2351
2352 if (!phi_hists_good.empty())
2353 {
2354 draw_comparison(phi_hists_good,
2355 comparison_colors_good,
2356 labels_good,
2357 "#phi [radian]",
2358 "N_{Uncorrected}^{INTT Doublets}/N_{Silicon seeds}^{good}",
2359 {0, (!is_simulation) ? 5 : 1.7},
2360 false,
2361 false,
2362 510,
2363 std::make_pair(true, 1.0),
2364 plotdir + "/" + outstem + "_ratio1D_nraw_over_seed_good_phi_comparison"
2365 );
2366
2367 draw_comparison(phi_hists_good_inverse,
2368 comparison_colors_good,
2369 labels_good,
2370 "#phi [radian]",
2371 "N_{Silicon seeds}^{good}/N_{Uncorrected}^{INTT Doublets}",
2372 {0, (!is_simulation) ? 5 : 2},
2373 false,
2374 false,
2375 510,
2376 std::make_pair(true, 1.0),
2377 plotdir + "/" + outstem + "_ratio1D_seed_good_over_nraw_phi_comparison"
2378 );
2379 }
2380
2381 for (auto &outputs : ratio_sets)
2382 {
2383 delete outputs.eta_ratio;
2384 delete outputs.phi_ratio;
2385 delete outputs.eta_ratio_good;
2386 delete outputs.phi_ratio_good;
2387 }
2388
2389 for (auto *hist : eta_hists_inverse)
2390 delete hist;
2391 for (auto *hist : phi_hists_inverse)
2392 delete hist;
2393 for (auto *hist : eta_hists_good_inverse)
2394 delete hist;
2395 for (auto *hist : phi_hists_good_inverse)
2396 delete hist;
2397 }
