Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:10 

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 ifstream f_luminosity("luminosity.txt");
0014 double luminosity;
0015 f_luminosity >> luminosity;
0016 double delta_eta = 1.4;
0017 int n_underflow_bin = 2;
0018 int n_overflow_bin = 1;
0019 std::cout << "Luminosity: " << luminosity << " mb^-1" << std::endl;
0020 std::cout << "Delta eta: " << delta_eta << std::endl;
0021 
0022 void draw_1D_multiple_plot_uncertainty(std::vector<TH1F*> h_input, float x_min, float x_max, float ratio_min, float ratio_max, std::string xtitle, std::string ytitle, std::vector<std::string> legend, float xstart_legend, float ystart_legend, std::vector<int> color, std::vector<int> markerstyle, std::string output_name);
0023 
0024 void get_finalspectrum() {
0025 
0026   std::vector<TH1F*> h_input;
0027   std::vector<TGraphAsymmErrors*> g_input;
0028   std::vector<int> color;
0029   std::vector<int> markerstyle;
0030   std::vector<std::string> text;
0031   std::vector<std::string> legend;
0032 
0033   TFile *f_data = new TFile("output_data.root", "READ");
0034   TH1D* h_recojet_pt_record = (TH1D*)f_data->Get("h_recojet_pt_record");
0035   h_recojet_pt_record->Rebin(10);
0036   TH1D* h_recojet_pt_record_dijet = (TH1D*)f_data->Get("h_recojet_pt_record_dijet");
0037   h_recojet_pt_record_dijet->Rebin(10);
0038   TH1D* h_recojet_pt_record_frac = (TH1D*)f_data->Get("h_recojet_pt_record_frac");
0039   h_recojet_pt_record_frac->Rebin(10);
0040   
0041   h_input.push_back((TH1F*)h_recojet_pt_record);
0042   h_input.push_back((TH1F*)h_recojet_pt_record_dijet);
0043   h_input.push_back((TH1F*)h_recojet_pt_record_frac);
0044   color.push_back(kBlack);
0045   color.push_back(kAzure+1);
0046   color.push_back(kRed+1);
0047   markerstyle.push_back(20);
0048   markerstyle.push_back(21);
0049   markerstyle.push_back(20);
0050   text.push_back("#bf{#it{sPHENIX}} Internal");
0051   text.push_back("p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}");
0052   text.push_back("anti-k_{T}#it{R} = 0.4");
0053   text.push_back("|#eta^{jet}| < 0.7");
0054   legend.push_back("Reco jet spectrum");
0055   legend.push_back("After dijet cut");
0056   legend.push_back("After fraction cut");
0057   draw_1D_multiple_plot(h_input, color, markerstyle,
0058                         false, 10, true,
0059                         true, 10, 100, false,
0060                         true, 0.3, 4e+06, true,
0061                         true, "p_{T}^{jet} [GeV]", "dN_{jet} / dp_{T} [1/GeV]",
0062                         true, text, 0.45, 0.9, 0.04,
0063                         true, legend, 0.48, 0.65, 0.04,
0064                         "figure/raw_reco_jet_spectrum.pdf");
0065   h_input.clear();
0066   color.clear();
0067   markerstyle.clear();
0068   text.clear();
0069   legend.clear();
0070 
0071   f_data->Close();
0072   std::cout << "check0" << std::endl;
0073 
0074   TFile* f_spectrum = new TFile("output_unfolded.root", "READ");
0075 
0076   TH1D* h_unfold_calib_dijet_reweighted_1 = new TH1D("h_spectra_calib_dijet_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_reweighted_1");
0077   TH1D* h_unfold_calib_dijet_effdown_reweighted_1 = new TH1D("h_spectra_calib_dijet_effdown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_effdown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_effdown_reweighted_1");
0078   TH1D* h_unfold_calib_dijet_effup_reweighted_1 = new TH1D("h_spectra_calib_dijet_effup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_effup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_effup_reweighted_1");
0079   TH1D* h_unfold_calib_dijet_jesdown_reweighted_1 = new TH1D("h_spectra_calib_dijet_jesdown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_jesdown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_jesdown_reweighted_1");
0080   TH1D* h_unfold_calib_dijet_jesup_reweighted_1 = new TH1D("h_spectra_calib_dijet_jesup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_jesup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_jesup_reweighted_1");
0081   TH1D* h_unfold_calib_dijet_jerdown_reweighted_1 = new TH1D("h_spectra_calib_dijet_jerdown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_jerdown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_jerdown_reweighted_1");
0082   TH1D* h_unfold_calib_dijet_jerup_reweighted_1 = new TH1D("h_spectra_calib_dijet_jerup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_jerup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_jerup_reweighted_1");
0083 
0084   TH1D* h_unfold_calib_frac_reweighted_1 = new TH1D("h_spectra_calib_frac_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_reweighted_1");
0085   TH1D* h_unfold_calib_frac_effdown_reweighted_1 = new TH1D("h_spectra_calib_frac_effdown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_effdown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_effdown_reweighted_1");
0086   TH1D* h_unfold_calib_frac_effup_reweighted_1 = new TH1D("h_spectra_calib_frac_effup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_effup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_effup_reweighted_1");
0087   TH1D* h_unfold_calib_frac_jesdown_reweighted_1 = new TH1D("h_spectra_calib_frac_jesdown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_jesdown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_jesdown_reweighted_1");
0088   TH1D* h_unfold_calib_frac_jesup_reweighted_1 = new TH1D("h_spectra_calib_frac_jesup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_jesup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_jesup_reweighted_1");
0089   TH1D* h_unfold_calib_frac_jerdown_reweighted_1 = new TH1D("h_spectra_calib_frac_jerdown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_jerdown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_jerdown_reweighted_1");
0090   TH1D* h_unfold_calib_frac_jerup_reweighted_1 = new TH1D("h_spectra_calib_frac_jerup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_jerup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_jerup_reweighted_1");
0091 
0092   TH1D* h_unfold_calib_dijet_reweighted_2 = new TH1D("h_spectra_calib_dijet_reweighted_2", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_reweighted_2_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_reweighted_2");
0093   TH1D* h_unfold_calib_dijet_reweighted_3 = new TH1D("h_spectra_calib_dijet_reweighted_3", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_reweighted_3_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_reweighted_3");
0094   TH1D* h_unfold_calib_dijet_mbddown_reweighted_1 = new TH1D("h_spectra_calib_dijet_mbddown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_mbddown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_mbddown_reweighted_1");
0095   TH1D* h_unfold_calib_dijet_mbdup_reweighted_1 = new TH1D("h_spectra_calib_dijet_mbdup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_mbdup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_mbdup_reweighted_1");
0096 
0097   TH1D* h_unfold_calib_frac_reweighted_2 = new TH1D("h_spectra_calib_frac_reweighted_2", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_reweighted_2_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_reweighted_2");
0098   TH1D* h_unfold_calib_frac_reweighted_3 = new TH1D("h_spectra_calib_frac_reweighted_3", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_reweighted_3_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_reweighted_3");
0099   TH1D* h_unfold_calib_frac_mbddown_reweighted_1 = new TH1D("h_spectra_calib_frac_mbddown_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_mbddown_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_mbddown_reweighted_1");
0100   TH1D* h_unfold_calib_frac_mbdup_reweighted_1 = new TH1D("h_spectra_calib_frac_mbdup_reweighted_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_mbdup_reweighted_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_mbdup_reweighted_1");
0101 
0102   TH1D* h_unfold_calib_dijet_1 = new TH1D("h_spectra_calib_dijet_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_dijet_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_dijet_1");
0103   TH1D* h_unfold_calib_frac_1 = new TH1D("h_spectra_calib_frac_1", "", calibnpt, calibptbins); TH1D* h_unfold_calib_frac_1_temp = (TH1D*)f_spectrum->Get("h_unfold_calib_frac_1");
0104 
0105   for (int ib = 1; ib <= h_unfold_calib_dijet_reweighted_1->GetNbinsX(); ++ib) {
0106     h_unfold_calib_dijet_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0107     h_unfold_calib_dijet_effdown_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_effdown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_effdown_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_effdown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0108     h_unfold_calib_dijet_effup_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_effup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_effup_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_effup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0109     h_unfold_calib_dijet_jesdown_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_jesdown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_jesdown_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_jesdown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0110     h_unfold_calib_dijet_jesup_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_jesup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_jesup_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_jesup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0111     h_unfold_calib_dijet_jerdown_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_jerdown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_jerdown_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_jerdown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0112     h_unfold_calib_dijet_jerup_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_jerup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_jerup_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_jerup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0113 
0114     h_unfold_calib_frac_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_reweighted_1->SetBinError(ib, h_unfold_calib_frac_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0115     h_unfold_calib_frac_effdown_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_effdown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_effdown_reweighted_1->SetBinError(ib, h_unfold_calib_frac_effdown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0116     h_unfold_calib_frac_effup_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_effup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_effup_reweighted_1->SetBinError(ib, h_unfold_calib_frac_effup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0117     h_unfold_calib_frac_jesdown_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_jesdown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_jesdown_reweighted_1->SetBinError(ib, h_unfold_calib_frac_jesdown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0118     h_unfold_calib_frac_jesup_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_jesup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_jesup_reweighted_1->SetBinError(ib, h_unfold_calib_frac_jesup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0119     h_unfold_calib_frac_jerdown_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_jerdown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_jerdown_reweighted_1->SetBinError(ib, h_unfold_calib_frac_jerdown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0120     h_unfold_calib_frac_jerup_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_jerup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_jerup_reweighted_1->SetBinError(ib, h_unfold_calib_frac_jerup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0121 
0122     h_unfold_calib_dijet_reweighted_2->SetBinContent(ib, h_unfold_calib_dijet_reweighted_2_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_reweighted_2->SetBinError(ib, h_unfold_calib_dijet_reweighted_2_temp->GetBinError(ib + n_underflow_bin));
0123     h_unfold_calib_dijet_reweighted_3->SetBinContent(ib, h_unfold_calib_dijet_reweighted_3_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_reweighted_3->SetBinError(ib, h_unfold_calib_dijet_reweighted_3_temp->GetBinError(ib + n_underflow_bin));
0124     h_unfold_calib_dijet_mbddown_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_mbddown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_mbddown_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_mbddown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0125     h_unfold_calib_dijet_mbdup_reweighted_1->SetBinContent(ib, h_unfold_calib_dijet_mbdup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_mbdup_reweighted_1->SetBinError(ib, h_unfold_calib_dijet_mbdup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0126 
0127     h_unfold_calib_frac_reweighted_2->SetBinContent(ib, h_unfold_calib_frac_reweighted_2_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_reweighted_2->SetBinError(ib, h_unfold_calib_frac_reweighted_2_temp->GetBinError(ib + n_underflow_bin));
0128     h_unfold_calib_frac_reweighted_3->SetBinContent(ib, h_unfold_calib_frac_reweighted_3_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_reweighted_3->SetBinError(ib, h_unfold_calib_frac_reweighted_3_temp->GetBinError(ib + n_underflow_bin));
0129     h_unfold_calib_frac_mbddown_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_mbddown_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_mbddown_reweighted_1->SetBinError(ib, h_unfold_calib_frac_mbddown_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0130     h_unfold_calib_frac_mbdup_reweighted_1->SetBinContent(ib, h_unfold_calib_frac_mbdup_reweighted_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_mbdup_reweighted_1->SetBinError(ib, h_unfold_calib_frac_mbdup_reweighted_1_temp->GetBinError(ib + n_underflow_bin));
0131 
0132     h_unfold_calib_dijet_1->SetBinContent(ib, h_unfold_calib_dijet_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_dijet_1->SetBinError(ib, h_unfold_calib_dijet_1_temp->GetBinError(ib + n_underflow_bin));
0133     h_unfold_calib_frac_1->SetBinContent(ib, h_unfold_calib_frac_1_temp->GetBinContent(ib + n_underflow_bin)); h_unfold_calib_frac_1->SetBinError(ib, h_unfold_calib_frac_1_temp->GetBinError(ib + n_underflow_bin));
0134   }
0135 
0136   // Get Uncertainty
0137   double uncertainty_unfold_upper[9] = {0};
0138   double uncertainty_unfold_lower[9] = {0};
0139   double uncertainty_cut_upper[9] = {0};
0140   double uncertainty_cut_lower[9] = {0};
0141   double uncertainty_mbd_upper[9] = {0}; // in ratio
0142   double uncertainty_mbd_lower[9] = {0}; // in ratio
0143   double uncertainty_other_upper[9] = {0};
0144   double uncertainty_other_lower[9] = {0};
0145   double uncertainty_final_upper[9] = {0};
0146   double uncertainty_final_lower[9] = {0};
0147 
0148   double uncertainty_unfold_upper_record[9] = {0}, uncertainty_unfold_lower_record[9] = {0};
0149   double uncertainty_cut_upper_record[9] = {0}, uncertainty_cut_lower_record[9] = {0};
0150   double uncertainty_jes_upper_record[9] = {0}, uncertainty_jes_lower_record[9] = {0};
0151   double uncertainty_jer_upper_record[9] = {0}, uncertainty_jer_lower_record[9] = {0};
0152   double uncertainty_jeteff_upper_record[9] = {0}, uncertainty_jeteff_lower_record[9] = {0};
0153   double uncertainty_mbdeff_upper_record[9] = {0}, uncertainty_mbdeff_lower_record[9] = {0};
0154   double uncertainty_mbdcs_upper_record[9] = {0}, uncertainty_mbdcs_lower_record[9] = {0};
0155   double uncertainty_stat_upper_record[9] = {0}, uncertainty_stat_lower_record[9] = {0};
0156   double uncertainty_total_upper_record[9] = {0}, uncertainty_total_lower_record[9] = {0};
0157 
0158   TH1D* h_unfold_calib_nominal = new TH1D("h_unfold_calib_nominal", "", calibnpt, calibptbins);
0159   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0160     h_unfold_calib_nominal->SetBinContent(ib, (h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) + h_unfold_calib_frac_reweighted_1->GetBinContent(ib)) / 2.);
0161     h_unfold_calib_nominal->SetBinError(ib, TMath::Sqrt(TMath::Power(h_unfold_calib_dijet_reweighted_1->GetBinError(ib), 2) + TMath::Power(h_unfold_calib_frac_reweighted_1->GetBinError(ib), 2)) / 2.);
0162   }
0163 
0164   TH1D* h_unfold_calib_nominal_unfold = new TH1D("h_unfold_calib_nominal_unfold", "", calibnpt, calibptbins);
0165   for (int ib = 1; ib <= h_unfold_calib_nominal_unfold->GetNbinsX(); ++ib) {
0166     h_unfold_calib_nominal_unfold->SetBinContent(ib, (h_unfold_calib_dijet_1->GetBinContent(ib) + h_unfold_calib_frac_1->GetBinContent(ib)) / 2.);
0167     h_unfold_calib_nominal_unfold->SetBinError(ib, TMath::Sqrt(TMath::Power(h_unfold_calib_dijet_1->GetBinError(ib), 2) + TMath::Power(h_unfold_calib_frac_1->GetBinError(ib), 2)) / 2.);
0168   }
0169 
0170   // MBD Cross Section Uncertainty
0171   for (int ib = 0; ib < 9; ++ib) {
0172     uncertainty_mbd_upper[ib] = 4.3 / 26.1;
0173     uncertainty_mbd_lower[ib] = 1.1 / 26.1;
0174     uncertainty_mbdcs_upper_record[ib] = uncertainty_mbd_upper[ib];
0175     uncertainty_mbdcs_lower_record[ib] = uncertainty_mbd_lower[ib];
0176   }
0177 
0178   // Unfolding Uncertainty
0179   h_input.push_back((TH1F*)h_unfold_calib_nominal);
0180   h_input.push_back((TH1F*)h_unfold_calib_nominal_unfold);
0181   color.push_back(kBlack);
0182   color.push_back(kRed);
0183   markerstyle.push_back(20);
0184   markerstyle.push_back(20);
0185   legend.push_back("Nominal distribution");
0186   legend.push_back("Nominal distribution without reweighting");
0187   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0.7, 1.3, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_unfold.pdf");
0188   h_input.clear();
0189   color.clear();
0190   markerstyle.clear();
0191   legend.clear();
0192   std::cout << std::endl << "Uncertainty from unfolding: " << std::endl << "up\t\tdown" << std::endl;
0193   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0194     double unc = h_unfold_calib_nominal_unfold->GetBinContent(ib) - h_unfold_calib_nominal->GetBinContent(ib);
0195     double unc_up = 0;
0196     double unc_down = 0;
0197     unc_up = TMath::Sqrt(unc_up * unc_up + unc * unc);
0198     unc_down = TMath::Sqrt(unc_down * unc_down + unc * unc);
0199 
0200     uncertainty_unfold_upper[ib-1] = TMath::Sqrt(uncertainty_unfold_upper[ib-1] * uncertainty_unfold_upper[ib-1] + unc * unc);
0201     uncertainty_unfold_lower[ib-1] = TMath::Sqrt(uncertainty_unfold_lower[ib-1] * uncertainty_unfold_lower[ib-1] + unc * unc);
0202 
0203     uncertainty_unfold_upper_record[ib-1] = unc_up / h_unfold_calib_nominal->GetBinContent(ib);
0204     uncertainty_unfold_lower_record[ib-1] = unc_down / h_unfold_calib_nominal->GetBinContent(ib);
0205     std::cout << unc_up / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%      \t" << unc_down / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%" << std::endl;
0206   }
0207 
0208   // Beam Background Cut Uncertainty
0209   h_input.push_back((TH1F*)h_unfold_calib_nominal);
0210   h_input.push_back((TH1F*)h_unfold_calib_dijet_reweighted_1);
0211   h_input.push_back((TH1F*)h_unfold_calib_frac_reweighted_1);
0212   color.push_back(kBlack);
0213   color.push_back(kAzure+1);
0214   color.push_back(kRed+1);
0215   markerstyle.push_back(20);
0216   markerstyle.push_back(20);
0217   markerstyle.push_back(20);
0218   legend.push_back("Nominal distribution");
0219   legend.push_back("With dijet cut");
0220   legend.push_back("With fraction cut");
0221   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0.8, 1.2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_cut.pdf");
0222   h_input.clear();
0223   color.clear();
0224   markerstyle.clear();
0225   legend.clear();
0226   std::cout << std::endl << "Uncertainty from cut: " << std::endl << "up\t\tdown" << std::endl;
0227   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0228     double unc1 = h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) - h_unfold_calib_nominal->GetBinContent(ib);
0229     double unc2 = h_unfold_calib_frac_reweighted_1->GetBinContent(ib) - h_unfold_calib_nominal->GetBinContent(ib);
0230     double unc_up = 0;
0231     double unc_down = 0;
0232     if (unc1 > 0) {
0233       uncertainty_cut_upper[ib-1] = TMath::Sqrt(uncertainty_cut_upper[ib-1] * uncertainty_cut_upper[ib-1] + unc1 * unc1);
0234       unc_up = TMath::Sqrt(unc_up * unc_up + unc1 * unc1);
0235     } else {
0236       uncertainty_cut_lower[ib-1] = TMath::Sqrt(uncertainty_cut_lower[ib-1] * uncertainty_cut_lower[ib-1] + unc1 * unc1);
0237       unc_down = TMath::Sqrt(unc_down * unc_down + unc1 * unc1);
0238     }
0239     if (unc2 > 0) {
0240       uncertainty_cut_upper[ib-1] = TMath::Sqrt(uncertainty_cut_upper[ib-1] * uncertainty_cut_upper[ib-1] + unc2 * unc2);
0241       unc_up = TMath::Sqrt(unc_up * unc_up + unc2 * unc2);
0242     } else {
0243       uncertainty_cut_lower[ib-1] = TMath::Sqrt(uncertainty_cut_lower[ib-1] * uncertainty_cut_lower[ib-1] + unc2 * unc2);
0244       unc_down = TMath::Sqrt(unc_down * unc_down + unc2 * unc2);
0245     }
0246     uncertainty_cut_upper_record[ib-1] = unc_up / h_unfold_calib_nominal->GetBinContent(ib);
0247     uncertainty_cut_lower_record[ib-1] = unc_down / h_unfold_calib_nominal->GetBinContent(ib);
0248     std::cout << unc_up / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%      \t" << unc_down / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%" << std::endl;
0249   }
0250 
0251   // Jet Trigger Efficiency Uncertainty
0252   h_input.push_back((TH1F*)h_unfold_calib_dijet_reweighted_1);
0253   h_input.push_back((TH1F*)h_unfold_calib_dijet_effdown_reweighted_1);
0254   h_input.push_back((TH1F*)h_unfold_calib_dijet_effup_reweighted_1);
0255   color.push_back(kBlack);
0256   color.push_back(kRed);
0257   color.push_back(kBlue);
0258   markerstyle.push_back(20);
0259   markerstyle.push_back(20);
0260   markerstyle.push_back(20);
0261   legend.push_back("Unfolded spectrum (Dijet cut)");
0262   legend.push_back("With jet trigger efficiency shift down");
0263   legend.push_back("With jet trigger efficiency shift up");
0264   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0.8, 1.2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_eff_dijet.pdf");
0265   h_input.clear();
0266   legend.clear();
0267   h_input.push_back((TH1F*)h_unfold_calib_frac_reweighted_1);
0268   h_input.push_back((TH1F*)h_unfold_calib_frac_effdown_reweighted_1);
0269   h_input.push_back((TH1F*)h_unfold_calib_frac_effup_reweighted_1);
0270   legend.push_back("Unfolded spectrum (Fraction cut)");
0271   legend.push_back("With jet trigger efficiency shift down");
0272   legend.push_back("With jet trigger efficiency shift up");
0273   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0.8, 1.2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_eff_frac.pdf");
0274   h_input.clear();
0275   color.clear();
0276   markerstyle.clear();
0277   legend.clear();
0278   std::cout << std::endl << "Uncertainty from jet trigger efficiency: " << std::endl << "up\t\t\tdown\t\t\tunc1\t\t\tunc2\t\t\tdijet1\t\t\tdijet2\t\t\tfrac1\t\t\tfrac2" << std::endl;
0279   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0280     double unc1_dijet = h_unfold_calib_dijet_effdown_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0281     double unc2_dijet = h_unfold_calib_dijet_effup_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0282     double unc1_frac = h_unfold_calib_frac_effdown_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0283     double unc2_frac = h_unfold_calib_frac_effup_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0284     double unc1 = (unc1_dijet + unc1_frac) / 2.;
0285     double unc2 = (unc2_dijet + unc2_frac) / 2.;
0286     if (unc1 > 0) {
0287       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc1 * unc1);
0288       uncertainty_jeteff_upper_record[ib-1] = TMath::Sqrt(uncertainty_jeteff_upper_record[ib-1] * uncertainty_jeteff_upper_record[ib-1] + unc1 * unc1);
0289     } else {
0290       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc1 * unc1);
0291       uncertainty_jeteff_lower_record[ib-1] = TMath::Sqrt(uncertainty_jeteff_lower_record[ib-1] * uncertainty_jeteff_lower_record[ib-1] + unc1 * unc1);
0292     }
0293     if (unc2 > 0) {
0294       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc2 * unc2);
0295       uncertainty_jeteff_upper_record[ib-1] = TMath::Sqrt(uncertainty_jeteff_upper_record[ib-1] * uncertainty_jeteff_upper_record[ib-1] + unc2 * unc2);
0296     } else {
0297       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc2 * unc2);
0298       uncertainty_jeteff_lower_record[ib-1] = TMath::Sqrt(uncertainty_jeteff_lower_record[ib-1] * uncertainty_jeteff_lower_record[ib-1] + unc2 * unc2);
0299     }
0300     uncertainty_jeteff_upper_record[ib-1] = uncertainty_jeteff_upper_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0301     uncertainty_jeteff_lower_record[ib-1] = uncertainty_jeteff_lower_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0302     std::cout << uncertainty_other_upper[ib-1] << " (" << uncertainty_other_upper[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << uncertainty_other_lower[ib-1] << " (" << uncertainty_other_lower[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0303     std::cout << unc1 << " (" << unc1 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << unc2 << " (" << unc2 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0304     std::cout << unc1_dijet << " (" << unc1_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_dijet << " (" << unc2_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t";
0305     std::cout << unc1_frac << " (" << unc1_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_frac << " (" << unc2_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%)" << std::endl;
0306   }
0307 
0308   // JES Uncertainty
0309   h_input.push_back((TH1F*)h_unfold_calib_dijet_reweighted_1);
0310   h_input.push_back((TH1F*)h_unfold_calib_dijet_jesdown_reweighted_1);
0311   h_input.push_back((TH1F*)h_unfold_calib_dijet_jesup_reweighted_1);
0312   color.push_back(kBlack);
0313   color.push_back(kRed);
0314   color.push_back(kBlue);
0315   markerstyle.push_back(20);
0316   markerstyle.push_back(20);
0317   markerstyle.push_back(20);
0318   legend.push_back("Unfolded spectrum (Dijet cut)");
0319   legend.push_back("With JES shift down");
0320   legend.push_back("With JES shift up");
0321   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0, 2.07, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_jes_dijet.pdf");
0322   h_input.clear();
0323   legend.clear();
0324   h_input.push_back((TH1F*)h_unfold_calib_frac_reweighted_1);
0325   h_input.push_back((TH1F*)h_unfold_calib_frac_jesdown_reweighted_1);
0326   h_input.push_back((TH1F*)h_unfold_calib_frac_jesup_reweighted_1);
0327   legend.push_back("Unfolded spectrum (Fraction cut)");
0328   legend.push_back("With JES shift down");
0329   legend.push_back("With JES shift up");
0330   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0, 2.07, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_jes_frac.pdf");
0331   h_input.clear();
0332   color.clear();
0333   markerstyle.clear();
0334   legend.clear();
0335   std::cout << std::endl << "Uncertainty from JES: " << std::endl << "up\t\t\tdown\t\t\tunc1\t\t\tunc2\t\t\tdijet1\t\t\tdijet2\t\t\tfrac1\t\t\tfrac2" << std::endl;
0336   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0337     double unc1_dijet = h_unfold_calib_dijet_jesdown_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0338     double unc2_dijet = h_unfold_calib_dijet_jesup_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0339     double unc1_frac = h_unfold_calib_frac_jesdown_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0340     double unc2_frac = h_unfold_calib_frac_jesup_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0341     double unc1 = (unc1_dijet + unc1_frac) / 2.;
0342     double unc2 = (unc2_dijet + unc2_frac) / 2.;
0343     if (unc1 > 0) {
0344       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc1 * unc1);
0345       uncertainty_jes_upper_record[ib-1] = TMath::Sqrt(uncertainty_jes_upper_record[ib-1] * uncertainty_jes_upper_record[ib-1] + unc1 * unc1);
0346     } else {
0347       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc1 * unc1);
0348       uncertainty_jes_lower_record[ib-1] = TMath::Sqrt(uncertainty_jes_lower_record[ib-1] * uncertainty_jes_lower_record[ib-1] + unc1 * unc1);
0349     }
0350     if (unc2 > 0) {
0351       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc2 * unc2);
0352       uncertainty_jes_upper_record[ib-1] = TMath::Sqrt(uncertainty_jes_upper_record[ib-1] * uncertainty_jes_upper_record[ib-1] + unc2 * unc2);
0353     } else {
0354       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc2 * unc2);
0355       uncertainty_jes_lower_record[ib-1] = TMath::Sqrt(uncertainty_jes_lower_record[ib-1] * uncertainty_jes_lower_record[ib-1] + unc2 * unc2);
0356     }
0357     uncertainty_jes_upper_record[ib-1] = uncertainty_jes_upper_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0358     uncertainty_jes_lower_record[ib-1] = uncertainty_jes_lower_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0359     std::cout << uncertainty_other_upper[ib-1] << " (" << uncertainty_other_upper[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << uncertainty_other_lower[ib-1] << " (" << uncertainty_other_lower[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0360     std::cout << unc1 << " (" << unc1 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << unc2 << " (" << unc2 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0361     std::cout << unc1_dijet << " (" << unc1_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_dijet << " (" << unc2_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t";
0362     std::cout << unc1_frac << " (" << unc1_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_frac << " (" << unc2_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%)" << std::endl;
0363   }
0364 
0365   // JER Uncertainty
0366   h_input.push_back((TH1F*)h_unfold_calib_dijet_reweighted_1);
0367   h_input.push_back((TH1F*)h_unfold_calib_dijet_jerdown_reweighted_1);
0368   h_input.push_back((TH1F*)h_unfold_calib_dijet_jerup_reweighted_1);
0369   color.push_back(kBlack);
0370   color.push_back(kRed);
0371   color.push_back(kBlue);
0372   markerstyle.push_back(20);
0373   markerstyle.push_back(20);
0374   markerstyle.push_back(20);
0375   legend.push_back("Unfolded spectrum (Dijet cut)");
0376   legend.push_back("With JER shift down");
0377   legend.push_back("With JER shift up");
0378   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0, 2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_jer_dijet.pdf");
0379   h_input.clear();
0380   legend.clear();
0381   h_input.push_back((TH1F*)h_unfold_calib_frac_reweighted_1);
0382   h_input.push_back((TH1F*)h_unfold_calib_frac_jerdown_reweighted_1);
0383   h_input.push_back((TH1F*)h_unfold_calib_frac_jerup_reweighted_1);
0384   legend.push_back("Unfolded spectrum (Fraction cut)");
0385   legend.push_back("With JER shift down");
0386   legend.push_back("With JER shift up");
0387   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0, 2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_jer_frac.pdf");
0388   h_input.clear();
0389   color.clear();
0390   markerstyle.clear();
0391   legend.clear();
0392   std::cout << std::endl << "Uncertainty from JER: " << std::endl << "up\t\t\tdown\t\t\tunc1\t\t\tunc2\t\t\tdijet1\t\t\tdijet2\t\t\tfrac1\t\t\tfrac2" << std::endl;
0393   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0394     double unc1_dijet = h_unfold_calib_dijet_jerdown_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0395     double unc2_dijet = h_unfold_calib_dijet_jerup_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0396     double unc1_frac = h_unfold_calib_frac_jerdown_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0397     double unc2_frac = h_unfold_calib_frac_jerup_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0398     double unc1 = (unc1_dijet + unc1_frac) / 2.;
0399     double unc2 = (unc2_dijet + unc2_frac) / 2.;
0400     if (unc1 > 0) {
0401       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc1 * unc1);
0402       uncertainty_jer_upper_record[ib-1] = TMath::Sqrt(uncertainty_jer_upper_record[ib-1] * uncertainty_jer_upper_record[ib-1] + unc1 * unc1);
0403     } else {
0404       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc1 * unc1);
0405       uncertainty_jer_lower_record[ib-1] = TMath::Sqrt(uncertainty_jer_lower_record[ib-1] * uncertainty_jer_lower_record[ib-1] + unc1 * unc1);
0406     }
0407     if (unc2 > 0) {
0408       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc2 * unc2);
0409       uncertainty_jer_upper_record[ib-1] = TMath::Sqrt(uncertainty_jer_upper_record[ib-1] * uncertainty_jer_upper_record[ib-1] + unc2 * unc2);
0410     } else {
0411       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc2 * unc2);
0412       uncertainty_jer_lower_record[ib-1] = TMath::Sqrt(uncertainty_jer_lower_record[ib-1] * uncertainty_jer_lower_record[ib-1] + unc2 * unc2);
0413     }
0414     uncertainty_jer_upper_record[ib-1] = uncertainty_jer_upper_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0415     uncertainty_jer_lower_record[ib-1] = uncertainty_jer_lower_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0416     std::cout << uncertainty_other_upper[ib-1] << " (" << uncertainty_other_upper[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << uncertainty_other_lower[ib-1] << " (" << uncertainty_other_lower[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0417     std::cout << unc1 << " (" << unc1 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << unc2 << " (" << unc2 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0418     std::cout << unc1_dijet << " (" << unc1_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_dijet << " (" << unc2_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t";
0419     std::cout << unc1_frac << " (" << unc1_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_frac << " (" << unc2_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%)" << std::endl;
0420   }
0421 
0422   // MBD Trigger Efficiency Uncertainty
0423   h_input.push_back((TH1F*)h_unfold_calib_dijet_reweighted_1);
0424   h_input.push_back((TH1F*)h_unfold_calib_dijet_mbddown_reweighted_1);
0425   h_input.push_back((TH1F*)h_unfold_calib_dijet_mbdup_reweighted_1);
0426   color.push_back(kBlack);
0427   color.push_back(kRed);
0428   color.push_back(kBlue);
0429   markerstyle.push_back(20);
0430   markerstyle.push_back(20);
0431   markerstyle.push_back(20);
0432   legend.push_back("Unfolded spectrum (Dijet cut)");
0433   legend.push_back("With MBD trigger efficiency shift down");
0434   legend.push_back("With MBD trigger efficiency shift up");
0435   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0.8, 1.2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_mbd_dijet.pdf");
0436   h_input.clear();
0437   legend.clear();
0438   h_input.push_back((TH1F*)h_unfold_calib_frac_reweighted_1);
0439   h_input.push_back((TH1F*)h_unfold_calib_frac_mbddown_reweighted_1);
0440   h_input.push_back((TH1F*)h_unfold_calib_frac_mbdup_reweighted_1);
0441   legend.push_back("Unfolded spectrum (Fraction cut)");
0442   legend.push_back("With MBD trigger efficiency shift down");
0443   legend.push_back("With MBD trigger efficiency shift up");
0444   draw_1D_multiple_plot_uncertainty(h_input, calibptbins[0], calibptbins[calibnpt], 0.8, 1.2, "p_{T}^{jet} [GeV]", "d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]", legend, 0.5, 0.7, color, markerstyle, "figure/specturm_uncertainty_mbd_frac.pdf");
0445   h_input.clear();
0446   color.clear();
0447   markerstyle.clear();
0448   legend.clear();
0449   std::cout << std::endl << "Uncertainty from MBD trigger efficiency: " << std::endl << "up\t\t\tdown\t\t\tunc1\t\t\tunc2\t\t\tdijet1\t\t\tdijet2\t\t\tfrac1\t\t\tfrac2" << std::endl;
0450   for (int ib = 1; ib <= h_unfold_calib_nominal->GetNbinsX(); ++ib) {
0451     double unc1_dijet = h_unfold_calib_dijet_mbddown_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0452     double unc2_dijet = h_unfold_calib_dijet_mbdup_reweighted_1->GetBinContent(ib) - h_unfold_calib_dijet_reweighted_1->GetBinContent(ib);
0453     double unc1_frac = h_unfold_calib_frac_mbddown_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0454     double unc2_frac = h_unfold_calib_frac_mbdup_reweighted_1->GetBinContent(ib) - h_unfold_calib_frac_reweighted_1->GetBinContent(ib);
0455     double unc1 = (unc1_dijet + unc1_frac) / 2.;
0456     double unc2 = (unc2_dijet + unc2_frac) / 2.;
0457     if (unc1 > 0) {
0458       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc1 * unc1);
0459       uncertainty_mbdeff_upper_record[ib-1] = TMath::Sqrt(uncertainty_mbdeff_upper_record[ib-1] * uncertainty_mbdeff_upper_record[ib-1] + unc1 * unc1);
0460     } else {
0461       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc1 * unc1);
0462       uncertainty_mbdeff_lower_record[ib-1] = TMath::Sqrt(uncertainty_mbdeff_lower_record[ib-1] * uncertainty_mbdeff_lower_record[ib-1] + unc1 * unc1);
0463     }
0464     if (unc2 > 0) {
0465       uncertainty_other_upper[ib-1] = TMath::Sqrt(uncertainty_other_upper[ib-1] * uncertainty_other_upper[ib-1] + unc2 * unc2);
0466       uncertainty_mbdeff_upper_record[ib-1] = TMath::Sqrt(uncertainty_mbdeff_upper_record[ib-1] * uncertainty_mbdeff_upper_record[ib-1] + unc2 * unc2);
0467     } else {
0468       uncertainty_other_lower[ib-1] = TMath::Sqrt(uncertainty_other_lower[ib-1] * uncertainty_other_lower[ib-1] + unc2 * unc2);
0469       uncertainty_mbdeff_lower_record[ib-1] = TMath::Sqrt(uncertainty_mbdeff_lower_record[ib-1] * uncertainty_mbdeff_lower_record[ib-1] + unc2 * unc2);
0470     }
0471     uncertainty_mbdeff_upper_record[ib-1] = uncertainty_mbdeff_upper_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0472     uncertainty_mbdeff_lower_record[ib-1] = uncertainty_mbdeff_lower_record[ib-1] / h_unfold_calib_nominal->GetBinContent(ib);
0473     std::cout << uncertainty_other_upper[ib-1] << " (" << uncertainty_other_upper[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << uncertainty_other_lower[ib-1] << " (" << uncertainty_other_lower[ib-1] / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0474     std::cout << unc1 << " (" << unc1 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t" << unc2 << " (" << unc2 / h_unfold_calib_nominal->GetBinContent(ib) * 100 << "\%) \t";
0475     std::cout << unc1_dijet << " (" << unc1_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_dijet << " (" << unc2_dijet / h_unfold_calib_dijet_reweighted_1->GetBinContent(ib) * 100 << "\%) \t";
0476     std::cout << unc1_frac << " (" << unc1_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%) \t" << unc2_frac << " (" << unc2_frac / h_unfold_calib_frac_reweighted_1->GetBinContent(ib) * 100 << "\%)" << std::endl;
0477   }
0478 
0479   // Get total uncertainty
0480   std::cout << std::endl << "Total uncertainty with background cut and unfolding uncertainty: " << std::endl << "nomiY\t\t\tnoup\t\t\tnodown" << std::endl;
0481   for (int ib = 0; ib < calibnpt; ++ib) {
0482     uncertainty_final_upper[ib] = TMath::Sqrt(uncertainty_other_upper[ib] * uncertainty_other_upper[ib] + uncertainty_cut_upper[ib] * uncertainty_cut_upper[ib] + uncertainty_unfold_upper[ib] * uncertainty_unfold_upper[ib]);
0483     uncertainty_final_lower[ib] = TMath::Sqrt(uncertainty_other_lower[ib] * uncertainty_other_lower[ib] + uncertainty_cut_lower[ib] * uncertainty_cut_lower[ib] + uncertainty_unfold_lower[ib] * uncertainty_unfold_lower[ib]);
0484     std::cout << h_unfold_calib_nominal->GetBinContent(ib+1) << "  \t\t" << uncertainty_final_upper[ib] << " (" << uncertainty_final_upper[ib] / h_unfold_calib_nominal->GetBinContent(ib+1) * 100 << "\%)\t" << uncertainty_final_lower[ib] << " (" << uncertainty_final_lower[ib] / h_unfold_calib_nominal->GetBinContent(ib+1) * 100 << "\%)" << std::endl;
0485   }
0486 
0487   std::cout << std::endl << std::endl;
0488   std::cout << "pT range\td^{2}#sigma/(d#eta dp_{T})\t#delta#sigma_{stat}\t\t#delta#sigma_{sys}^{low}\t#delta#sigma_{sys}^{high}" << std::endl;
0489   std::cout << "[GeV]\t\t[pb/GeV]\t\t\t[pb/GeV]\t\t\t[pb/GeV]\t\t\t[pb/GeV]" << std::endl;
0490   double sampled_N = luminosity * 26.1;
0491   for (int ib = 0; ib < calibnpt; ++ib) {
0492     h_unfold_calib_nominal->SetBinContent(ib+1, h_unfold_calib_nominal->GetBinContent(ib+1) / (float)(sampled_N * delta_eta * h_unfold_calib_nominal->GetBinWidth(ib+1)));
0493     h_unfold_calib_nominal->SetBinError(ib+1, h_unfold_calib_nominal->GetBinError(ib+1) / (float)(sampled_N * delta_eta * h_unfold_calib_nominal->GetBinWidth(ib+1)));
0494     uncertainty_final_upper[ib] = uncertainty_final_upper[ib] / (float)(sampled_N * delta_eta * h_unfold_calib_nominal->GetBinWidth(ib+1)) / (float)h_unfold_calib_nominal->GetBinContent(ib+1);
0495     uncertainty_final_lower[ib] = uncertainty_final_lower[ib] / (float)(sampled_N * delta_eta * h_unfold_calib_nominal->GetBinWidth(ib+1)) / (float)h_unfold_calib_nominal->GetBinContent(ib+1);
0496     h_unfold_calib_nominal->SetBinContent(ib+1, h_unfold_calib_nominal->GetBinContent(ib+1) * 26.1);
0497     h_unfold_calib_nominal->SetBinError(ib+1, h_unfold_calib_nominal->GetBinError(ib+1) * 26.1);
0498     uncertainty_final_upper[ib] = h_unfold_calib_nominal->GetBinContent(ib+1) * TMath::Sqrt(uncertainty_final_upper[ib] * uncertainty_final_upper[ib] + uncertainty_mbd_upper[ib] * uncertainty_mbd_upper[ib]);
0499     uncertainty_final_lower[ib] = h_unfold_calib_nominal->GetBinContent(ib+1) * TMath::Sqrt(uncertainty_final_lower[ib] * uncertainty_final_lower[ib] + uncertainty_mbd_lower[ib] * uncertainty_mbd_lower[ib]);
0500     std::cout << calibptbins[ib] << " - " << calibptbins[ib+1] << "  \t" << h_unfold_calib_nominal->GetBinContent(ib+1) << "  \t\t\t" << h_unfold_calib_nominal->GetBinError(ib+1) << " (" << h_unfold_calib_nominal->GetBinError(ib+1) / h_unfold_calib_nominal->GetBinContent(ib+1) * 100 << "\%)\t\t" << uncertainty_final_lower[ib] << " (" << uncertainty_final_lower[ib] / h_unfold_calib_nominal->GetBinContent(ib+1) * 100 << "\%)\t\t" << uncertainty_final_upper[ib] << " (" << uncertainty_final_upper[ib] / h_unfold_calib_nominal->GetBinContent(ib+1) * 100 << "\%)" << std::endl;
0501   }
0502 
0503   // Draw the final spectrum
0504   // sPHENIX result
0505   double g_sphenix_x[calibnpt], g_sphenix_xerrdown[calibnpt], g_sphenix_xerrup[calibnpt], g_sphenix_y[calibnpt], g_sphenix_yerrdown[calibnpt], g_sphenix_yerrup[calibnpt], g_sphenix_yerrstat[calibnpt];
0506   for (int ib = 0; ib < calibnpt; ++ ib) {
0507     g_sphenix_x[ib] = (calibptbins[ib] + calibptbins[ib+1]) / 2.;
0508     g_sphenix_xerrdown[ib] = (calibptbins[ib+1] - calibptbins[ib]) / 2.;
0509     g_sphenix_xerrup[ib] = (calibptbins[ib+1] - calibptbins[ib]) / 2.;
0510     g_sphenix_y[ib] = h_unfold_calib_nominal->GetBinContent(ib+1);
0511     g_sphenix_yerrstat[ib] = h_unfold_calib_nominal->GetBinError(ib+1);
0512     g_sphenix_yerrdown[ib] = TMath::Sqrt(g_sphenix_yerrstat[ib]*g_sphenix_yerrstat[ib] + uncertainty_final_lower[ib]*uncertainty_final_lower[ib]);
0513     g_sphenix_yerrup[ib] = TMath::Sqrt(g_sphenix_yerrstat[ib]*g_sphenix_yerrstat[ib] + uncertainty_final_upper[ib]*uncertainty_final_upper[ib]);
0514 
0515     uncertainty_stat_upper_record[ib] = g_sphenix_yerrstat[ib] / g_sphenix_y[ib];
0516     uncertainty_stat_lower_record[ib] = g_sphenix_yerrstat[ib] / g_sphenix_y[ib];
0517     uncertainty_total_upper_record[ib] = g_sphenix_yerrup[ib] / g_sphenix_y[ib];
0518     uncertainty_total_lower_record[ib] = g_sphenix_yerrdown[ib] / g_sphenix_y[ib];
0519   }
0520   TGraphAsymmErrors* g_sphenix_result = new TGraphAsymmErrors(calibnpt, g_sphenix_x, g_sphenix_y, g_sphenix_xerrdown, g_sphenix_xerrup, g_sphenix_yerrdown, g_sphenix_yerrup);
0521   g_sphenix_result->SetMarkerStyle(20);
0522   g_sphenix_result->SetMarkerColor(kAzure+2);
0523   g_sphenix_result->SetLineColor(kAzure+2);
0524   g_sphenix_result->SetFillColorAlpha(kAzure + 1, 0.5);
0525   TGraphErrors* g_sphenix_result_point = new TGraphErrors(calibnpt, g_sphenix_x, g_sphenix_y, g_sphenix_xerrdown, g_sphenix_yerrstat);
0526   g_sphenix_result_point->SetMarkerStyle(20);
0527   g_sphenix_result_point->SetMarkerColor(kAzure+2);
0528   g_sphenix_result_point->SetLineColor(kAzure+2);
0529   g_sphenix_result_point->SetLineWidth(2);
0530 
0531   // Pythia result
0532   TFile *f_sim_specturm = new TFile("output_sim.root", "READ");
0533   TH1D* h_pythia_spectrum_dijet = (TH1D*)f_sim_specturm->Get("h_truth_calib_dijet");
0534   TH1D* h_pythia_spectrum_frac = (TH1D*)f_sim_specturm->Get("h_truth_calib_frac");
0535   TH1D* h_pythia_spectrum = (TH1D*)h_pythia_spectrum_dijet->Clone("h_pythia_spectrum");
0536   h_pythia_spectrum->Add(h_pythia_spectrum_frac);
0537   h_pythia_spectrum->Scale(1. / 2.);
0538   h_pythia_spectrum->Scale(1e+06); // N_event is scaled down by 1e+06 in analysis sim.
0539   for (int ib = 1; ib <= h_pythia_spectrum->GetNbinsX(); ++ib) {
0540     h_pythia_spectrum->SetBinContent(ib, h_pythia_spectrum->GetBinContent(ib)/ (float)(delta_eta * h_pythia_spectrum->GetBinWidth(ib)));
0541     h_pythia_spectrum->SetBinError(ib, h_pythia_spectrum->GetBinError(ib)/ (float)(delta_eta * h_pythia_spectrum->GetBinWidth(ib)));
0542   }
0543   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];
0544   for (int ib = 0; ib < truthnpt; ++ib) {
0545     g_pythia_x[ib] = (truthptbins[ib] + truthptbins[ib+1]) / 2;
0546     g_pythia_xerrdown[ib] = (truthptbins[ib+1] - truthptbins[ib]) / 2;
0547     g_pythia_xerrup[ib] = (truthptbins[ib+1] - truthptbins[ib]) / 2;
0548     g_pythia_y[ib] = h_pythia_spectrum->GetBinContent(ib+1);
0549     g_pythia_yerrdown[ib] = h_pythia_spectrum->GetBinError(ib+1);
0550     g_pythia_yerrup[ib] = h_pythia_spectrum->GetBinError(ib+1);
0551   }
0552   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);
0553   g_pythia_result->SetMarkerStyle(20);
0554   g_pythia_result->SetMarkerColor(kOrange);
0555   g_pythia_result->SetLineColor(kOrange);
0556 
0557   // Pythia ratio
0558   TH1D* h_pythia_ratio = (TH1D*)h_unfold_calib_nominal->Clone("h_pythia_ratio");
0559   for (int ib = 1; ib <= h_pythia_ratio->GetNbinsX(); ++ib) {
0560     h_pythia_ratio->SetBinContent(ib, h_pythia_spectrum->GetBinContent(ib+n_underflow_bin) / h_pythia_ratio->GetBinContent(ib));
0561     h_pythia_ratio->SetBinError(ib, h_pythia_ratio->GetBinContent(ib) * TMath::Sqrt(TMath::Power(h_pythia_spectrum->GetBinError(ib+n_underflow_bin) / h_pythia_spectrum->GetBinContent(ib+n_underflow_bin), 2) + TMath::Power(h_unfold_calib_nominal->GetBinError(ib) / h_unfold_calib_nominal->GetBinContent(ib), 2)));
0562   }
0563 
0564   // STAR old result
0565   double starpointxold[] = {8.3, 10.3, 12.6, 15.5, 19.0, 23.4, 28.7, 35.3, 43.3};
0566   const static int starnpointold = sizeof(starpointxold) / sizeof(starpointxold[0]);
0567   double starpointxerrdownold[starnpointold] = {0.7, 1.0, 1.2, 1.4, 1.7, 2.1, 2.5, 3.1, 3.7};
0568   double starpointxerruppld[starnpointold] = {1.0, 1.1, 1.5, 1.8, 2.3, 2.8, 3.5, 4.3, 5.4};
0569   double starpointyold[] = {6.40e+05*2*TMath::Pi(), 2.40e+05*2*TMath::Pi(), 5.90e+04*2*TMath::Pi(), 1.19e+04*2*TMath::Pi(), 3.50e+03*2*TMath::Pi(), 7.80e+02*2*TMath::Pi(), 1.41e+02*2*TMath::Pi(), 2.20e+01*2*TMath::Pi(), 2.60e+00*2*TMath::Pi()};
0570   double starpointyerrdownold[] = {0.90e+05, 0.30e+05, 0.60e+04, 0.10e+04, 0.20e+03, 0.40e+02, 0.09e+02, 0.30e+01, 0.60e+00};
0571   double starpointyerrupold[] = {3.10e+05, 1.10e+05, 2.80e+04, 0.57e+04, 1.70e+03, 0.37e+02, 0.68e+02, 1.10e+01, 1.30e+00};
0572   TGraphAsymmErrors* g_star_result_old = new TGraphAsymmErrors(starnpointold, starpointxold, starpointyold, starpointxerrdownold, starpointxerruppld, starpointyerrdownold, starpointyerrupold);
0573   g_star_result_old->SetMarkerStyle(20);
0574   g_star_result_old->SetMarkerColor(kRed);
0575   g_star_result_old->SetLineColor(kRed);
0576 
0577   // STAR result
0578   double starpointx[] = {7.8, 10.4, 13.85, 18.45, 24.6, 32.8, 43.75};
0579   const static int starnpoint = sizeof(starpointx) / sizeof(starpointx[0]);
0580   double starpointxerrdown[starnpoint] = {1.1, 1.5, 1.95, 2.65, 3.5, 4.7, 6.25};
0581   double starpointxerrup[starnpoint] = {1.1, 1.5, 1.95, 2.65, 3.5, 4.7, 6.25};
0582   double starpointy[starnpoint] = {5.655596e+06, 9.684281e+05, 1.603961e+05, 2.641261e+04, 3.391784e+03, 2.871360e+02, 2.178426e+01};
0583   double starpointyerrdown[starnpoint] = {7.270251e+05, 1.409256e+05, 1.897396e+04, 3.829357e+03, 5.537599e+02, 5.478192e+01, 5.218297e+00};
0584   double starpointyerrup[starnpoint] = {7.486261e+05, 1.452015e+05, 1.897113e+04, 4.021414e+03, 5.756504e+02, 6.336199e+01, 5.765260e+00};
0585   TGraphAsymmErrors* g_star_result = new TGraphAsymmErrors(starnpoint, starpointx, starpointy, starpointxerrdown, starpointxerrup, starpointyerrdown, starpointyerrup);
0586   g_star_result->SetMarkerStyle(20);
0587   g_star_result->SetMarkerColor(kRed);
0588   g_star_result->SetLineColor(kRed);
0589 
0590   // PHENIX result
0591   double phenixpointx[] = {8.45355, 9.45319, 10.8252, 13.0123, 15.7055, 18.739, 22.0908, 26.2576, 31.1771, 37.4739};
0592   const static int phenixnpoint = sizeof(phenixpointx) / sizeof(phenixpointx[0]);
0593   double phenixpointxerrdown[phenixnpoint] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0594   double phenixpointxerrup[phenixnpoint] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0595   double phenixpointy[phenixnpoint] = {5.86972e+05, 2.99829e+05, 1.3335e+05, 4.07571e+04, 1.24647e+04, 4.28294e+03, 1.39962e+03, 4.11469e+02, 1.01631e+02, 2.49519e+01};
0596   double phenixpointysysdown[phenixnpoint] = {2.32513e+05, 1.08448e+05, 4.80741e+04, 1.44942e+04, 4.19428e+03, 1.21361e+03, 3.87209e+02, 9.41715e+01, 1.20866e+01, 1.06009};
0597   double phenixpointysysup[phenixnpoint] = {1.90074e+05, 9.64104e+04, 4.43763e+04, 1.42428e+04, 4.28854e+03, 1.30758e+03, 4.68572e+02, 1.15733e+02, 1.5672e+01, 2.62781};
0598   double phenixpointystat[phenixnpoint] = {3.5247e+03, 2.16461e+03, 1.15278e+03, 4.6706e+02, 1.86848e+02, 1.04529e+02, 4.38652e+01, 2.17658e+01, 9.3723, 4.46448};
0599   double phenixpointyerrdown[phenixnpoint];
0600   double phenixpointyerrup[phenixnpoint];
0601   for (int i = 0; i < phenixnpoint; ++i) {
0602     phenixpointyerrdown[i] = TMath::Sqrt(phenixpointystat[i]*phenixpointystat[i] + phenixpointysysdown[i]*phenixpointysysdown[i]);
0603     phenixpointyerrup[i] = TMath::Sqrt(phenixpointystat[i]*phenixpointystat[i] + phenixpointysysup[i]*phenixpointysysup[i]);
0604   }
0605   TGraphAsymmErrors* g_phenix_result = new TGraphAsymmErrors(phenixnpoint, phenixpointx, phenixpointy, phenixpointxerrdown, phenixpointxerrup, phenixpointyerrdown, phenixpointyerrup);
0606   g_phenix_result->SetMarkerStyle(20);
0607   g_phenix_result->SetMarkerColor(kBlue);
0608   g_phenix_result->SetLineColor(kBlue);
0609 
0610   // NLO result
0611   ifstream myfile1, myfile05, myfile2;
0612   TGraph *tjet = new TGraph();
0613   TGraph *tjet05 = new TGraph();
0614   TGraph *tjet2 = new TGraph();
0615   double scalefactor = 2.0 * 3.14159;
0616   myfile1.open("sphenix_nlo/sPhenix-sc1.dat");
0617   myfile05.open("sphenix_nlo/sPhenix-sc05.dat");
0618   myfile2.open("sphenix_nlo/sPhenix-sc2.dat");
0619   double drop1, drop2, drop3, pt1, yield1, pt05, yield05, pt2, yield2;
0620   for (int index=0; index<66; index++) {
0621     myfile1 >> drop1 >> drop2 >> drop3 >> pt1 >> yield1;
0622     yield1 = yield1 * scalefactor * pt1;
0623     tjet->SetPoint(index, pt1, yield1);
0624 
0625     myfile05 >> drop1 >> drop2 >> drop3 >> pt05 >> yield05;
0626     yield05 = yield05 * scalefactor * pt05;
0627     tjet05->SetPoint(index, pt05, yield05);
0628 
0629     myfile2 >> drop1 >> drop2 >> drop3 >> pt2 >> yield2;
0630     yield2 = yield2 * scalefactor * pt2;
0631     tjet2->SetPoint(index, pt2, yield2);
0632   }
0633   myfile1.close();
0634   myfile05.close();
0635   myfile2.close();
0636   tjet->SetMarkerStyle(21);
0637   tjet->SetMarkerColor(2);
0638   tjet->SetLineColor(2);
0639   tjet05->SetMarkerStyle(21);
0640   tjet05->SetMarkerColor(2);
0641   tjet05->SetLineColor(2);
0642   tjet2->SetMarkerStyle(21);
0643   tjet2->SetMarkerColor(2);
0644   tjet2->SetLineColor(2);
0645 
0646   TCanvas *can_nlo = new TCanvas("can_nlo", "", 850, 800);
0647   gStyle->SetPalette(57);
0648   can_nlo->SetTopMargin(0.03);
0649   can_nlo->SetLeftMargin(0.15);
0650   can_nlo->SetBottomMargin(0.12);
0651   can_nlo->SetRightMargin(0.05);
0652   can_nlo->SetLogy();
0653   TH2F *frame = new TH2F("frame", "", 10, 15, 72, 10, 1e-03, 2e+05);
0654   frame->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0655   frame->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0656   frame->GetXaxis()->SetTitleOffset(0.98);
0657   frame->GetYaxis()->SetTitleOffset(1.15);
0658   frame->GetXaxis()->SetLabelSize(0.045);
0659   frame->GetYaxis()->SetLabelSize(0.045);
0660   frame->GetXaxis()->CenterTitle();
0661   frame->GetYaxis()->CenterTitle();
0662   frame->Draw();
0663   g_sphenix_result->SetMarkerStyle(20);
0664   g_sphenix_result->SetMarkerColor(kAzure+2);
0665   g_sphenix_result->SetLineColor(kAzure+2);
0666   g_sphenix_result->SetFillColorAlpha(kAzure + 1, 0.5);
0667   g_sphenix_result->Draw("2 same");
0668   g_sphenix_result_point->SetMarkerStyle(20);
0669   g_sphenix_result_point->SetMarkerColor(kAzure+2);
0670   g_sphenix_result_point->SetLineColor(kAzure+2);
0671   g_sphenix_result_point->SetLineWidth(2);
0672   g_sphenix_result_point->Draw("P same");
0673   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0674   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0675   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0676   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0677   TLegend *leg = new TLegend(0.2, 0.2, 0.55, 0.28);
0678   leg->SetBorderSize(0);
0679   leg->SetFillStyle(0);
0680   leg->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0681   leg->AddEntry(tjet, "NLO pQCD", "pl");
0682   leg->Draw();
0683   tjet->Draw("pl same");
0684   tjet05->Draw("l same");
0685   tjet2->Draw("l same");
0686   can_nlo->SaveAs("figure/spectrum_calib_nlo.pdf");
0687 
0688   TCanvas *can_sphenix = new TCanvas("can_sphenix", "", 850, 800);
0689   gStyle->SetPalette(57);
0690   can_sphenix->SetTopMargin(0.03);
0691   can_sphenix->SetLeftMargin(0.15);
0692   can_sphenix->SetBottomMargin(0.12);
0693   can_sphenix->SetRightMargin(0.05);
0694   can_sphenix->SetLogy();
0695   TH2F *frame_sphenix = new TH2F("frame_sphenix", "", 10, 15, 72, 10, 1e-03, 2e+05);
0696   frame_sphenix->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0697   frame_sphenix->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0698   frame_sphenix->GetXaxis()->SetTitleOffset(0.98);
0699   frame_sphenix->GetYaxis()->SetTitleOffset(1.15);
0700   frame_sphenix->GetXaxis()->SetLabelSize(0.045);
0701   frame_sphenix->GetYaxis()->SetLabelSize(0.045);
0702   frame_sphenix->GetXaxis()->CenterTitle();
0703   frame_sphenix->GetYaxis()->CenterTitle();
0704   frame_sphenix->Draw();
0705   g_sphenix_result->Draw("2 same");
0706   g_sphenix_result_point->Draw("P same");
0707   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0708   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0709   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0710   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0711   TLegend *leg_sphenix = new TLegend(0.2, 0.2, 0.7, 0.26);
0712   leg_sphenix->SetBorderSize(0);
0713   leg_sphenix->SetFillStyle(0);
0714   leg_sphenix->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0715   leg_sphenix->Draw();
0716   can_sphenix->SaveAs("figure/spectrum_calib_sphenix.pdf");
0717 
0718   TCanvas *can_data = new TCanvas("can_data", "", 850, 800);
0719   gStyle->SetPalette(57);
0720   can_data->SetTopMargin(0.03);
0721   can_data->SetLeftMargin(0.15);
0722   can_data->SetBottomMargin(0.12);
0723   can_data->SetRightMargin(0.05);
0724   can_data->SetLogy();
0725   TH2F *frame_data = new TH2F("frame_data", "", 10, 15, 72, 10, 1e-03, 2e+05);
0726   frame_data->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0727   frame_data->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0728   frame_data->GetXaxis()->SetTitleOffset(0.98);
0729   frame_data->GetYaxis()->SetTitleOffset(1.15);
0730   frame_data->GetXaxis()->SetLabelSize(0.045);
0731   frame_data->GetYaxis()->SetLabelSize(0.045);
0732   frame_data->GetXaxis()->CenterTitle();
0733   frame_data->GetYaxis()->CenterTitle();
0734   frame_data->Draw();
0735   g_sphenix_result->Draw("2 same");
0736   g_sphenix_result_point->Draw("P same");
0737   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0738   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0739   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0740   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0741   TLegend *leg_data = new TLegend(0.2, 0.2, 0.7, 0.32);
0742   leg_data->SetBorderSize(0);
0743   leg_data->SetFillStyle(0);
0744   leg_data->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0745   leg_data->AddEntry(g_star_result, "STAR Run 2012 Preliminary, anti-k_{T} R=0.6", "pl");
0746   leg_data->AddEntry(g_phenix_result, "PHENIX Run12+15, anti-k_{T} R=0.3", "pl");
0747   leg_data->Draw();
0748   g_star_result->Draw("p same");
0749   g_phenix_result->Draw("p same");
0750   can_data->SaveAs("figure/spectrum_calib_datacompare.pdf");
0751 
0752   TCanvas *can_dataold = new TCanvas("can_dataold", "", 850, 800);
0753   gStyle->SetPalette(57);
0754   can_dataold->SetTopMargin(0.03);
0755   can_dataold->SetLeftMargin(0.15);
0756   can_dataold->SetBottomMargin(0.12);
0757   can_dataold->SetRightMargin(0.05);
0758   can_dataold->SetLogy();
0759   TH2F *frame_dataold = new TH2F("frame_dataold", "", 10, 15, 72, 10, 1e-03, 2e+05);
0760   frame_dataold->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0761   frame_dataold->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0762   frame_dataold->GetXaxis()->SetTitleOffset(0.98);
0763   frame_dataold->GetYaxis()->SetTitleOffset(1.15);
0764   frame_dataold->GetXaxis()->SetLabelSize(0.045);
0765   frame_dataold->GetYaxis()->SetLabelSize(0.045);
0766   frame_dataold->GetXaxis()->CenterTitle();
0767   frame_dataold->GetYaxis()->CenterTitle();
0768   frame_dataold->Draw();
0769   g_sphenix_result->Draw("2 same");
0770   g_sphenix_result_point->Draw("P same");
0771   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0772   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0773   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0774   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0775   TLegend *leg_dataold = new TLegend(0.2, 0.2, 0.7, 0.32);
0776   leg_dataold->SetBorderSize(0);
0777   leg_dataold->SetFillStyle(0);
0778   leg_dataold->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0779   leg_dataold->AddEntry(g_star_result_old, "STAR Midpoint Cone#it{R}=0.4 (Combined HT)", "pl");
0780   leg_dataold->AddEntry(g_phenix_result, "PHENIX Run12+15, anti-k_{T} R=0.3", "pl");
0781   leg_dataold->Draw();
0782   g_star_result_old->Draw("p same");
0783   g_phenix_result->Draw("p same");
0784   can_dataold->SaveAs("figure/spectrum_calib_datacompareold.pdf");
0785 
0786   TCanvas *can_sim = new TCanvas("can_sim", "", 850, 800);
0787   gStyle->SetPalette(57);
0788   can_sim->SetTopMargin(0.03);
0789   can_sim->SetLeftMargin(0.15);
0790   can_sim->SetBottomMargin(0.12);
0791   can_sim->SetRightMargin(0.05);
0792   can_sim->SetLogy();
0793   TH2F *frame_sim = new TH2F("frame_sim", "", 10, 15, 72, 10, 1e-03, 2e+05);
0794   frame_sim->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0795   frame_sim->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0796   frame_sim->GetXaxis()->SetTitleOffset(0.98);
0797   frame_sim->GetYaxis()->SetTitleOffset(1.15);
0798   frame_sim->GetXaxis()->SetLabelSize(0.045);
0799   frame_sim->GetYaxis()->SetLabelSize(0.045);
0800   frame_sim->GetXaxis()->CenterTitle();
0801   frame_sim->GetYaxis()->CenterTitle();
0802   frame_sim->Draw();
0803   g_sphenix_result->Draw("2 same");
0804   g_sphenix_result_point->Draw("P same");
0805   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0806   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0807   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0808   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0809   TLegend *leg_sim = new TLegend(0.2, 0.2, 0.55, 0.28);
0810   leg_sim->SetBorderSize(0);
0811   leg_sim->SetFillStyle(0);
0812   leg_sim->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0813   leg_sim->AddEntry(g_pythia_result, "PYTHIA8 truth jet spectrum", "pl");
0814   leg_sim->Draw();
0815   g_pythia_result->Draw("p same");
0816   can_sim->SaveAs("figure/spectrum_calib_simcompare.pdf");
0817 
0818   TCanvas *can_simnlo = new TCanvas("can_simnlo", "", 850, 800);
0819   gStyle->SetPalette(57);
0820   can_simnlo->SetTopMargin(0.03);
0821   can_simnlo->SetLeftMargin(0.15);
0822   can_simnlo->SetBottomMargin(0.12);
0823   can_simnlo->SetRightMargin(0.05);
0824   can_simnlo->SetLogy();
0825   TH2F *frame_simnlo = new TH2F("frame_simnlo", "", 10, 15, 72, 10, 1e-03, 2e+05);
0826   frame_simnlo->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0827   frame_simnlo->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0828   frame_simnlo->GetXaxis()->SetTitleOffset(0.98);
0829   frame_simnlo->GetYaxis()->SetTitleOffset(1.15);
0830   frame_simnlo->GetXaxis()->SetLabelSize(0.045);
0831   frame_simnlo->GetYaxis()->SetLabelSize(0.045);
0832   frame_simnlo->GetXaxis()->CenterTitle();
0833   frame_simnlo->GetYaxis()->CenterTitle();
0834   frame_simnlo->Draw();
0835   g_sphenix_result->Draw("2 same");
0836   g_sphenix_result_point->Draw("P same");
0837   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0838   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0839   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0840   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0841   TLegend *leg_simnlo = new TLegend(0.2, 0.2, 0.55, 0.28);
0842   leg_simnlo->SetBorderSize(0);
0843   leg_simnlo->SetFillStyle(0);
0844   leg_simnlo->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0845   leg_simnlo->AddEntry(g_pythia_result, "PYTHIA8 truth jet spectrum", "pl");
0846   leg_simnlo->AddEntry(tjet, "NLO pQCD", "pl");
0847   leg_simnlo->Draw();
0848   g_pythia_result->Draw("p same");
0849   tjet->Draw("pl same");
0850   tjet05->Draw("l same");
0851   tjet2->Draw("l same");
0852   can_simnlo->SaveAs("figure/spectrum_calib_simnlocompare.pdf");
0853 
0854   TCanvas *can_simratio = new TCanvas("can_simratio", "", 850, 1200);
0855   can_simratio->Divide(1, 2);
0856   gStyle->SetPalette(57);
0857   TPad *pad_1 = (TPad*)can_simratio->cd(1);
0858   pad_1->SetPad(0, 0.33333, 1, 1);
0859   pad_1->SetTopMargin(0.03);
0860   pad_1->SetLeftMargin(0.13);
0861   pad_1->SetBottomMargin(0.035);
0862   pad_1->SetRightMargin(0.08);
0863   pad_1->SetLogy();
0864   TH2F *frame_simratio = new TH2F("frame_simratio", "", 10, 15, 72, 10, 1e-03, 2e+05);
0865   frame_simratio->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0866   frame_simratio->GetYaxis()->SetTitle("d^{2}#sigma/(d#eta dp_{T}) [pb/GeV]");
0867   frame_simratio->GetXaxis()->SetTitleOffset(0.98);
0868   frame_simratio->GetYaxis()->SetTitleOffset(1.15);
0869   frame_simratio->GetXaxis()->SetLabelSize(0.045);
0870   frame_simratio->GetYaxis()->SetLabelSize(0.045);
0871   frame_simratio->GetXaxis()->SetLabelOffset(2);
0872   frame_simratio->GetXaxis()->CenterTitle();
0873   frame_simratio->GetYaxis()->CenterTitle();
0874   frame_simratio->Draw();
0875   g_sphenix_result->Draw("2 same");
0876   g_sphenix_result_point->Draw("P same");
0877   myText(0.4, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0878   myText(0.4, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0879   myText(0.4, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0880   myText(0.4, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
0881   TLegend *leg_simratio = new TLegend(0.2, 0.2, 0.55, 0.28);
0882   leg_simratio->SetBorderSize(0);
0883   leg_simratio->SetFillStyle(0);
0884   leg_simratio->AddEntry(g_sphenix_result, "sPHENIX Run 2024 data", "pf");
0885   leg_simratio->AddEntry(g_pythia_result, "PYTHIA8 truth jet spectrum", "pl");
0886   leg_simratio->Draw();
0887   g_pythia_result->Draw("p same");
0888   TPad *pad_2 = (TPad*)can_simratio->cd(2);
0889   pad_2->SetPad(0, 0, 1, 0.333333);
0890   pad_2->SetTopMargin(0.02);
0891   pad_2->SetLeftMargin(0.13);
0892   pad_2->SetBottomMargin(0.25);
0893   pad_2->SetRightMargin(0.08);
0894   TH2F *frame_ratio = new TH2F("frame_ratio", "", 10, 15, 72, 2, 0.6, 1.7);
0895   frame_ratio->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0896   frame_ratio->GetYaxis()->SetTitle("PYTHIA8/Data");
0897   frame_ratio->GetXaxis()->SetTitleOffset(0.98);
0898   frame_ratio->GetYaxis()->SetTitleOffset(0.5);
0899   frame_ratio->GetXaxis()->SetLabelSize(0.045);
0900   frame_ratio->GetYaxis()->SetLabelSize(0.045);
0901   frame_ratio->GetXaxis()->CenterTitle();
0902   frame_ratio->GetYaxis()->CenterTitle();
0903   frame_ratio->GetYaxis()->SetTitleOffset(frame_simratio->GetYaxis()->GetTitleOffset()*1/2.);
0904   frame_ratio->GetYaxis()->SetLabelOffset(frame_simratio->GetYaxis()->GetLabelOffset()*1/2.);
0905   frame_ratio->GetXaxis()->SetLabelSize(frame_simratio->GetXaxis()->GetLabelSize()*2);
0906   frame_ratio->GetYaxis()->SetLabelSize(frame_simratio->GetYaxis()->GetLabelSize()*2);
0907   frame_ratio->GetXaxis()->SetTitleSize(frame_simratio->GetXaxis()->GetTitleSize()*2);
0908   frame_ratio->GetYaxis()->SetTitleSize(frame_simratio->GetYaxis()->GetTitleSize()*2);
0909   frame_ratio->Draw();
0910   h_pythia_ratio->SetMarkerStyle(20);
0911   h_pythia_ratio->SetMarkerColor(kBlack);
0912   h_pythia_ratio->Draw("P same");
0913   can_simratio->SaveAs("figure/spectrum_calib_simratio.pdf");
0914 
0915   // Draw the uncertainty
0916   TH1D* h_uncertainty_unfold_upper = new TH1D("h_uncertainty_unfold_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_unfold_lower = new TH1D("h_uncertainty_unfold_lower", "", calibnpt, calibptbins);
0917   TH1D* h_uncertainty_cut_upper = new TH1D("h_uncertainty_cut_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_cut_lower = new TH1D("h_uncertainty_cut_lower", "", calibnpt, calibptbins);
0918   TH1D* h_uncertainty_jes_upper = new TH1D("h_uncertainty_jes_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_jes_lower = new TH1D("h_uncertainty_jes_lower", "", calibnpt, calibptbins);
0919   TH1D* h_uncertainty_jer_upper = new TH1D("h_uncertainty_jer_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_jer_lower = new TH1D("h_uncertainty_jer_lower", "", calibnpt, calibptbins);
0920   TH1D* h_uncertainty_jeteff_upper = new TH1D("h_uncertainty_jeteff_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_jeteff_lower = new TH1D("h_uncertainty_jeteff_lower", "", calibnpt, calibptbins);
0921   TH1D* h_uncertainty_mbdeff_upper = new TH1D("h_uncertainty_mbdeff_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_mbdeff_lower = new TH1D("h_uncertainty_mbdeff_lower", "", calibnpt, calibptbins);
0922   TH1D* h_uncertainty_mbdcs_upper = new TH1D("h_uncertainty_mbdcs_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_mbdcs_lower = new TH1D("h_uncertainty_mbdcs_lower", "", calibnpt, calibptbins);
0923   TH1D* h_uncertainty_stat_upper = new TH1D("h_uncertainty_stat_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_stat_lower = new TH1D("h_uncertainty_stat_lower", "", calibnpt, calibptbins);
0924   TH1D* h_uncertainty_total_upper = new TH1D("h_uncertainty_total_upper", "", calibnpt, calibptbins); TH1D* h_uncertainty_total_lower = new TH1D("h_uncertainty_total_lower", "", calibnpt, calibptbins);
0925   for (int ib = 1; ib <= h_uncertainty_cut_upper->GetNbinsX(); ++ib) {
0926     h_uncertainty_unfold_upper->SetBinContent(ib, uncertainty_unfold_upper_record[ib-1]);
0927     h_uncertainty_unfold_lower->SetBinContent(ib, -uncertainty_unfold_lower_record[ib-1]);
0928     h_uncertainty_cut_upper->SetBinContent(ib, uncertainty_cut_upper_record[ib-1]);
0929     h_uncertainty_cut_lower->SetBinContent(ib, -uncertainty_cut_lower_record[ib-1]);
0930     h_uncertainty_jes_upper->SetBinContent(ib, uncertainty_jes_upper_record[ib-1]);
0931     h_uncertainty_jes_lower->SetBinContent(ib, -uncertainty_jes_lower_record[ib-1]);
0932     h_uncertainty_jer_upper->SetBinContent(ib, uncertainty_jer_upper_record[ib-1]);
0933     h_uncertainty_jer_lower->SetBinContent(ib, -uncertainty_jer_lower_record[ib-1]);
0934     h_uncertainty_jeteff_upper->SetBinContent(ib, uncertainty_jeteff_upper_record[ib-1]);
0935     h_uncertainty_jeteff_lower->SetBinContent(ib, -uncertainty_jeteff_lower_record[ib-1]);
0936     h_uncertainty_mbdeff_upper->SetBinContent(ib, uncertainty_mbdeff_upper_record[ib-1]);
0937     h_uncertainty_mbdeff_lower->SetBinContent(ib, -uncertainty_mbdeff_lower_record[ib-1]);
0938     h_uncertainty_mbdcs_upper->SetBinContent(ib, uncertainty_mbdcs_upper_record[ib-1]);
0939     h_uncertainty_mbdcs_lower->SetBinContent(ib, -uncertainty_mbdcs_lower_record[ib-1]);
0940     h_uncertainty_stat_upper->SetBinContent(ib, uncertainty_stat_upper_record[ib-1]);
0941     h_uncertainty_stat_lower->SetBinContent(ib, -uncertainty_stat_lower_record[ib-1]);
0942     h_uncertainty_total_upper->SetBinContent(ib, uncertainty_total_upper_record[ib-1]);
0943     h_uncertainty_total_lower->SetBinContent(ib, -uncertainty_total_lower_record[ib-1]);
0944   }
0945   TCanvas *can_unc = new TCanvas("can_unc", "", 850, 800);
0946   gStyle->SetPalette(57);
0947   can_unc->SetTopMargin(0.03);
0948   can_unc->SetLeftMargin(0.15);
0949   can_unc->SetBottomMargin(0.12);
0950   can_unc->SetRightMargin(0.05);
0951   TH2F *frame_unc = new TH2F("frame_unc", "", 10, 15, 72, 2, -1.5, 1.5);
0952   frame_unc->GetXaxis()->SetTitle("p_{T}^{jet} [GeV]");
0953   frame_unc->GetYaxis()->SetTitle("Relative Uncertainty");
0954   frame_unc->GetXaxis()->SetTitleOffset(0.98);
0955   frame_unc->GetYaxis()->SetTitleOffset(1.17);
0956   frame_unc->GetXaxis()->SetLabelSize(0.045);
0957   frame_unc->GetYaxis()->SetLabelSize(0.045);
0958   frame_unc->GetXaxis()->CenterTitle();
0959   frame_unc->GetYaxis()->CenterTitle();
0960   frame_unc->Draw();
0961   h_uncertainty_unfold_upper->SetLineColor(kYellow+1); h_uncertainty_unfold_lower->SetLineColor(kYellow+1); h_uncertainty_unfold_upper->Draw("hist same"); h_uncertainty_unfold_lower->Draw("hist same");
0962   h_uncertainty_cut_upper->SetLineColor(kAzure+2); h_uncertainty_cut_lower->SetLineColor(kAzure+2); h_uncertainty_cut_upper->Draw("hist same"); h_uncertainty_cut_lower->Draw("hist same");
0963   h_uncertainty_jes_upper->SetLineColor(kGreen+2); h_uncertainty_jes_lower->SetLineColor(kGreen+2); h_uncertainty_jes_upper->Draw("hist same"); h_uncertainty_jes_lower->Draw("hist same");
0964   h_uncertainty_jer_upper->SetLineColor(kMagenta); h_uncertainty_jer_lower->SetLineColor(kMagenta); h_uncertainty_jer_upper->Draw("hist same"); h_uncertainty_jer_lower->Draw("hist same");
0965   h_uncertainty_jeteff_upper->SetLineColor(kOrange+1); h_uncertainty_jeteff_lower->SetLineColor(kOrange+1); h_uncertainty_jeteff_upper->Draw("hist same"); h_uncertainty_jeteff_lower->Draw("hist same");
0966   h_uncertainty_mbdeff_upper->SetLineColor(kRed); h_uncertainty_mbdeff_lower->SetLineColor(kRed); h_uncertainty_mbdeff_upper->Draw("hist same"); h_uncertainty_mbdeff_lower->Draw("hist same");
0967   h_uncertainty_mbdcs_upper->SetLineColor(kCyan); h_uncertainty_mbdcs_lower->SetLineColor(kCyan); h_uncertainty_mbdcs_upper->Draw("hist same"); h_uncertainty_mbdcs_lower->Draw("hist same");
0968   h_uncertainty_stat_upper->SetLineColor(kViolet-1); h_uncertainty_stat_lower->SetLineColor(kViolet-1); h_uncertainty_stat_upper->Draw("hist same"); h_uncertainty_stat_lower->Draw("hist same");
0969   h_uncertainty_total_upper->SetLineColor(kBlack); h_uncertainty_total_lower->SetLineColor(kBlack);
0970   h_uncertainty_total_upper->SetMarkerStyle(20); h_uncertainty_total_lower->SetMarkerStyle(20); 
0971   h_uncertainty_total_upper->Draw("p same"); h_uncertainty_total_lower->Draw("p same");
0972   myText(0.17, 0.92, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
0973   myText(0.17, 0.87, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
0974   myText(0.17, 0.82, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
0975   myText(0.17, 0.77, 1, "|#eta^{jet}| < 0.7", 0.04);
0976   TLegend *leg_unc1 = new TLegend(0.155, 0.15, 0.5, 0.34);
0977   leg_unc1->SetBorderSize(0);
0978   leg_unc1->SetFillStyle(0);
0979   leg_unc1->AddEntry(h_uncertainty_unfold_upper, "Unfolding", "l");
0980   leg_unc1->AddEntry(h_uncertainty_cut_upper, "Background Cut", "l");
0981   leg_unc1->AddEntry(h_uncertainty_jes_upper, "Jet Energy Scale", "l");
0982   leg_unc1->AddEntry(h_uncertainty_jer_upper, "Jet Energy Resolution", "l");
0983   leg_unc1->AddEntry(h_uncertainty_jeteff_upper, "Jet Efficiency", "l");
0984   leg_unc1->Draw();
0985   TLegend *leg_unc2 = new TLegend(0.505, 0.15, 0.85, 0.30);
0986   leg_unc2->SetBorderSize(0);
0987   leg_unc2->SetFillStyle(0);
0988   leg_unc2->AddEntry(h_uncertainty_mbdeff_upper, "MBD Efficiency", "l");
0989   leg_unc2->AddEntry(h_uncertainty_mbdcs_upper, "MBD Cross Section", "l");
0990   leg_unc2->AddEntry(h_uncertainty_stat_upper, "Statistical", "l");
0991   leg_unc2->AddEntry(h_uncertainty_total_upper, "Total Systematic", "p");
0992   leg_unc2->Draw();
0993   TLine *line = new TLine(15, 0, 72, 0);
0994   line->SetLineStyle(10);
0995   line->Draw("same");
0996   can_unc->SaveAs("figure/spectrum_uncertainty.pdf");
0997 }
0998 
0999 void draw_1D_multiple_plot_uncertainty(std::vector<TH1F*> h_input, float x_min, float x_max, float ratio_min, float ratio_max, std::string xtitle, std::string ytitle, std::vector<std::string> legend, float xstart_legend, float ystart_legend, std::vector<int> color, std::vector<int> markerstyle, std::string output_name) {
1000   if (h_input.size() <= 1) {
1001     std::cout << "Error: no enough input histograms for draw_1D_multiple_plot" << std::endl;
1002     return;
1003   }
1004   std::vector<TH1F*> h;
1005   for (int i = 0; i < h_input.size(); ++i) {
1006     h.push_back((TH1F*)h_input[i]->Clone(Form("h_%d", i)));
1007   }
1008   TCanvas *can = new TCanvas("can", "", 800, 889);
1009   can->Divide(1, 2);
1010   gStyle->SetPalette(57);
1011   TPad *pad_1 = (TPad*)can->cd(1);
1012   pad_1->SetPad(0, 0.4, 1, 1);
1013   pad_1->SetTopMargin(0.03);
1014   pad_1->SetLeftMargin(0.13);
1015   pad_1->SetBottomMargin(0.035);
1016   pad_1->SetRightMargin(0.08);
1017   pad_1->SetLogy();
1018   for (int ih = 0; ih < h.size(); ++ih) {
1019     h.at(ih)->SetMarkerStyle(markerstyle.at(ih));
1020     h.at(ih)->SetMarkerColor(color.at(ih));
1021     h.at(ih)->SetLineColor(color.at(ih));
1022     for (int ib = 1; ib <= h.at(1)->GetNbinsX(); ++ib) {
1023       h.at(ih)->SetBinContent(ib, h.at(ih)->GetBinContent(ib)/ (float)(/*2*TMath::Pi() * */luminosity * delta_eta * h.at(ih)->GetBinWidth(ib)));
1024       h.at(ih)->SetBinError(ib, h.at(ih)->GetBinError(ib)/ (float)(/*2*TMath::Pi() * */luminosity * delta_eta * h.at(ih)->GetBinWidth(ib)));
1025     }
1026   }
1027   h.at(0)->GetXaxis()->SetTitle(xtitle.c_str());
1028   h.at(0)->GetYaxis()->SetTitle(ytitle.c_str());
1029   h.at(0)->GetXaxis()->SetRangeUser(x_min, x_max);
1030   h.at(0)->GetXaxis()->SetTitleOffset(0.98);
1031   h.at(0)->GetYaxis()->SetTitleOffset(1.15);
1032   h.at(0)->GetXaxis()->SetLabelSize(0.045);
1033   h.at(0)->GetYaxis()->SetLabelSize(0.045);
1034   h.at(0)->GetXaxis()->SetLabelOffset(2);
1035   h.at(0)->GetXaxis()->CenterTitle();
1036   h.at(0)->GetYaxis()->CenterTitle();
1037   h.at(0)->Draw();
1038   for (int i = 0; i < h.size(); ++i) {
1039     if (i == 0) continue;
1040     h.at(i)->Draw("same");
1041   }
1042 
1043   myText(0.47, 0.9, 1, "#bf{#it{sPHENIX}} Internal", 0.04);
1044   myText(0.47, 0.85, 1, "p+p#sqrt{s} = 200 GeV,#scale[0.5]{#int}#kern[-0.1]{#it{L}}dt = 14.93pb^{-1}", 0.04);
1045   myText(0.47, 0.8, 1, "anti-k_{T}#it{R} = 0.4", 0.04);
1046   myText(0.47, 0.75, 1, "|#eta^{jet}| < 0.7", 0.04);
1047   for (int i = 0; i < legend.size(); i++) {
1048     myMarkerLineText(xstart_legend, ystart_legend-i*0.05, 1, color.at(i), markerstyle.at(i), color.at(i), 1, legend.at(i).c_str(), 0.04, true);
1049   }
1050 
1051   TPad *pad_2 = (TPad*)can->cd(2);
1052   pad_2->SetPad(0, 0, 1, 0.4);
1053   pad_2->SetTopMargin(0.02);
1054   pad_2->SetLeftMargin(0.13);
1055   pad_2->SetBottomMargin(0.25);
1056   pad_2->SetRightMargin(0.08);
1057   std::vector<TH1F*> h_ratio;
1058   for (int i = 0; i < h.size()-1; ++i) {
1059     TH1F *h_temp = (TH1F*)h.at(i+1)->Clone(Form("h_ratio_%d", i));
1060     h_temp->Divide(h.at(0));
1061     h_ratio.push_back(h_temp);
1062     h_ratio.at(i)->SetMarkerStyle(markerstyle.at(i+1));
1063     h_ratio.at(i)->SetMarkerColor(color.at(i+1));
1064     h_ratio.at(i)->SetLineColor(color.at(i+1));
1065   }
1066   h_ratio.at(0)->GetXaxis()->SetTitle(xtitle.c_str());
1067   h_ratio.at(0)->GetYaxis()->SetTitle("Ratio");
1068   h_ratio.at(0)->GetXaxis()->SetRangeUser(x_min, x_max);
1069   h_ratio.at(0)->GetYaxis()->SetRangeUser(ratio_min, ratio_max);
1070   h_ratio.at(0)->GetXaxis()->CenterTitle();
1071   h_ratio.at(0)->GetYaxis()->CenterTitle();
1072   h_ratio.at(0)->GetYaxis()->SetTitleOffset(h.at(0)->GetYaxis()->GetTitleOffset()*4/6.);
1073   h_ratio.at(0)->GetYaxis()->SetLabelOffset(h.at(0)->GetYaxis()->GetLabelOffset()*4/6.);
1074   h_ratio.at(0)->GetXaxis()->SetLabelSize(h.at(0)->GetXaxis()->GetLabelSize()*6/4.);
1075   h_ratio.at(0)->GetYaxis()->SetLabelSize(h.at(0)->GetYaxis()->GetLabelSize()*6/4.);
1076   h_ratio.at(0)->GetXaxis()->SetTitleSize(h.at(0)->GetXaxis()->GetTitleSize()*6/4.);
1077   h_ratio.at(0)->GetYaxis()->SetTitleSize(h.at(0)->GetYaxis()->GetTitleSize()*6/4.);
1078   h_ratio.at(0)->Draw();
1079   for (int i = 0; i < h_ratio.size(); ++i) {
1080     if (i == 0) continue;
1081     h_ratio.at(i)->Draw("same");
1082   }
1083 
1084   TLine *line;
1085   line = new TLine(x_min, 1, x_max, 1);
1086   line->SetLineColor(kBlack);
1087   line->SetLineStyle(3);
1088   line->Draw("same");
1089 
1090   can->SaveAs(output_name.c_str());
1091   delete can;
1092 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!