Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
calLinearity_ShowerCalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file calLinearity_ShowerCalib.C
1 #include "TSystem.h"
2 #include "TFile.h"
3 #include "TH1F.h"
4 #include "TF1.h"
5 #include "TMath.h"
6 #include "TCanvas.h"
7 #include "TGraphAsymmErrors.h"
8 #include "TLine.h"
9 #include "TLegend.h"
10 #include "TString.h"
11 
12 #include <iostream>
13 #include <string>
14 #include <cstdio>
15 #include <cstring>
16 
17 using namespace std;
18 
19 float ErrorAdd(float x, float y)
20 {
21  return sqrt(x*x+y*y);
22 }
23 
24 float ErrTimes(float x, float y, float dx, float dy)
25 {
26  return x*y*ErrorAdd(dx/x,dy/y);
27 }
28 
29 float ErrDiv(float x, float y, float dx, float dy)
30 {
31  return x/y*ErrorAdd(dx/x,dy/y);
32 }
33 
34 void calLinearity_ShowerCalib()
35 {
36  TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
37  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
38  float mBeamEnergy[6] = {4.0,8.0,12.0,16.0,24.0,30.0};
39  string runID[6] = {"1241","0422","0571","1683","1666","1684"};
40 
41  TFile *File_InPut[6];
42  TH2F *h_mAsymmEnergy_pion_ShowerCalib[6];
43  TH1F *h_mEnergy_pion_ShowerCalib[6];
44 
45  for(int i_energy = 0; i_energy < 6; ++i_energy)
46  {
47  string inputfile = Form("/sphenix/user/xusun/software/data/beam/ShowerCalibAna/Proto4ShowerInfoRAW_%s.root",runID[i_energy].c_str());
48  File_InPut[i_energy] = TFile::Open(inputfile.c_str());
49  string title = Form("%1.0f GeV & #pi^{-}",mBeamEnergy[i_energy]);
50 
51  h_mAsymmEnergy_pion_ShowerCalib[i_energy] = (TH2F*)File_InPut[i_energy]->Get("h_mAsymmEnergy_pion_ShowerCalib")->Clone();
52  h_mAsymmEnergy_pion_ShowerCalib[i_energy]->SetTitle(title.c_str());
53  h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->SetTitle("Energy (GeV)");
54  h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->CenterTitle();
55  h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->SetNdivisions(505);
56  h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->SetTitle("E_{Asymm}");
57  h_mAsymmEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->CenterTitle();
58 
59  h_mEnergy_pion_ShowerCalib[i_energy] = (TH1F*)h_mAsymmEnergy_pion_ShowerCalib[i_energy]->ProjectionY()->Clone();
60  h_mEnergy_pion_ShowerCalib[i_energy]->SetStats(0);
61  h_mEnergy_pion_ShowerCalib[i_energy]->SetTitle(title.c_str());
62  h_mEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->SetTitle("Energy (GeV)");
63  h_mEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->CenterTitle();
64  h_mEnergy_pion_ShowerCalib[i_energy]->GetXaxis()->SetNdivisions(505);
65  // h_mEnergy_pion_ShowerCalib[i_energy]->SetLineColor(2);
66  h_mEnergy_pion_ShowerCalib[i_energy]->GetYaxis()->SetRangeUser(0,1.2*h_mEnergy_pion_ShowerCalib[i_energy]->GetMaximum());
67  }
68 
69  TCanvas *c_energy_asymm = new TCanvas("c_energy_asymm","c_energy_asymm",800,1200);
70  c_energy_asymm->Divide(2,3);
71  for(int i_pad = 0; i_pad < 6; ++i_pad)
72  {
73  c_energy_asymm->cd(i_pad+1);
74  c_energy_asymm->cd(i_pad+1)->SetLeftMargin(0.15);
75  c_energy_asymm->cd(i_pad+1)->SetBottomMargin(0.15);
76  c_energy_asymm->cd(i_pad+1)->SetTicks(1,1);
77  c_energy_asymm->cd(i_pad+1)->SetGrid(0,0);
78  h_mAsymmEnergy_pion_ShowerCalib[i_pad]->Draw("colz");
79  }
80  c_energy_asymm->SaveAs("./figures/sPHENIX_CollMeeting/c_energy_asymm_ShowerCalib.eps");
81 
82  TF1 *f_landau[6];
83  float mpv[6], err_mpv[6];
84  float sigma[6], err_sigma[6];
85  TCanvas *c_showercalib = new TCanvas("c_showercalib","c_showercalib",800,1200);
86  c_showercalib->Divide(2,3);
87  for(int i_pad = 0; i_pad < 6; ++i_pad)
88  {
89  c_showercalib->cd(i_pad+1);
90  c_showercalib->cd(i_pad+1)->SetLeftMargin(0.15);
91  c_showercalib->cd(i_pad+1)->SetBottomMargin(0.15);
92  c_showercalib->cd(i_pad+1)->SetTicks(1,1);
93  c_showercalib->cd(i_pad+1)->SetGrid(0,0);
94  h_mEnergy_pion_ShowerCalib[i_pad]->Draw();
95  string FuncName = Form("f_landau_%1.0fGeV",mBeamEnergy[i_pad]);
96  f_landau[i_pad] = new TF1(FuncName.c_str(),"landau",0,100);
97  f_landau[i_pad]->SetParameter(0,1.0);
98  f_landau[i_pad]->SetParameter(1,mBeamEnergy[i_pad]);
99  f_landau[i_pad]->SetParameter(2,1.0);
100  f_landau[i_pad]->SetRange(0.5*mBeamEnergy[i_pad],mBeamEnergy[i_pad]*2.0);
101  if(i_pad > 3) f_landau[i_pad]->SetRange(0.75*mBeamEnergy[i_pad],mBeamEnergy[i_pad]*1.9);
102  h_mEnergy_pion_ShowerCalib[i_pad]->Fit(f_landau[i_pad],"R");
103 
104  mpv[i_pad] = f_landau[i_pad]->GetParameter(1);
105  err_mpv[i_pad] = f_landau[i_pad]->GetParError(1);
106  g_linearity->SetPoint(i_pad,mBeamEnergy[i_pad],mpv[i_pad]);
107  g_linearity->SetPointError(i_pad,0.0,0.0,err_mpv[i_pad],err_mpv[i_pad]);
108 
109  sigma[i_pad] = f_landau[i_pad]->GetParameter(2);
110  err_sigma[i_pad] = f_landau[i_pad]->GetParError(2);
111  float val_resolution = sigma[i_pad]/mpv[i_pad];
112  float err_resolution = ErrDiv(sigma[i_pad],mpv[i_pad],err_sigma[i_pad],err_mpv[i_pad]);
113  g_resolution->SetPoint(i_pad,mBeamEnergy[i_pad],val_resolution);
114  g_resolution->SetPointError(i_pad,0.0,0.0,err_resolution,err_resolution);
115 
116  cout << "input energy is " << mBeamEnergy[i_pad] << ", calibrated energy is " << mpv[i_pad] << " +/- " << err_mpv[i_pad] << ", sigma is " << sigma[i_pad] << endl;
117  }
118  c_showercalib->SaveAs("./figures/sPHENIX_CollMeeting/c_showercalib.eps");
119 
120  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
121  c_Linearity->cd();
122  c_Linearity->cd()->SetLeftMargin(0.15);
123  c_Linearity->cd()->SetBottomMargin(0.15);
124  c_Linearity->cd()->SetTicks(1,1);
125  c_Linearity->cd()->SetGrid(0,0);
126 
127  TH1F *h_play = new TH1F("h_play","h_play",100,-0.5,99.5);
128  for(int i_bin = 0; i_bin < 100; ++i_bin)
129  {
130  h_play->SetBinContent(i_bin+1,-100.0);
131  h_play->SetBinError(i_bin+1,1.0);
132  }
133  h_play->SetTitle("");
134  h_play->SetStats(0);
135 
136  h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
137  h_play->GetXaxis()->CenterTitle();
138  h_play->GetXaxis()->SetRangeUser(0,40);
139 
140  h_play->GetYaxis()->SetTitle("Calibrated Energy (GeV)");
141  h_play->GetYaxis()->CenterTitle();
142  h_play->GetYaxis()->SetRangeUser(0,40);
143  h_play->DrawCopy("pE");
144 
145  TGraphAsymmErrors *g_linearity_old = new TGraphAsymmErrors();
146  string inputlinearity = "linearity.txt";
147  FILE *flinear = fopen(inputlinearity.c_str(), "r");
148  if(flinear == NULL)
149  {
150  perror("Error opening file!");
151  }
152  else
153  {
154  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
155  char line[1024];
156  int line_counter = 0;
157  while(fgets(line,1024,flinear))
158  {
159  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
160  // 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;
161  g_linearity_old->SetPoint(line_counter,x,y);
162  g_linearity_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
163  line_counter++;
164  }
165  }
166 
167  g_linearity_old->SetMarkerStyle(21);
168  g_linearity_old->SetMarkerColor(2);
169  g_linearity_old->SetMarkerSize(1.8);
170  g_linearity_old->Draw("pE same");
171 
172  g_linearity->SetMarkerStyle(20);
173  g_linearity->SetMarkerColor(kGray+2);
174  g_linearity->SetMarkerSize(2.0);
175  g_linearity->Draw("pE same");
176 
177  TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
178  l_unity->SetLineColor(4);
179  l_unity->SetLineStyle(2);
180  l_unity->SetLineWidth(2);
181  l_unity->Draw("l same");
182 
183  TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
184  leg_linear->SetBorderSize(0);
185  leg_linear->SetFillColor(0);
186  leg_linear->AddEntry(g_linearity_old,"#pi- T1044-2017","p");
187  leg_linear->AddEntry(g_linearity,"#pi- T1044-2018","p");
188  leg_linear->AddEntry(l_unity,"unity","l");
189  leg_linear->Draw("same");
190 
191  c_Linearity->SaveAs("./figures/sPHENIX_CollMeeting/c_Linearity.eps");
192 
193  TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
194  c_Resolution->cd();
195  c_Resolution->cd()->SetLeftMargin(0.15);
196  c_Resolution->cd()->SetBottomMargin(0.15);
197  c_Resolution->cd()->SetTicks(1,1);
198  c_Resolution->cd()->SetGrid(0,0);
199 
200  h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
201  h_play->GetXaxis()->CenterTitle();
202  h_play->GetXaxis()->SetRangeUser(0,40);
203 
204  h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
205  h_play->GetYaxis()->CenterTitle();
206  h_play->GetYaxis()->SetRangeUser(0,0.7);
207  h_play->DrawCopy("pE");
208 
209  TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
210  string inputresolution = "resolution.txt";
211  FILE *fres = fopen(inputresolution.c_str(), "r");
212  if(fres == NULL)
213  {
214  perror("Error opening file!");
215  }
216  else
217  {
218  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
219  char line[1024];
220  int line_counter = 0;
221  while(fgets(line,1024,fres))
222  {
223  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
224  // 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;
225  g_resolution_old->SetPoint(line_counter,x,y);
226  g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
227  line_counter++;
228  }
229  }
230 
231  g_resolution_old->SetMarkerStyle(21);
232  g_resolution_old->SetMarkerColor(2);
233  g_resolution_old->SetMarkerSize(1.8);
234  g_resolution_old->Draw("pE same");
235 
236  g_resolution->SetMarkerStyle(20);
237  g_resolution->SetMarkerColor(kGray+2);
238  g_resolution->SetMarkerSize(2.0);
239  g_resolution->Draw("pE same");
240 
241  TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
242  leg_res->SetBorderSize(0);
243  leg_res->SetFillColor(0);
244  leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
245  leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
246  leg_res->Draw("same");
247 
248  c_Resolution->SaveAs("./figures/sPHENIX_CollMeeting/c_Resolution.eps");
249 
250 }