PPG09_JetSpectrum_Analysis/offline/get_purityefficiency.C

0001 #include <iostream>
0002 #include <TFile.h>
0003 #include <TH1D.h>
0004 #include "/sphenix/user/hanpuj/CaloDataAna24_skimmed/src/draw_template.C"
0005 #include "unfold_Def.h"
0006
0007 void get_purityefficiency_hist(TH1D *&h_purity, TH1D *&h_efficiency, string surfix, TFile *f_in, TFile *f_out, int jet_radius_index) {
0008   // Get histograms from input file
0009   TH1D* h_truth = (TH1D*)f_in->Get(("h_truth"+surfix).c_str());
0010   TH1D* h_measure = (TH1D*)f_in->Get(("h_measure"+surfix).c_str());
0011   TH2D* h_respmatrix = (TH2D*)f_in->Get(("h_respmatrix"+surfix).c_str());
0012   TH1D* h_fake = (TH1D*)f_in->Get(("h_fake"+surfix).c_str());
0013   TH1D* h_miss = (TH1D*)f_in->Get(("h_miss"+surfix).c_str());
0014   if (!h_truth || !h_measure || !h_respmatrix || !h_fake || !h_miss) {
0015     std::cerr << "Error: One or more histograms not found in input file for surfix: " << surfix << std::endl;
0016     return;
0017   }
0018   // Get purity and efficiency histograms
0019   h_purity = (TH1D*)h_respmatrix->ProjectionX(("h_purity"+surfix).c_str());
0020   h_purity->Divide(h_purity, h_measure, 1, 1, "B");
0021   h_efficiency = (TH1D*)h_respmatrix->ProjectionY(("h_efficiency"+surfix).c_str());
0022   h_efficiency->Divide(h_efficiency, h_truth, 1, 1, "B");
0023   // Check purity and efficiency
0024   TH1D* h_purity_check = (TH1D*)h_respmatrix->ProjectionX(("h_purity"+surfix+"_check").c_str());
0025   h_purity_check->Add(h_fake);
0026   TH1D* h_efficiency_check = (TH1D*)h_respmatrix->ProjectionY(("h_efficiency"+surfix+"_check").c_str());
0027   h_efficiency_check->Add(h_miss);
0028   for (int i = 1; i <= h_purity_check->GetNbinsX(); i++) {
0029     float purity_check = h_purity_check->GetBinContent(i) / (double)h_measure->GetBinContent(i);
0030     if (purity_check != 1.0) std::cout << "Case: " << surfix << " in purity check failed at bin " << i << " with ratio = " << purity_check << " with purity content = " << h_purity_check->GetBinContent(i) << " and measurement content = " << h_measure->GetBinContent(i) << std::endl;
0031   }
0032   for (int i = 1; i <= h_efficiency_check->GetNbinsX(); i++) {
0033     float efficiency_check = h_efficiency_check->GetBinContent(i) / (double)h_truth->GetBinContent(i);
0034     if (efficiency_check != 1.0) std::cout << "Case: " << surfix << " in efficiency check failed at bin " << i << " with ratio = " << efficiency_check << " with purity content = " << h_efficiency_check->GetBinContent(i) << " and measurement content = " << h_truth->GetBinContent(i) << std::endl;
0035   }
0036   // Write histograms to output file
0037   f_out->cd();
0038   h_purity->Write();
0039   h_efficiency->Write();
0040
0041   TH1D* h_purity_proj = (TH1D*)h_respmatrix->ProjectionX("h_purity_proj");
0042   TH1D* h_efficiency_proj = (TH1D*)h_respmatrix->ProjectionY("h_efficiency_proj");
0043   std::vector<TH1F*> h_input;
0044   std::vector<int> color;
0045   std::vector<int> markerstyle;
0046   std::vector<std::string> text;
0047   std::vector<std::string> legend;
0048
0049   h_input.push_back((TH1F*)h_measure);
0050   h_input.push_back((TH1F*)h_purity_proj);
0051   color.push_back(kBlack);
0052   color.push_back(kRed);
0053   markerstyle.push_back(20);
0054   markerstyle.push_back(20);
0055   text.push_back("#bf{#it{sPHENIX}} Simulation");
0056   text.push_back("PYTHIA8 p+p#sqrt{s} = 200 GeV");
0057   text.push_back(("anti-k_{t} #kern[-0.5]{#it{R}} = 0." + std::to_string(jet_radius_index)).c_str());
0058   text.push_back(("|#eta^{jet}| < 0." + std::to_string((int)(11-jet_radius_index))).c_str());
0059   text.push_back("No z-vertex cut");
0060   legend.push_back("Reco jet spectrum");
0061   legend.push_back("Matched reco jet spectrum");
0062   draw_1D_multiple_plot_ratio(h_input, color, markerstyle,
0063                               false, 10, true,
0064                               true, 15, 75, false,
0065                               false, 0, 0.5, true,
0066                               false, 0.75, 1.1,
0067                               true, "p_{T}^{reco jet} [GeV]", "Arbitrary Unit", "Purity", 1,
0068                               true, text, 0.45, 0.9, 0.05,
0069                               true, legend, 0.25, 0.2, 0.05,
0070                               Form("figure_purityefficiency/purity%s_r0%d.png", surfix.c_str(), jet_radius_index));
0071   h_input.clear();
0072   legend.clear();
0073
0074   h_input.push_back((TH1F*)h_truth);
0075   h_input.push_back((TH1F*)h_efficiency_proj);
0076   legend.push_back("Truth jet spectrum");
0077   legend.push_back("Matched truth jet spectrum");
0078   draw_1D_multiple_plot_ratio(h_input, color, markerstyle,
0079                               false, 10, true,
0080                               true, 15, 75, false,
0081                               false, 1e-9, 1e-1, true,
0082                               false, 0., 1.2,
0083                               true, "p_{T}^{truth jet} [GeV]", "Arbitrary Unit", "Efficiency", 1,
0084                               true, text, 0.45, 0.9, 0.05,
0085                               true, legend, 0.25, 0.2, 0.05,
0086                               Form("figure_purityefficiency/efficiency%s_r0%d.png", surfix.c_str(), jet_radius_index));
0087   h_input.clear();
0088   color.clear();
0089   markerstyle.clear();
0090   text.clear();
0091   legend.clear();
0092 }
0093
0094 void get_purityefficiency_clonehist(TH1D *&h_purity, TH1D *&h_efficiency, string surfix, TH1D *h_purity_forclone, TH1D *h_efficiency_forclone, TFile *f_out) {
0095   h_purity = dynamic_cast<TH1D*>(h_purity_forclone->Clone(("h_purity"+surfix).c_str()));
0096   h_efficiency = dynamic_cast<TH1D*>((TH1D*)h_efficiency_forclone->Clone(("h_efficiency"+surfix).c_str()));
0097   f_out->cd();
0098   h_purity->Write();
0099   h_efficiency->Write();
0100 }
0101
0102 void do_puritycorr(TH1D* h_result, TH1D* h_purity, string result_surfix, string data_surfix, TFile* f_data, TFile* f_dataout) {
0103   TH1D* h_calibjet_pt = (TH1D*)f_data->Get(("h_calibjet_pt_scaled" + data_surfix).c_str());
0104   h_result = (TH1D*)h_calibjet_pt->Clone(("h_calibjet_pt_puritycorr" + result_surfix).c_str());
0105   h_result->Multiply(h_purity);
0106   f_dataout->cd();
0107   h_result->Write();
0108 }
0109
0110
0111 void get_purityefficiency(int radius_index = 4) {
0112   // Read Files
0113   TFile *f_in = new TFile(Form("output_sim_r0%d.root", radius_index), "READ");
0114   if (!f_in) {
0115     std::cout << "Error: cannot open output_sim_r0" << radius_index << ".root" << std::endl;
0116     return;
0117   }
0118   TFile *f_out = new TFile(Form("output_purityefficiency_r0%d.root", radius_index), "RECREATE");
0119   TFile* f_data = new TFile(Form("output_data_r0%d.root", radius_index), "READ");
0120   TFile* f_dataout = new TFile(Form("output_data_puritycorr_r0%d.root", radius_index), "RECREATE");
0121
0122   // Get purity and efficiency histograms
0123   TH1D* h_purity_nominal, *h_efficiency_nominal; get_purityefficiency_hist(h_purity_nominal, h_efficiency_nominal, "_nominal", f_in, f_out, radius_index);
0124   TH1D* h_purity_jesup, *h_efficiency_jesup; get_purityefficiency_hist(h_purity_jesup, h_efficiency_jesup, "_jesup", f_in, f_out, radius_index);
0125   TH1D* h_purity_jesdown, *h_efficiency_jesdown; get_purityefficiency_hist(h_purity_jesdown, h_efficiency_jesdown, "_jesdown", f_in, f_out, radius_index);
0126   TH1D* h_purity_jerup, *h_efficiency_jerup; get_purityefficiency_hist(h_purity_jerup, h_efficiency_jerup, "_jerup", f_in, f_out, radius_index);
0127   TH1D* h_purity_jerdown, *h_efficiency_jerdown; get_purityefficiency_hist(h_purity_jerdown, h_efficiency_jerdown, "_jerdown", f_in, f_out, radius_index);
0128   TH1D* h_purity_jettrigup, *h_efficiency_jettrigup; get_purityefficiency_clonehist(h_purity_jettrigup, h_efficiency_jettrigup, "_jettrigup", h_purity_nominal, h_efficiency_nominal, f_out);
0129   TH1D* h_purity_jettrigdown, *h_efficiency_jettrigdown; get_purityefficiency_clonehist(h_purity_jettrigdown, h_efficiency_jettrigdown, "_jettrigdown", h_purity_nominal, h_efficiency_nominal, f_out);
0130   TH1D* h_purity_jettimingup, *h_efficiency_jettimingup; get_purityefficiency_clonehist(h_purity_jettimingup, h_efficiency_jettimingup, "_jettimingup", h_purity_nominal, h_efficiency_nominal, f_out);
0131   TH1D* h_purity_jettimingdown, *h_efficiency_jettimingdown; get_purityefficiency_clonehist(h_purity_jettimingdown, h_efficiency_jettimingdown, "_jettimingdown", h_purity_nominal, h_efficiency_nominal, f_out);
0132
0133   // Do purity correction
0134   TH1D* h_calibjet_pt_puritycorr_nominal; do_puritycorr(h_calibjet_pt_puritycorr_nominal, h_purity_nominal, "_nominal", "_nominal", f_data, f_dataout);
0135   TH1D* h_calibjet_pt_puritycorr_jesup; do_puritycorr(h_calibjet_pt_puritycorr_jesup, h_purity_jesup, "_jesup", "_nominal", f_data, f_dataout);
0136   TH1D* h_calibjet_pt_puritycorr_jesdown; do_puritycorr(h_calibjet_pt_puritycorr_jesdown, h_purity_jesdown, "_jesdown", "_nominal", f_data, f_dataout);
0137   TH1D* h_calibjet_pt_puritycorr_jerup; do_puritycorr(h_calibjet_pt_puritycorr_jerup, h_purity_jerup, "_jerup", "_nominal", f_data, f_dataout);
0138   TH1D* h_calibjet_pt_puritycorr_jerdown; do_puritycorr(h_calibjet_pt_puritycorr_jerdown, h_purity_jerdown, "_jerdown", "_nominal", f_data, f_dataout);
0139   TH1D* h_calibjet_pt_puritycorr_jettrigup; do_puritycorr(h_calibjet_pt_puritycorr_jettrigup, h_purity_jettrigup, "_jettrigup", "_jettrigup", f_data, f_dataout);
0140   TH1D* h_calibjet_pt_puritycorr_jettrigdown; do_puritycorr(h_calibjet_pt_puritycorr_jettrigdown, h_purity_jettrigdown, "_jettrigdown", "_jettrigdown", f_data, f_dataout);
0141   TH1D* h_calibjet_pt_puritycorr_jettimingup; do_puritycorr(h_calibjet_pt_puritycorr_jettimingup, h_purity_jettimingup, "_jettimingup", "_jettimingup", f_data, f_dataout);
0142   TH1D* h_calibjet_pt_puritycorr_jettimingdown; do_puritycorr(h_calibjet_pt_puritycorr_jettimingdown, h_purity_jettimingdown, "_jettimingdown", "_jettimingdown", f_data, f_dataout);
0143 }