Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawPrototype2EMCalTower_ResolutionRecalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawPrototype2EMCalTower_ResolutionRecalib.C
1 // $Id: $
2 
11 #include <cmath>
12 #include <TFile.h>
13 #include <TString.h>
14 #include <TLine.h>
15 #include <TTree.h>
16 #include <TLatex.h>
17 #include <TGraphErrors.h>
18 #include <cassert>
19 #include "SaveCanvas.C"
20 #include "SetOKStyle.C"
21 using namespace std;
22 
23 bool mask = true;
24 
25 void
26 DrawPrototype2EMCalTower_ResolutionRecalib(
27 //
28  const TString base =
29  "/phenix/u/jinhuang/tmp/miliped_work/Production_0414/calirbated_")
30 {
31  SetOKStyle();
32  gStyle->SetOptStat(0);
33  gStyle->SetGridStyle(0);
34  if (mask)
35  gStyle->SetOptFit(0);
36  else
37  gStyle->SetOptFit(0);
38  TVirtualFitter::SetDefaultFitter("Minuit2");
39 
40  const int N = 6;
41  const double Es[] =
42  { 2, 3, 4, 8, 12, 16 };
43  const double runs[] =
44  { 2042, 2040, 2039, 2038, 2067, 2063 };
45 
46  double mean[1000];
47  double dmean[1000];
48  double res[1000];
49  double dres[1000];
50 
51  TText * t;
52  TCanvas *c1 = new TCanvas(
53  "DrawPrototype2EMCalTower_ResolutionRecalib_RunSummary"
54  + TString(mask ? "Mask" : ""),
55  "DrawPrototype2EMCalTower_ResolutionRecalib_RunSummary", 1800, 900);
56  c1->Divide(3, 2);
57  int idx = 1;
58  TPad * p;
59  for (int i = 0; i < N; ++i)
60  {
61 
62  p = (TPad *) c1->cd(idx++);
63  c1->Update();
64 
65  p->SetGridx(0);
66  p->SetGridy(0);
67 
68  const double E = Es[i];
69 
70  vector<double> v(4);
71 
72 // v = GetOnePlot(base, runs[i], cut);
73  v = RecalibratorMakeOnePlot(base, runs[i], E);
74 
75  mean[i] = v[0];
76  dmean[i] = v[1];
77  res[i] = v[2] / v[0];
78  dres[i] = v[3] / v[0];
79 
80  cout << "mean[i] = " << mean[i] << ", " //
81  << "dmean[i] = " << dmean[i] << ", " //
82  << "res[i] = " << res[i] << ", " //
83  << "dres[i] = " << dres[i] << endl;
84  }
85 
86  SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
87 
88  TGraphErrors * ge_linear = new TGraphErrors(N, Es, mean, 0, dmean);
89  TGraphErrors * ge_res = new TGraphErrors(N, Es, res, 0, dres);
90 
91  ge_linear->SetLineColor(kBlue + 3);
92  ge_linear->SetMarkerColor(kBlue + 3);
93  ge_linear->SetLineWidth(2);
94  ge_linear->SetMarkerStyle(kFullCircle);
95  ge_linear->SetMarkerSize(2);
96 
97  ge_res->SetLineColor(kBlue + 3);
98  ge_res->SetMarkerColor(kBlue + 3);
99  ge_res->SetLineWidth(2);
100  ge_res->SetMarkerStyle(kFullCircle);
101  ge_res->SetMarkerSize(2);
102 
103  TF1 * f_calo_r = new TF1("f_calo_r", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5, 25);
104  TF1 * f_calo_l = new TF1("f_calo_l", "pol2", 0.5, 25);
105 
106  f_calo_r->SetLineColor(kBlue + 3);
107  f_calo_r->SetLineWidth(3);
108 
109  TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
110  25);
111  f_calo_sim->SetParameters(2.4, 11.8);
112  f_calo_sim->SetLineWidth(3);
113  f_calo_sim->SetLineColor(kGreen + 2);
114 
115  TF1 * f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
116 // f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
117  f_calo_l_sim->SetParameters(-0., 1, -0.);
118 // f_calo_l_sim->SetLineWidth(1);
119  f_calo_l_sim->SetLineColor(kGreen + 2);
120  f_calo_l_sim->SetLineWidth(3);
121 
122  TFile * f_ref =
123  new TFile(
124  "/gpfs/mnt/gpfs02/sphenix/user/jinhuang/Prototype_2016/Production_0417_CEMC_MIP_set2_v3/DrawPrototype2EMCalTower_Resolution_Run2042_col1_row2_5x5_Valid_HODO_center_col1_row2.root");
125  assert(f_ref->IsOpen());
126 
127  TGraphErrors * ge_res_ref = (TGraphErrors *) f_ref->GetObjectChecked("Graph",
128  "TGraphErrors");
129  assert(ge_res_ref);
130  ge_res_ref->SetObjectStat(0);
131 
132  TF1 * f_calo_r_ref = (TF1 *) f_ref->GetObjectChecked("f_calo_r", "TF1");
133  assert(f_calo_r_ref);
134 
135  ge_res_ref->SetLineColor(kBlue);
136  ge_res_ref->SetMarkerColor(kBlue);
137  ge_res_ref->SetLineWidth(2);
138  ge_res_ref->SetMarkerStyle(kOpenCircle);
139  ge_res_ref->SetMarkerSize(2);
140 
141  f_calo_r_ref->SetLineColor(kBlue);
142  f_calo_r_ref->SetMarkerColor(kBlue);
143  f_calo_r_ref->SetLineWidth(2);
144 
145  TText * t;
146  TCanvas *c1 = new TCanvas(
147  Form("DrawPrototype2EMCalTower_ResolutionRecalib%.0f_", runs[0]),
148  Form("DrawPrototype2EMCalTower_ResolutionRecalib%.0f_", runs[0]), 1300,
149  600);
150  c1->Divide(2, 1);
151  int idx = 1;
152  TPad * p;
153 
154  p = (TPad *) c1->cd(idx++);
155  c1->Update();
156 // p->SetLogy();
157  p->SetGridx(0);
158  p->SetGridy(0);
159 
160  TLegend* leg = new TLegend(.15, .7, .6, .85);
161 
162  p->DrawFrame(0, 0, 25, 25,
163  Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
164  TLine * l = new TLine(0, 0, 25, 25);
165  l->SetLineColor(kGray);
166 // l->Draw();
167 
168  f_calo_l_sim->Draw("same");
169  ge_linear->Draw("p");
170 // ge_linear->Fit(f_calo_l, "RM0");
171 // f_calo_l->Draw("same");
172 
173  leg->AddEntry(ge_linear, "T-1044 data", "ep");
174  leg->AddEntry(f_calo_l_sim, "Unity", "l");
175  leg->Draw();
176 
177  p = (TPad *) c1->cd(idx++);
178  c1->Update();
179 // p->SetLogy();
180  p->SetGridx(0);
181  p->SetGridy(0);
182 
183  TH1 * hframe = p->DrawFrame(0, 0, 25, 0.2,
184  Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
185 
186  ge_res->Fit(f_calo_r, "RM0Q");
187 
188  f_calo_r_ref->Draw("same");
189  ge_res_ref->Draw("p");
190  f_calo_r->Draw("same");
191  f_calo_sim->Draw("same");
192  ge_res->Draw("p");
193 
194  TLegend* leg = new TLegend(.2, .6, .85, .9);
195  leg->AddEntry(ge_res_ref, "T-1044, MIP calib.", "ep");
196  leg->AddEntry(f_calo_r_ref,
197  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
198  f_calo_r_ref->GetParameter(0), f_calo_r_ref->GetParameter(1)), "l");
199 // leg->AddEntry(new TH1(), "", "l");
200  leg->AddEntry((TObject*)0, " ", "");
201  leg->AddEntry(ge_res, "T-1044, e-shower calib.", "ep");
202  leg->AddEntry(f_calo_r,
203  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
204  f_calo_r->GetParameter(0), f_calo_r->GetParameter(1)), "l");
205  leg->Draw();
206 
207  TLegend* leg = new TLegend(.2, .15, .85, .3);
208 
209  leg->AddEntry(f_calo_sim,
210  Form("Simulation, #DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
211  f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)), "l");
212  leg->Draw();
213 
214  hframe->SetTitle("Electron Resolution");
215 
216  SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
217 }
218 
219 vector<double>
220 RecalibratorMakeOnePlot(const TString base, const int run, const double E)
221 {
222 // /phenix/u/jinhuang/tmp/miliped_work/Production_0414/calirbated_2040.dat
223 
224  TString fname = base + Form("%d.dat", run);
225 // Prototype_e-_2_Seg0_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUMall_event.root
226 
227  cout << "Process " << fname << endl;
228 
229  fstream f(fname, ios_base::in);
230  assert(f.is_open());
231 
232  TH1 * EnergySum_LG_production = new TH1F(
233  Form("EnergySum_LG_production_%d", run),
234  Form(";Run %d: 5x5 EMCal Energy Sum (GeV);Count / bin", run), 400, 0, 30);
235  TH1 * EnergySum_LG_recalib = new TH1F(Form("EnergySum_LG_recalib_%d", run),
236  Form(";Run %d: 5x5 EMCal Energy Sum (GeV);Count / bin", run), 400, 0, 30);
237 
238  int count = 0;
239  while (!f.eof())
240  {
241  string line;
242  getline(f, line);
243 
244  if (line.length() == 0)
245  continue;
246 
247  double old_E = -1, new_E = -1;
248  char tab;
249 
250  stringstream sline(line);
251  sline >> old_E >> tab >> new_E;
252 
253 // cout <<"line -> old_E = "<<old_E<<", new_E = "<<new_E<<endl;
254  assert(old_E > 0);
255  assert(new_E > 0);
256 
257  EnergySum_LG_production->Fill(old_E);
258  EnergySum_LG_recalib->Fill(new_E);
259 
260  }
261 
262  EnergySum_LG_production->SetLineColor(kBlue + 3);
263  EnergySum_LG_production->SetLineWidth(2);
264  EnergySum_LG_production->SetMarkerStyle(kFullCircle);
265  EnergySum_LG_production->SetMarkerColor(kBlue + 3);
266  EnergySum_LG_recalib->SetLineColor(kRed + 3);
267  EnergySum_LG_recalib->SetLineWidth(2);
268  EnergySum_LG_recalib->SetMarkerStyle(kFullCircle);
269  EnergySum_LG_recalib->SetMarkerColor(kRed + 3);
270 
271  EnergySum_LG_recalib->Sumw2();
272 
273  TH1 * h = (TH1 *) EnergySum_LG_recalib->DrawClone();
274  if (!mask)
275  EnergySum_LG_production->DrawClone("same");
276 
277  TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus", h->GetMean() - 1 * h->GetRMS(),
278  h->GetMean() + 4 * h->GetRMS());
279  fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
280  h->GetMean() + 2 * h->GetRMS());
281  h->Fit(fgaus_g, "MR0N");
282 
283  TF1 * fgaus = new TF1("fgaus_LG", "gaus",
284  fgaus_g->GetParameter(1) - 2 * fgaus_g->GetParameter(2),
285  fgaus_g->GetParameter(1) + 2 * fgaus_g->GetParameter(2));
286  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
287  fgaus_g->GetParameter(2));
288  fgaus->SetLineColor(kRed);
289  fgaus->SetLineWidth(4);
290  h->Fit(fgaus, "MR");
291 
292 // h->Sumw2();
293  h->GetXaxis()->SetRangeUser(h->GetMean() - 5 * h->GetRMS(),
294  h->GetMean() + 5 * h->GetRMS());
295  if (!mask)
296  h->SetTitle(
297  Form("#DeltaE/<E> = %.1f%% @ Beam = %.0f GeV",
298  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1), E));
299  else
300  h->SetTitle(Form("Beam = %.0f GeV", E));
301 
302  assert(fgaus);
303 
304  vector<double> v(4);
305  v[0] = fgaus->GetParameter(1);
306  v[1] = fgaus->GetParError(1);
307  v[2] = fgaus->GetParameter(2);
308  v[3] = fgaus->GetParError(2);
309 
310  return v;
311 }