PPG09_JetSpectrum_Analysis/offline/do_unfolding_iter.C

0001 #include <TChain.h>
0002 #include <TFile.h>
0003 #include <TH1F.h>
0004 #include <TH2F.h>
0005 #include <TF1.h>
0006 #include <TMath.h>
0007 #include <iostream>
0008 #include <vector>
0009 #include <string>
0010 #include <cmath>
0011 #include "RooUnfold.h"
0012 #include "RooUnfoldResponse.h"
0013 #include "RooUnfoldBayes.h"
0014
0015 void do_unfolding_iteration_check(TH2D* h_respmatrix, TH1D* h_spectrum, string figure_name) {
0016   TH1D* h_meas = (TH1D*)h_respmatrix->ProjectionX("h_meas");
0017   TH1D* h_truth = (TH1D*)h_respmatrix->ProjectionY("h_truth");
0018   const int niter = 20;
0019   TH1D* h_unfolded[niter];
0020   for (int i = 0; i < niter; ++i) {
0021     RooUnfoldResponse* response = new RooUnfoldResponse(h_meas, h_truth, h_respmatrix, "response", "");
0022     RooUnfoldBayes unfold(response, h_spectrum, i + 1);
0023     h_unfolded[i] = (TH1D*)unfold.Hunfold(RooUnfold::kErrors);
0024     h_unfolded[i]->SetName(Form("h_unfolded_iter_%d", i + 1));
0025   }
0026
0027   TGraph* g_var_iter = new TGraph(niter-1);
0028   TGraph* g_var_error = new TGraph(niter-1);
0029   TGraph* g_var_total = new TGraph(niter-1);
0030   for (int i = 1; i < niter; ++i) {
0031     double var_iter = 0;
0032     double var_error = 0;
0033     double var_total = 0;
0034     for (int j = 3; j <= h_unfolded[i]->GetNbinsX()-1; ++j) {
0035       double diff = h_unfolded[i]->GetBinContent(j) - h_unfolded[i-1]->GetBinContent(j);
0036       double cont = (h_unfolded[i]->GetBinContent(j) + h_unfolded[i-1]->GetBinContent(j)) / 2.0;
0037       var_iter += diff * diff / (double)(cont * cont);
0038       var_error += h_unfolded[i]->GetBinError(j) * h_unfolded[i]->GetBinError(j) / (double)(h_unfolded[i]->GetBinContent(j) * h_unfolded[i]->GetBinContent(j));
0039     }
0040     var_iter = std::sqrt(var_iter / (double)(h_unfolded[i]->GetNbinsX()-3));
0041     var_error = std::sqrt(var_error / (double)(h_unfolded[i]->GetNbinsX()-3));
0042     var_total = std::sqrt(var_iter * var_iter + var_error * var_error);
0043
0044     g_var_iter->SetPoint(i-1, i+1, var_iter);
0045     g_var_error->SetPoint(i-1, i+1, var_error);
0046     g_var_total->SetPoint(i-1, i+1, var_total);
0047   }
0048   TCanvas* can = new TCanvas("can", "Unfolding Iteration Check", 800, 600);
0049   g_var_iter->SetMarkerStyle(20);
0050   g_var_iter->SetMarkerColor(kBlue);
0051   g_var_iter->SetLineColor(kBlue);
0052   g_var_error->SetMarkerStyle(20);
0053   g_var_error->SetMarkerColor(kRed);
0054   g_var_error->SetLineColor(kRed);
0055   g_var_total->SetMarkerStyle(20);
0056   g_var_total->SetMarkerColor(kBlack);
0057   g_var_total->SetLineColor(kBlack);
0058   g_var_total->SetMinimum(0);
0059   g_var_total->Draw("AP");
0060   g_var_iter->Draw("P same");
0061   g_var_error->Draw("P same");
0062   g_var_iter->GetXaxis()->SetTitle("Iteration");
0063   can->SaveAs((figure_name + ".png").c_str());
0064 }
0065
0066 void get_unfolded_spectrum(TH2D* h_respmatrix, TH1D* h_spectrum, TH1D*& h_unfolded, int niter, std::string hist_name, TH1D* h_eff) {
0067   TH1D* h_meas = (TH1D*)h_respmatrix->ProjectionX("h_meas");
0068   TH1D* h_truth = (TH1D*)h_respmatrix->ProjectionY("h_truth");
0069   RooUnfoldResponse* response = new RooUnfoldResponse(h_meas, h_truth, h_respmatrix, "response", "");
0070   RooUnfoldBayes unfold(response, h_spectrum, niter);
0071   h_unfolded = (TH1D*)unfold.Hunfold(RooUnfold::kErrors);
0072   h_unfolded->SetName(hist_name.c_str());
0073   h_unfolded->Divide(h_eff);
0074 }
0075
0076 ////////////////////////////////////////// Main Function //////////////////////////////////////////
0077 void do_unfolding_iter(int radius_index = 4) {
0078   //********** General Set up **********//
0079   const float PI = TMath::Pi();
0080   float jet_radius = 0.1 * radius_index;
0081
0082   //********** Files **********//
0083   TFile *f_out = new TFile(Form("output_unfolded_r0%d.root", radius_index), "RECREATE");
0084   TFile *f_data = new TFile(Form("output_data_puritycorr_r0%d.root", radius_index), "READ");
0085   TFile *f_in_rm = new TFile(Form("output_sim_r0%d.root", radius_index), "READ");
0086   TFile *f_efficiency = new TFile(Form("output_purityefficiency_r0%d.root", radius_index), "READ");
0087
0088   // Get data spectrum
0089   TH1D* h_calibjet_pt_puritycorr_nominal = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_nominal");
0090   TH1D* h_calibjet_pt_puritycorr_jesup = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jesup");
0091   TH1D* h_calibjet_pt_puritycorr_jesdown = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jesdown");
0092   TH1D* h_calibjet_pt_puritycorr_jerup = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jerup");
0093   TH1D* h_calibjet_pt_puritycorr_jerdown = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jerdown");
0094   TH1D* h_calibjet_pt_puritycorr_jettrigup = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jettrigup");
0095   TH1D* h_calibjet_pt_puritycorr_jettrigdown = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jettrigdown");
0096   TH1D* h_calibjet_pt_puritycorr_jettimingup = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jettimingup");
0097   TH1D* h_calibjet_pt_puritycorr_jettimingdown = (TH1D*)f_data->Get("h_calibjet_pt_puritycorr_jettimingdown");
0098
0099   // Response matrix
0100   TH2D* h_respmatrix_nominal = (TH2D*)f_in_rm->Get("h_respmatrix_nominal");
0101   TH2D* h_respmatrix_jesup = (TH2D*)f_in_rm->Get("h_respmatrix_jesup");
0102   TH2D* h_respmatrix_jesdown = (TH2D*)f_in_rm->Get("h_respmatrix_jesdown");
0103   TH2D* h_respmatrix_jerup = (TH2D*)f_in_rm->Get("h_respmatrix_jerup");
0104   TH2D* h_respmatrix_jerdown = (TH2D*)f_in_rm->Get("h_respmatrix_jerdown");
0105   TH2D* h_respmatrix_jettrigup = (TH2D*)h_respmatrix_nominal->Clone("h_respmatrix_jettrigup");
0106   TH2D* h_respmatrix_jettrigdown = (TH2D*)h_respmatrix_nominal->Clone("h_respmatrix_jettrigdown");
0107   TH2D* h_respmatrix_jettimingup = (TH2D*)h_respmatrix_nominal->Clone("h_respmatrix_jettimingup");
0108   TH2D* h_respmatrix_jettimingdown = (TH2D*)h_respmatrix_nominal->Clone("h_respmatrix_jettimingdown");
0109
0110   // Get efficiency
0111   TH1D* h_efficiency_nominal = (TH1D*)f_efficiency->Get("h_efficiency_nominal");
0112   TH1D* h_efficiency_jesup = (TH1D*)f_efficiency->Get("h_efficiency_jesup");
0113   TH1D* h_efficiency_jesdown = (TH1D*)f_efficiency->Get("h_efficiency_jesdown");
0114   TH1D* h_efficiency_jerup = (TH1D*)f_efficiency->Get("h_efficiency_jerup");
0115   TH1D* h_efficiency_jerdown = (TH1D*)f_efficiency->Get("h_efficiency_jerdown");
0116   TH1D* h_efficiency_jettrigup = (TH1D*)f_efficiency->Get("h_efficiency_jettrigup");
0117   TH1D* h_efficiency_jettrigdown = (TH1D*)f_efficiency->Get("h_efficiency_jettrigdown");
0118   TH1D* h_efficiency_jettimingup = (TH1D*)f_efficiency->Get("h_efficiency_jettimingup");
0119   TH1D* h_efficiency_jettimingdown = (TH1D*)f_efficiency->Get("h_efficiency_jettimingdown");
0120
0121   //********** Response Matrix **********//
0122   do_unfolding_iteration_check(h_respmatrix_nominal, h_calibjet_pt_puritycorr_nominal, Form("unfolding_iteration_check_nominal_r0%d", radius_index));
0123   int iteration = 4;
0124
0125   TH1D* h_unfold_nominal; get_unfolded_spectrum(h_respmatrix_nominal, h_calibjet_pt_puritycorr_nominal, h_unfold_nominal, iteration, "h_unfold_nominal", h_efficiency_nominal);
0126   TH1D* h_unfold_jesup; get_unfolded_spectrum(h_respmatrix_jesup, h_calibjet_pt_puritycorr_jesup, h_unfold_jesup, iteration, "h_unfold_jesup", h_efficiency_jesup);
0127   TH1D* h_unfold_jesdown; get_unfolded_spectrum(h_respmatrix_jesdown, h_calibjet_pt_puritycorr_jesdown, h_unfold_jesdown, iteration, "h_unfold_jesdown", h_efficiency_jesdown);
0128   TH1D* h_unfold_jerup; get_unfolded_spectrum(h_respmatrix_jerup, h_calibjet_pt_puritycorr_jerup, h_unfold_jerup, iteration, "h_unfold_jerup", h_efficiency_jerup);
0129   TH1D* h_unfold_jerdown; get_unfolded_spectrum(h_respmatrix_jerdown, h_calibjet_pt_puritycorr_jerdown, h_unfold_jerdown, iteration, "h_unfold_jerdown", h_efficiency_jerdown);
0130   TH1D* h_unfold_jettrigup; get_unfolded_spectrum(h_respmatrix_jettrigup, h_calibjet_pt_puritycorr_jettrigup, h_unfold_jettrigup, iteration, "h_unfold_jettrigup", h_efficiency_jettrigup);
0131   TH1D* h_unfold_jettrigdown; get_unfolded_spectrum(h_respmatrix_jettrigdown, h_calibjet_pt_puritycorr_jettrigdown, h_unfold_jettrigdown, iteration, "h_unfold_jettrigdown", h_efficiency_jettrigdown);
0132   TH1D* h_unfold_jettimingup; get_unfolded_spectrum(h_respmatrix_jettimingup, h_calibjet_pt_puritycorr_jettimingup, h_unfold_jettimingup, iteration, "h_unfold_jettimingup", h_efficiency_jettimingup);
0133   TH1D* h_unfold_jettimingdown; get_unfolded_spectrum(h_respmatrix_jettimingdown, h_calibjet_pt_puritycorr_jettimingdown, h_unfold_jettimingdown, iteration, "h_unfold_jettimingdown", h_efficiency_jettimingdown);
0134   TH1D* h_unfold_unfolditerup; get_unfolded_spectrum(h_respmatrix_nominal, h_calibjet_pt_puritycorr_nominal, h_unfold_unfolditerup, iteration+1, "h_unfold_unfolditerup", h_efficiency_nominal);
0135   TH1D* h_unfold_unfolditerdown; get_unfolded_spectrum(h_respmatrix_nominal, h_calibjet_pt_puritycorr_nominal, h_unfold_unfolditerdown, iteration-1, "h_unfold_unfolditerdown", h_efficiency_nominal);
0136
0137   //********** Writing **********//
0138   std::cout << "Writing histograms..." << std::endl;
0139   f_out->cd();
0140
0141   h_unfold_nominal->Write();
0142   h_unfold_jesup->Write();
0143   h_unfold_jesdown->Write();
0144   h_unfold_jerup->Write();
0145   h_unfold_jerdown->Write();
0146   h_unfold_jettrigup->Write();
0147   h_unfold_jettrigdown->Write();
0148   h_unfold_jettimingup->Write();
0149   h_unfold_jettimingdown->Write();
0150   h_unfold_unfolditerup->Write();
0151   h_unfold_unfolditerdown->Write();
0152
0153   f_out->Close();
0154   std::cout << "All done!" << std::endl;
0155 }