Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
extractLevelingEnergyLinearity.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file extractLevelingEnergyLinearity.C
1 #include <string>
2 #include <TString.h>
3 #include <TFile.h>
4 #include <TH2F.h>
5 #include <TH1F.h>
6 #include <TF1.h>
7 #include <TCanvas.h>
8 #include <TGraphAsymmErrors.h>
9 #include <TProfile.h>
10 #include <TLine.h>
11 
12 float ErrorAdd(float x, float y)
13 {
14  return sqrt(x*x+y*y);
15 }
16 
17 float ErrTimes(float x, float y, float dx, float dy)
18 {
19  return x*y*ErrorAdd(dx/x,dy/y);
20 }
21 
22 float ErrDiv(float x, float y, float dx, float dy)
23 {
24  return x/y*ErrorAdd(dx/x,dy/y);
25 }
26 
27 void extractLevelingEnergyLinearity()
28 {
29  string input_5GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_1087.root";
30  TFile *File_5GeV = TFile::Open(input_5GeV.c_str());
31  // TH2F *h_mAsymmEnergy_pion_5GeV = (TH2F*)File_5GeV->Get("h_mAsymmEnergy_electron");
32  TH2F *h_mAsymmEnergy_pion_5GeV = (TH2F*)File_5GeV->Get("h_mAsymmEnergy_pion_leveling");
33  TH1F *h_mEnergy_pion_5GeV = (TH1F*)h_mAsymmEnergy_pion_5GeV->ProjectionY()->Clone();
34 
35  string input_8GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0422.root";
36  TFile *File_8GeV = TFile::Open(input_8GeV.c_str());
37  // TH2F *h_mAsymmEnergy_pion_8GeV = (TH2F*)File_8GeV->Get("h_mAsymmEnergy_electron");
38  TH2F *h_mAsymmEnergy_pion_8GeV = (TH2F*)File_8GeV->Get("h_mAsymmEnergy_pion_leveling");
39  TH1F *h_mEnergy_pion_8GeV = (TH1F*)h_mAsymmEnergy_pion_8GeV->ProjectionY()->Clone();
40 
41  string input_12GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0571.root";
42  TFile *File_12GeV = TFile::Open(input_12GeV.c_str());
43  // TH2F *h_mAsymmEnergy_pion_12GeV = (TH2F*)File_12GeV->Get("h_mAsymmEnergy_electron");
44  TH2F *h_mAsymmEnergy_pion_12GeV = (TH2F*)File_12GeV->Get("h_mAsymmEnergy_pion_leveling");
45  TH1F *h_mEnergy_pion_12GeV = (TH1F*)h_mAsymmEnergy_pion_12GeV->ProjectionY()->Clone();
46 
47  string input_60GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0821.root";
48  TFile *File_60GeV = TFile::Open(input_60GeV.c_str());
49  // TH2F *h_mAsymmEnergy_pion_60GeV = (TH2F*)File_60GeV->Get("h_mAsymmEnergy_electron");
50  TH2F *h_mAsymmEnergy_pion_60GeV = (TH2F*)File_60GeV->Get("h_mAsymmEnergy_pion");
51  TH1F *h_mEnergy_pion_60GeV = (TH1F*)h_mAsymmEnergy_pion_60GeV->ProjectionY()->Clone();
52 
53  TCanvas *c_AsymmEnergy = new TCanvas("c_AsymmEnergy","c_AsymmEnergy",1600,400);
54  c_AsymmEnergy->Divide(4,1);
55  for(int i_pad = 0; i_pad < 4; ++i_pad)
56  {
57  c_AsymmEnergy->cd(i_pad+1);
58  c_AsymmEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
59  c_AsymmEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
60  c_AsymmEnergy->cd(i_pad+1)->SetTicks(1,1);
61  c_AsymmEnergy->cd(i_pad+1)->SetGrid(0,0);
62  }
63  c_AsymmEnergy->cd(1);
64  h_mAsymmEnergy_pion_5GeV->Draw("colz");
65  h_mAsymmEnergy_pion_5GeV->ProfileX("p_mAsymmEnergy_pion_5GeV",1,-1,"i");
66  p_mAsymmEnergy_pion_5GeV->SetMarkerStyle(20);
67  p_mAsymmEnergy_pion_5GeV->SetMarkerColor(1);
68  p_mAsymmEnergy_pion_5GeV->SetMarkerSize(1.0);
69  p_mAsymmEnergy_pion_5GeV->Draw("pE same");
70 
71  c_AsymmEnergy->cd(2);
72  h_mAsymmEnergy_pion_8GeV->Draw("colz");
73  h_mAsymmEnergy_pion_8GeV->ProfileX("p_mAsymmEnergy_pion_8GeV",1,-1,"i");
74  p_mAsymmEnergy_pion_8GeV->SetMarkerStyle(20);
75  p_mAsymmEnergy_pion_8GeV->SetMarkerColor(1);
76  p_mAsymmEnergy_pion_8GeV->SetMarkerSize(1.0);
77  p_mAsymmEnergy_pion_8GeV->Draw("pE same");
78 
79  c_AsymmEnergy->cd(3);
80  h_mAsymmEnergy_pion_12GeV->Draw("colz");
81  h_mAsymmEnergy_pion_12GeV->ProfileX("p_mAsymmEnergy_pion_12GeV",1,-1,"i");
82  p_mAsymmEnergy_pion_12GeV->SetMarkerStyle(20);
83  p_mAsymmEnergy_pion_12GeV->SetMarkerColor(1);
84  p_mAsymmEnergy_pion_12GeV->SetMarkerSize(1.0);
85  p_mAsymmEnergy_pion_12GeV->Draw("pE same");
86 
87  c_AsymmEnergy->cd(4);
88  h_mAsymmEnergy_pion_60GeV->Draw("colz");
89  h_mAsymmEnergy_pion_60GeV->ProfileX("p_mAsymmEnergy_pion_60GeV",1,-1,"i");
90  p_mAsymmEnergy_pion_60GeV->SetMarkerStyle(20);
91  p_mAsymmEnergy_pion_60GeV->SetMarkerColor(1);
92  p_mAsymmEnergy_pion_60GeV->SetMarkerSize(1.0);
93  p_mAsymmEnergy_pion_60GeV->Draw("pE same");
94 
95  float mean_energy[4] = {5.0,8.0,12.0,60.0};
96  float val_mean[4];
97  float err_mean[4];
98  float val_sigma[3];
99  float err_sigma[3];
100  float val_resolution[3];
101  float err_resolution[3];
102 
103  TCanvas *c_Energy = new TCanvas("c_Energy","c_Energy",1600,400);
104  c_Energy->Divide(4,1);
105  for(int i_pad = 0; i_pad < 4; ++i_pad)
106  {
107  c_Energy->cd(i_pad+1);
108  c_Energy->cd(i_pad+1)->SetLeftMargin(0.15);
109  c_Energy->cd(i_pad+1)->SetBottomMargin(0.15);
110  c_Energy->cd(i_pad+1)->SetTicks(1,1);
111  c_Energy->cd(i_pad+1)->SetGrid(0,0);
112  }
113  c_Energy->cd(1);
114  h_mEnergy_pion_5GeV->Draw("hE");
115  TF1 *f_gaus_5GeV = new TF1("f_gaus_5GeV","gaus",0,100);
116  f_gaus_5GeV->SetParameter(0,1.0);
117  f_gaus_5GeV->SetParameter(1,h_mEnergy_pion_5GeV->GetMean());
118  f_gaus_5GeV->SetParameter(2,1.0);
119  f_gaus_5GeV->SetRange(1.2,2.8);
120  h_mEnergy_pion_5GeV->Fit(f_gaus_5GeV,"NR");
121  f_gaus_5GeV->SetLineColor(2);
122  f_gaus_5GeV->SetLineStyle(2);
123  f_gaus_5GeV->SetLineWidth(2);
124  f_gaus_5GeV->Draw("l same");
125  val_mean[0] = f_gaus_5GeV->GetParameter(1);
126  err_mean[0] = f_gaus_5GeV->GetParError(1);
127  val_sigma[0] = f_gaus_5GeV->GetParameter(2); // extract calibrated energy
128  err_sigma[0] = f_gaus_5GeV->GetParError(2);
129  val_resolution[0] = val_sigma[0]/val_mean[0];
130  err_resolution[0] = ErrDiv(val_sigma[0],val_mean[0],err_sigma[0],err_mean[0]);
131 
132  c_Energy->cd(2);
133  h_mEnergy_pion_8GeV->Draw("hE");
134  TF1 *f_gaus_8GeV = new TF1("f_gaus_8GeV","gaus",0,100);
135  f_gaus_8GeV->SetParameter(0,1.0);
136  f_gaus_8GeV->SetParameter(1,h_mEnergy_pion_8GeV->GetMean());
137  f_gaus_8GeV->SetParameter(2,1.0);
138  f_gaus_8GeV->SetRange(1.0,4.2);
139  h_mEnergy_pion_8GeV->Fit(f_gaus_8GeV,"NR");
140  f_gaus_8GeV->SetLineColor(2);
141  f_gaus_8GeV->SetLineStyle(2);
142  f_gaus_8GeV->SetLineWidth(2);
143  f_gaus_8GeV->Draw("l same");
144  val_mean[1] = f_gaus_8GeV->GetParameter(1);
145  err_mean[1] = f_gaus_8GeV->GetParError(1);
146  val_sigma[1] = f_gaus_8GeV->GetParameter(2);
147  err_sigma[1] = f_gaus_8GeV->GetParError(2);
148  val_resolution[1] = val_sigma[1]/val_mean[1];
149  err_resolution[1] = ErrDiv(val_sigma[1],val_mean[1],err_sigma[1],err_mean[1]);
150 
151  c_Energy->cd(3);
152  h_mEnergy_pion_12GeV->Draw("hE");
153  TF1 *f_gaus_12GeV = new TF1("f_gaus_12GeV","gaus",0,100);
154  f_gaus_12GeV->SetParameter(0,1.0);
155  f_gaus_12GeV->SetParameter(1,h_mEnergy_pion_12GeV->GetMean());
156  f_gaus_12GeV->SetParameter(2,1.0);
157  f_gaus_12GeV->SetRange(1.5,7.0);
158  h_mEnergy_pion_12GeV->Fit(f_gaus_12GeV,"NR");
159  f_gaus_12GeV->SetLineColor(2);
160  f_gaus_12GeV->SetLineStyle(2);
161  f_gaus_12GeV->SetLineWidth(2);
162  f_gaus_12GeV->Draw("l same");
163  val_mean[2] = f_gaus_12GeV->GetParameter(1);
164  err_mean[2] = f_gaus_12GeV->GetParError(1);
165  val_sigma[2] = f_gaus_12GeV->GetParameter(2);
166  err_sigma[2] = f_gaus_12GeV->GetParError(2);
167  val_resolution[2] = val_sigma[2]/val_mean[2];
168  err_resolution[2] = ErrDiv(val_sigma[2],val_mean[2],err_sigma[2],err_mean[2]);
169 
170  c_Energy->cd(4);
171  h_mEnergy_pion_60GeV->Draw("hE");
172  TF1 *f_gaus_60GeV = new TF1("f_gaus_60GeV","gaus",0,100);
173  f_gaus_60GeV->SetParameter(0,1.0);
174  f_gaus_60GeV->SetParameter(1,h_mEnergy_pion_60GeV->GetMean());
175  f_gaus_60GeV->SetParameter(2,1.0);
176  f_gaus_60GeV->SetRange(10.2,30.2);
177  h_mEnergy_pion_60GeV->Fit(f_gaus_60GeV,"NR");
178  f_gaus_60GeV->SetLineColor(2);
179  f_gaus_60GeV->SetLineStyle(2);
180  f_gaus_60GeV->SetLineWidth(2);
181  f_gaus_60GeV->Draw("l same");
182  val_mean[3] = f_gaus_60GeV->GetParameter(1);
183  err_mean[3] = f_gaus_60GeV->GetParError(1);
184 
185  c_Energy->SaveAs("../figures/HCAL_ShowerCalib/c_Energy_LevelingCorr.eps");
186 
187  TGraphAsymmErrors *g_lieanrity = new TGraphAsymmErrors();
188  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
189  for(int i_point = 0; i_point < 4; ++i_point)
190  {
191  g_lieanrity->SetPoint(i_point,mean_energy[i_point],val_mean[i_point]);
192  g_lieanrity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
193 
194  g_resolution->SetPoint(i_point,mean_energy[i_point],val_resolution[i_point]);
195  g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
196  }
197 
198  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
199  c_Linearity->cd();
200  c_Linearity->cd()->SetLeftMargin(0.15);
201  c_Linearity->cd()->SetBottomMargin(0.15);
202  c_Linearity->cd()->SetTicks(1,1);
203  c_Linearity->cd()->SetGrid(0,0);
204 
205  TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
206  for(int i_bin = 0; i_bin < 100; ++i_bin)
207  {
208  h_play->SetBinContent(i_bin+1,-10.0);
209  h_play->SetBinError(i_bin+1,1.0);
210  }
211  h_play->SetTitle("");
212  h_play->SetStats(0);
213  h_play->GetXaxis()->SetTitle("input Energy (GeV)");
214  h_play->GetXaxis()->CenterTitle();
215  h_play->GetXaxis()->SetNdivisions(505);
216  h_play->GetXaxis()->SetRangeUser(0.0,80.0);
217 
218  h_play->GetYaxis()->SetTitle("Tower Calibrated Energy (GeV)");
219  h_play->GetYaxis()->CenterTitle();
220  h_play->GetYaxis()->SetRangeUser(0.0,40.0);
221  h_play->DrawCopy("pE");
222 
223  g_lieanrity->SetMarkerStyle(24);
224  g_lieanrity->SetMarkerColor(2);
225  g_lieanrity->SetMarkerSize(1.2);
226  g_lieanrity->Draw("PE same");
227 
228  TF1 *f_pol = new TF1("f_pol","pol1",0.0,20.0);
229  f_pol->SetParameter(0,0.0);
230  f_pol->SetParameter(1,1.0);
231  f_pol->SetRange(0.0,80.0);
232  g_lieanrity->Fit(f_pol,"NR");
233  f_pol->SetLineColor(2);
234  f_pol->SetLineStyle(2);
235  f_pol->SetLineWidth(2);
236  f_pol->Draw("l same");
237 
238  c_Linearity->SaveAs("../figures/HCAL_ShowerCalib/c_Linearity_LevelingCorr.eps");
239 
240  /*
241  TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
242  l_unity->SetLineColor(4);
243  l_unity->SetLineStyle(2);
244  l_unity->SetLineWidth(2);
245  l_unity->Draw("l same");
246  */
247 
248  TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
249  c_Resolution->cd();
250  c_Resolution->cd()->SetLeftMargin(0.15);
251  c_Resolution->cd()->SetBottomMargin(0.15);
252  c_Resolution->cd()->SetTicks(1,1);
253  c_Resolution->cd()->SetGrid(0,0);
254 
255  h_play->GetYaxis()->SetTitle("$DeltaE/<E>");
256  h_play->GetYaxis()->CenterTitle();
257  h_play->GetYaxis()->SetRangeUser(0.0,0.8);
258  h_play->DrawCopy("pE");
259 
260  g_resolution->SetMarkerStyle(24);
261  g_resolution->SetMarkerColor(2);
262  g_resolution->SetMarkerSize(1.2);
263  g_resolution->Draw("PE same");
264 }