Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
extractEnergyLinearity_OHCal.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file extractEnergyLinearity_OHCal.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_OHCal()
28 {
29  string input_5GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_1087.root";
30  TFile *File_5GeV = TFile::Open(input_5GeV.c_str());
31  TH1F *h_mEnergy_pion_5GeV = (TH1F*)File_5GeV->Get("h_mEnergyOut_pion_ShowerCalib");
32 
33  string input_8GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0422.root";
34  TFile *File_8GeV = TFile::Open(input_8GeV.c_str());
35  TH1F *h_mEnergy_pion_8GeV = (TH1F*)File_8GeV->Get("h_mEnergyOut_pion_ShowerCalib");
36 
37  string input_12GeV = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_0571.root";
38  TFile *File_12GeV = TFile::Open(input_12GeV.c_str());
39  TH1F *h_mEnergy_pion_12GeV = (TH1F*)File_12GeV->Get("h_mEnergyOut_pion_ShowerCalib");
40 
41  float mean_energy[4] = {5.0,8.0,12.0,60.0};
42  float val_mean[4];
43  float err_mean[4];
44  float val_sigma[4];
45  float err_sigma[4];
46  float val_resolution[4];
47  float err_resolution[4];
48 
49  TCanvas *c_Energy = new TCanvas("c_Energy","c_Energy",1500,500);
50  c_Energy->Divide(3,1);
51  for(int i_pad = 0; i_pad < 3; ++i_pad)
52  {
53  c_Energy->cd(i_pad+1);
54  c_Energy->cd(i_pad+1)->SetLeftMargin(0.15);
55  c_Energy->cd(i_pad+1)->SetBottomMargin(0.15);
56  c_Energy->cd(i_pad+1)->SetTicks(1,1);
57  c_Energy->cd(i_pad+1)->SetGrid(0,0);
58  }
59  c_Energy->cd(1);
60  h_mEnergy_pion_5GeV->SetTitle("5 GeV");
61  h_mEnergy_pion_5GeV->GetXaxis()->SetTitle("Shower Calibrated Energy (GeV)");
62  h_mEnergy_pion_5GeV->Draw("hE");
63  TF1 *f_gaus_5GeV = new TF1("f_gaus_5GeV","gaus",0,100);
64  f_gaus_5GeV->SetParameter(0,1.0);
65  f_gaus_5GeV->SetParameter(1,h_mEnergy_pion_5GeV->GetMean());
66  f_gaus_5GeV->SetParameter(2,1.0);
67  f_gaus_5GeV->SetRange(3.5,10.0);
68  h_mEnergy_pion_5GeV->Fit(f_gaus_5GeV,"NR");
69  f_gaus_5GeV->SetLineColor(2);
70  f_gaus_5GeV->SetLineStyle(2);
71  f_gaus_5GeV->SetLineWidth(2);
72  f_gaus_5GeV->Draw("l same");
73  val_mean[0] = f_gaus_5GeV->GetParameter(1);
74  err_mean[0] = f_gaus_5GeV->GetParError(1);
75  val_sigma[0] = f_gaus_5GeV->GetParameter(2); // extract calibrated energy
76  err_sigma[0] = f_gaus_5GeV->GetParError(2);
77  val_resolution[0] = val_sigma[0]/val_mean[0];
78  err_resolution[0] = ErrDiv(val_sigma[0],val_mean[0],err_sigma[0],err_mean[0]);
79 
80  c_Energy->cd(2);
81  h_mEnergy_pion_8GeV->SetTitle("8 GeV");
82  h_mEnergy_pion_8GeV->GetXaxis()->SetTitle("Shower Calibrated Energy (GeV)");
83  h_mEnergy_pion_8GeV->Draw("hE");
84  TF1 *f_gaus_8GeV = new TF1("f_gaus_8GeV","gaus",0,100);
85  f_gaus_8GeV->SetParameter(0,1.0);
86  f_gaus_8GeV->SetParameter(1,h_mEnergy_pion_8GeV->GetMean());
87  f_gaus_8GeV->SetParameter(2,1.0);
88  f_gaus_8GeV->SetRange(3.0,14.2);
89  h_mEnergy_pion_8GeV->Fit(f_gaus_8GeV,"NR");
90  f_gaus_8GeV->SetLineColor(2);
91  f_gaus_8GeV->SetLineStyle(2);
92  f_gaus_8GeV->SetLineWidth(2);
93  f_gaus_8GeV->Draw("l same");
94  val_mean[1] = f_gaus_8GeV->GetParameter(1);
95  err_mean[1] = f_gaus_8GeV->GetParError(1);
96  val_sigma[1] = f_gaus_8GeV->GetParameter(2);
97  err_sigma[1] = f_gaus_8GeV->GetParError(2);
98  val_resolution[1] = val_sigma[1]/val_mean[1];
99  err_resolution[1] = ErrDiv(val_sigma[1],val_mean[1],err_sigma[1],err_mean[1]);
100 
101  c_Energy->cd(3);
102  h_mEnergy_pion_12GeV->SetTitle("12 GeV");
103  h_mEnergy_pion_12GeV->GetXaxis()->SetTitle("Shower Calibrated Energy (GeV)");
104  h_mEnergy_pion_12GeV->Draw("hE");
105  TF1 *f_gaus_12GeV = new TF1("f_gaus_12GeV","gaus",0,100);
106  f_gaus_12GeV->SetParameter(0,1.0);
107  f_gaus_12GeV->SetParameter(1,h_mEnergy_pion_12GeV->GetMean());
108  f_gaus_12GeV->SetParameter(2,1.0);
109  f_gaus_12GeV->SetRange(5.0,20.0);
110  h_mEnergy_pion_12GeV->Fit(f_gaus_12GeV,"NR");
111  f_gaus_12GeV->SetLineColor(2);
112  f_gaus_12GeV->SetLineStyle(2);
113  f_gaus_12GeV->SetLineWidth(2);
114  f_gaus_12GeV->Draw("l same");
115  val_mean[2] = f_gaus_12GeV->GetParameter(1);
116  err_mean[2] = f_gaus_12GeV->GetParError(1);
117  val_sigma[2] = f_gaus_12GeV->GetParameter(2);
118  err_sigma[2] = f_gaus_12GeV->GetParError(2);
119  val_resolution[2] = val_sigma[2]/val_mean[2];
120  err_resolution[2] = ErrDiv(val_sigma[2],val_mean[2],err_sigma[2],err_mean[2]);
121 
122  c_Energy->SaveAs("../figures/HCAL_ShowerCalib/c_Energy_ShowerCalib_OHCAL.eps");
123 
124  //------------------------linearity------------------------
125  TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
126  for(int i_point = 0; i_point < 3; ++i_point)
127  {
128  g_linearity->SetPoint(i_point,mean_energy[i_point],val_mean[i_point]);
129  g_linearity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
130  }
131 
132  TGraphAsymmErrors *g_linearity_old = new TGraphAsymmErrors();
133  string inputlinearity = "linearity.txt";
134  FILE *flinear = fopen(inputlinearity.c_str(), "r");
135  if(flinear == NULL)
136  {
137  perror("Error opening file!");
138  }
139  else
140  {
141  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
142  char line[1024];
143  int line_counter = 0;
144  while(fgets(line,1024,flinear))
145  {
146  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
147  // 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;
148  g_linearity_old->SetPoint(line_counter,x,y);
149  g_linearity_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
150  line_counter++;
151  }
152  }
153 
154  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
155  c_Linearity->cd();
156  c_Linearity->cd()->SetLeftMargin(0.15);
157  c_Linearity->cd()->SetBottomMargin(0.15);
158  c_Linearity->cd()->SetTicks(1,1);
159  c_Linearity->cd()->SetGrid(0,0);
160 
161  TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
162  for(int i_bin = 0; i_bin < 100; ++i_bin)
163  {
164  h_play->SetBinContent(i_bin+1,-10.0);
165  h_play->SetBinError(i_bin+1,1.0);
166  }
167  h_play->SetTitle("");
168  h_play->SetStats(0);
169  h_play->GetXaxis()->SetTitle("input Energy (GeV)");
170  h_play->GetXaxis()->CenterTitle();
171  h_play->GetXaxis()->SetNdivisions(505);
172  h_play->GetXaxis()->SetRangeUser(0.0,40.0);
173 
174  h_play->GetYaxis()->SetTitle("Shower Calibrated Energy (GeV)");
175  h_play->GetYaxis()->CenterTitle();
176  h_play->GetYaxis()->SetRangeUser(0.0,40.0);
177  h_play->DrawCopy("pE");
178 
179  g_linearity->SetMarkerStyle(20);
180  g_linearity->SetMarkerColor(kGray+2);
181  g_linearity->SetMarkerSize(2.0);
182  g_linearity->Draw("pE same");
183 
184  g_linearity_old->SetMarkerStyle(21);
185  g_linearity_old->SetMarkerColor(2);
186  g_linearity_old->SetMarkerSize(1.8);
187  g_linearity_old->Draw("pE same");
188 
189  /*
190  TF1 *f_pol = new TF1("f_pol","pol1",0.0,20.0);
191  f_pol->FixParameter(0,0.0);
192  f_pol->SetParameter(1,1.0);
193  f_pol->SetRange(0.0,80.0);
194  g_lieanrity->Fit(f_pol,"NR");
195  f_pol->SetLineColor(2);
196  f_pol->SetLineStyle(2);
197  f_pol->SetLineWidth(2);
198  f_pol->Draw("l same");
199  */
200 
201  TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
202  l_unity->SetLineColor(4);
203  l_unity->SetLineStyle(2);
204  l_unity->SetLineWidth(2);
205  l_unity->Draw("l same");
206 
207  TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
208  leg_linear->SetBorderSize(0);
209  leg_linear->SetFillColor(0);
210  leg_linear->AddEntry(g_linearity_old,"#pi- T1044-2017","p");
211  leg_linear->AddEntry(g_linearity,"#pi- T1044-2018","p");
212  leg_linear->AddEntry(l_unity,"unity","l");
213  leg_linear->Draw("same");
214 
215  c_Linearity->SaveAs("../figures/HCAL_ShowerCalib/c_Linearity_ShowerCalib_OHCAL.eps");
216 
217  //------------------------linearity------------------------
218 
219  //------------------------resolution-----------------------
220  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
221  for(int i_point = 0; i_point < 3; ++i_point)
222  {
223  g_resolution->SetPoint(i_point,mean_energy[i_point],val_resolution[i_point]);
224  g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
225  }
226 
227  TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
228  string inputresolution = "resolution.txt";
229  FILE *fres = fopen(inputresolution.c_str(), "r");
230  if(fres == NULL)
231  {
232  perror("Error opening file!");
233  }
234  else
235  {
236  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
237  char line[1024];
238  int line_counter = 0;
239  while(fgets(line,1024,fres))
240  {
241  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
242  // 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;
243  g_resolution_old->SetPoint(line_counter,x,y);
244  g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
245  line_counter++;
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(20);
261  g_resolution->SetMarkerColor(kGray+2);
262  g_resolution->SetMarkerSize(2.0);
263  g_resolution->Draw("pE same");
264 
265  g_resolution_old->SetMarkerStyle(21);
266  g_resolution_old->SetMarkerColor(2);
267  g_resolution_old->SetMarkerSize(1.8);
268  g_resolution_old->Draw("pE same");
269 
270  TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
271  leg_res->SetBorderSize(0);
272  leg_res->SetFillColor(0);
273  leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
274  leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
275  leg_res->Draw("same");
276 
277  c_Resolution->SaveAs("../figures/HCAL_ShowerCalib/c_Resolution_ShowerCalib_OHCAL.eps");
278  //------------------------resolution-----------------------
279 
280  string outputfile;
281  outputfile = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/T1044_2018_pion_hcalout.root";
282  TFile *File_OutPut = new TFile(outputfile.c_str(),"RECREATE");
283  File_OutPut->cd();
284  g_linearity->SetName("g_linearity_2018_pion_hcalout");
285  g_linearity->Write();
286  g_resolution->SetName("g_resolution_2018_pion_hcalout");
287  g_resolution->Write();
288  File_OutPut->Close();
289 }