Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
extractMIPEnergy.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file extractMIPEnergy.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 <iostream>
10 #include <fstream>
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 extractMIPEnergy()
28 {
29  string inputenergy[4] = {"5GeV","8GeV","12GeV","60GeV"};
30  string runId[4] = {"1087","0422","0571","0821"};
31  float mean_energy[4] = {5.0,8.0,12.0,60.0};
32 
33  TH2F *h_mAsymmEnergy_pion[4];
34  TH1F *h_mMIPEnergy_pion[4];
35  TFile *File_InPut[4];
36 
37  for(int i_energy = 0; i_energy < 4; ++i_energy)
38  {
39  string inputfile = Form("/sphenix/user/xusun/TestBeam/ShowerCalibAna/Proto4ShowerInfoRAW_%s.root",runId[i_energy].c_str());
40  File_InPut[i_energy] = TFile::Open(inputfile.c_str());
41  h_mAsymmEnergy_pion[i_energy] = (TH2F*)File_InPut[i_energy]->Get("h_mAsymmEnergy_pion");
42  h_mMIPEnergy_pion[i_energy] = (TH1F*)h_mAsymmEnergy_pion[i_energy]->ProjectionY()->Clone();
43  }
44 
45  TF1 *f_gaus_MIP[4];
46  float MIP_mean[4];
47  float MIP_width[4];
48  TCanvas *c_MIPEnergy = new TCanvas("c_MIPEnergy","c_MIPEnergy",1500,500);
49  c_MIPEnergy->Divide(3,1);
50  for(int i_pad = 0; i_pad < 3; ++i_pad)
51  {
52  c_MIPEnergy->cd(i_pad+1);
53  c_MIPEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
54  c_MIPEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
55  c_MIPEnergy->cd(i_pad+1)->SetTicks(1,1);
56  c_MIPEnergy->cd(i_pad+1)->SetGrid(0,0);
57  h_mMIPEnergy_pion[i_pad]->SetTitle(inputenergy[i_pad].c_str());
58  h_mMIPEnergy_pion[i_pad]->GetXaxis()->SetTitle("Tower Calibrated Energy (GeV)");
59  h_mMIPEnergy_pion[i_pad]->Draw("hE");
60  string FuncName = Form("f_gaus_MIP_%d",i_pad);
61  f_gaus_MIP[i_pad] = new TF1(FuncName.c_str(),"gaus",0,5);
62  f_gaus_MIP[i_pad]->SetParameter(0,1.0);
63  f_gaus_MIP[i_pad]->SetParameter(1,1.0);
64  f_gaus_MIP[i_pad]->SetParameter(2,0.1);
65  f_gaus_MIP[i_pad]->SetRange(0.4,0.9);
66  h_mMIPEnergy_pion[i_pad]->Fit(f_gaus_MIP[i_pad],"NQR");
67  f_gaus_MIP[i_pad]->SetLineColor(2);
68  f_gaus_MIP[i_pad]->SetLineStyle(2);
69  f_gaus_MIP[i_pad]->SetLineWidth(2);
70  f_gaus_MIP[i_pad]->Draw("l same");
71  MIP_mean[i_pad] = f_gaus_MIP[i_pad]->GetParameter(1);
72  MIP_width[i_pad] = f_gaus_MIP[i_pad]->GetParameter(2);
73  cout << "MIP for " << inputenergy[i_pad] << " is: " << MIP_mean[i_pad] << " +/- " << MIP_width[i_pad] << endl;
74  }
75 
76  int i_MIP = 1; // use 8 GeV as defualt MIP cut
77  cout << "default MIP for all energy chosen to be " << inputenergy[i_MIP] << " is: " << MIP_mean[i_MIP] << " +/- " << MIP_width[i_MIP] << endl;
78 
79  ofstream File_OutPut("MIPEnergy.txt");
80  File_OutPut << "default MIP at " << inputenergy[i_MIP] << " is: " << MIP_mean[i_MIP] << " +/- " << MIP_width[i_MIP] << endl;
81  File_OutPut.close();
82 
83  c_MIPEnergy->SaveAs("../figures/HCAL_ShowerCalib/c_MIPEnergy.eps");
84 
85 #if 0
86  TF1 *f_gaus_Reco[4];
87  float Reco_mean[4];
88  float err_mean[4];
89  float Reco_width[4];
90  float err_width[4];
91  float Reco_resolution[4];
92  float err_resolution[4];
93  TCanvas *c_RecoEnergy = new TCanvas("c_RecoEnergy","c_RecoEnergy",1500,500);
94  c_RecoEnergy->Divide(3,1);
95  for(int i_pad = 0; i_pad < 3; ++i_pad)
96  {
97  c_RecoEnergy->cd(i_pad+1);
98  c_RecoEnergy->cd(i_pad+1)->SetLeftMargin(0.15);
99  c_RecoEnergy->cd(i_pad+1)->SetBottomMargin(0.15);
100  c_RecoEnergy->cd(i_pad+1)->SetTicks(1,1);
101  c_RecoEnergy->cd(i_pad+1)->SetGrid(0,0);
102  h_mMIPEnergy_pion[i_pad]->Draw("hE");
103 
104  string FuncName = Form("f_gaus_Reco_%d",i_pad);
105  f_gaus_Reco[i_pad] = new TF1(FuncName.c_str(),"gaus",0,100);
106  f_gaus_Reco[i_pad]->SetParameter(0,1.0);
107  f_gaus_Reco[i_pad]->SetParameter(1,mean_energy[i_pad]/3.0);
108  f_gaus_Reco[i_pad]->SetParameter(2,1.0);
109  f_gaus_Reco[i_pad]->SetRange(MIP_mean[i_MIP]+2.0*MIP_width[i_MIP],0.5*mean_energy[i_pad]);
110  h_mMIPEnergy_pion[i_pad]->Fit(f_gaus_Reco[i_pad],"NQR");
111  f_gaus_Reco[i_pad]->SetLineColor(2);
112  f_gaus_Reco[i_pad]->SetLineStyle(2);
113  f_gaus_Reco[i_pad]->SetLineWidth(2);
114  f_gaus_Reco[i_pad]->Draw("l same");
115 
116  Reco_mean[i_pad] = f_gaus_Reco[i_pad]->GetParameter(1);
117  err_mean[i_pad] = f_gaus_Reco[i_pad]->GetParError(1);
118  Reco_width[i_pad] = f_gaus_Reco[i_pad]->GetParameter(2);
119  err_width[i_pad] = f_gaus_Reco[i_pad]->GetParError(2);
120  cout << "Reco for " << inputenergy[i_pad] << " is: " << Reco_mean[i_pad] << " +/- " << Reco_width[i_pad] << endl;
121 
122  Reco_resolution[i_pad] = Reco_width[i_pad]/Reco_mean[i_pad];
123  err_resolution[i_pad] = ErrDiv(Reco_width[i_pad],Reco_mean[i_pad],err_width[i_pad],err_mean[i_pad]);
124  }
125 
126  TGraphAsymmErrors *g_lieanrity = new TGraphAsymmErrors();
127  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
128  for(int i_point = 0; i_point < 3; ++i_point)
129  {
130  g_lieanrity->SetPoint(i_point,mean_energy[i_point],Reco_mean[i_point]);
131  g_lieanrity->SetPointError(i_point,0.0,0.0,err_mean[i_point],err_mean[i_point]);
132 
133  g_resolution->SetPoint(i_point,mean_energy[i_point],Reco_resolution[i_point]);
134  g_resolution->SetPointError(i_point,0.0,0.0,err_resolution[i_point],err_resolution[i_point]);
135  }
136 
137  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
138  c_Linearity->cd();
139  c_Linearity->cd()->SetLeftMargin(0.15);
140  c_Linearity->cd()->SetBottomMargin(0.15);
141  c_Linearity->cd()->SetTicks(1,1);
142  c_Linearity->cd()->SetGrid(0,0);
143 
144  TH1F *h_play = new TH1F("h_play","h_play",100,0.0,100.0);
145  for(int i_bin = 0; i_bin < 100; ++i_bin)
146  {
147  h_play->SetBinContent(i_bin+1,-10.0);
148  h_play->SetBinError(i_bin+1,1.0);
149  }
150  h_play->SetTitle("");
151  h_play->SetStats(0);
152  h_play->GetXaxis()->SetTitle("input Energy (GeV)");
153  h_play->GetXaxis()->CenterTitle();
154  h_play->GetXaxis()->SetNdivisions(505);
155  h_play->GetXaxis()->SetRangeUser(0.0,20.0);
156 
157  h_play->GetYaxis()->SetTitle("Tower Calibrated Energy (GeV)");
158  h_play->GetYaxis()->CenterTitle();
159  h_play->GetYaxis()->SetRangeUser(0.0,10.0);
160  h_play->DrawCopy("pE");
161 
162  g_lieanrity->SetMarkerStyle(24);
163  g_lieanrity->SetMarkerColor(2);
164  g_lieanrity->SetMarkerSize(1.2);
165  g_lieanrity->Draw("PE same");
166 
167  TF1 *f_pol = new TF1("f_pol","pol1",0.0,20.0);
168  f_pol->SetParameter(0,0.0);
169  f_pol->SetParameter(1,1.0);
170  f_pol->SetRange(0.0,20.0);
171  g_lieanrity->Fit(f_pol,"NR");
172  f_pol->SetLineColor(2);
173  f_pol->SetLineStyle(2);
174  f_pol->SetLineWidth(2);
175  f_pol->Draw("l same");
176 
177  TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
178  string inputresolution = "resolution.txt";
179  FILE *fres = fopen(inputresolution.c_str(), "r");
180  if(fres == NULL)
181  {
182  perror("Error opening file!");
183  }
184  else
185  {
186  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
187  char line[1024];
188  int line_counter = 0;
189  while(fgets(line,1024,fres))
190  {
191  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
192  // 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;
193  g_resolution_old->SetPoint(line_counter,x,y);
194  g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
195  line_counter++;
196  }
197  }
198  TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
199  c_Resolution->cd();
200  c_Resolution->cd()->SetLeftMargin(0.15);
201  c_Resolution->cd()->SetBottomMargin(0.15);
202  c_Resolution->cd()->SetTicks(1,1);
203  c_Resolution->cd()->SetGrid(0,0);
204 
205  h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
206  h_play->GetYaxis()->CenterTitle();
207  h_play->GetYaxis()->SetRangeUser(0.0,0.8);
208  h_play->DrawCopy("pE");
209 
210  g_resolution->SetMarkerStyle(20);
211  g_resolution->SetMarkerColor(kGray+2);
212  g_resolution->SetMarkerSize(2.0);
213  g_resolution->Draw("pE same");
214 
215  g_resolution_old->SetMarkerStyle(21);
216  g_resolution_old->SetMarkerColor(2);
217  g_resolution_old->SetMarkerSize(1.8);
218  g_resolution_old->Draw("pE same");
219 
220  TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
221  leg_res->SetBorderSize(0);
222  leg_res->SetFillColor(0);
223  leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
224  leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
225  leg_res->Draw("same");
226 #endif
227 }