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