Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_resolution.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_resolution.C
1 static const int N_ENERGIES = 8;
2 static const int N_CUTLEVELS = 8;
3 #include "load_files.C"
4 void plot_resolution() {
5 
6  gStyle->SetOptFit(1111);
7  gStyle->SetOptTitle(1);
8  gStyle->SetStatW(0.3);
9 
10  int scan = 2;
11  int location = 1;
12  int energies[N_ENERGIES];
13 
14  string filenames[N_ENERGIES];
15 
16  load_files(scan,location,energies,filenames);
17 
18  TFile *fin[N_ENERGIES];
19  TTree *trees[N_ENERGIES];
20 
21  float Etotal_t, Emax_t, E3by3_t, E5by5_t, C1_t, C2_inner_t, C2_outer_t, C2_inner_new_t, C2_outer_new_t;
22  float Horz_HODO_R0_t, Horz_HODO_R1_t, Horz_HODO_R2_t, Horz_HODO_R3_t, Horz_HODO_R4_t;
23  float Horz_HODO_R5_t, Horz_HODO_R6_t, Horz_HODO_R7_t, Vert_HODO_R0_t, Vert_HODO_R1_t;
24  float Vert_HODO_R2_t, Vert_HODO_R3_t, Vert_HODO_R4_t, Vert_HODO_R5_t, Vert_HODO_R6_t, Vert_HODO_R7_t;
25  float Veto1_t, Veto2_t, Veto3_t, Veto4_t, TowerID_t,Tower_column_t,Tower_row_t;
26  float TowerE_column_0_t;
27  float TowerE_column_1_t;
28  float TowerE_column_2_t;
29  float TowerE_column_3_t;
30  float TowerE_column_4_t;
31  float TowerE_column_5_t;
32  float TowerE_column_6_t;
33  float TowerE_column_7_t;
34  float TowerE_row_0_t;
35  float TowerE_row_1_t;
36  float TowerE_row_2_t;
37  float TowerE_row_3_t;
38  float TowerE_row_4_t;
39  float TowerE_row_5_t;
40  float TowerE_row_6_t;
41  float TowerE_row_7_t;
42  float HorzTowerID_t;
43  float VertTowerID_t;
44  float SaveHoriz_TowerID0_t;
45  float SaveHoriz_TowerID1_t;
46  float SaveHoriz_TowerID2_t;
47  float SaveHoriz_TowerID3_t;
48  float SaveHoriz_TowerID4_t;
49  float SaveHoriz_TowerID5_t;
50  float SaveHoriz_TowerID6_t;
51  float SaveHoriz_TowerID7_t;
52 
53  float SaveVert_TowerID0_t;
54  float SaveVert_TowerID1_t;
55  float SaveVert_TowerID2_t;
56  float SaveVert_TowerID3_t;
57  float SaveVert_TowerID4_t;
58  float SaveVert_TowerID5_t;
59  float SaveVert_TowerID6_t;
60  float SaveVert_TowerID7_t;
61 
62  for(int i=0; i<N_ENERGIES; i++){
63  fin[i] = new TFile(filenames[i].c_str());
64 
65  ostringstream name;
66  name << "test1";
67  trees[i] = (TTree *)fin[i]->Get(name.str().c_str());
68 
69  name.str("");
70  name.clear();
71 
72  name << "tree_" << energies[i];
73  trees[i]->SetName(name.str().c_str());
74  trees[i]->SetBranchAddress("Etotal_t",&Etotal_t);
75  trees[i]->SetBranchAddress("E3by3_t",&E3by3_t);
76  trees[i]->SetBranchAddress("Emax_t",&Emax_t);
77 
78  trees[i]->SetBranchAddress("E5by5_t",&E5by5_t);
79  trees[i]->SetBranchAddress("C1_t",&C1_t);
80  trees[i]->SetBranchAddress("C2_inner_t",&C2_inner_t);
81  trees[i]->SetBranchAddress("C2_outer_t",&C2_outer_t);
82 
83  if (scan==3) {
84  trees[i]->SetBranchAddress("C2_inner_new_t",&C2_inner_new_t);
85  trees[i]->SetBranchAddress("C2_outer_new_t",&C2_outer_new_t);
86  }
87 
88  trees[i]->SetBranchAddress("Horz_HODO_R0_t",&Horz_HODO_R0_t);
89  trees[i]->SetBranchAddress("Horz_HODO_R1_t",&Horz_HODO_R1_t);
90  trees[i]->SetBranchAddress("Horz_HODO_R2_t",&Horz_HODO_R2_t);
91  trees[i]->SetBranchAddress("Horz_HODO_R3_t",&Horz_HODO_R3_t);
92  trees[i]->SetBranchAddress("Horz_HODO_R4_t",&Horz_HODO_R4_t);
93  trees[i]->SetBranchAddress("Horz_HODO_R5_t",&Horz_HODO_R5_t);
94  trees[i]->SetBranchAddress("Horz_HODO_R6_t",&Horz_HODO_R6_t);
95  trees[i]->SetBranchAddress("Horz_HODO_R7_t",&Horz_HODO_R7_t);
96  trees[i]->SetBranchAddress("Vert_HODO_R0_t",&Vert_HODO_R0_t);
97  trees[i]->SetBranchAddress("Vert_HODO_R1_t",&Vert_HODO_R1_t);
98  trees[i]->SetBranchAddress("Vert_HODO_R2_t",&Vert_HODO_R2_t);
99  trees[i]->SetBranchAddress("Vert_HODO_R3_t",&Vert_HODO_R3_t);
100  trees[i]->SetBranchAddress("Vert_HODO_R4_t",&Vert_HODO_R4_t);
101  trees[i]->SetBranchAddress("Vert_HODO_R5_t",&Vert_HODO_R5_t);
102  trees[i]->SetBranchAddress("Vert_HODO_R6_t",&Vert_HODO_R6_t);
103  trees[i]->SetBranchAddress("Vert_HODO_R7_t",&Vert_HODO_R7_t);
104  trees[i]->SetBranchAddress("Veto1_t",&Veto1_t);
105  trees[i]->SetBranchAddress("Veto2_t",&Veto2_t);
106  trees[i]->SetBranchAddress("Veto3_t",&Veto3_t);
107  trees[i]->SetBranchAddress("Veto4_t",&Veto4_t);
108  trees[i]->SetBranchAddress("TowerID_t",&TowerID_t);
109  trees[i]->SetBranchAddress("Tower_column_t",&Tower_column_t);
110  trees[i]->SetBranchAddress("Tower_row_t",&Tower_row_t);
111  trees[i]->SetBranchAddress("TowerE_column_0_t",&TowerE_column_0_t);
112  trees[i]->SetBranchAddress("TowerE_column_1_t",&TowerE_column_1_t);
113  trees[i]->SetBranchAddress("TowerE_column_2_t",&TowerE_column_2_t);
114  trees[i]->SetBranchAddress("TowerE_column_3_t",&TowerE_column_3_t);
115  trees[i]->SetBranchAddress("TowerE_column_4_t",&TowerE_column_4_t);
116  trees[i]->SetBranchAddress("TowerE_column_5_t",&TowerE_column_5_t);
117  trees[i]->SetBranchAddress("TowerE_column_6_t",&TowerE_column_6_t);
118  trees[i]->SetBranchAddress("TowerE_column_7_t",&TowerE_column_7_t);
119  trees[i]->SetBranchAddress("TowerE_row_0_t",&TowerE_row_0_t);
120  trees[i]->SetBranchAddress("TowerE_row_1_t",&TowerE_row_1_t);
121  trees[i]->SetBranchAddress("TowerE_row_2_t",&TowerE_row_2_t);
122  trees[i]->SetBranchAddress("TowerE_row_3_t",&TowerE_row_3_t);
123  trees[i]->SetBranchAddress("TowerE_row_4_t",&TowerE_row_4_t);
124  trees[i]->SetBranchAddress("TowerE_row_5_t",&TowerE_row_5_t);
125  trees[i]->SetBranchAddress("TowerE_row_6_t",&TowerE_row_6_t);
126  trees[i]->SetBranchAddress("TowerE_row_7_t",&TowerE_row_7_t);
127  trees[i]->SetBranchAddress("HorzTowerID_t",&HorzTowerID_t);
128  trees[i]->SetBranchAddress("VertTowerID_t",&VertTowerID_t);
129 
130  trees[i]->SetBranchAddress("SaveHoriz_TowerID0_t",&SaveHoriz_TowerID0_t);
131  trees[i]->SetBranchAddress("SaveHoriz_TowerID1_t",&SaveHoriz_TowerID1_t);
132  trees[i]->SetBranchAddress("SaveHoriz_TowerID2_t",&SaveHoriz_TowerID2_t);
133  trees[i]->SetBranchAddress("SaveHoriz_TowerID3_t",&SaveHoriz_TowerID3_t);
134  trees[i]->SetBranchAddress("SaveHoriz_TowerID4_t",&SaveHoriz_TowerID4_t);
135  trees[i]->SetBranchAddress("SaveHoriz_TowerID5_t",&SaveHoriz_TowerID5_t);
136  trees[i]->SetBranchAddress("SaveHoriz_TowerID6_t",&SaveHoriz_TowerID6_t);
137  trees[i]->SetBranchAddress("SaveHoriz_TowerID7_t",&SaveHoriz_TowerID7_t);
138 
139  trees[i]->SetBranchAddress("SaveVert_TowerID0_t",&SaveVert_TowerID0_t);
140  trees[i]->SetBranchAddress("SaveVert_TowerID1_t",&SaveVert_TowerID1_t);
141  trees[i]->SetBranchAddress("SaveVert_TowerID2_t",&SaveVert_TowerID2_t);
142  trees[i]->SetBranchAddress("SaveVert_TowerID3_t",&SaveVert_TowerID3_t);
143  trees[i]->SetBranchAddress("SaveVert_TowerID4_t",&SaveVert_TowerID4_t);
144  trees[i]->SetBranchAddress("SaveVert_TowerID5_t",&SaveVert_TowerID5_t);
145  trees[i]->SetBranchAddress("SaveVert_TowerID6_t",&SaveVert_TowerID6_t);
146  trees[i]->SetBranchAddress("SaveVert_TowerID7_t",&SaveVert_TowerID7_t);
147  }
148 
149  TF1 *fgauss[N_CUTLEVELS][N_ENERGIES]; //for the fit arrays
150  TH1D *hE5x5[N_CUTLEVELS][N_ENERGIES]; //for the histogram arrays
151 
152  double range_histogram[2][N_ENERGIES];
153  range_histogram[0][0] = 0;
154  range_histogram[1][0] = 5;
155  range_histogram[0][1] = 0;
156  range_histogram[1][1] = 15;
157  range_histogram[0][2] = 0;
158  range_histogram[1][2] = 15;
159  range_histogram[0][3] = 0;
160  range_histogram[1][3] = 18;
161  range_histogram[0][4] = 1;
162  range_histogram[1][4] = 22;
163  range_histogram[0][5] = 1;
164  range_histogram[1][5] = 30;
165  range_histogram[0][6] = 1;
166  range_histogram[1][6] = 22;
167  range_histogram[0][7] = 1;
168  range_histogram[1][7] = 25;
169 
170  double range_fit[2][N_ENERGIES];
171  //scan 3
172  if (scan==3) {
173  range_fit[0][0] = 1.7;
174  range_fit[1][0] = 2.7;
175  range_fit[0][1] = 3.5;
176  range_fit[1][1] = 5.5;
177  range_fit[0][2] = 5.7;
178  range_fit[1][2] = 7.8;
179  range_fit[0][3] = 8.0;
180  range_fit[1][3] = 10.8;
181  range_fit[0][4] = 11.5;
182  range_fit[1][4] = 15.8;
183  range_fit[0][5] = 15.0;
184  range_fit[1][5] = 20;
185  range_fit[0][6] = 10.2;
186  range_fit[1][6] = 13.4;
187  range_fit[0][7] = 11.6;
188  range_fit[1][7] = 16.0;
189  }
190 
191  if (scan==2) {
192  if (location==1) {
193  range_fit[0][0] = 0.7;
194  range_fit[1][0] = 1.4;
195  range_fit[0][1] = 1.5;
196  range_fit[1][1] = 2.7;
197  range_fit[0][2] = 2.5;
198  range_fit[1][2] = 4.0;
199  range_fit[0][3] = 3.3;
200  range_fit[1][3] = 4.9;
201  range_fit[0][4] = 5.2;
202  range_fit[1][4] = 7.5;
203  range_fit[0][5] = 7.3;
204  range_fit[1][5] = 9.9;
205  range_fit[0][6] = 10.9;
206  range_fit[1][6] = 15.0;
207  range_fit[0][7] = 14.7;
208  range_fit[1][7] = 19.2;
209 
210  }
211  if (location==3) {
212  range_fit[0][0] = 1.8;
213  range_fit[1][0] = 3.5;
214  range_fit[0][1] = 1.8;
215  range_fit[1][1] = 2.8;
216  range_fit[0][2] = 2.8;
217  range_fit[1][2] = 4.5;
218  range_fit[0][3] = 4.0;
219  range_fit[1][3] = 5.8;
220  range_fit[0][4] = 6.5;
221  range_fit[1][4] = 8.0;
222  range_fit[0][5] = 8.0;
223  range_fit[1][5] = 11.5;
224  range_fit[0][6] = 13.5;
225  range_fit[1][6] = 16.5;
226  range_fit[0][7] = 17.0;
227  range_fit[1][7] = 21.5;
228  }
229  }
230 
231  /*
232  cutlevel 0: no cuts
233  cutlevel 1: cherenkov
234  cutlevel 2: veto
235  cutlevel 3: hodo
236  cutlevel 4: cherenkov + veto
237  cutlevel 5: cherenkov + hodo
238  cutlevel 6: veto + hodo
239  cutlevel 7: cherenkov + veto + hodo
240  */
241  string cutlevel[N_CUTLEVELS];
242  cutlevel[0] = "no_cuts";
243  cutlevel[1] = "cherenkov";
244  cutlevel[2] = "veto";
245  cutlevel[3] = "hodoscope";
246  cutlevel[4] = "cherenkov+veto";
247  cutlevel[5] = "cherenkov+hodo";
248  cutlevel[6] = "veto+hodo";
249  cutlevel[7] = "cherenkov+veto+hodo";
250 
251  for(int i=0; i<N_ENERGIES; i++){
252 
253  cout << "starting: " << energies[i] << " GeV" << endl;
254 
255  //setup the energy distribution plots
256  int entries = trees[i]->GetEntries();
257  for(int j=0; j<N_CUTLEVELS; j++){
258  //name the plots
259  ostringstream name;
260  name << "hE5x5_" << cutlevel[j] << "_" << energies[i]<< " GeV";
261  hE5x5[j][i] = new TH1D(name.str().c_str(),name.str().c_str(),200,range_histogram[0][i],range_histogram[1][i]);
262  hE5x5[j][i]->Sumw2();
263 
264  name.str("");
265  name.clear();
266  name << "fgauss_" << j << "_" << i;
267  fgauss[j][i] = new TF1(name.str().c_str(),"gaus",range_fit[0][i],range_fit[1][i]);
268  }
269 
270  for(j = 0; j<entries; j++){//loop to fill histograms
271  trees[i]->GetEntry(j);
272 
273  bool passCherenkov = false;
274  bool passVeto = false;
275  //no cuts
276  hE5x5[0][i]->Fill(E5by5_t);
277 
278  //Cherenkov cut
279  if (scan==3) {
280  if(fabs(C2_inner_new_t + C2_outer_new_t) > 100){
281  passCherenkov = true;
282  hE5x5[1][i]->Fill(E5by5_t);
283  }
284  }
285  if (scan==2) {
286  if(fabs(C2_inner_t + C2_outer_t) > 100){
287  passCherenkov = true;
288  hE5x5[1][i]->Fill(E5by5_t);
289  }
290  }
291  //veto cut
292  if (Veto1_t < 15 || Veto2_t < 15 || Veto3_t < 15 || Veto4_t < 15) {
293  passVeto = true;
294  hE5x5[2][i]->Fill(E5by5_t);
295  }
296 
297  //vertical hodoscope cut
298  bool passVertical = false;
299  if (scan==2) {
300  if (location==1) {
301  if (fabs(Vert_HODO_R2_t) > 30||fabs(Vert_HODO_R3_t)>30) {
302  passVertical = true;
303  }
304  }
305  if (location==2) {
306  if (fabs(Vert_HODO_R2_t) > 30||fabs(Vert_HODO_R3_t)>30) {
307  passVertical = true;
308  }
309  }
310  if (location==3) {
311  if (fabs(Vert_HODO_R4_t) > 30) {
312  passVertical = true;
313  }
314  }
315  }
316  if (scan==3) {
317  if (fabs(Vert_HODO_R5_t) > 30||fabs(Vert_HODO_R6_t)>30) {
318  passVertical = true;
319  }
320  }
321  //horizontal hodoscope cut
322  bool passHorizontal = false;
323  if (scan==2) {
324  if (location==1) {
325  if (fabs(Horz_HODO_R5_t) > 30) {
326  passHorizontal = true;
327  }
328  }
329  if (location==2) {
330  if (fabs(Horz_HODO_R3_t) > 30)||fabs(Horz_HODO_R4_t) > 30) {
331  passHorizontal = true;
332  }
333  }
334  if (location==3) {
335  if (fabs(Horz_HODO_R3_t) > 30)||fabs(Horz_HODO_R4_t) > 30) {
336  passHorizontal = true;
337  }
338  }
339  }
340  if (scan==3) {
341  if (fabs(Horz_HODO_R4_t) > 30) {
342  passHorizontal = true;
343  }
344  }
345 
346  if (passHorizontal && passVertical) {
347  hE5x5[3][i]->Fill(E5by5_t);
348  }
349 
350  //veto+chrenkov cut
351  if(passCherenkov && passVeto){
352  hE5x5[4][i]->Fill(E5by5_t);
353  }
354  //cherenkov+hodo
355  if (passVertical && passHorizontal && passCherenkov) {
356  hE5x5[5][i]->Fill(E5by5_t);
357  }
358  //veto+hodo
359  if (passVeto && passVertical && passHorizontal) {
360  hE5x5[6][i]->Fill(E5by5_t);
361  }
362  //veto + hodo +cherenkov
363  if(passCherenkov && passVeto && passVertical && passHorizontal) {
364  hE5x5[7][i]->Fill(E5by5_t);
365  }
366  }
367  }
368 
369  TCanvas *c2[N_CUTLEVELS]; //canvas with energy resolutions
370 
371  //~~~
372  for(int i=0; i<N_CUTLEVELS; i++){
373  ostringstream name;
374  if (i==0)
375  name << "no cuts_";
376  if (i==1)
377  name<< "cherenkov cuts ";
378  if (i==2)
379  name<< "veto cuts ";
380  if (i==3)
381  name<<" hodoscope cuts";
382  if (i==4)
383  name<< "chrenkov+veto";
384  if (i==5)
385  name<<"cherenkov+hodoscopes";
386  if (i==6)
387  name<<"veto+hodoscopes";
388  if (i==7)
389  name<<"cherenkov + veto + hodoscopes";
390 
391  c2[i] = new TCanvas(name.str().c_str(),name.str().c_str(),1500,700);//for the histograms of energies, all the cuts
392  c2[i]->Divide(4,2);
393  //define the naming convaention on the overall canvas
394  for(int j =0; j<N_ENERGIES; j++)//loop through the enrgies on the same plot
395  {
396  c2[i]->cd(j+1);
397  gPad->SetLogy();
398  hE5x5[i][j]->GetXaxis()->SetTitle("energy in 5x5 (GeV) ");
399  hE5x5[i][j]->Draw();
400  hE5x5[i][j]->Fit(fgauss[i][j],"R");
401 
402  }//end j loop through 8 energy plots
403 
404  }//end of i loop
405 
406  //==========================================================================================
407  TCanvas *c3 = new TCanvas("c3","EMCAL Energy Resolution",1500,700);
408  c3->Divide(4,2);
409 
410  Float_t *mean [N_ENERGIES];
411 
412  for (int j=0; j<N_CUTLEVELS; j++)
413  {
414  Float_t getMEAN[N_ENERGIES];
415  Float_t getSIGMA[N_ENERGIES];
416  Float_t getSIGMAerror[N_ENERGIES];
417  Float_t getMEANerror[N_ENERGIES];
418 
419  Float_t getAllerror[N_ENERGIES];
420  Float_t y_1[N_ENERGIES]; //stores the ratio of sigma/mean
421  Float_t y_2[N_ENERGIES]; //stores the mean
422 
423  cout<<cutlevel[j]<<":"<<endl;
424  for (int i=0;i<N_ENERGIES;i++)
425  {
426  getMEAN[i]=fgauss[j][i]->GetParameter(1);
427  getSIGMA[i]=fgauss[j][i]->GetParameter(2);
428  getMEANerror[i]=fgauss[j][i]->GetParError(1);
429  getSIGMAerror[i]=fgauss[j][i]->GetParError(2);
430 
431  getAllerror[i]=(getSIGMA[i]/getMEAN[i])*TMath::Sqrt( (getSIGMAerror[i]/getSIGMA[i])*(getSIGMAerror[i]/getSIGMA[i])+(getMEANerror[i]/getMEAN[i])*(getMEANerror[i]/getMEAN[i]) );
432  y_1[i]=getSIGMA[i]/getMEAN[i];
433  y_2[i]=getMEAN[i];
434  cout<<"All error "<<i<<" is :"<<getAllerror[i]<<endl;
435  cout<<"All y "<<i<<" is :"<<y_1[i]<<endl;
436 
437  } //end of energy loop
438  c3->cd(j+1);
439  c3->SetFillColor(0);
440  c3->SetGrid();
441  c3->GetFrame()->SetBorderSize(12);
442 
443  const Int_t n1=8;
444  TF1 *f2=new TF1("f2","TMath::Sqrt([0]*[0]+([1]*[1])*TMath::Power(x,-1))",0.5,28); //convolved
445  if (scan==3) {
446  Float_t x1[n1] = {2,4,6,8,12,16,24,28};
447  }
448  if (scan==2) {
449  Float_t x1[n1] = {1,2,3,4,6,8,12,16};
450  }
451  Float_t ex1[n1] = {0,0,0,0,0,0,0,0};
452  TGraphErrors *grEtotal = new TGraphErrors(n1,x1,y_1,ex1,getAllerror);
453  grEtotal->Fit("f2","R");//fit the Etotal convolved
454  grEtotal->SetMarkerColor(kBlack);
455  grEtotal->SetMarkerStyle(20);
456  ostringstream name;
457  name<<cutlevel[j];
458  grEtotal->SetTitle(name.str().c_str());
459  grEtotal->GetXaxis()->SetTitle("Energy (GeV)");
460  grEtotal->GetXaxis()->SetRangeUser(0,20);
461  grEtotal->GetYaxis()->SetTitle("#sigma(E)/<E>");
462  grEtotal->GetHistogram()->SetMaximum(0.25); // along
463  grEtotal->GetHistogram()->SetMinimum(0.02); // Y
464  grEtotal->Draw("AP");
465  }//end of N_cutlevels
466  //c3->SaveAs("EMCal_Energy_Resolution_scan3.pdf");
467 
468  TCanvas *c6 = new TCanvas("c6","EMCAL Linearity & Energy Resolution ",1500,700);
469 
470  c6->Divide(2,1);
471  c6->cd(2);
472  const Int_t n1=8;
473  TF1 *f2=new TF1("f2","TMath::Sqrt([0]*[0]+([1]*[1])*TMath::Power(x,-1))",0.5,28); //convolved
474  if (scan==3) {
475  Float_t x1[n1] = {2,4,6,8,12,16,24,28};
476  }
477  if (scan==2) {
478  Float_t x1[n1] = {1,2,3,4,6,8,12,16};
479  }
480  Float_t ex1[n1] = {0,0,0,0,0,0,0,0};
481  TGraphErrors *grEtotal = new TGraphErrors(n1,x1,y_1,ex1,getAllerror);
482  grEtotal->Fit("f2","R");//fit the Etotal convolved
483  grEtotal->SetMarkerColor(kBlack);
484  grEtotal->SetMarkerStyle(20);
485  ostringstream name;
486  name<<cutlevel[j];
487  grEtotal->SetTitle("cherenkov+veto+hodo");
488  grEtotal->GetXaxis()->SetTitle("Energy (GeV)");
489  grEtotal->GetXaxis()->SetRangeUser(0,30);
490  grEtotal->GetYaxis()->SetTitle("#sigma(E)/<E>");
491  grEtotal->GetHistogram()->SetMaximum(0.3); // along
492  grEtotal->GetHistogram()->SetMinimum(-0.1); // Y
493  grEtotal->Draw("AP");
494  c6->cd(1);
495  //make a linearity plot
496  gStyle->SetStatX(0.9);
497  gStyle->SetStatY(0.35);
498 
499  TF1 * f_calo_l = new TF1("f_calo_l", "pol1", 0.5, 20);
500  TGraphErrors *grLinearity = new TGraphErrors(n1,x1,y_2,ex1,getAllerror);
501  grLinearity->Fit("f_calo_l","R");//fit the Etotal convolved
502  grLinearity->SetMarkerColor(kBlack);
503  grLinearity->SetMarkerStyle(20);
504  grLinearity->SetTitle("Linearity 3rd scan");
505  grLinearity->GetXaxis()->SetTitle("Input Energy (GeV)");
506  grLinearity->GetXaxis()->SetRangeUser(1,30);
507  grLinearity->GetYaxis()->SetTitle("Measured Energy (GeV)");
508  grLinearity->Draw("AP");
509  TGraphErrors *y_x = new TGraphErrors(n1,x1,x1,ex1,ex1);
510  y_x->SetLineColor(kGray);
511  y_x->Draw("same");
512 
513 }
514 
515