Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawPrototype2EMCalTower_Resolution.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawPrototype2EMCalTower_Resolution.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 void
24 DrawPrototype2EMCalTower_Resolution(
25  //
26  const TString base =
27  "/gpfs/mnt/gpfs02/sphenix/user/jinhuang/Prototype_2016/Production_0417_CEMC_MIP_set2_v3/",
28  TString cut = "col1_row2_5x5_Valid_HODO_center_col1_row2")
29 {
30  SetOKStyle();
31  gStyle->SetOptStat(0);
32  gStyle->SetOptFit(0);
33  TVirtualFitter::SetDefaultFitter("Minuit2");
34 
35  const int N = 6;
36  const double Es[] =
37  { 2, 3, 4, 8, 12, 16 };
38  const double runs[] =
39  { 2042, 2040, 2039, 2038, 2067, 2063 };
40 
41  double mean[1000];
42  double dmean[1000];
43  double res[1000];
44  double dres[1000];
45 
46  for (int i = 0; i < N; ++i)
47  {
48  const double E = Es[i];
49 
50  vector<double> v(4);
51 
52  v = GetOnePlot(base, runs[i], cut);
53 
54  mean[i] = v[0];
55  dmean[i] = v[1];
56  res[i] = v[2] / v[0];
57  dres[i] = v[3] / v[0];
58 
59  cout << "mean[i] = " << mean[i] << ", " //
60  << "dmean[i] = " << dmean[i] << ", " //
61  << "res[i] = " << res[i] << ", " //
62  << "dres[i] = " << dres[i] << endl;
63  }
64 
65  TGraphErrors * ge_linear = new TGraphErrors(N, Es, mean, 0, dmean);
66  TGraphErrors * ge_res = new TGraphErrors(N, Es, res, 0, dres);
67 
68  ge_linear->SetLineColor(kBlue + 3);
69  ge_linear->SetMarkerColor(kBlue + 3);
70  ge_linear->SetLineWidth(2);
71  ge_linear->SetMarkerStyle(kFullCircle);
72  ge_linear->SetMarkerSize(1.5);
73 
74  ge_res->SetLineColor(kBlue + 3);
75  ge_res->SetMarkerColor(kBlue + 3);
76  ge_res->SetLineWidth(2);
77  ge_res->SetMarkerStyle(kFullCircle);
78  ge_res->SetMarkerSize(1.5);
79 
80  TF1 * f_calo_r = new TF1("f_calo_r", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5, 25);
81  TF1 * f_calo_l = new TF1("f_calo_l", "pol2", 0.5, 25);
82 
83  TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
84  25);
85  f_calo_sim->SetParameters(2.4, 11.8);
86  f_calo_sim->SetLineWidth(1);
87  f_calo_sim->SetLineColor(kGreen + 2);
88 
89  TF1 * f_calo_l_sim = new TF1("f_calo_l_sim", "pol2", 0.5, 25);
90  f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
91  f_calo_l_sim->SetLineWidth(1);
92  f_calo_l_sim->SetLineColor(kGreen + 2);
93 
94  TText * t;
95  TCanvas *c1 = new TCanvas(
96  Form("DrawPrototype2EMCalTower_Resolution_Run%.0f_", runs[0]) + cut,
97  Form("DrawPrototype2EMCalTower_Resolution_Run%.0f_", runs[0]) + cut, 1300,
98  600);
99  c1->Divide(2, 1);
100  int idx = 1;
101  TPad * p;
102 
103  p = (TPad *) c1->cd(idx++);
104  c1->Update();
105 // p->SetLogy();
106  p->SetGridx(0);
107  p->SetGridy(0);
108 
109  p->DrawFrame(0, 0, 25, 30,
110  Form("Linearity;Input energy (GeV);Measured Energy (GeV)", E));
111  TLine * l = new TLine(0, 0, 25, 25);
112  l->SetLineColor(kGray);
113 // l->Draw();
114 
115  ge_linear->Draw("p");
116  ge_linear->Fit(f_calo_l, "RM0");
117  f_calo_l->Draw("same");
118  f_calo_l_sim->Draw("same");
119 
120  p = (TPad *) c1->cd(idx++);
121  c1->Update();
122 // p->SetLogy();
123  p->SetGridx(0);
124  p->SetGridy(0);
125 
126  TH1 * hframe = p->DrawFrame(0, 0, 25, 0.2,
127  Form("Resolution;Input energy (GeV);#DeltaE/<E>", E));
128 
129  ge_res->Draw("p");
130  ge_res->Fit(f_calo_r, "RM0");
131  f_calo_r->Draw("same");
132  f_calo_sim->Draw("same");
133 
134  hframe->SetTitle(
135  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
136  f_calo_r->GetParameter(0), f_calo_r->GetParameter(1)));
137 
138  SaveCanvas(c1, base + TString(c1->GetName()), kTRUE);
139 }
140 
141 vector<double>
142 GetOnePlot(const TString base, const int run, TString cut)
143 {
144 // beam_00002063-0000_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUM_Energy_Sum_col1_row2_5x5_tower_center_col1_row2.png
145 // TString fname = base + TString("Prototype_") + particle + "_" + sE
146 // + "_SegALL_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUM_all_event.root";
147  TString fname =
148  base
149  + Form(
150  "/beam_0000%d-0000_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUM_Energy_Sum_",
151  run) + cut + ".root";
152 // Prototype_e-_2_Seg0_DSTReader.root_DrawPrototype2EMCalTower_EMCDistribution_SUMall_event.root
153 
154  cout << "Process " << fname << endl;
155  TFile * f = new TFile(fname);
156  assert(f->IsOpen());
157 
158  TH1 * hEnergySum = (TH1 *) f->GetObjectChecked("EnergySum_LG", "TH1");
159  assert(hEnergySum);
160  new TCanvas();
161  TH1 * h = hEnergySum->DrawClone();
162 
163  hEnergySum->Scale(1. / hEnergySum->Integral(1, -1));
164 // TF1 * fgaus = hEnergySum->GetFunction("fgaus_LG");
165 // assert(fgaus);
166 
167  TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus", h->GetMean() - 1 * h->GetRMS(),
168  h->GetMean() + 4 * h->GetRMS());
169  fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
170  h->GetMean() + 2 * h->GetRMS());
171  h->Fit(fgaus_g, "MR0N");
172 
173  TF1 * fgaus = new TF1("fgaus_LG", "gaus",
174  fgaus_g->GetParameter(1) - 2 * fgaus_g->GetParameter(2),
175  fgaus_g->GetParameter(1) + 2 * fgaus_g->GetParameter(2));
176  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
177  fgaus_g->GetParameter(2));
178  fgaus->SetLineColor(kRed);
179  fgaus->SetLineWidth(4);
180  h->Fit(fgaus, "MR");
181 
182  vector<double> v(4);
183  v[0] = fgaus->GetParameter(1);
184  v[1] = fgaus->GetParError(1);
185  v[2] = fgaus->GetParameter(2);
186  v[3] = fgaus->GetParError(2);
187 
188  return v;
189 }
190