Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
extractEnergyLinearity_Electron.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file extractEnergyLinearity_Electron.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_Electron()
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_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_electron_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_electron_ShowerCalib");
42  TH1F *h_mEnergy_pion_12GeV = (TH1F*)h_mAsymmEnergy_pion_12GeV->ProjectionY()->Clone();
43 
44  TCanvas *c_AsymmEnergy = new TCanvas("c_AsymmEnergy","c_AsymmEnergy",1500,500);
45  c_AsymmEnergy->Divide(3,1);
46  for(int i_pad = 0; i_pad < 3; ++i_pad)
47  {
48  c_AsymmEnergy->cd(i_pad+1);
49  c_AsymmEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
50  c_AsymmEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
51  c_AsymmEnergy->cd(i_pad+1)->SetTicks(1,1);
52  c_AsymmEnergy->cd(i_pad+1)->SetGrid(0,0);
53  }
54  c_AsymmEnergy->cd(1);
55  h_mAsymmEnergy_pion_5GeV->Draw("colz");
56  h_mAsymmEnergy_pion_5GeV->ProfileX("p_mAsymmEnergy_pion_5GeV",1,-1,"i");
57  p_mAsymmEnergy_pion_5GeV->SetMarkerStyle(20);
58  p_mAsymmEnergy_pion_5GeV->SetMarkerColor(1);
59  p_mAsymmEnergy_pion_5GeV->SetMarkerSize(1.0);
60  // p_mAsymmEnergy_pion_5GeV->Draw("pE same");
61 
62  c_AsymmEnergy->cd(2);
63  h_mAsymmEnergy_pion_8GeV->Draw("colz");
64  h_mAsymmEnergy_pion_8GeV->ProfileX("p_mAsymmEnergy_pion_8GeV",1,-1,"i");
65  p_mAsymmEnergy_pion_8GeV->SetMarkerStyle(20);
66  p_mAsymmEnergy_pion_8GeV->SetMarkerColor(1);
67  p_mAsymmEnergy_pion_8GeV->SetMarkerSize(1.0);
68  // p_mAsymmEnergy_pion_8GeV->Draw("pE same");
69 
70  c_AsymmEnergy->cd(3);
71  h_mAsymmEnergy_pion_12GeV->Draw("colz");
72  h_mAsymmEnergy_pion_12GeV->ProfileX("p_mAsymmEnergy_pion_12GeV",1,-1,"i");
73  p_mAsymmEnergy_pion_12GeV->SetMarkerStyle(20);
74  p_mAsymmEnergy_pion_12GeV->SetMarkerColor(1);
75  p_mAsymmEnergy_pion_12GeV->SetMarkerSize(1.0);
76  // p_mAsymmEnergy_pion_12GeV->Draw("pE same");
77 
78  c_AsymmEnergy->SaveAs("../figures/HCAL_ShowerCalib/c_AsymmEnergy_ShowerCalib.eps");
79 
80  float mean_energy[4] = {5.0,8.0,12.0,60.0};
81  float val_mean[4];
82  float err_mean[4];
83  float val_sigma[4];
84  float err_sigma[4];
85  float val_resolution[4];
86  float err_resolution[4];
87 
88  TCanvas *c_Energy = new TCanvas("c_Energy","c_Energy",1500,500);
89  c_Energy->Divide(3,1);
90  for(int i_pad = 0; i_pad < 3; ++i_pad)
91  {
92  c_Energy->cd(i_pad+1);
93  c_Energy->cd(i_pad+1)->SetLeftMargin(0.15);
94  c_Energy->cd(i_pad+1)->SetBottomMargin(0.15);
95  c_Energy->cd(i_pad+1)->SetTicks(1,1);
96  c_Energy->cd(i_pad+1)->SetGrid(0,0);
97  }
98  c_Energy->cd(1);
99  h_mEnergy_pion_5GeV->Draw("hE");
100  TF1 *f_gaus_5GeV = new TF1("f_gaus_5GeV","gaus",0,100);
101  f_gaus_5GeV->SetParameter(0,1.0);
102  f_gaus_5GeV->SetParameter(1,h_mEnergy_pion_5GeV->GetMean());
103  f_gaus_5GeV->SetParameter(2,1.0);
104  f_gaus_5GeV->SetRange(5.0,10.0);
105  h_mEnergy_pion_5GeV->Fit(f_gaus_5GeV,"NR");
106  f_gaus_5GeV->SetLineColor(2);
107  f_gaus_5GeV->SetLineStyle(2);
108  f_gaus_5GeV->SetLineWidth(2);
109  f_gaus_5GeV->Draw("l same");
110  val_mean[0] = f_gaus_5GeV->GetParameter(1);
111  err_mean[0] = f_gaus_5GeV->GetParError(1);
112  val_sigma[0] = f_gaus_5GeV->GetParameter(2); // extract calibrated energy
113  err_sigma[0] = f_gaus_5GeV->GetParError(2);
114  val_resolution[0] = val_sigma[0]/val_mean[0];
115  err_resolution[0] = ErrDiv(val_sigma[0],val_mean[0],err_sigma[0],err_mean[0]);
116 
117  c_Energy->cd(2);
118  h_mEnergy_pion_8GeV->Draw("hE");
119  TF1 *f_gaus_8GeV = new TF1("f_gaus_8GeV","gaus",0,100);
120  f_gaus_8GeV->SetParameter(0,1.0);
121  f_gaus_8GeV->SetParameter(1,h_mEnergy_pion_8GeV->GetMean());
122  f_gaus_8GeV->SetParameter(2,1.0);
123  f_gaus_8GeV->SetRange(7.0,15.0);
124  h_mEnergy_pion_8GeV->Fit(f_gaus_8GeV,"NR");
125  f_gaus_8GeV->SetLineColor(2);
126  f_gaus_8GeV->SetLineStyle(2);
127  f_gaus_8GeV->SetLineWidth(2);
128  f_gaus_8GeV->Draw("l same");
129  val_mean[1] = f_gaus_8GeV->GetParameter(1);
130  err_mean[1] = f_gaus_8GeV->GetParError(1);
131  val_sigma[1] = f_gaus_8GeV->GetParameter(2);
132  err_sigma[1] = f_gaus_8GeV->GetParError(2);
133  val_resolution[1] = val_sigma[1]/val_mean[1];
134  err_resolution[1] = ErrDiv(val_sigma[1],val_mean[1],err_sigma[1],err_mean[1]);
135 
136  c_Energy->cd(3);
137  h_mEnergy_pion_12GeV->Draw("hE");
138  TF1 *f_gaus_12GeV = new TF1("f_gaus_12GeV","gaus",0,100);
139  f_gaus_12GeV->SetParameter(0,1.0);
140  f_gaus_12GeV->SetParameter(1,h_mEnergy_pion_12GeV->GetMean());
141  f_gaus_12GeV->SetParameter(2,1.0);
142  f_gaus_12GeV->SetRange(10.0,24.0);
143  h_mEnergy_pion_12GeV->Fit(f_gaus_12GeV,"NR");
144  f_gaus_12GeV->SetLineColor(2);
145  f_gaus_12GeV->SetLineStyle(2);
146  f_gaus_12GeV->SetLineWidth(2);
147  f_gaus_12GeV->Draw("l same");
148  val_mean[2] = f_gaus_12GeV->GetParameter(1);
149  err_mean[2] = f_gaus_12GeV->GetParError(1);
150  val_sigma[2] = f_gaus_12GeV->GetParameter(2);
151  err_sigma[2] = f_gaus_12GeV->GetParError(2);
152  val_resolution[2] = val_sigma[2]/val_mean[2];
153  err_resolution[2] = ErrDiv(val_sigma[2],val_mean[2],err_sigma[2],err_mean[2]);
154 
155  c_Energy->SaveAs("../figures/HCAL_ShowerCalib/c_Energy_ShowerCalib.eps");
156 
157  //------------------------linearity------------------------
158  TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
159  for(int i_point = 0; i_point < 3; ++i_point)
160  {
161  g_linearity->SetPoint(i_point,mean_energy[i_point],val_mean[i_point]);
162  g_linearity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
163  }
164 
165  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
166  c_Linearity->cd();
167  c_Linearity->cd()->SetLeftMargin(0.15);
168  c_Linearity->cd()->SetBottomMargin(0.15);
169  c_Linearity->cd()->SetTicks(1,1);
170  c_Linearity->cd()->SetGrid(0,0);
171 
172  TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
173  for(int i_bin = 0; i_bin < 100; ++i_bin)
174  {
175  h_play->SetBinContent(i_bin+1,-10.0);
176  h_play->SetBinError(i_bin+1,1.0);
177  }
178  h_play->SetTitle("");
179  h_play->SetStats(0);
180  h_play->GetXaxis()->SetTitle("input Energy (GeV)");
181  h_play->GetXaxis()->CenterTitle();
182  h_play->GetXaxis()->SetNdivisions(505);
183  h_play->GetXaxis()->SetRangeUser(0.0,20.0);
184 
185  h_play->GetYaxis()->SetTitle("Tower Calibrated Energy (GeV)");
186  h_play->GetYaxis()->CenterTitle();
187  h_play->GetYaxis()->SetRangeUser(0.0,20.0);
188  h_play->DrawCopy("pE");
189 
190  g_linearity->SetMarkerStyle(20);
191  g_linearity->SetMarkerColor(kGray+2);
192  g_linearity->SetMarkerSize(2.0);
193  g_linearity->Draw("pE same");
194 
195  TLine *l_unity = new TLine(1.0,1.0,79.0,79.0);
196  l_unity->SetLineColor(4);
197  l_unity->SetLineStyle(2);
198  l_unity->SetLineWidth(2);
199  l_unity->Draw("l same");
200 
201  TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
202  leg_linear->SetBorderSize(0);
203  leg_linear->SetFillColor(0);
204  leg_linear->AddEntry(g_linearity,"e- T1044-2018","p");
205  leg_linear->AddEntry(l_unity,"unity","l");
206  leg_linear->Draw("same");
207 
208  c_Linearity->SaveAs("../figures/HCAL_ShowerCalib/c_Linearity_ShowerCalib.eps");
209 
210  //------------------------linearity------------------------
211 
212  //------------------------resolution-----------------------
213  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
214  for(int i_point = 0; i_point < 3; ++i_point)
215  {
216  g_resolution->SetPoint(i_point,mean_energy[i_point],val_resolution[i_point]);
217  g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
218  }
219 
220  TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
221  c_Resolution->cd();
222  c_Resolution->cd()->SetLeftMargin(0.15);
223  c_Resolution->cd()->SetBottomMargin(0.15);
224  c_Resolution->cd()->SetTicks(1,1);
225  c_Resolution->cd()->SetGrid(0,0);
226 
227  h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
228  h_play->GetYaxis()->CenterTitle();
229  h_play->GetYaxis()->SetRangeUser(0.0,0.8);
230  h_play->DrawCopy("pE");
231 
232  g_resolution->SetMarkerStyle(20);
233  g_resolution->SetMarkerColor(kGray+2);
234  g_resolution->SetMarkerSize(2.0);
235  g_resolution->Draw("pE same");
236 
237  TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
238  leg_res->SetBorderSize(0);
239  leg_res->SetFillColor(0);
240  leg_res->AddEntry(g_resolution,"e- T1044-2018","p");
241  leg_res->Draw("same");
242 
243  c_Resolution->SaveAs("../figures/HCAL_ShowerCalib/c_Resolution_ShowerCalib.eps");
244  //------------------------resolution-----------------------
245 
246  string outputfile;
247  outputfile = "/sphenix/user/xusun/TestBeam/ShowerCalibAna/T1044_2018_electron.root";
248  TFile *File_OutPut = new TFile(outputfile.c_str(),"RECREATE");
249  File_OutPut->cd();
250  g_linearity->SetName("g_linearity_2018_electron");
251  g_linearity->Write();
252  g_resolution->SetName("g_resolution_2018_electron");
253  g_resolution->Write();
254  File_OutPut->Close();
255 }