Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
extractEnergyLinearity_Pion.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file extractEnergyLinearity_Pion.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 extractEnergyLinearity_Pion()
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_pion_ShowerCalib");
32  TH1F *h_mEnergy_pion_5GeV = (TH1F*)h_mAsymmEnergy_pion_5GeV->ProjectionY()->Clone();
33 
34  string input_8GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0422.root";
35  TFile *File_8GeV = TFile::Open(input_8GeV.c_str());
36  TH2F *h_mAsymmEnergy_pion_8GeV = (TH2F*)File_8GeV->Get("h_mAsymmEnergy_pion_ShowerCalib");
37  TH1F *h_mEnergy_pion_8GeV = (TH1F*)h_mAsymmEnergy_pion_8GeV->ProjectionY()->Clone();
38 
39  string input_12GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0571.root";
40  TFile *File_12GeV = TFile::Open(input_12GeV.c_str());
41  TH2F *h_mAsymmEnergy_pion_12GeV = (TH2F*)File_12GeV->Get("h_mAsymmEnergy_pion_ShowerCalib");
42  TH1F *h_mEnergy_pion_12GeV = (TH1F*)h_mAsymmEnergy_pion_12GeV->ProjectionY()->Clone();
43 
44  string input_60GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0821.root";
45  TFile *File_60GeV = TFile::Open(input_60GeV.c_str());
46  TH2F *h_mAsymmEnergy_pion_60GeV = (TH2F*)File_60GeV->Get("h_mAsymmEnergy_pion_ShowerCalib");
47  TH1F *h_mEnergy_pion_60GeV = (TH1F*)h_mAsymmEnergy_pion_60GeV->ProjectionY()->Clone();
48 
49  TCanvas *c_AsymmEnergy = new TCanvas("c_AsymmEnergy","c_AsymmEnergy",1600,400);
50  c_AsymmEnergy->Divide(4,1);
51  for(int i_pad = 0; i_pad < 4; ++i_pad)
52  {
53  c_AsymmEnergy->cd(i_pad+1);
54  c_AsymmEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
55  c_AsymmEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
56  c_AsymmEnergy->cd(i_pad+1)->SetTicks(1,1);
57  c_AsymmEnergy->cd(i_pad+1)->SetGrid(0,0);
58  }
59  c_AsymmEnergy->cd(1);
60  h_mAsymmEnergy_pion_5GeV->Draw("colz");
61  h_mAsymmEnergy_pion_5GeV->ProfileX("p_mAsymmEnergy_pion_5GeV",1,-1,"i");
62  p_mAsymmEnergy_pion_5GeV->SetMarkerStyle(20);
63  p_mAsymmEnergy_pion_5GeV->SetMarkerColor(1);
64  p_mAsymmEnergy_pion_5GeV->SetMarkerSize(1.0);
65  // p_mAsymmEnergy_pion_5GeV->Draw("pE same");
66 
67  c_AsymmEnergy->cd(2);
68  h_mAsymmEnergy_pion_8GeV->Draw("colz");
69  h_mAsymmEnergy_pion_8GeV->ProfileX("p_mAsymmEnergy_pion_8GeV",1,-1,"i");
70  p_mAsymmEnergy_pion_8GeV->SetMarkerStyle(20);
71  p_mAsymmEnergy_pion_8GeV->SetMarkerColor(1);
72  p_mAsymmEnergy_pion_8GeV->SetMarkerSize(1.0);
73  // p_mAsymmEnergy_pion_8GeV->Draw("pE same");
74 
75  c_AsymmEnergy->cd(3);
76  h_mAsymmEnergy_pion_12GeV->Draw("colz");
77  h_mAsymmEnergy_pion_12GeV->ProfileX("p_mAsymmEnergy_pion_12GeV",1,-1,"i");
78  p_mAsymmEnergy_pion_12GeV->SetMarkerStyle(20);
79  p_mAsymmEnergy_pion_12GeV->SetMarkerColor(1);
80  p_mAsymmEnergy_pion_12GeV->SetMarkerSize(1.0);
81  // p_mAsymmEnergy_pion_12GeV->Draw("pE same");
82 
83  c_AsymmEnergy->cd(4);
84  h_mAsymmEnergy_pion_60GeV->Draw("colz");
85  h_mAsymmEnergy_pion_60GeV->ProfileX("p_mAsymmEnergy_pion_60GeV",1,-1,"i");
86  p_mAsymmEnergy_pion_60GeV->SetMarkerStyle(20);
87  p_mAsymmEnergy_pion_60GeV->SetMarkerColor(1);
88  p_mAsymmEnergy_pion_60GeV->SetMarkerSize(1.0);
89  // p_mAsymmEnergy_pion_60GeV->Draw("pE same");
90 
91  c_AsymmEnergy->SaveAs("../figures/HCAL_ShowerCalib/c_AsymmEnergy_ShowerCalib.eps");
92 
93  float mean_energy[4] = {5.0,8.0,12.0,60.0};
94  float val_mean[4];
95  float err_mean[4];
96  float val_sigma[4];
97  float err_sigma[4];
98  float val_resolution[4];
99  float err_resolution[4];
100 
101  TCanvas *c_Energy = new TCanvas("c_Energy","c_Energy",1600,400);
102  c_Energy->Divide(4,1);
103  for(int i_pad = 0; i_pad < 4; ++i_pad)
104  {
105  c_Energy->cd(i_pad+1);
106  c_Energy->cd(i_pad+1)->SetLeftMargin(0.15);
107  c_Energy->cd(i_pad+1)->SetBottomMargin(0.15);
108  c_Energy->cd(i_pad+1)->SetTicks(1,1);
109  c_Energy->cd(i_pad+1)->SetGrid(0,0);
110  }
111  c_Energy->cd(1);
112  h_mEnergy_pion_5GeV->Draw("hE");
113  TF1 *f_gaus_5GeV = new TF1("f_gaus_5GeV","gaus",0,100);
114  f_gaus_5GeV->SetParameter(0,1.0);
115  f_gaus_5GeV->SetParameter(1,h_mEnergy_pion_5GeV->GetMean());
116  f_gaus_5GeV->SetParameter(2,1.0);
117  f_gaus_5GeV->SetRange(3.5,10.0);
118  h_mEnergy_pion_5GeV->Fit(f_gaus_5GeV,"NR");
119  f_gaus_5GeV->SetLineColor(2);
120  f_gaus_5GeV->SetLineStyle(2);
121  f_gaus_5GeV->SetLineWidth(2);
122  f_gaus_5GeV->Draw("l same");
123  val_mean[0] = f_gaus_5GeV->GetParameter(1);
124  err_mean[0] = f_gaus_5GeV->GetParError(1);
125  val_sigma[0] = f_gaus_5GeV->GetParameter(2); // extract calibrated energy
126  err_sigma[0] = f_gaus_5GeV->GetParError(2);
127  val_resolution[0] = val_sigma[0]/val_mean[0];
128  err_resolution[0] = ErrDiv(val_sigma[0],val_mean[0],err_sigma[0],err_mean[0]);
129 
130  c_Energy->cd(2);
131  h_mEnergy_pion_8GeV->Draw("hE");
132  TF1 *f_gaus_8GeV = new TF1("f_gaus_8GeV","gaus",0,100);
133  f_gaus_8GeV->SetParameter(0,1.0);
134  f_gaus_8GeV->SetParameter(1,h_mEnergy_pion_8GeV->GetMean());
135  f_gaus_8GeV->SetParameter(2,1.0);
136  f_gaus_8GeV->SetRange(3.0,14.2);
137  h_mEnergy_pion_8GeV->Fit(f_gaus_8GeV,"NR");
138  f_gaus_8GeV->SetLineColor(2);
139  f_gaus_8GeV->SetLineStyle(2);
140  f_gaus_8GeV->SetLineWidth(2);
141  f_gaus_8GeV->Draw("l same");
142  val_mean[1] = f_gaus_8GeV->GetParameter(1);
143  err_mean[1] = f_gaus_8GeV->GetParError(1);
144  val_sigma[1] = f_gaus_8GeV->GetParameter(2);
145  err_sigma[1] = f_gaus_8GeV->GetParError(2);
146  val_resolution[1] = val_sigma[1]/val_mean[1];
147  err_resolution[1] = ErrDiv(val_sigma[1],val_mean[1],err_sigma[1],err_mean[1]);
148 
149  c_Energy->cd(3);
150  h_mEnergy_pion_12GeV->Draw("hE");
151  TF1 *f_gaus_12GeV = new TF1("f_gaus_12GeV","gaus",0,100);
152  f_gaus_12GeV->SetParameter(0,1.0);
153  f_gaus_12GeV->SetParameter(1,h_mEnergy_pion_12GeV->GetMean());
154  f_gaus_12GeV->SetParameter(2,1.0);
155  f_gaus_12GeV->SetRange(5.0,20.0);
156  h_mEnergy_pion_12GeV->Fit(f_gaus_12GeV,"NR");
157  f_gaus_12GeV->SetLineColor(2);
158  f_gaus_12GeV->SetLineStyle(2);
159  f_gaus_12GeV->SetLineWidth(2);
160  f_gaus_12GeV->Draw("l same");
161  val_mean[2] = f_gaus_12GeV->GetParameter(1);
162  err_mean[2] = f_gaus_12GeV->GetParError(1);
163  val_sigma[2] = f_gaus_12GeV->GetParameter(2);
164  err_sigma[2] = f_gaus_12GeV->GetParError(2);
165  val_resolution[2] = val_sigma[2]/val_mean[2];
166  err_resolution[2] = ErrDiv(val_sigma[2],val_mean[2],err_sigma[2],err_mean[2]);
167 
168  c_Energy->cd(4);
169  h_mEnergy_pion_60GeV->Draw("hE");
170  TF1 *f_gaus_60GeV = new TF1("f_gaus_60GeV","gaus",0,100);
171  f_gaus_60GeV->SetParameter(0,1.0);
172  f_gaus_60GeV->SetParameter(1,h_mEnergy_pion_60GeV->GetMean());
173  f_gaus_60GeV->SetParameter(2,1.0);
174  f_gaus_60GeV->SetRange(20.2,80.2);
175  h_mEnergy_pion_60GeV->Fit(f_gaus_60GeV,"NR");
176  f_gaus_60GeV->SetLineColor(2);
177  f_gaus_60GeV->SetLineStyle(2);
178  f_gaus_60GeV->SetLineWidth(2);
179  f_gaus_60GeV->Draw("l same");
180  val_mean[3] = f_gaus_60GeV->GetParameter(1);
181  err_mean[3] = f_gaus_60GeV->GetParError(1);
182  val_sigma[3] = f_gaus_60GeV->GetParameter(2);
183  err_sigma[3] = f_gaus_60GeV->GetParError(2);
184  val_resolution[3] = val_sigma[3]/val_mean[3];
185  err_resolution[3] = ErrDiv(val_sigma[3],val_mean[3],err_sigma[3],err_mean[3]);
186 
187  c_Energy->SaveAs("../figures/HCAL_ShowerCalib/c_Energy_ShowerCalib.eps");
188 
189  //------------------------linearity------------------------
190  TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
191  for(int i_point = 0; i_point < 4; ++i_point)
192  {
193  g_linearity->SetPoint(i_point,mean_energy[i_point],val_mean[i_point]);
194  g_linearity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
195  }
196 
197  TGraphAsymmErrors *g_linearity_old = new TGraphAsymmErrors();
198  string inputlinearity = "linearity.txt";
199  FILE *flinear = fopen(inputlinearity.c_str(), "r");
200  if(flinear == NULL)
201  {
202  perror("Error opening file!");
203  }
204  else
205  {
206  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
207  char line[1024];
208  int line_counter = 0;
209  while(fgets(line,1024,flinear))
210  {
211  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
212  // cout << "x = " << x << ", y = " << y << ", err_x_low = " << err_x_low << ", err_x_high = " << err_x_high << ", err_y_low = " << err_y_low << ", err_y_high = " << err_y_high << endl;
213  g_linearity_old->SetPoint(line_counter,x,y);
214  g_linearity_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
215  line_counter++;
216  }
217  }
218 
219  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
220  c_Linearity->cd();
221  c_Linearity->cd()->SetLeftMargin(0.15);
222  c_Linearity->cd()->SetBottomMargin(0.15);
223  c_Linearity->cd()->SetTicks(1,1);
224  c_Linearity->cd()->SetGrid(0,0);
225 
226  TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
227  for(int i_bin = 0; i_bin < 100; ++i_bin)
228  {
229  h_play->SetBinContent(i_bin+1,-10.0);
230  h_play->SetBinError(i_bin+1,1.0);
231  }
232  h_play->SetTitle("");
233  h_play->SetStats(0);
234  h_play->GetXaxis()->SetTitle("input Energy (GeV)");
235  h_play->GetXaxis()->CenterTitle();
236  h_play->GetXaxis()->SetNdivisions(505);
237  h_play->GetXaxis()->SetRangeUser(0.0,80.0);
238 
239  h_play->GetYaxis()->SetTitle("Tower Calibrated Energy (GeV)");
240  h_play->GetYaxis()->CenterTitle();
241  h_play->GetYaxis()->SetRangeUser(0.0,80.0);
242  h_play->DrawCopy("pE");
243 
244  g_linearity->SetMarkerStyle(20);
245  g_linearity->SetMarkerColor(kGray+2);
246  g_linearity->SetMarkerSize(2.0);
247  g_linearity->Draw("pE same");
248 
249  g_linearity_old->SetMarkerStyle(21);
250  g_linearity_old->SetMarkerColor(2);
251  g_linearity_old->SetMarkerSize(1.8);
252  g_linearity_old->Draw("pE same");
253 
254  /*
255  TF1 *f_pol = new TF1("f_pol","pol1",0.0,20.0);
256  f_pol->FixParameter(0,0.0);
257  f_pol->SetParameter(1,1.0);
258  f_pol->SetRange(0.0,80.0);
259  g_lieanrity->Fit(f_pol,"NR");
260  f_pol->SetLineColor(2);
261  f_pol->SetLineStyle(2);
262  f_pol->SetLineWidth(2);
263  f_pol->Draw("l same");
264  */
265 
266  TLine *l_unity = new TLine(1.0,1.0,79.0,79.0);
267  l_unity->SetLineColor(4);
268  l_unity->SetLineStyle(2);
269  l_unity->SetLineWidth(2);
270  l_unity->Draw("l same");
271 
272  TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
273  leg_linear->SetBorderSize(0);
274  leg_linear->SetFillColor(0);
275  leg_linear->AddEntry(g_linearity_old,"#pi- T1044-2017","p");
276  leg_linear->AddEntry(g_linearity,"#pi- T1044-2018","p");
277  leg_linear->AddEntry(l_unity,"unity","l");
278  leg_linear->Draw("same");
279 
280  c_Linearity->SaveAs("../figures/HCAL_ShowerCalib/c_Linearity_ShowerCalib.eps");
281 
282  //------------------------linearity------------------------
283 
284  //------------------------resolution-----------------------
285  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
286  for(int i_point = 0; i_point < 4; ++i_point)
287  {
288  g_resolution->SetPoint(i_point,mean_energy[i_point],val_resolution[i_point]);
289  g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
290  }
291 
292  TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
293  string inputresolution = "resolution.txt";
294  FILE *fres = fopen(inputresolution.c_str(), "r");
295  if(fres == NULL)
296  {
297  perror("Error opening file!");
298  }
299  else
300  {
301  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
302  char line[1024];
303  int line_counter = 0;
304  while(fgets(line,1024,fres))
305  {
306  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
307  // cout << "x = " << x << ", y = " << y << ", err_x_low = " << err_x_low << ", err_x_high = " << err_x_high << ", err_y_low = " << err_y_low << ", err_y_high = " << err_y_high << endl;
308  g_resolution_old->SetPoint(line_counter,x,y);
309  g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
310  line_counter++;
311  }
312  }
313  TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
314  c_Resolution->cd();
315  c_Resolution->cd()->SetLeftMargin(0.15);
316  c_Resolution->cd()->SetBottomMargin(0.15);
317  c_Resolution->cd()->SetTicks(1,1);
318  c_Resolution->cd()->SetGrid(0,0);
319 
320  h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
321  h_play->GetYaxis()->CenterTitle();
322  h_play->GetYaxis()->SetRangeUser(0.0,0.8);
323  h_play->DrawCopy("pE");
324 
325  g_resolution->SetMarkerStyle(20);
326  g_resolution->SetMarkerColor(kGray+2);
327  g_resolution->SetMarkerSize(2.0);
328  g_resolution->Draw("pE same");
329 
330  g_resolution_old->SetMarkerStyle(21);
331  g_resolution_old->SetMarkerColor(2);
332  g_resolution_old->SetMarkerSize(1.8);
333  g_resolution_old->Draw("pE same");
334 
335  TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
336  leg_res->SetBorderSize(0);
337  leg_res->SetFillColor(0);
338  leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
339  leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
340  leg_res->Draw("same");
341 
342  c_Resolution->SaveAs("../figures/HCAL_ShowerCalib/c_Resolution_ShowerCalib.eps");
343  //------------------------resolution-----------------------
344 
345  string outputfile;
346  outputfile = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/T1044_2018_pion.root";
347  TFile *File_OutPut = new TFile(outputfile.c_str(),"RECREATE");
348  File_OutPut->cd();
349  g_linearity->SetName("g_linearity_2018_pion");
350  g_linearity->Write();
351  g_resolution->SetName("g_resolution_2018_pion");
352  g_resolution->Write();
353  File_OutPut->Close();
354 
355  outputfile = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/T1044_2017_pion.root";
356  TFile *File_OutPut_old = new TFile(outputfile.c_str(),"RECREATE");
357  File_OutPut_old->cd();
358  g_linearity_old->SetName("g_linearity_2017_pion");
359  g_linearity_old->Write();
360  g_resolution_old->SetName("g_resolution_2017_pion");
361  g_resolution_old->Write();
362  File_OutPut_old->Close();
363 }