File indexing completed on 2025-12-16 09:17:31
0001 #include <iostream>
0002 #include <fstream>
0003 #include <TFile.h>
0004 #include <TH1D.h>
0005 #include <TH2D.h>
0006 #include <TMath.h>
0007 #include <vector>
0008 #include "TLatex.h"
0009 #include <string>
0010 #include "unfold_Def.h"
0011 #include "/sphenix/user/hanpuj/CaloDataAna24_skimmed/src/draw_template.C"
0012
0013 void do_normalization(TH1D* h_unfolded, double luminosity, TH1D* h_etacorrection, float jet_radius);
0014 TH1D* get_statistical_uncertainty(TH1D* h_unfold, string name);
0015 void get_uncertainty(TH1D* &h_uncertainty_up, TH1D* &h_uncertainty_down, TH1D* h_nominal, TH1D* h_shiftup, TH1D* h_shiftdown, string name);
0016 void get_unfold_uncertainty(TH1D* &h_uncertainty_up, TH1D* &h_uncertainty_down, TH1D* h_nominal, TH1D* h_shift, string name);
0017 void get_total_uncertainty(TH1D* &h_total_up, TH1D* &h_total_down,TH1D* h_statistical, std::vector<TH1D*> h_systematic_up, std::vector<TH1D*> h_systematic_down);
0018 void draw_total_uncertainty(TH1D* h_total_up, TH1D* h_total_down, TH1D* h_statistical, std::vector<TH1D*> h_systematic_up, std::vector<TH1D*> h_systematic_down, std::vector<int> color, std::vector<string> legend, float jet_radius);
0019
0020 void get_finalspectrum(float jet_radius = 0.4) {
0021
0022 ifstream f_luminosity("luminosity.txt");
0023 if (!f_luminosity.is_open()) {
0024 std::cout << "Failed to open luminosity file!" << std::endl;
0025 return;
0026 }
0027 double luminosity22;
0028 f_luminosity >> luminosity22;
0029 std::cout << "Luminosity: " << luminosity22 << std::endl;
0030
0031 TFile* f_spectrum = new TFile(Form("output_unfolded_r0%d.root", (int)(10*jet_radius)), "READ");
0032 if (!f_spectrum->IsOpen()) {
0033 std::cout << "Failed to open unfolded spectrum file!" << std::endl;
0034 return;
0035 }
0036 TFile* f_etacorrection = new TFile("input_etacorrection.root", "READ");
0037 if (!f_etacorrection->IsOpen()) {
0038 std::cout << "Failed to open eta correction file!" << std::endl;
0039 return;
0040 }
0041 TH1D* h_etacorrection = (TH1D*)f_etacorrection->Get("h_etacorrection");
0042 if (!h_etacorrection) {
0043 std::cout << "Failed to get eta correction histogram!" << std::endl;
0044 return;
0045 }
0046
0047 TH1D* h_unfold_nominal = (TH1D*)f_spectrum->Get("h_unfold_nominal"); do_normalization(h_unfold_nominal, luminosity22, h_etacorrection, jet_radius);
0048 TH1D* h_unfold_jesup = (TH1D*)f_spectrum->Get("h_unfold_jesup"); do_normalization(h_unfold_jesup, luminosity22, h_etacorrection, jet_radius);
0049 TH1D* h_unfold_jesdown = (TH1D*)f_spectrum->Get("h_unfold_jesdown"); do_normalization(h_unfold_jesdown, luminosity22, h_etacorrection, jet_radius);
0050 TH1D* h_unfold_jerup = (TH1D*)f_spectrum->Get("h_unfold_jerup"); do_normalization(h_unfold_jerup, luminosity22, h_etacorrection, jet_radius);
0051 TH1D* h_unfold_jerdown = (TH1D*)f_spectrum->Get("h_unfold_jerdown"); do_normalization(h_unfold_jerdown, luminosity22, h_etacorrection, jet_radius);
0052 TH1D* h_unfold_jettrigup = (TH1D*)f_spectrum->Get("h_unfold_jettrigup"); do_normalization(h_unfold_jettrigup, luminosity22, h_etacorrection, jet_radius);
0053 TH1D* h_unfold_jettrigdown = (TH1D*)f_spectrum->Get("h_unfold_jettrigdown"); do_normalization(h_unfold_jettrigdown, luminosity22, h_etacorrection, jet_radius);
0054 TH1D* h_unfold_jettimingup = (TH1D*)f_spectrum->Get("h_unfold_jettimingup"); do_normalization(h_unfold_jettimingup, luminosity22, h_etacorrection, jet_radius);
0055 TH1D* h_unfold_jettimingdown = (TH1D*)f_spectrum->Get("h_unfold_jettimingdown"); do_normalization(h_unfold_jettimingdown, luminosity22, h_etacorrection, jet_radius);
0056 TH1D* h_unfold_unfolditerup = (TH1D*)f_spectrum->Get("h_unfold_unfolditerup"); do_normalization(h_unfold_unfolditerup, luminosity22, h_etacorrection, jet_radius);
0057 TH1D* h_unfold_unfolditerdown = (TH1D*)f_spectrum->Get("h_unfold_unfolditerdown"); do_normalization(h_unfold_unfolditerdown, luminosity22, h_etacorrection, jet_radius);
0058 TH1D* h_uncertainty_stat = get_statistical_uncertainty(h_unfold_nominal, "h_uncertainty_stat");
0059 TH1D* h_uncertainty_jesup,* h_uncertainty_jesdown; get_uncertainty(h_uncertainty_jesup, h_uncertainty_jesdown, h_unfold_nominal, h_unfold_jesup, h_unfold_jesdown, "h_uncertainty_jes");
0060 TH1D* h_uncertainty_jerup,* h_uncertainty_jerdown; get_uncertainty(h_uncertainty_jerup, h_uncertainty_jerdown, h_unfold_nominal, h_unfold_jerup, h_unfold_jerdown, "h_uncertainty_jer");
0061 TH1D* h_uncertainty_jettrigup,* h_uncertainty_jettrigdown; get_uncertainty(h_uncertainty_jettrigup, h_uncertainty_jettrigdown, h_unfold_nominal, h_unfold_jettrigup, h_unfold_jettrigdown, "h_uncertainty_jettrig");
0062 TH1D* h_uncertainty_jettimingup,* h_uncertainty_jettimingdown; get_uncertainty(h_uncertainty_jettimingup, h_uncertainty_jettimingdown, h_unfold_nominal, h_unfold_jettimingup, h_unfold_jettimingdown, "h_uncertainty_jettiming");
0063 TH1D* h_uncertainty_unfolditerup,* h_uncertainty_unfolditerdown; get_uncertainty(h_uncertainty_unfolditerup, h_uncertainty_unfolditerdown, h_unfold_nominal, h_unfold_unfolditerup, h_unfold_unfolditerdown, "h_uncertainty_unfolditer");
0064 std::vector<TH1D*> h_uncertainty_up = {h_uncertainty_jesup, h_uncertainty_jerup, h_uncertainty_jettrigup, h_uncertainty_jettimingup, h_uncertainty_unfolditerup};
0065 std::vector<TH1D*> h_uncertainty_down = {h_uncertainty_jesdown, h_uncertainty_jerdown, h_uncertainty_jettrigdown, h_uncertainty_jettimingdown, h_uncertainty_unfolditerdown};
0066 TH1D* h_uncertainty_total_up, *h_uncertainty_total_down; get_total_uncertainty(h_uncertainty_total_up, h_uncertainty_total_down, h_uncertainty_stat, h_uncertainty_up, h_uncertainty_down);
0067
0068 std::vector<int> color;
0069 std::vector<string> legend;
0070 color = {kRed, kYellow+1, kAzure+2, kGreen+2, kOrange-2, kMagenta};
0071 legend = {"Jet Energy Scale", "Jet Energy Resolution", "Jet Trigger Efficiency", "Jet Timing Cut Efficiency", "Unfolding Iterations"};
0072 draw_total_uncertainty(h_uncertainty_total_up, h_uncertainty_total_down, h_uncertainty_stat, h_uncertainty_up, h_uncertainty_down, color, legend, jet_radius);
0073 color.clear();
0074 legend.clear();
0075
0076
0077 double g_sphenix_x[truthnpt], g_sphenix_xerrdown[truthnpt], g_sphenix_xerrup[truthnpt], g_sphenix_y[truthnpt], g_sphenix_yerrdown[truthnpt], g_sphenix_yerrup[truthnpt], g_sphenix_stat[truthnpt];
0078 for (int ib = 0; ib < truthnpt; ++ib) {
0079 g_sphenix_x[ib] = (truthptbins[ib] + truthptbins[ib+1]) / 2.;
0080 g_sphenix_xerrdown[ib] = (truthptbins[ib+1] - truthptbins[ib]) / 2.;
0081 g_sphenix_xerrup[ib] = (truthptbins[ib+1] - truthptbins[ib]) / 2.;
0082 g_sphenix_y[ib] = h_unfold_nominal->GetBinContent(ib + 1);
0083 g_sphenix_yerrdown[ib] = -1 * h_uncertainty_total_down->GetBinContent(ib + 1) * h_unfold_nominal->GetBinContent(ib + 1);
0084 g_sphenix_yerrup[ib] = h_uncertainty_total_up->GetBinContent(ib + 1) * h_unfold_nominal->GetBinContent(ib + 1);
0085 g_sphenix_stat[ib] = h_uncertainty_stat->GetBinContent(ib + 1) * h_unfold_nominal->GetBinContent(ib + 1);
0086 std::cout << "sPHENIX pt bin: " << g_sphenix_x[ib] << " " << g_sphenix_y[ib] << " " << g_sphenix_yerrdown[ib] << " " << g_sphenix_yerrup[ib] << " " << g_sphenix_stat[ib] << std::endl;
0087 }
0088 TGraphAsymmErrors* g_sphenix_result = new TGraphAsymmErrors(truthnpt, g_sphenix_x, g_sphenix_y, g_sphenix_xerrdown, g_sphenix_xerrup, g_sphenix_yerrdown, g_sphenix_yerrup);
0089 g_sphenix_result->SetMarkerStyle(20);
0090 g_sphenix_result->SetMarkerSize(1.5);
0091 g_sphenix_result->SetMarkerColor(kAzure+2);
0092 g_sphenix_result->SetLineColor(kAzure+2);
0093 g_sphenix_result->SetLineWidth(0);
0094 g_sphenix_result->SetFillColorAlpha(kAzure + 1, 0.5);
0095 TGraphErrors* g_sphenix_result_point = new TGraphErrors(truthnpt, g_sphenix_x, g_sphenix_y, g_sphenix_xerrdown, g_sphenix_stat);
0096 g_sphenix_result_point->SetMarkerStyle(20);
0097 g_sphenix_result_point->SetMarkerSize(1.5);
0098 g_sphenix_result_point->SetMarkerColor(kAzure+2);
0099 g_sphenix_result_point->SetLineColor(kAzure+2);
0100 g_sphenix_result_point->SetLineWidth(2);
0101
0102
0103 double g_sphenix_qm_x[9], g_sphenix_qm_xerrdown[9], g_sphenix_qm_xerrup[9], g_sphenix_qm_y[9], g_sphenix_qm_yerrdown[9], g_sphenix_qm_yerrup[9], g_sphenix_qm_stat[9];
0104 ifstream file_qm("sphenix_qm_result_point.txt");
0105 for (int ib = 0; ib < 9; ++ib) {
0106 file_qm >> g_sphenix_qm_x[ib] >> g_sphenix_qm_xerrdown[ib] >> g_sphenix_qm_xerrup[ib] >> g_sphenix_qm_y[ib] >> g_sphenix_qm_yerrdown[ib] >> g_sphenix_qm_yerrup[ib] >> g_sphenix_qm_stat[ib];
0107 }
0108 TGraphAsymmErrors* g_sphenix_qm_result = new TGraphAsymmErrors(9, g_sphenix_qm_x, g_sphenix_qm_y, g_sphenix_qm_xerrdown, g_sphenix_qm_xerrup, g_sphenix_qm_yerrdown, g_sphenix_qm_yerrup);
0109 g_sphenix_qm_result->SetMarkerStyle(20);
0110 g_sphenix_qm_result->SetMarkerSize(1.5);
0111 g_sphenix_qm_result->SetMarkerColor(kRed+1);
0112 g_sphenix_qm_result->SetLineColor(kRed+1);
0113 g_sphenix_qm_result->SetLineWidth(0);
0114 g_sphenix_qm_result->SetFillColorAlpha(kRed-4, 0.5);
0115 TGraphErrors* g_sphenix_qm_result_point = new TGraphErrors(9, g_sphenix_qm_x, g_sphenix_qm_y, g_sphenix_qm_xerrdown, g_sphenix_qm_stat);
0116 g_sphenix_qm_result_point->SetMarkerStyle(20);
0117 g_sphenix_qm_result_point->SetMarkerSize(1.5);
0118 g_sphenix_qm_result_point->SetMarkerColor(kRed+1);
0119 g_sphenix_qm_result_point->SetLineColor(kRed+1);
0120 g_sphenix_qm_result_point->SetLineWidth(2);
0121
0122
0123 TFile *f_sim_specturm = new TFile("output_sim_r04.root", "READ");
0124 TH1D* h_pythia_spectrum = (TH1D*)f_sim_specturm->Get("h_truth_nominal");
0125 h_pythia_spectrum->Multiply(h_etacorrection);
0126 for (int ib = 1; ib <= h_pythia_spectrum->GetNbinsX(); ++ib) {
0127 h_pythia_spectrum->SetBinContent(ib, h_pythia_spectrum->GetBinContent(ib)/ (float)(1.4 * h_pythia_spectrum->GetBinWidth(ib)));
0128 h_pythia_spectrum->SetBinError(ib, h_pythia_spectrum->GetBinError(ib)/ (float)(1.4 * h_pythia_spectrum->GetBinWidth(ib)));
0129 }
0130 double g_pythia_x[truthnpt], g_pythia_xerrdown[truthnpt], g_pythia_xerrup[truthnpt], g_pythia_y[truthnpt], g_pythia_yerrdown[truthnpt], g_pythia_yerrup[truthnpt];
0131 for (int ib = 0; ib < truthnpt; ++ib) {
0132 g_pythia_x[ib] = (truthptbins[ib] + truthptbins[ib+1]) / 2;
0133 g_pythia_xerrdown[ib] = (truthptbins[ib+1] - truthptbins[ib]) / 2;
0134 g_pythia_xerrup[ib] = (truthptbins[ib+1] - truthptbins[ib]) / 2;
0135 g_pythia_y[ib] = h_pythia_spectrum->GetBinContent(ib+1);
0136 g_pythia_yerrdown[ib] = h_pythia_spectrum->GetBinError(ib+1);
0137 g_pythia_yerrup[ib] = h_pythia_spectrum->GetBinError(ib+1);
0138 }
0139 TGraphAsymmErrors* g_pythia_result = new TGraphAsymmErrors(truthnpt, g_pythia_x, g_pythia_y, g_pythia_xerrdown, g_pythia_xerrup, g_pythia_yerrdown, g_pythia_yerrup);
0140 g_pythia_result->SetMarkerStyle(21);
0141 g_pythia_result->SetMarkerSize(1.5);
0142 g_pythia_result->SetMarkerColor(kMagenta+2);
0143 g_pythia_result->SetLineColor(kMagenta+2);
0144
0145
0146 TFile *f_sim_specturm_herwig = new TFile("herwigXsections.root", "READ");
0147 TH1D* h_herwig_spectrum_04 = (TH1D*)f_sim_specturm_herwig->Get("Herwigxsections");
0148 double g_herwig_x_04[11], g_herwig_y_04[11], g_herwig_xerrdown_04[11], g_herwig_xerrup_04[11], g_herwig_yerrdown_04[11], g_herwig_yerrup_04[11];
0149 for (int ib = 0; ib < 11; ++ib) {
0150 g_herwig_x_04[ib] = h_herwig_spectrum_04->GetBinCenter(ib+1);
0151 g_herwig_xerrdown_04[ib] = h_herwig_spectrum_04->GetBinWidth(ib+1) / 2;
0152 g_herwig_xerrup_04[ib] = h_herwig_spectrum_04->GetBinWidth(ib+1) / 2;
0153 g_herwig_y_04[ib] = h_herwig_spectrum_04->GetBinContent(ib+1);
0154 g_herwig_yerrdown_04[ib] = h_herwig_spectrum_04->GetBinError(ib+1);
0155 g_herwig_yerrup_04[ib] = h_herwig_spectrum_04->GetBinError(ib+1);
0156 }
0157 TGraphAsymmErrors* g_herwig_result_04 = new TGraphAsymmErrors(11, g_herwig_x_04, g_herwig_y_04, g_herwig_xerrdown_04, g_herwig_xerrup_04, g_herwig_yerrdown_04, g_herwig_yerrup_04);
0158 g_herwig_result_04->SetMarkerStyle(33);
0159 g_herwig_result_04->SetMarkerSize(1.5);
0160 g_herwig_result_04->SetMarkerColor(kGreen+2);
0161 g_herwig_result_04->SetLineColor(kGreen+2);
0162
0163
0164 ifstream myfile1, myfile05, myfile2;
0165 TGraph *tjet = new TGraph();
0166 TGraph *tjet05 = new TGraph();
0167 TGraph *tjet2 = new TGraph();
0168 double scalefactor = 2.0 * 3.14159;
0169 myfile1.open("sphenix_nlo/sPhenix-sc1.dat");
0170 myfile05.open("sphenix_nlo/sPhenix-sc05.dat");
0171 myfile2.open("sphenix_nlo/sPhenix-sc2.dat");
0172 double drop1, drop2, drop3, pt1, yield1, pt05, yield05, pt2, yield2;
0173 std::vector<double> xpoint, ypoint05, ypoint2;
0174 for (int index=0; index<66; index++) {
0175 myfile1 >> drop1 >> drop2 >> drop3 >> pt1 >> yield1;
0176 yield1 = yield1 * scalefactor * pt1;
0177 tjet->SetPoint(index, pt1, yield1);
0178 xpoint.push_back(pt1);
0179
0180 myfile05 >> drop1 >> drop2 >> drop3 >> pt05 >> yield05;
0181 yield05 = yield05 * scalefactor * pt05;
0182 tjet05->SetPoint(index, pt05, yield05);
0183 ypoint05.push_back(yield05);
0184
0185 myfile2 >> drop1 >> drop2 >> drop3 >> pt2 >> yield2;
0186 yield2 = yield2 * scalefactor * pt2;
0187 tjet2->SetPoint(index, pt2, yield2);
0188 ypoint2.push_back(yield2);
0189 }
0190 std::vector <double> xpoint_shadow, ypoint_shadow;
0191 for (int i = 0; i < xpoint.size(); i++) {
0192 xpoint_shadow.push_back(xpoint[i]);
0193 ypoint_shadow.push_back(ypoint05[i]);
0194 }
0195 for (int i = xpoint.size(); i > 0; i--) {
0196 xpoint_shadow.push_back(xpoint[i-1]);
0197 ypoint_shadow.push_back(ypoint2[i-1]);
0198 }
0199 TGraph *tjet_shadow = new TGraph(2*xpoint.size(), &xpoint_shadow[0], &ypoint_shadow[0]);
0200 tjet_shadow->SetMarkerStyle(20);
0201 tjet_shadow->SetMarkerSize(1.5);
0202 tjet_shadow->SetMarkerColor(kRed-4);
0203 tjet_shadow->SetLineColor(kRed);
0204 tjet_shadow->SetFillColorAlpha(kRed-4, 0.5);
0205
0206 myfile1.close();
0207 myfile05.close();
0208 myfile2.close();
0209 tjet->SetMarkerStyle(21);
0210 tjet->SetMarkerColor(2);
0211 tjet->SetLineColor(2);
0212 tjet05->SetMarkerStyle(21);
0213 tjet05->SetMarkerColor(2);
0214 tjet05->SetLineColor(2);
0215 tjet2->SetMarkerStyle(21);
0216 tjet2->SetMarkerColor(2);
0217 tjet2->SetLineColor(2);
0218
0219 TCanvas *can_cross_section = new TCanvas("can_cross_section", "", 2550, 2400);
0220 gStyle->SetPalette(57);
0221 can_cross_section->SetTopMargin(0.03);
0222 can_cross_section->SetLeftMargin(0.15);
0223 can_cross_section->SetBottomMargin(0.12);
0224 can_cross_section->SetRightMargin(0.05);
0225 can_cross_section->SetLogy();
0226 TH2F *frame_cross_section = new TH2F("frame_cross_section", "", 10, 15, 75, 10, 1e-05, 1e+05);
0227 frame_cross_section->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0228 frame_cross_section->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0229 frame_cross_section->GetXaxis()->SetTitleOffset(0.98);
0230 frame_cross_section->GetYaxis()->SetTitleOffset(1.15);
0231 frame_cross_section->GetXaxis()->SetLabelSize(0.045);
0232 frame_cross_section->GetYaxis()->SetLabelSize(0.045);
0233 frame_cross_section->GetXaxis()->CenterTitle();
0234 frame_cross_section->GetYaxis()->CenterTitle();
0235 frame_cross_section->Draw();
0236 g_sphenix_result->Draw("2 same");
0237 myText(0.6, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0238 myText(0.6, 0.8, 1, "anti-k_{t} #kern[-0.3]{#it{R}} = 0.4", 0.04);
0239 myText(0.6, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0240 TLegend *leg_cross_section = new TLegend(0.15, 0.18, 0.9, 0.33);
0241 leg_cross_section->SetBorderSize(0);
0242 leg_cross_section->SetFillStyle(0);
0243 leg_cross_section->SetTextSize(0.03);
0244 leg_cross_section->SetTextFont(42);
0245 g_pythia_result->SetMarkerStyle(54);
0246 g_herwig_result_04->SetMarkerStyle(56);
0247 g_pythia_result->SetMarkerSize(4.5);
0248 g_herwig_result_04->SetMarkerSize(6.5);
0249 g_sphenix_result->SetMarkerSize(4.5);
0250 g_sphenix_result_point->SetMarkerSize(4.5);
0251 leg_cross_section->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0252 leg_cross_section->AddEntry(g_pythia_result, "PYTHIA8 truth jet spectrum", "pl");
0253 leg_cross_section->AddEntry(g_herwig_result_04, "HERWIG truth jet spectrum", "pl");
0254 leg_cross_section->AddEntry(tjet_shadow, "NLO pQCD (No Hadronization)", "f");
0255 leg_cross_section->AddEntry(tjet_shadow, "Werner Vogelsang", "");
0256 leg_cross_section->Draw();
0257 g_sphenix_result->Draw("P same");
0258 g_sphenix_result_point->Draw("P same");
0259 g_pythia_result->Draw("p same");
0260 g_herwig_result_04->Draw("p same");
0261 tjet_shadow->Draw("F same");
0262 can_cross_section->SaveAs(Form("figure/cross_section_jet0%d.png", (int)(10*jet_radius)));
0263 }
0264
0265 void do_normalization(TH1D* h_unfolded, double luminosity, TH1D* h_etacorrection, float jet_radius) {
0266 h_unfolded->Multiply(h_etacorrection);
0267 double delta_eta = 2 * (1.1 - jet_radius);
0268 for (int i = 1; i <= h_unfolded->GetNbinsX(); ++i) {
0269 double binwidth = h_unfolded->GetBinWidth(i);
0270 h_unfolded->SetBinContent(i, h_unfolded->GetBinContent(i) / (double)(luminosity * binwidth * delta_eta));
0271 h_unfolded->SetBinError(i, h_unfolded->GetBinError(i) / (double)(luminosity * binwidth * delta_eta));
0272 }
0273 }
0274
0275 TH1D* get_statistical_uncertainty(TH1D* h_unfold, string name) {
0276 TH1D* h_uncertainty = (TH1D*)h_unfold->Clone(name.c_str());
0277 h_uncertainty->Reset();
0278 for (int i = 1; i <= h_unfold->GetNbinsX(); ++i) {
0279 double content = h_unfold->GetBinContent(i);
0280 if (content == 0) continue;
0281 double error = h_unfold->GetBinError(i);
0282 double uncertainty = error / content;
0283 h_uncertainty->SetBinContent(i, uncertainty);
0284 h_uncertainty->SetBinError(i, 0);
0285 }
0286 return h_uncertainty;
0287 }
0288
0289 void get_uncertainty(TH1D* &h_uncertainty_up, TH1D* &h_uncertainty_down, TH1D* h_nominal, TH1D* h_shiftup, TH1D* h_shiftdown, string name) {
0290 h_uncertainty_up = (TH1D*)h_nominal->Clone((name + "_up").c_str());
0291 h_uncertainty_down = (TH1D*)h_nominal->Clone((name + "_down").c_str());
0292 h_uncertainty_up->Reset();
0293 h_uncertainty_down->Reset();
0294 for (int i = 1; i <= h_nominal->GetNbinsX(); ++i) {
0295 double up = 0, down = 0;
0296 if (h_nominal->GetBinContent(i) != 0) {
0297 double diff_up = (h_shiftup->GetBinContent(i) - h_nominal->GetBinContent(i)) / (double)h_nominal->GetBinContent(i);
0298 if (diff_up >= 0) up = TMath::Sqrt(TMath::Power(diff_up, 2) + TMath::Power(up, 2));
0299 else down = TMath::Sqrt(TMath::Power(diff_up, 2) + TMath::Power(down, 2));
0300 double diff_down = (h_shiftdown->GetBinContent(i) - h_nominal->GetBinContent(i)) / (double)h_nominal->GetBinContent(i);
0301 if (diff_down >= 0) up = TMath::Sqrt(TMath::Power(diff_down, 2) + TMath::Power(up, 2));
0302 else down = TMath::Sqrt(TMath::Power(diff_down, 2) + TMath::Power(down, 2));
0303 }
0304 h_uncertainty_up->SetBinContent(i, up);
0305 h_uncertainty_up->SetBinError(i, 0);
0306 h_uncertainty_down->SetBinContent(i, down);
0307 h_uncertainty_down->SetBinError(i, 0);
0308 }
0309 }
0310
0311 void get_unfold_uncertainty(TH1D* &h_uncertainty_up, TH1D* &h_uncertainty_down, TH1D* h_nominal, TH1D* h_shift, string name) {
0312 h_uncertainty_up = (TH1D*)h_nominal->Clone((name + "_up").c_str());
0313 h_uncertainty_down = (TH1D*)h_nominal->Clone((name + "_down").c_str());
0314 h_uncertainty_up->Reset();
0315 h_uncertainty_down->Reset();
0316 for (int i = 1; i <= h_nominal->GetNbinsX(); ++i) {
0317 double up = 0, down = 0;
0318 if (h_nominal->GetBinContent(i) == 0) continue;
0319
0320 double diff_up = (h_shift->GetBinContent(i) - h_nominal->GetBinContent(i)) / (double)h_nominal->GetBinContent(i);
0321 if (diff_up >= 0) up = TMath::Sqrt(TMath::Power(diff_up, 2) + TMath::Power(up, 2));
0322 else down = TMath::Sqrt(TMath::Power(diff_up, 2) + TMath::Power(down, 2));
0323 double diff_down = -1*(h_shift->GetBinContent(i) - h_nominal->GetBinContent(i)) / (double)h_nominal->GetBinContent(i);
0324 if (diff_down >= 0) up = TMath::Sqrt(TMath::Power(diff_down, 2) + TMath::Power(up, 2));
0325 else down = TMath::Sqrt(TMath::Power(diff_down, 2) + TMath::Power(down, 2));
0326
0327 h_uncertainty_up->SetBinContent(i, up);
0328 h_uncertainty_up->SetBinError(i, 0);
0329 h_uncertainty_down->SetBinContent(i, down);
0330 h_uncertainty_down->SetBinError(i, 0);
0331 }
0332 }
0333
0334 void get_total_uncertainty(TH1D* &h_total_up, TH1D* &h_total_down,TH1D* h_statistical, std::vector<TH1D*> h_systematic_up, std::vector<TH1D*> h_systematic_down) {
0335 h_total_up = (TH1D*)h_statistical->Clone("h_total_up");
0336 h_total_down = (TH1D*)h_statistical->Clone("h_total_down");
0337 h_total_up->Reset();
0338 h_total_down->Reset();
0339 for (int i = 1; i <= h_statistical->GetNbinsX(); ++i) {
0340 double stat_error = h_statistical->GetBinContent(i);
0341
0342 double sys_error_up = 0;
0343 for (auto& h_up : h_systematic_up) {
0344 if (h_up->GetBinContent(i) > 0) {
0345 sys_error_up += h_up->GetBinContent(i) * h_up->GetBinContent(i);
0346 }
0347 }
0348 double total_error_up = TMath::Sqrt(stat_error * stat_error + sys_error_up);
0349 h_total_up->SetBinContent(i, total_error_up);
0350 h_total_up->SetBinError(i, 0);
0351
0352 double sys_error_down = 0;
0353 for (auto& h_down : h_systematic_down) {
0354 if (h_down->GetBinContent(i) > 0) {
0355 sys_error_down += h_down->GetBinContent(i) * h_down->GetBinContent(i);
0356 }
0357 }
0358 double total_error_down = TMath::Sqrt(stat_error * stat_error + sys_error_down);
0359 h_total_down->SetBinContent(i, total_error_down);
0360 h_total_down->SetBinError(i, 0);
0361 }
0362 }
0363
0364 void draw_total_uncertainty(TH1D* h_total_up, TH1D* h_total_down, TH1D* h_statistical, std::vector<TH1D*> h_systematic_up, std::vector<TH1D*> h_systematic_down, std::vector<int> color, std::vector<string> legend, float jet_radius) {
0365 TH1D* h_statistical_up = (TH1D*)h_statistical->Clone("h_statistical_up");
0366 TH1D* h_statistical_down = (TH1D*)h_statistical->Clone("h_statistical_down");
0367 h_statistical_down->Scale(-1.0);
0368 h_statistical_up->SetLineColor(kViolet-1);
0369 h_statistical_down->SetLineColor(kViolet-1);
0370 for (int i = 0; i < h_systematic_down.size(); ++i) {
0371 h_systematic_up[i]->SetLineColor(color[i]);
0372 h_systematic_down[i]->SetLineColor(color[i]);
0373 h_systematic_down[i]->Scale(-1.0);
0374 }
0375 h_total_up->SetMarkerStyle(20);
0376 h_total_down->SetMarkerStyle(20);
0377 h_total_up->SetLineColor(kBlack);
0378 h_total_down->SetLineColor(kBlack);
0379 h_total_up->SetMarkerColor(kBlack);
0380 h_total_down->SetMarkerColor(kBlack);
0381 h_total_down->Scale(-1.0);
0382
0383 TCanvas *can_unc = new TCanvas("can_unc", "", 850, 800);
0384 gStyle->SetPalette(57);
0385 can_unc->SetTopMargin(0.03);
0386 can_unc->SetLeftMargin(0.15);
0387 can_unc->SetBottomMargin(0.12);
0388 can_unc->SetRightMargin(0.05);
0389 TH2F *frame_unc = new TH2F("frame_unc", "", 10, 15, 75, 2, -1.8, 2.5);
0390 frame_unc->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0391 frame_unc->GetYaxis()->SetTitle("Relative Uncertainty");
0392 frame_unc->GetXaxis()->SetTitleOffset(0.98);
0393 frame_unc->GetYaxis()->SetTitleOffset(1.17);
0394 frame_unc->GetXaxis()->SetLabelSize(0.045);
0395 frame_unc->GetYaxis()->SetLabelSize(0.045);
0396 frame_unc->GetXaxis()->CenterTitle();
0397 frame_unc->GetYaxis()->CenterTitle();
0398 frame_unc->Draw();
0399 h_statistical_up->Draw("hist same");
0400 h_statistical_down->Draw("hist same");
0401 for (int i = 0; i < h_systematic_up.size(); ++i) {
0402 h_systematic_up[i]->Draw("hist same");
0403 h_systematic_down[i]->Draw("hist same");
0404 }
0405 h_total_up->Draw("p same");
0406 h_total_down->Draw("p same");
0407
0408 myText(0.17, 0.92, 1, "#bf{#it{sPHENIX}} Internal", 0.05);
0409 myText(0.17, 0.86, 1, "p+p#sqrt{s} = 200 GeV", 0.05);
0410 myText(0.17, 0.8, 1, "anti-k_{t} #kern[-0.3]{#it{R}} = 0.4", 0.05);
0411 myText(0.17, 0.74, 1, "|#eta^{jet}| < 0.7", 0.05);
0412 TLegend *leg_unc1 = new TLegend(0.15, 0.13, 0.75, 0.25);
0413 leg_unc1->SetNColumns(2);
0414 leg_unc1->SetBorderSize(0);
0415 leg_unc1->SetFillStyle(0);
0416 leg_unc1->SetTextFont(42);
0417 leg_unc1->AddEntry(h_statistical_up, "Statistical Uncertainty", "l");
0418 for (int i = 0; i < h_systematic_up.size(); ++i) {
0419 leg_unc1->AddEntry(h_systematic_up[i], legend[i].c_str(), "l");
0420 }
0421 leg_unc1->AddEntry(h_total_up, "Total Uncertainty", "p");
0422 leg_unc1->Draw();
0423 TLine *line = new TLine(15, 0, 75, 0);
0424 line->SetLineStyle(10);
0425 line->Draw("same");
0426 can_unc->SaveAs(Form("figure/uncertainty_jet0%d.png", (int)(10*jet_radius)));
0427
0428 delete can_unc;
0429 delete frame_unc;
0430 }
