Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawPrototype3ShowerCalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawPrototype3ShowerCalib.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 //#include "Prototype3_DSTReader.h"
24 
25 TCut event_sel;
26 TString cuts;
27 TFile * _file0 = NULL;
28 TTree * T = NULL;
29 
30 class lin_res
31 {
32 public:
33  TString name;
34  TGraphErrors * linearity;
35  TGraphErrors * resolution;
36  TF1 * f_res;
37 };
38 
39 void
40 DrawPrototype3ShowerCalib( //
41  const TString infile =
42  "/sphenix/user/jinhuang/Prototype_2017/ShowerCalib/JointEnergyScan1_Neg.lst_EMCalCalib.root"//
43 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/ShowerCalib/2nd_tower21.lst_EMCalCalib.root"//
44 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/ShowerCalib/2nd_tower45.lst_EMCalCalib.root"//
45 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/ShowerCalib/test.lst_EMCalCalib.root"//
46 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/ShowerCalib/JointEnergyScan2_tower45_Neg.lst_EMCalCalib.root"//
47 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/ShowerCalib/JointEnergyScan3_tower21_NormBias.lst_EMCalCalib.root"//
48 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/ShowerCalib/JointEnergyScan3_tower21_LowBias.lst_EMCalCalib.root"//
49  )
50 {
51 
52  SetOKStyle();
53  gStyle->SetOptStat(0);
54  gStyle->SetOptFit(1111);
55  gStyle->SetPadGridX(0);
56  gStyle->SetPadGridY(0);
57  TVirtualFitter::SetDefaultFitter("Minuit2");
58 
59  gSystem->Load("libPrototype3.so");
60  gSystem->Load("libProto3ShowCalib.so");
61 
62 // gROOT->LoadMacro("Prototype3_DSTReader.C+");
63 
64  if (!_file0)
65  {
66  TString chian_str = infile;
67  chian_str.ReplaceAll("ALL", "*");
68 
69  TChain * t = new TChain("T");
70  const int n = t->Add(chian_str);
71 
72  cout << "Loaded " << n << " root files with " << chian_str << endl;
73  assert(n > 0);
74 
75  T = t;
76 
77  _file0 = new TFile;
78  _file0->SetName(infile);
79  }
80 
81  assert(_file0);
82 
83 // event_sel = "1";
84 // cuts = "_all_data";
85 // event_sel = "abs(truth_z )<0.25 && abs(truth_y )<0.25";
86 // cuts = "_0DegreeRot_h1_v1";
87 // event_sel = "abs(truth_z )<1.25 && abs(truth_y )<1.25";
88 // cuts = "_0DegreeRot_h5_v5";
89 // event_sel = "abs(truth_z )<1.25 && abs(truth_y )<1.25 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
90 // cuts = "_0DegreeRot_h5_v5_col2_row2";
91 // event_sel = "abs(truth_z )<.25 && abs(truth_y )<.25 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
92 // cuts = "_0DegreeRot_h1_v1_col2_row2";
93 // event_sel = "abs(truth_z + 15.9)<0.25 && abs(truth_y - .1)<0.25";
94 // cuts = "_10DegreeRot_h1_v1";
95 // event_sel = "abs(truth_z + 15.9)<0.75 && abs(truth_y - .1)<0.75";
96 // cuts = "_10DegreeRot_h3_v3";
97 // event_sel = "abs(truth_z + 37)<0.25 && abs(truth_y )<0.25";
98 // cuts = "_45DegreeRot_h1_v1";
99 // event_sel = "good_e";
100 // cuts = "_good_e";
101 // event_sel = "info.beam_mom == -8 && good_e";
102 // cuts = "_8GeV_good_e";
103 // event_sel = "info.beam_mom == -12 && good_e";
104 // cuts = "_12GeV_good_e";
105 // event_sel = "info.beam_mom == -16 && good_e";
106 // cuts = "_16GeV_good_e";
107 // event_sel = "info.beam_mom == -6";
108 // cuts = "_Neg6GeV";
109 // event_sel = "good_e && info.hodo_h==2 && info.hodo_v==3"; // Tower 21
110 // cuts = "_good_data_h2_v3";
111 // event_sel = "good_e && info.hodo_h>=2 && info.hodo_h<=3 && info.hodo_v>=1 && info.hodo_v<=3"; // Tower 21
112 // cuts = "_good_data_h23_v123";
113 // event_sel = "good_e && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=0 && info.hodo_v<=4 && abs(info.beam_mom)<16"; // Tower 21
114 // cuts = "_good_data_h12345_v01234";
115  // event_sel = "good_e && info.hodo_h==3 && info.hodo_v==2"; // Tower 45
116  // cuts = "_good_data_h3_v2";
117  // event_sel = "good_e && info.hodo_h>=2 && info.hodo_h<=3 && info.hodo_v>=1 && info.hodo_v<=3"; // Tower 45
118  // cuts = "_good_data_h23_v123";
119  // event_sel = "good_e && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=1 && info.hodo_v<=5 && abs(info.beam_mom)<16"; // Tower 45
120  // cuts = "_good_data_h12345_v12345";
121  event_sel = "good_e && info.hodo_h==3 && info.hodo_v==6"; // JointEnergyScan2_tower45_Neg
122  cuts = "_good_data_h3_v6";
123 // event_sel = "good_e && info.hodo_h==2 && info.hodo_v==2"; // JointEnergyScan2_tower45_Neg
124 // cuts = "_good_data_h2_v2";
125 // event_sel = "good_e && info.hodo_h>=1 && info.hodo_h<=2 && info.hodo_v>=1 && info.hodo_v<=3"; // JointEnergyScan2_tower45_Neg
126 // cuts = "_good_data_h12_v123";
127 // event_sel = "good_e && info.hodo_h>=0 && info.hodo_h<=4 && info.hodo_v>=0 && info.hodo_v<=4"; // JointEnergyScan2_tower45_Neg
128 // cuts = "_good_data_h01234_v01234";
129 // event_sel = "good_e && info.hodo_h==4 && info.hodo_v==4"; // JointEnergyScan3_tower21_LowBias
130 // cuts = "_good_data_h4_v4";
131  // event_sel = "good_e && info.hodo_h>=4 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=5 "; // JointEnergyScan3_tower21_LowBias
132  // cuts = "_good_data_h45_v345";
133 // event_sel = "good_e && info.hodo_h>=4 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=4 "; // JointEnergyScan3_tower21_LowBias
134 // cuts = "_good_data_h45_v34";
135 // event_sel = "good_e && info.hodo_h>=2 && info.hodo_h<=6 && info.hodo_v>=2 && info.hodo_v<=6"; // JointEnergyScan3_tower21_LowBias
136 // cuts = "_good_data_h23456_v23456";
137 // event_sel = "good_e && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=1 && info.hodo_v<=5"; // JointEnergyScan3_tower21_LowBias
138 // cuts = "_good_data_h12345_v12345_InBlock";
139 // event_sel = "good_data && info.hodo_h>=2 && info.hodo_h<=4 && info.hodo_v>=4 && info.hodo_v<=6";
140 // cuts = "_good_data_h234_v456";
141 // event_sel = "good_data && info.hodo_h>=3 && info.hodo_h<=3 && info.hodo_v>=5 && info.hodo_v<=5";
142 // cuts = "_good_data_h3_v5";
143 // event_sel =
144 // "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7";
145 // cuts = "_Neg8GeV_quality_h12345_v34567";
146 // event_sel =
147 // "info.beam_mom == -4 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
148 // cuts = "_Neg4GeV_quality_h12345_v34567_col2_row2";
149 // event_sel =
150 // "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
151 // cuts = "_Neg8GeV_quality_h12345_v34567_col2_row2";
152 // event_sel =
153 // "info.beam_mom == -12 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
154 // cuts = "_Neg12GeV_quality_h12345_v34567_col2_row2";
155 // event_sel =
156 // "info.beam_mom == -4 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
157 // cuts = "_Neg4GeV_quality_h3_v5_col2_row2";
158 // event_sel =
159 // "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
160 // cuts = "_Neg8GeV_quality_h3_v5_col2_row2";
161 // event_sel =
162 // "info.beam_mom == -12 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
163 // cuts = "_Neg12GeV_quality_h3_v5_col2_row2";
164 // event_sel = // UpTilt5
165 // "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==4 && info.hodo_v==6 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
166 // cuts = "_Neg8GeV_quality_h4_v6_col2_row2";
167 // event_sel = "good_data && info.hodo_h>=4 && info.hodo_h<=6 && info.hodo_v>=2 && info.hodo_v<=4";
168 // cuts = "_good_data_h456_v234";
169 // event_sel = "good_data && info.hodo_h==5 && info.hodo_v==3";
170 // cuts = "_good_data_h5_v3";
171 // event_sel = "info.beam_mom == -8 && good_data && info.hodo_h==5 && info.hodo_v==3";
172 // cuts = "_Neg8GeV_good_data_h5_v3";
173 // event_sel = "good_data && info.hodo_h==5 && (info.hodo_v==4 || info.hodo_v==5)";
174 // cuts = "_good_data_h5_v54";
175 
176 // event_sel = "good_data && info.hodo_h==5 && info.hodo_v==5";
177 // cuts = "_good_data_h5_v5";
178 // event_sel = "good_data && info.hodo_v==5";
179 // cuts = "_good_data_hall_v5";
180 // event_sel = "good_data && ( info.hodo_v==5 || info.hodo_v==4)";
181 // cuts = "_good_data_hall_v45";
182 
183  T->SetAlias("SimEnergyScale", "1*1");
184 // // based on /phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/UIUC21.lst_EMCalCalib.root_DrawPrototype3ShowerCalib_LineShapeData_Neg8GeV_good_data_h5_v3.svg
185 // T->SetAlias("SimEnergyScale","8.74635e+00/7.60551");
186  // // based on /phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root_DrawPrototype3ShowerCalib_LineShapeData_Neg8GeV_quality_h3_v5_col2_row2.root_DrawPrototype3ShowerCalib_SumLineShapeCompare_Electron_8GeV_QGSP_BERT_HP.root
187 // T->SetAlias("SimEnergyScale","8.88178/8.16125e+00");
188  // Tilt0.lst <-> QGSP_BERT_HP Birk 0.151
189 // T->SetAlias("SimEnergyScale","8.88178/8.01845");
190  // Tilt0.lst <-> QGSP_BERT_HP Birk 0.18
191 // T->SetAlias("SimEnergyScale","8.88178/7.94838");
192 
193  // Dummy
194 // T->SetAlias("SimEnergyScale","1*1");
195 
196  cout << "Build event selection of " << (const char *) event_sel << endl;
197 
198  T->Draw(">>EventList", event_sel);
199  TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
200  cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
201  << endl;
202 
203  T->SetEventList(elist);
204 
205 // // data stuff
206  PositionDependenceData("clus_5x5_prod.sum_E");
207 // PositionDependenceData("clus_5x5_recalib.sum_E");
208  HodoscopeCheck();
209 // LineShapeData("abs(info.C2_sum)<100", "(info.C2_sum)>600 && (info.C2_sum)<1300"); // 4 GeV
210 // LineShapeData("abs(info.C2_sum)<100", "(info.C2_sum)>500 && (info.C2_sum)<1300"); // 8 GeV
211 // LineShapeData("abs(info.C2_sum)<100", "(info.C2_sum)>200 && (info.C2_sum)<1300"); // 12 GeV
212 
213  Get_Res_linear_Summmary("sum_E");
214 // Get_Res_linear_Summmary("sum_E*2");
215 
216  // simulation stuff
217 // SimPositionCheck(-0); // 0 degree tilted
218 // LineShapeSim();
219 
220 // PositionDependenceSim("clus_5x5_prod.sum_E", -0, 5); // 0 degree tilted
221 // SimPositionCheck(-15); // 10 degree tilted
222 // PositionDependenceSim("clus_5x5_prod.sum_E", -15, 5); // 10 degree tilted
223 // SimPositionCheck(-40+3); // 45 degree tilted
224 // Get_Res_linear_Summmary_Sim();
225 }
226 
227 void
228 PositionDependenceData(TString sTOWER = "clus_5x5_prod.sum_E",
229  const double z_shift = 0, const int n_div = 1)
230 {
231  TH3F * EnergySum_LG3 =
232  new TH3F("EnergySum_LG3",
233  ";Horizontal Hodoscope (5 mm);Vertical Hodoscope (5 mm);5x5 Cluster Energy (GeV)", //
234  8, -.5, 7.5, //
235  8, -.5, 7.5, //
236  200, 1, 20);
237 
238  T->Draw(sTOWER + ":7-hodo_v:hodo_h>>EnergySum_LG3", "", "goff");
239 
240  TProfile2D * EnergySum_LG3_prof_xy = EnergySum_LG3->Project3DProfile("yx");
241  TH2 * EnergySum_LG3_yx = EnergySum_LG3->Project3D("yx");
242  TH2 * EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
243  TH2 * EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
244 
245  TGraphErrors * ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
246  TGraphErrors * ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
247 
248  TText * t;
249  TCanvas *c1 = new TCanvas(
250  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
251  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
252  1000, 960);
253  c1->Divide(2, 2);
254  int idx = 1;
255  TPad * p;
256 
257  p = (TPad *) c1->cd(idx++);
258  c1->Update();
259 // p->SetLogy();
260  p->SetGridx(0);
261  p->SetGridy(0);
262 
263  EnergySum_LG3_prof_xy->SetMinimum(0);
264 
265  EnergySum_LG3_prof_xy->Draw("colz");
266  EnergySum_LG3_prof_xy->SetTitle(
267  "Energy response;Horizontal Hodoscope (5 mm);7 - Vertical Hodoscope (5 mm)");
268 
269  p = (TPad *) c1->cd(idx++);
270  c1->Update();
271 // p->SetLogy();
272  p->SetGridx(0);
273  p->SetGridy(0);
274 
275  EnergySum_LG3_yx->SetMinimum(0);
276  EnergySum_LG3_yx->Draw("colz");
277  EnergySum_LG3_yx->SetTitle(
278  "Event counts;Horizontal Hodoscope (5 mm);7 - Vertical Hodoscope (5 mm)");
279 
280  p = (TPad *) c1->cd(idx++);
281  c1->Update();
282 // p->SetLogy();
283  p->SetGridx(0);
284  p->SetGridy(0);
285 
286  EnergySum_LG3_zx->Draw("colz");
287  EnergySum_LG3_zx->SetTitle(
288  "Position scan;Horizontal Hodoscope (5 mm);5x5 Cluster Energy (GeV)");
289 
290  ge_EnergySum_LG3_zx->SetLineWidth(2);
291  ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
292  ge_EnergySum_LG3_zx->Draw("pe");
293 
294  p = (TPad *) c1->cd(idx++);
295  c1->Update();
296 // p->SetLogy();
297  p->SetGridx(0);
298  p->SetGridy(0);
299 
300  EnergySum_LG3_zy->Draw("colz");
301  EnergySum_LG3_zy->SetTitle(
302  "Position scan;7 - Vertical Hodoscope (5 mm);5x5 Cluster Energy (GeV)");
303 
304  ge_EnergySum_LG3_zy->SetLineWidth(2);
305  ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
306  ge_EnergySum_LG3_zy->Draw("pe");
307 
308 // p = (TPad *) c1->cd(idx++);
309 // c1->Update();
311 // p->SetGridx(0);
312 // p->SetGridy(0);
313 //
314 // TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
315 //
316 // TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
317 // fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
318 // h->GetMean() + 2 * h->GetRMS());
319 // h->Fit(fgaus, "M");
320 //
321 // h->Sumw2();
322 // h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
323 // h->GetMean() + 4 * h->GetRMS());
324 // EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
325 // h->GetMean() + 4 * h->GetRMS());
326 // EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
327 // h->GetMean() + 4 * h->GetRMS());
328 //
329 // h->SetLineWidth(2);
330 // h->SetMarkerStyle(kFullCircle);
331 //
332 // h->SetTitle(
333 // Form("#DeltaE/<E> = %.1f%%",
334 // 100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
335 
336  SaveCanvas(c1,
337  TString(_file0->GetName()) + TString("_DrawPrototype3ShowerCalib_")
338  + TString(c1->GetName()), kTRUE);
339 }
340 
341 void
342 PositionDependenceSim(TString sTOWER = "clus_5x5_prod.sum_E",
343  const double z_shift = 0, const int n_div = 1)
344 {
345  TH3F * EnergySum_LG3 =
346  new TH3F("EnergySum_LG3",
347  ";Beam Horizontal Pos (cm);Beam Vertical Pos (cm);5x5 Cluster Energy (GeV)", //
348  20 * n_div, z_shift - 5, z_shift + 5, //
349  20 * n_div, -5, 5, //
350  200, 1, 20);
351 
352  T->Draw(sTOWER + ":info.truth_y:info.truth_z>>EnergySum_LG3", "", "goff");
353 
354  TProfile2D * EnergySum_LG3_xy = EnergySum_LG3->Project3DProfile("yx");
355  TH2 * EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
356  TH2 * EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
357 
358  TGraphErrors * ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
359  TGraphErrors * ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
360 
361  TText * t;
362  TCanvas *c1 = new TCanvas(
363  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
364  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
365  1000, 960);
366  c1->Divide(2, 2);
367  int idx = 1;
368  TPad * p;
369 
370  p = (TPad *) c1->cd(idx++);
371  c1->Update();
372 // p->SetLogy();
373  p->SetGridx(0);
374  p->SetGridy(0);
375 
376  EnergySum_LG3_xy->Draw("colz");
377  EnergySum_LG3_xy->SetTitle(
378  "Position scan;Beam Horizontal Pos (cm);Beam Vertical Pos (cm)");
379 
380  p = (TPad *) c1->cd(idx++);
381  c1->Update();
382 // p->SetLogy();
383  p->SetGridx(0);
384  p->SetGridy(0);
385 
386  EnergySum_LG3_zx->Draw("colz");
387  EnergySum_LG3_zx->SetTitle(
388  "Position scan;Beam Horizontal Pos (cm);5x5 Cluster Energy (GeV)");
389 
390  ge_EnergySum_LG3_zx->SetLineWidth(2);
391  ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
392  ge_EnergySum_LG3_zx->Draw("pe");
393 
394  p = (TPad *) c1->cd(idx++);
395  c1->Update();
396 // p->SetLogy();
397  p->SetGridx(0);
398  p->SetGridy(0);
399 
400  EnergySum_LG3_zy->Draw("colz");
401  EnergySum_LG3_zy->SetTitle(
402  "Position scan;Beam Vertical Pos (cm);5x5 Cluster Energy (GeV)");
403 
404  ge_EnergySum_LG3_zy->SetLineWidth(2);
405  ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
406  ge_EnergySum_LG3_zy->Draw("pe");
407 
408  p = (TPad *) c1->cd(idx++);
409  c1->Update();
410 // p->SetLogy();
411  p->SetGridx(0);
412  p->SetGridy(0);
413 
414  TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
415 
416  TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
417  fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
418  h->GetMean() + 2 * h->GetRMS());
419  h->Fit(fgaus, "M");
420 
421  h->Sumw2();
422  h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
423  h->GetMean() + 4 * h->GetRMS());
424  EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
425  h->GetMean() + 4 * h->GetRMS());
426  EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
427  h->GetMean() + 4 * h->GetRMS());
428 
429  h->SetLineWidth(2);
430  h->SetMarkerStyle(kFullCircle);
431 
432  h->SetTitle(
433  Form("#DeltaE/<E> = %.1f%%",
434  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
435 
436  SaveCanvas(c1,
437  TString(_file0->GetName()) + TString("_DrawPrototype3ShowerCalib_")
438  + TString(c1->GetName()), kTRUE);
439 }
440 
441 void
442 LineShapeData(TCut c2_h = "abs(info.C2_sum)<100", TCut c2_e =
443  "(info.C2_sum)>500 && (info.C2_sum)<1300")
444 {
445 
446  vector<double> mom;
447 
448  TText * t;
449  TCanvas *c1 = new TCanvas("LineShapeData" + cuts, "LineShapeData" + cuts,
450  1800, 950);
451 
452  c1->Divide(2, 2);
453  int idx = 1;
454  TPad * p;
455 
456  p = (TPad *) c1->cd(idx++);
457  c1->Update();
458 // p->SetGridx(1);
459 // p->SetGridy(1);
460  p->SetLogy();
461 
462  T->Draw("info.C2_sum>>h_c2_sum(300,-500,2500)");
463  h_c2_sum->SetTitle("Cherenkov Checks;Sum C2 (ADC)");
464  T->Draw("info.C2_sum>>h_c2_h(300,-500,2500)", c2_h, "same");
465  T->Draw("info.C2_sum>>h_c2_e(300,-500,2500)", c2_e, "same");
466 
467  h_c2_h->SetLineColor(kBlue);
468  h_c2_e->SetLineColor(kRed);
469 
470  h_c2_sum->SetLineWidth(2);
471  h_c2_h->SetLineWidth(2);
472  h_c2_e->SetLineWidth(2);
473 
474  p = (TPad *) c1->cd(idx++);
475  c1->Update();
476  p->SetLogy();
477 
478  T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_sum(170,-1,16)");
479  h_5x5sum_c2_sum->SetTitle(
480  "Cluster spectrum decomposition;5x5 cluster energy (GeV)");
481  T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_h(170,-1,16)", c2_h, "same");
482  T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_rej_h(170,-1,16)", !c2_h,
483  "same");
484  T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_e(170,-1,16)", c2_e, "same");
485 
486  h_5x5sum_c2_h->SetLineColor(kBlue);
487  h_5x5sum_c2_rej_h->SetLineColor(kMagenta);
488  h_5x5sum_c2_e->SetLineColor(kRed);
489 
490  h_5x5sum_c2_sum->SetLineWidth(2);
491  h_5x5sum_c2_h->SetLineWidth(2);
492  h_5x5sum_c2_rej_h->SetLineWidth(2);
493  h_5x5sum_c2_e->SetLineWidth(2);
494 
495  p = (TPad *) c1->cd(idx++);
496  c1->Update();
497  p->SetLogy();
498 
499  TH1 * h_5x5sum_c2_h2 = h_5x5sum_c2_h->DrawClone();
500  h_5x5sum_c2_h2->Sumw2();
501  h_5x5sum_c2_h2->SetMarkerColor(kBlue);
502  h_5x5sum_c2_h2->SetMarkerStyle(kFullCircle);
503  h_5x5sum_c2_h2->SetTitle(";5x5 cluster energy (GeV)");
504 
505  TH1 * h_5x5sum_c2_e2 = h_5x5sum_c2_e->DrawClone("same");
506  h_5x5sum_c2_e2->Sumw2();
507  h_5x5sum_c2_e2->SetMarkerColor(kRed);
508  h_5x5sum_c2_e2->SetMarkerStyle(kFullCircle);
509  h_5x5sum_c2_e2->SetTitle(";5x5 cluster energy (GeV)");
510 
511  p = (TPad *) c1->cd(idx++);
512  c1->Update();
513  p->SetLogy();
514 
515  TH1F * h_5x5sum_c2_h3 = h_5x5sum_c2_h->DrawClone();
516  h_5x5sum_c2_h3->SetName("h_5x5sum_c2_h3");
517  h_5x5sum_c2_h3->Sumw2();
518  h_5x5sum_c2_h3->Scale(1. / h_5x5sum_c2_h3->GetSum());
519  h_5x5sum_c2_h3->SetMarkerColor(kBlue);
520  h_5x5sum_c2_h3->SetMarkerStyle(kFullCircle);
521  h_5x5sum_c2_h3->SetTitle(";5x5 cluster energy (GeV);Probability / bin");
522 
523  SaveCanvas(c1,
524  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
525  + TString(c1->GetName()), kTRUE);
526 }
527 
528 void
529 LineShapeSim()
530 {
531  vector<double> mom;
532 
533  TText * t;
534  TCanvas *c1 = new TCanvas("LineShapeSim" + cuts, "LineShapeSim" + cuts, 1000,
535  650);
536 
537 // c1->Divide(2,2);
538  int idx = 1;
539  TPad * p;
540 
541  p = (TPad *) c1->cd(idx++);
542  c1->Update();
543  p->SetLogy();
544 
545  T->Draw("clus_5x5_prod.sum_E*SimEnergyScale>>h_5x5sum_c2_sum(170,-1,16)");
546  h_5x5sum_c2_sum->SetTitle(
547  "Cluster spectrum decomposition;5x5 cluster energy (GeV)");
548  h_5x5sum_c2_sum->SetLineWidth(2);
549 
550  SaveCanvas(c1,
551  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
552  + TString(c1->GetName()), kTRUE);
553 
554 }
555 
556 void
557 HodoscopeCheck()
558 {
559  vector<double> mom;
560 
561  TText * t;
562  TCanvas *c1 = new TCanvas("HodoscopeCheck" + cuts, "HodoscopeCheck" + cuts,
563  1300, 950);
564 
565  c1->Divide(2, 1);
566  int idx = 1;
567  TPad * p;
568 
569  p = (TPad *) c1->cd(idx++);
570  c1->Update();
571  p->SetGridx(1);
572  p->SetGridy(1);
573  p->SetLogz();
574 
575  T->Draw("clus_5x5_prod.average_col:hodo_h>>h2_h(8,-.5,7.5,160,-.5,7.5)",
576  "good_e", "colz");
577  h2_h->SetTitle(
578  "Horizontal hodoscope check;Horizontal Hodoscope;5x5 cluster mean col");
579 
580  p = (TPad *) c1->cd(idx++);
581  c1->Update();
582  p->SetGridx(1);
583  p->SetGridy(1);
584  p->SetLogz();
585 
586  T->Draw("clus_5x5_prod.average_row:hodo_v>>h2_v(8,-.5,7.5,160,-.5,7.5)",
587  "good_e", "colz");
588  h2_v->SetTitle(
589  "Vertical hodoscope check;Vertical Hodoscope;5x5 cluster mean row");
590 
591  SaveCanvas(c1,
592  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
593  + TString(c1->GetName()), kTRUE);
594 
595  return mom;
596 }
597 
598 void
599 SimPositionCheck(const double shift_z = 0)
600 {
601  vector<double> mom;
602 
603  TText * t;
604  TCanvas *c1 = new TCanvas("SimPositionCheck" + cuts,
605  "SimPositionCheck" + cuts, 1300, 950);
606 
607  c1->Divide(2, 1);
608  int idx = 1;
609  TPad * p;
610 
611  p = (TPad *) c1->cd(idx++);
612  c1->Update();
613  p->SetGridx(1);
614  p->SetGridy(1);
615  p->SetLogz();
616 
617  T->Draw(
618  Form("clus_5x5_prod.average_col:truth_z>>h2_h(30,%f,%f,160,-.5,7.5)",
619  shift_z - 1.5, shift_z + 1.5), "1", "colz");
620  h2_h->SetTitle(
621  "Horizontal hodoscope check;Horizontal beam pos;5x5 cluster mean col");
622 
623  p = (TPad *) c1->cd(idx++);
624  c1->Update();
625  p->SetGridx(1);
626  p->SetGridy(1);
627  p->SetLogz();
628 
629  T->Draw("clus_5x5_prod.average_row:truth_y>>h2_v(30,-1.5,1.5,160,-.5,7.5)",
630  "1", "colz");
631  h2_v->SetTitle(
632  "Vertical hodoscope check;Vertical beam pos;5x5 cluster mean row");
633 
634  SaveCanvas(c1,
635  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
636  + TString(c1->GetName()), kTRUE);
637 
638  return;
639 }
640 
641 void
642 Get_Res_linear_Summmary(TString e_sum = "sum_E")
643 {
644 
645  vector<double> beam_mom(GetBeamMom());
646 
647 // return;
648 
649  lin_res ges_clus_5x5_prod = GetResolution("clus_5x5_prod", beam_mom, kBlue,e_sum);
650 // lin_res ges_clus_3x3_prod = GetResolution("clus_3x3_prod", beam_mom,
651 // kBlue + 3,e_sum);
652 // lin_res ges_clus_5x5_temp = GetResolution("clus_5x5_temp", beam_mom,
653 // kRed - 2,e_sum);
654 // lin_res ges_clus_5x5_recalib = GetResolution("clus_5x5_recalib", beam_mom,
655 // kRed + 3,e_sum);
656 
657  TCanvas *c1 = new TCanvas(Form("Res_linear") + cuts,
658  Form("Res_linear") + cuts, 1300, 600);
659  c1->Divide(2, 1);
660  int idx = 1;
661  TPad * p;
662 
663  p = (TPad *) c1->cd(idx++);
664  c1->Update();
665 // p->SetLogy();
666 
667  TLegend* leg = new TLegend(.15, .7, .6, .85);
668 
669  p->DrawFrame(0, 0, 25, 25,
670  Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
671  TLine * l = new TLine(0, 0, 25, 25);
672  l->SetLineColor(kGray);
673 // l->Draw();
674 
675  TF1 * f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
676 // f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
677  f_calo_l_sim->SetParameters(-0., 1, -0.);
678 // f_calo_l_sim->SetLineWidth(1);
679  f_calo_l_sim->SetLineColor(kGreen + 2);
680  f_calo_l_sim->SetLineWidth(3);
681 
682  f_calo_l_sim->Draw("same");
683  ges_clus_5x5_prod.linearity->Draw("p");
684 // ges_clus_3x3_prod.linearity->Draw("p");
685 // ges_clus_5x5_temp.linearity->Draw("p");
686 // ges_clus_5x5_recalib.linearity->Draw("p");
687 // ge_linear->Fit(f_calo_l, "RM0");
688 // f_calo_l->Draw("same");
689 
690  leg->AddEntry(ges_clus_5x5_prod.linearity, ges_clus_5x5_prod.name, "ep");
691 // leg->AddEntry(ges_clus_3x3_prod.linearity, ges_clus_3x3_prod.name, "ep");
692 // leg->AddEntry(ges_clus_5x5_temp.linearity, ges_clus_5x5_temp.name, "ep");
693 // leg->AddEntry(ges_clus_5x5_recalib.linearity, "clus_5x5_recalib", "ep");
694  leg->AddEntry(f_calo_l_sim, "Unity", "l");
695  leg->Draw();
696 
697  p = (TPad *) c1->cd(idx++);
698  c1->Update();
699 // p->SetLogy();
700  p->SetGridx(0);
701  p->SetGridy(0);
702 
703  TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
704  25);
705  f_calo_sim->SetParameters(3.7, 12.8);
706  f_calo_sim->SetLineWidth(3);
707  f_calo_sim->SetLineColor(kGreen + 2);
708 
709  TH1 * hframe = p->DrawFrame(0, 0, 25, 0.2,
710  Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
711 
712  TLegend* leg = new TLegend(.2, .6, .85, .9);
713 
714  ges_clus_5x5_prod.f_res->Draw("same");
715  ges_clus_5x5_prod.resolution->Draw("ep");
716 // ges_clus_3x3_prod.f_res->Draw("same");
717 // ges_clus_3x3_prod.resolution->Draw("ep");
718 // ges_clus_5x5_temp.f_res->Draw("same");
719 // ges_clus_5x5_temp.resolution->Draw("ep");
720 // ges_clus_5x5_recalib.f_res->Draw("same");
721 // ges_clus_5x5_recalib.resolution->Draw("ep");
722  f_calo_sim->Draw("same");
723 
724  leg->AddEntry(ges_clus_5x5_prod.resolution, ges_clus_5x5_prod.name, "ep");
725  leg->AddEntry(ges_clus_5x5_prod.f_res,
726  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
727  ges_clus_5x5_prod.f_res->GetParameter(0),
728  ges_clus_5x5_prod.f_res->GetParameter(1)), "l");
729 
730 // leg->AddEntry(ges_clus_3x3_prod.resolution, ges_clus_3x3_prod.name, "ep");
731 // leg->AddEntry(ges_clus_3x3_prod.f_res,
732 // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
733 // ges_clus_3x3_prod.f_res->GetParameter(0),
734 // ges_clus_3x3_prod.f_res->GetParameter(1)), "l");
735 //
736 // leg->AddEntry(ges_clus_5x5_temp.resolution, ges_clus_5x5_temp.name, "ep");
737 // leg->AddEntry(ges_clus_5x5_temp.f_res,
738 // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
739 // ges_clus_5x5_temp.f_res->GetParameter(0),
740 // ges_clus_5x5_temp.f_res->GetParameter(1)), "l");
741 //
742 // leg->AddEntry(ges_clus_5x5_recalib.resolution, "clus_5x5_recalib", "ep");
743 // leg->AddEntry(ges_clus_5x5_recalib.f_res,
744 // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
745 // ges_clus_5x5_recalib.f_res->GetParameter(0),
746 // ges_clus_5x5_recalib.f_res->GetParameter(1)), "l");
747 // leg->AddEntry(new TH1(), "", "l");
748 // leg->AddEntry((TObject*) 0, " ", "");
749 
750  leg->Draw();
751 
752  TLegend* leg = new TLegend(.1, .15, .85, .25);
753 
754  leg->AddEntry(f_calo_sim,
755  Form(
756  "#splitline{Simulation w/ flat light collection}{#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}}",
757  f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)), "l");
758  leg->Draw();
759 
760  hframe->SetTitle("Electron Resolution");
761 
762  SaveCanvas(c1,
763  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
764  + TString(c1->GetName()), kTRUE);
765 }
766 
767 void
768 Get_Res_linear_Summmary_Sim()
769 {
770 
771  vector<double> beam_mom(GetBeamMom());
772 
773 // return;
774 
775  lin_res ges_clus_5x5_prod = GetResolution("clus_5x5_prod", beam_mom, kBlue,"sum_E");
776  lin_res ges_clus_3x3_prod = GetResolution("clus_3x3_prod", beam_mom,
777  kBlue + 3,"sum_E");
778 
779  TCanvas *c1 = new TCanvas(Form("Res_linear") + cuts,
780  Form("Res_linear") + cuts, 1300, 600);
781  c1->Divide(2, 1);
782  int idx = 1;
783  TPad * p;
784 
785  p = (TPad *) c1->cd(idx++);
786  c1->Update();
787 // p->SetLogy();
788 
789  TLegend* leg = new TLegend(.15, .7, .6, .85);
790 
791  p->DrawFrame(0, 0, 25, 25,
792  Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
793  TLine * l = new TLine(0, 0, 25, 25);
794  l->SetLineColor(kGray);
795 // l->Draw();
796 
797  TF1 * f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
798 // f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
799  f_calo_l_sim->SetParameters(-0., 1, -0.);
800 // f_calo_l_sim->SetLineWidth(1);
801  f_calo_l_sim->SetLineColor(kGreen + 2);
802  f_calo_l_sim->SetLineWidth(3);
803 
804  f_calo_l_sim->Draw("same");
805  ges_clus_5x5_prod.linearity->Draw("p");
806  ges_clus_3x3_prod.linearity->Draw("p");
807 // ge_linear->Fit(f_calo_l, "RM0");
808 // f_calo_l->Draw("same");
809 
810  leg->AddEntry(ges_clus_5x5_prod.linearity, ges_clus_5x5_prod.name, "ep");
811  leg->AddEntry(ges_clus_3x3_prod.linearity, ges_clus_3x3_prod.name, "ep");
812  leg->AddEntry(f_calo_l_sim, "Unity", "l");
813  leg->Draw();
814 
815  p = (TPad *) c1->cd(idx++);
816  c1->Update();
817 // p->SetLogy();
818  p->SetGridx(0);
819  p->SetGridy(0);
820 
821  TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
822  25);
823  f_calo_sim->SetParameters(2.4, 11.8);
824  f_calo_sim->SetLineWidth(3);
825  f_calo_sim->SetLineColor(kGreen + 2);
826 
827  TH1 * hframe = p->DrawFrame(0, 0, 25, 0.3,
828  Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
829 
830  TLegend* leg = new TLegend(.2, .6, .85, .9);
831 
832  ges_clus_5x5_prod.f_res->Draw("same");
833  ges_clus_5x5_prod.resolution->Draw("ep");
834  ges_clus_3x3_prod.f_res->Draw("same");
835  ges_clus_3x3_prod.resolution->Draw("ep");
836  f_calo_sim->Draw("same");
837 
838  leg->AddEntry(ges_clus_5x5_prod.resolution, ges_clus_5x5_prod.name, "ep");
839  leg->AddEntry(ges_clus_5x5_prod.f_res,
840  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
841  ges_clus_5x5_prod.f_res->GetParameter(0),
842  ges_clus_5x5_prod.f_res->GetParameter(1)), "l");
843 
844  leg->AddEntry(ges_clus_3x3_prod.resolution, ges_clus_3x3_prod.name, "ep");
845  leg->AddEntry(ges_clus_3x3_prod.f_res,
846  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
847  ges_clus_3x3_prod.f_res->GetParameter(0),
848  ges_clus_3x3_prod.f_res->GetParameter(1)), "l");
849 
850 // leg->AddEntry(new TH1(), "", "l");
851 // leg->AddEntry((TObject*) 0, " ", "");
852 
853  leg->Draw();
854 
855  TLegend* leg = new TLegend(.2, .1, .85, .3);
856 
857  leg->AddEntry(f_calo_sim,
858  Form("Prelim. Sim., #DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
859  f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)), "l");
860  leg->Draw();
861 
862  hframe->SetTitle("Electron Resolution");
863 
864  SaveCanvas(c1,
865  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
866  + TString(c1->GetName()), kTRUE);
867 }
868 
869 vector<double>
870 GetBeamMom()
871 {
872  vector<double> mom;
873 
874  TH1F * hbeam_mom = new TH1F("hbeam_mom", ";beam momentum (GeV)", 32, .5,
875  32.5);
876 
877  TText * t;
878  TCanvas *c1 = new TCanvas("GetBeamMom" + cuts, "GetBeamMom" + cuts, 1800,
879  900);
880 
881  T->Draw("abs(info.beam_mom)>>hbeam_mom");
882 
883  for (int bin = 1; bin < hbeam_mom->GetNbinsX(); bin++)
884  {
885  if (hbeam_mom->GetBinContent(bin) > 40)
886  {
887  const double momentum = hbeam_mom->GetBinCenter(bin);
888 
889  if (momentum == 1 || momentum == 2 || momentum == 3 || momentum == 4
890  || momentum == 6 || momentum == 8 || momentum == 12
891  || momentum == 16 || momentum == 24 || momentum == 32)
892  {
893  mom.push_back(momentum);
894 
895  cout << "GetBeamMom - " << momentum << " GeV for "
896  << hbeam_mom->GetBinContent(bin) << " event" << endl;
897  }
898  }
899  }
900 
901  SaveCanvas(c1,
902  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
903  + TString(c1->GetName()), kTRUE);
904 
905  return mom;
906 }
907 
908 lin_res
909 GetResolution(TString cluster_name, vector<double> beam_mom, Color_t col,TString e_sum = "sum_E")
910 {
911  vector<double> mean;
912  vector<double> mean_err;
913  vector<double> res;
914  vector<double> res_err;
915 
916  TCanvas *c1 = new TCanvas("GetResolution_LineShape_" + cluster_name + cuts,
917  "GetResolution_LineShape_" + cluster_name + cuts, 1800, 900);
918 
919  c1->Divide(4, 2);
920  int idx = 1;
921  TPad * p;
922 
923  for (int i = 0; i < beam_mom.size(); ++i)
924  {
925  p = (TPad *) c1->cd(idx++);
926  c1->Update();
927 
928  const double momemtum = beam_mom[i];
929  const TString histname = Form("hLineShape%.0fGeV_", momemtum)
930  + cluster_name;
931 
932  TH1F * h = new TH1F(histname, histname + ";Observed energy (GeV)",
933  (momemtum <= 8 ? 25 : 50), momemtum / 2, momemtum * 1.5);
934  T->Draw(cluster_name + "."+e_sum+">>" + histname,
935  Form("abs(abs(info.beam_mom)-%f)/%f<.06", momemtum, momemtum));
936 
937  h->Sumw2();
938 
939  TF1 * fgaus_g = new TF1("fgaus_LG_g_" + cluster_name, "gaus",
940  h->GetMean() - 1 * h->GetRMS(), h->GetMean() + 4 * h->GetRMS());
941  fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
942  h->GetMean() + 2 * h->GetRMS());
943  h->Fit(fgaus_g, "MR0N");
944 
945  TF1 * fgaus = new TF1("fgaus_LG_" + cluster_name, "gaus",
946  fgaus_g->GetParameter(1) - 3 * fgaus_g->GetParameter(2),
947  fgaus_g->GetParameter(1) + 3 * fgaus_g->GetParameter(2));
948  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
949  fgaus_g->GetParameter(2));
950  h->Fit(fgaus, "MR");
951 
952  h->SetLineWidth(2);
953  h->SetMarkerStyle(kFullCircle);
954  fgaus->SetLineWidth(2);
955 
956  h->SetTitle(
957  Form("%.0f GeV/c: #DeltaE/<E> = %.1f%%", momemtum,
958  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
959 
960  mean.push_back(fgaus->GetParameter(1));
961  mean_err.push_back(fgaus->GetParError(1));
962  res.push_back(fgaus->GetParameter(2) / fgaus->GetParameter(1));
963  res_err.push_back(fgaus->GetParError(2) / fgaus->GetParameter(1));
964 
965  }
966 
967  SaveCanvas(c1,
968  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_"
969  + TString(c1->GetName()), kTRUE);
970 
971  TGraphErrors * ge_linear = new TGraphErrors(beam_mom.size(), &beam_mom[0],
972  &mean[0], 0, &mean_err[0]);
973 
974  TGraphErrors * ge_res = new TGraphErrors(beam_mom.size(), &beam_mom[0],
975  &res[0], 0, &res_err[0]);
976  ge_res->GetHistogram()->SetStats(0);
977  ge_res->Print();
978 
979  lin_res ret;
980  ret.name = cluster_name;
981  ret.linearity = ge_linear;
982  ret.resolution = ge_res;
983  ret.f_res = new TF1("f_calo_r_" + cluster_name, "sqrt([0]*[0]+[1]*[1]/x)/100",
984  0.5, 25);
985  ge_res->Fit(ret.f_res, "RM0QN");
986 
987  ge_linear->SetLineColor(col);
988  ge_linear->SetMarkerColor(col);
989  ge_linear->SetLineWidth(2);
990  ge_linear->SetMarkerStyle(kFullCircle);
991  ge_linear->SetMarkerSize(1.5);
992 
993  ge_res->SetLineColor(col);
994  ge_res->SetMarkerColor(col);
995  ge_res->SetLineWidth(2);
996  ge_res->SetMarkerStyle(kFullCircle);
997  ge_res->SetMarkerSize(1.5);
998 
999  ge_res->GetHistogram()->SetStats(0);
1000 
1001  ret.f_res->SetLineColor(col);
1002  ret.f_res->SetLineWidth(3);
1003 
1004  return ret;
1005 }
1006 
1007 TGraphErrors *
1008 FitProfile(const TH2 * h2)
1009 {
1010 
1011  TProfile * p2 = h2->ProfileX();
1012 
1013  int n = 0;
1014  double x[1000];
1015  double ex[1000];
1016  double y[1000];
1017  double ey[1000];
1018 
1019  for (int i = 1; i <= h2->GetNbinsX(); i++)
1020  {
1021  TH1D * h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
1022 
1023  if (h1->GetSum() < 30)
1024  {
1025  cout << "FitProfile - ignore bin " << i << endl;
1026  continue;
1027  }
1028  else
1029  {
1030  cout << "FitProfile - fit bin " << i << endl;
1031  }
1032 
1033  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
1034  p2->GetBinError(i) * 4);
1035 
1036  TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
1037  -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
1038 
1039  fgaus.SetParameter(1, p2->GetBinContent(i));
1040  fgaus.SetParameter(2, p2->GetBinError(i));
1041 
1042  h1->Fit(&fgaus, "MQ");
1043 
1044  f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
1045  fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
1046  fgaus.GetParameter(2) / 4, 0);
1047 
1048  h1->Fit(&f2, "MQ");
1049 
1050 // new TCanvas;
1051 // h1->Draw();
1052 // fgaus.Draw("same");
1053 // break;
1054 
1055  x[n] = p2->GetBinCenter(i);
1056  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
1057  y[n] = fgaus.GetParameter(1);
1058  ey[n] = fgaus.GetParError(1);
1059 
1060 // p2->SetBinContent(i, fgaus.GetParameter(1));
1061 // p2->SetBinError(i, fgaus.GetParameter(2));
1062 
1063  n++;
1064  delete h1;
1065  }
1066 
1067  return new TGraphErrors(n, x, y, ex, ey);
1068 }
1069