Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
calLinearity.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file calLinearity.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 
11 #include <iostream>
12 #include <string>
13 #include <cstdio>
14 #include <cstring>
15 
16 using namespace std;
17 
18 float ErrorAdd(float x, float y)
19 {
20  return sqrt(x*x+y*y);
21 }
22 
23 float ErrTimes(float x, float y, float dx, float dy)
24 {
25  return x*y*ErrorAdd(dx/x,dy/y);
26 }
27 
28 float ErrDiv(float x, float y, float dx, float dy)
29 {
30  return x/y*ErrorAdd(dx/x,dy/y);
31 }
32 
33 void calLinearity()
34 {
35  gSystem->Load("libPrototype4.so");
36  gSystem->Load("libProto4_HCalShowerCalib.so");
37 
38  TGraphAsymmErrors *g_linearity = new TGraphAsymmErrors();
39  TGraphAsymmErrors *g_resolution = new TGraphAsymmErrors();
40  float mBeamEnergy[5] = {8.0,12.0,16.0,24.0,30.0};
41 
42  string inputlist = "EnergyScan.list";
43  FILE *fp = fopen(inputlist.c_str(), "r");
44  if(fp == NULL)
45  {
46  perror("Error opening file!");
47  }
48  else
49  {
50  char line[1024];
51  int line_counter = 0;
52  while(fgets(line,1024,fp))
53  {
54  line[strcspn(line, "\n")] = '\0';
55  string inputfile = line;
56  cout << "read in: " << inputfile.c_str() << endl;
57  TFile *File_HCAL = TFile::Open(inputfile.c_str());
58  TH1F *h_BeamMom = (TH1F*)File_HCAL->Get("hBeam_Mom")->Clone();
59  // float inputEnergy = TMath::Abs(h_BeamMom->GetMean()); // get input energy
60  float inputEnergy = mBeamEnergy[line_counter]; // get input energy
61 
62  TH1F *h_calib = (TH1F*)File_HCAL->Get("h_hcal_total_calib")->Clone();
63  TF1 *f_gaus = new TF1("f_gaus","gaus",0,100);
64  f_gaus->SetParameter(0,1.0);
65  f_gaus->SetParameter(1,inputEnergy);
66  f_gaus->SetParameter(2,1.0);
67  h_calib->Fit(f_gaus,"NQ");
68  float val_calib = f_gaus->GetParameter(1); // extract calibrated energy
69  float err_calib = f_gaus->GetParError(1);
70  g_linearity->SetPoint(line_counter,inputEnergy,val_calib);
71  g_linearity->SetPointError(line_counter,0.0,0.0,err_calib,err_calib);
72 
73  float val_sigma = f_gaus->GetParameter(2); // extract calibrated energy
74  float err_sigma = f_gaus->GetParError(2);
75  float val_resolution = val_sigma/val_calib;
76  float err_resolution = ErrDiv(val_sigma,val_calib,err_sigma,err_calib);
77  g_resolution->SetPoint(line_counter,inputEnergy,val_resolution);
78  g_resolution->SetPointError(line_counter,0.0,0.0,err_resolution,err_resolution);
79 
80  line_counter++;
81  File_HCAL->Close();
82  cout << "input energy is " << inputEnergy << ", calibrated energy is " << val_calib << " +/- " << err_calib << ", sigma is " << val_sigma << endl;
83  }
84  }
85 
86  TCanvas *c_Linearity = new TCanvas("c_Linearity","c_Linearity",10,10,800,800);
87  c_Linearity->cd();
88  c_Linearity->cd()->SetLeftMargin(0.15);
89  c_Linearity->cd()->SetBottomMargin(0.15);
90  c_Linearity->cd()->SetTicks(1,1);
91  c_Linearity->cd()->SetGrid(0,0);
92 
93  TH1F *h_play = new TH1F("h_play","h_play",100,-0.5,99.5);
94  for(int i_bin = 0; i_bin < 100; ++i_bin)
95  {
96  h_play->SetBinContent(i_bin+1,-100.0);
97  h_play->SetBinError(i_bin+1,1.0);
98  }
99  h_play->SetTitle("");
100  h_play->SetStats(0);
101 
102  h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
103  h_play->GetXaxis()->CenterTitle();
104  h_play->GetXaxis()->SetRangeUser(0,40);
105 
106  h_play->GetYaxis()->SetTitle("Calibrated Energy (GeV)");
107  h_play->GetYaxis()->CenterTitle();
108  h_play->GetYaxis()->SetRangeUser(0,40);
109  h_play->DrawCopy("pE");
110 
111  TGraphAsymmErrors *g_linearity_old = new TGraphAsymmErrors();
112  string inputlinearity = "linearity.txt";
113  FILE *flinear = fopen(inputlinearity.c_str(), "r");
114  if(flinear == NULL)
115  {
116  perror("Error opening file!");
117  }
118  else
119  {
120  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
121  char line[1024];
122  int line_counter = 0;
123  while(fgets(line,1024,flinear))
124  {
125  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
126  // 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;
127  g_linearity_old->SetPoint(line_counter,x,y);
128  g_linearity_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
129  line_counter++;
130  }
131  }
132 
133  g_linearity_old->SetMarkerStyle(21);
134  g_linearity_old->SetMarkerColor(2);
135  g_linearity_old->SetMarkerSize(1.8);
136  g_linearity_old->Draw("pE same");
137 
138  g_linearity->SetMarkerStyle(20);
139  g_linearity->SetMarkerColor(kGray+2);
140  g_linearity->SetMarkerSize(2.0);
141  g_linearity->Draw("pE same");
142 
143  TLine *l_unity = new TLine(1.0,1.0,39.0,39.0);
144  l_unity->SetLineColor(4);
145  l_unity->SetLineStyle(2);
146  l_unity->SetLineWidth(2);
147  l_unity->Draw("l same");
148 
149  TLegend *leg_linear = new TLegend(0.2,0.6,0.5,0.8);
150  leg_linear->SetBorderSize(0);
151  leg_linear->SetFillColor(0);
152  leg_linear->AddEntry(g_linearity_old,"#pi- T1044-2017","p");
153  leg_linear->AddEntry(g_linearity,"#pi- T1044-2018","p");
154  leg_linear->AddEntry(l_unity,"unity","l");
155  leg_linear->Draw("same");
156 
157  c_Linearity->SaveAs("/gpfs/mnt/gpfs04/sphenix/user/xusun/HCAL_BeamTest/figures/c_Linearity.eps");
158 
159  TCanvas *c_Resolution = new TCanvas("c_Resolution","c_Resolution",10,10,800,800);
160  c_Resolution->cd();
161  c_Resolution->cd()->SetLeftMargin(0.15);
162  c_Resolution->cd()->SetBottomMargin(0.15);
163  c_Resolution->cd()->SetTicks(1,1);
164  c_Resolution->cd()->SetGrid(0,0);
165 
166  h_play->GetXaxis()->SetTitle("Input Energy (GeV)");
167  h_play->GetXaxis()->CenterTitle();
168  h_play->GetXaxis()->SetRangeUser(0,40);
169 
170  h_play->GetYaxis()->SetTitle("#DeltaE/<E>");
171  h_play->GetYaxis()->CenterTitle();
172  h_play->GetYaxis()->SetRangeUser(0,0.7);
173  h_play->DrawCopy("pE");
174 
175  TGraphAsymmErrors *g_resolution_old = new TGraphAsymmErrors();
176  string inputresolution = "resolution.txt";
177  FILE *fres = fopen(inputresolution.c_str(), "r");
178  if(fres == NULL)
179  {
180  perror("Error opening file!");
181  }
182  else
183  {
184  float x, y, err_x_low, err_x_high, err_y_low, err_y_high;
185  char line[1024];
186  int line_counter = 0;
187  while(fgets(line,1024,fres))
188  {
189  sscanf(&line[0], "%f %f %f %f %f %f", &x, &y, &err_x_low, &err_x_high, &err_y_low, &err_y_high);
190  // 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;
191  g_resolution_old->SetPoint(line_counter,x,y);
192  g_resolution_old->SetPointError(line_counter,err_x_low,err_x_high,err_y_low,err_y_high);
193  line_counter++;
194  }
195  }
196 
197  g_resolution_old->SetMarkerStyle(21);
198  g_resolution_old->SetMarkerColor(2);
199  g_resolution_old->SetMarkerSize(1.8);
200  g_resolution_old->Draw("pE same");
201 
202  g_resolution->SetMarkerStyle(20);
203  g_resolution->SetMarkerColor(kGray+2);
204  g_resolution->SetMarkerSize(2.0);
205  g_resolution->Draw("pE same");
206 
207  TLegend *leg_res = new TLegend(0.55,0.6,0.85,0.8);
208  leg_res->SetBorderSize(0);
209  leg_res->SetFillColor(0);
210  leg_res->AddEntry(g_resolution_old,"#pi- T1044-2017","p");
211  leg_res->AddEntry(g_resolution,"#pi- T1044-2018","p");
212  leg_res->Draw("same");
213 
214  c_Resolution->SaveAs("/gpfs/mnt/gpfs04/sphenix/user/xusun/HCAL_BeamTest/figures/c_Resolution.eps");
215 }