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