EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LoopEvalPortableCircularCut.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LoopEvalPortableCircularCut.C
1 /*
2 > LoopEvalPortableCircularCut.C(TString detector, int print = 0, int mips = 1, int debug = 0, Double_t energyCutAggregate = 0.0, Double_t energyCut = 0.0)
3 - Creates analysis plots for the individual detector passed as an argument
4 - Processing - Eta Cuts, Manual Clustering (or elliptical cuts based on difference between generated and detected azimuth and polar angle), Recalibration, Tower energy cuts on individual towers as well as tower energy aggregated over an event
5 - Arguments
6  # detector - CEMC, EEMC, FEMC, FHCAL, HCALIN, HCALOUT
7  # print - saves plots as .png files if not 0
8  # mips - enable MIPs characterization (used with muons as particles; returns an equation that describes energy deposition by MIPs as a function of the polar angle, which can be used for a theta-dependent energy cut that eliminates MIPs)
9  # debug - prints debugging messages if not 0; Pipe the output to a file while using this
10  # energyCutAggregate - specify the value for the tower energy cut on FEMC+FHCAL
11  # energyCut - specify the value for the tower energy cut on individual towers
12 - Output file - Energy_verification_EtaCut_CircularCut_<detector>.root, and the .png plots if generated
13 */
14 
15 /*
16  Authors - Siddhant Rathi (me190003061@iiti.ac.in)
17  Sagar Joshi (ee190002054@iiti.ac.in)
18 
19  version - 2.0
20 */
21 
22 #include <iostream>
23 #include <cmath>
24 #include <eicqa_modules/EvalRootTTree.h>
25 #include <eicqa_modules/EvalHit.h>
26 #include "TMath.h"
27 #include "TStyle.h"
28 
29 R__LOAD_LIBRARY(libeicqa_modules.so)
30 
31 void LoopEvalPortableCircularCut(TString detector, int print = 0, int mips = 1, int debug = 0, Double_t energyCutAggregate = 0.0, Double_t energyCut = 0.0){
32 
33  TFile *f1 = new TFile("merged_Eval_" + detector + ".root","READ");
34 
35  TTree* T1 = (TTree*)f1->Get("T");
36  EvalRootTTree *evaltree1 = nullptr;
37 
38  gStyle->SetCanvasPreferGL(kTRUE);
39  TCanvas *c = new TCanvas();
40  c->SetTickx();
41  c->SetTicky();
42 
43  // Modifying the default plotting style
44  gStyle->SetOptTitle(0);
45  gStyle->SetOptFit(102);
46  gStyle->SetTitleXOffset(1);
47  gStyle->SetTitleYOffset(1);
48  gStyle->SetLabelSize(0.05);
49  gStyle->SetTitleXSize(0.05);
50  gStyle->SetTitleYSize(0.05);
51 
52  long double total_te = 0;
53  long double total_te_CircularCut = 0;
54  long double total_ge = 0;
55 
56  int nSlicesx; // Number of ge-axis slices in sigma_e vs ge plot
57  int nSlicesy; // Number of sigma_e-axis slices in sigma_e vs ge plot
58  Double_t eta_min, eta_max; // Eta range of detectors
59  Double_t x_radius, y_radius; // Length of semi-minor and semi-major axes of circular (elliptical) cut
60  Double_t fit_min, fit_max; // Fit range of the first slice of [(te-ge)/ge vs ge] plot
61  Double_t ttheta_gtheta_min, ttheta_gtheta_max; // Range of both axes in ttheta vs gtheta plot
62  Double_t towerCounts_max, towerCounts_CircularCut_max; // Maximum value of x axis in 1d distribution showing number of towers
63  Double_t sigma_min, sigma_max; // Range of Y-axis in sigma_e vs ge plot
64  Double_t mean_min, mean_max; // Range of Y-axis in mean_e vs ge plot
65  Double_t chi2_min, chi2_max; // Range of Y-axis in chi2_e vs ge plot
66  Double_t recalibration_firstSlice; // Recalibration factor of first slice (needed to be done manually because of low statistics)
67  Double_t te_minus_ge_by_ge_ge_min, te_minus_ge_by_ge_ge_max; // Range of y-axis in [(te-ge)/ge vs ge] plot
68  TString fitting_function; // Function fitted to slices of [(te-ge)/ge vs ge] plot
69  TString cut_text, eRes, eRes1;
70 
71  if (detector == "FEMC"){
72  x_radius = 100;//0.13
73  y_radius = 100;//0.35
74  fit_min = -0.35;
75  fit_max = 0.2;
76  eta_min = 1.3;//1.3
77  eta_max = 3.3;//3.3
78  sigma_min = 0.;
79  sigma_max = 0.5;
80  mean_min = -0.1;
81  mean_max = 0.1;
82  chi2_min = 0;
83  chi2_max = 4.00;
84  ttheta_gtheta_min = 0;
85  ttheta_gtheta_max = 0.6;
86  towerCounts_max = 400;
87  towerCounts_CircularCut_max = 200;
88  te_minus_ge_by_ge_ge_min = -0.99;
89  te_minus_ge_by_ge_ge_max = 1;
90  recalibration_firstSlice = 0.8200;
91  fitting_function = "gaus";
92  cut_text = " {1.3 < true_eta < 3.3} ";
93  eRes = "0.02 + 0.08/sqrt(x) + 0.02/x";
94  eRes1 = "0.02 + 0.08/#sqrt{ge} + 0.02/ge";
95  nSlicesx = 10;
96  nSlicesy = 1000;
97 
98  if(debug==1){
99  std::cout<<"FEMC eta cut & circular cut applied"<<"\n\n";
100  }
101  }
102 
103  else if(detector == "EEMC"){
104  x_radius = 100;//0.10
105  y_radius = 100;//0.40
106  fit_min = 0.14;
107  fit_max = 0.21;
108  eta_min = -3.5;//-3.45
109  eta_max = -1.7;//-1.75
110  sigma_min = 0;
111  sigma_max = 0.3;
112  mean_min = -0.1;
113  mean_max = 0.1;
114  chi2_min = 0;
115  chi2_max = 3.1;
116  ttheta_gtheta_min = 2.7;
117  ttheta_gtheta_max = 3.1;
118  towerCounts_max = 400;
119  towerCounts_CircularCut_max = 200;
120  te_minus_ge_by_ge_ge_min = -0.58;
121  te_minus_ge_by_ge_ge_max = 0.45;
122  recalibration_firstSlice = 0.75;
123  fitting_function = "gaus";
124  cut_text = " {-3.5 < true_eta < -1.7} ";
125  eRes = "0.01 + 0.025/sqrt(x) + 0.01/x";
126  eRes1 = "0.01 + 0.025/#sqrt{ge} + 0.01/ge";
127  nSlicesx = 10;
128  nSlicesy = 350;
129 
130  if(debug==1){
131  std::cout<<"EEMC eta cut & circular cut applied"<<"\n\n";
132  }
133  }
134 
135  else if(detector == "CEMC"){
136  x_radius = 100;//0.10
137  y_radius = 100;//0.20
138  fit_min = -0.4;
139  fit_max = 0.4;
140  eta_min = -1.5;//-1.5
141  eta_max = 1.2;//1.2
142  sigma_min = 0;
143  sigma_max = 0.5;
144  mean_min = -0.1;
145  mean_max = 0.1;
146  chi2_min = 0;
147  chi2_max = 2.7;
148  ttheta_gtheta_min = 0.45;
149  ttheta_gtheta_max = 2.6;
150  towerCounts_max = 400;
151  towerCounts_CircularCut_max = 200;
152  te_minus_ge_by_ge_ge_min = -0.99;
153  te_minus_ge_by_ge_ge_max = 1;
154  recalibration_firstSlice = 0.96;
155  fitting_function = "gaus";
156  cut_text = " {-1.5 < true_eta < 1.2} ";
157  eRes = "0.025 + 0.13/sqrt(x) + 0.02/x";
158  eRes1 = "0.025 + 0.13/#sqrt{ge} + 0.02/ge";
159  nSlicesx = 15;
160  nSlicesy = 1000;
161 
162  if(debug==1){
163  std::cout<<"CEMC eta cut & circular cut applied"<<"\n\n";
164  }
165  }
166 
167  else if(detector == "FHCAL"){
168  x_radius = 100;//0.15
169  y_radius = 100;//0.45
170  fit_min = -0.35;
171  fit_max = -0.05;
172  eta_min = 1.2;//1.2
173  eta_max = 3.5;//3.5
174  sigma_min = 0;
175  sigma_max = 0.5;
176  mean_min = 0.5;
177  mean_max =0.95;
178  chi2_min = 0;
179  chi2_max = 1.82;
180  ttheta_gtheta_min = 0;
181  ttheta_gtheta_max = 0.6;
182  towerCounts_max = 400;
183  towerCounts_CircularCut_max = 200;
184  te_minus_ge_by_ge_ge_min = -0.99;
185  te_minus_ge_by_ge_ge_max = 1;
186  recalibration_firstSlice = 1.414;
187  fitting_function = "gaus";
188  cut_text = " {1.2 < true_eta < 3.5} ";
189  eRes = "1";
190  eRes1 = "1";
191  nSlicesx = 15;
192  nSlicesy = 1000;
193 
194  if(debug==1){
195  std::cout<<"FHCAL eta cut & circular cut applied"<<"\n\n";
196  }
197  }
198 
199  else if(detector == "HCALIN"){
200  x_radius = 100;//0.15
201  y_radius = 100;//0.25
202  fit_min = -0.35;
203  fit_max = -0.05;
204  eta_min = -1.1;//-1.1
205  eta_max = 1.1;//1.1
206  sigma_min = 0;
207  sigma_max = 0.5;
208  mean_min = 0.5;
209  mean_max = 0.95;
210  chi2_min = 0;
211  chi2_max = 1.82;
212  ttheta_gtheta_min = 0;
213  ttheta_gtheta_max = 0.6;
214  towerCounts_max = 400;
215  towerCounts_CircularCut_max = 200;
216  te_minus_ge_by_ge_ge_min = -0.99;
217  te_minus_ge_by_ge_ge_max = 1;
218  recalibration_firstSlice = 1.414;
219  fitting_function = "gaus";
220  cut_text = " {-1.1 < true_eta < 1.1} ";
221  eRes = "1";
222  eRes1 = "1";
223  nSlicesx = 15;
224  nSlicesy = 1000;
225 
226  if(debug==1){
227  std::cout<<"HCALIN eta cut & circular cut applied"<<"\n\n";
228  }
229  }
230 
231  else if(detector == "HCALOUT"){
232  x_radius = 100;//0.2
233  y_radius = 100;//0.3
234  fit_min = -0.35;
235  fit_max = -0.05;
236  eta_min = -1.1;//-1.1
237  eta_max = 1.1;//1.1
238  sigma_min = 0;
239  sigma_max = 0.5;
240  mean_min = 0.5;
241  mean_max = 0.95;
242  chi2_min = 0;
243  chi2_max = 1.82;
244  ttheta_gtheta_min = 0;
245  ttheta_gtheta_max = 0.6;
246  towerCounts_max = 400;
247  towerCounts_CircularCut_max = 200;
248  te_minus_ge_by_ge_ge_min = -0.99;
249  te_minus_ge_by_ge_ge_max = 1;
250  recalibration_firstSlice = 1.414;
251  fitting_function = "gaus";
252  cut_text = " {-1.1 < true_eta < 1.1} ";
253  eRes = "1";
254  eRes1 = "1";
255  nSlicesx = 15;
256  nSlicesy = 1000;
257 
258  if(debug==1){
259  std::cout<<"HCALOUT eta cut & circular cut applied"<<"\n\n";
260  }
261  }
262 
263  else{
264  std::cout<<"Please try again.";
265  return 1;
266  }
267 
268  //Initialising histogram variables
269 
270  TH2D *te_minus_ge_by_ge_ge_EtaCut = new TH2D("te_minus_ge_by_ge_ge_EtaCut","#frac{#Delta e_{agg}}{truth e} vs truth e",200,0,30,200,-1.5,1);
271  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut = new TH2D("te_minus_ge_by_ge_ge_EtaCut_CircularCut","#frac{#Delta e_{agg}}{truth e} vs truth e",200,0,30,200,-1.5,2);
272  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated = new TH2D("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated","#frac{#Delta e_{agg}}{truth e} vs truth e",200,0,30,200,-1.5,1);
273  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp = new TH2D("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp","#frac{#Delta e_{agg}}{truth e} vs truth e",nSlicesx,0,30,nSlicesy,te_minus_ge_by_ge_ge_min,te_minus_ge_by_ge_ge_max); // histogram from which mean vs ge, sigma vs ge, and reduced_chi2 vs ge plots are derived
274  TH2D *te_aggregate_by_geta = new TH2D("te_aggregate_by_geta","",200,-4,4,200,0,30);
275  TH2D *te_aggregate_by_ge = new TH2D("te_aggregate_by_ge","",200,0,30,200,0,30);
276 
277  TH2D *te_by_ge_ge_EtaCut = new TH2D("te_by_ge_ge_EtaCut","te_{agg}/ge vs ge",200,0,30,200,-0.5,3.5);
278  TH2D *te_by_ge_ge_EtaCut_CircularCut = new TH2D("te_by_ge_ge_EtaCut_CircularCut","te_{agg}/ge vs ge",200,0,30,200,-0.5,3.5);
279  auto *mean_te_by_ge_ge_EtaCut_CircularCut = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut","Mean_{te/ge}",nSlicesx,0,30,-0.5,3.5);
280 
281  TH2D *tphi_gphi_EtaCut = new TH2D("tphi_gphi_EtaCut","tphi vs gphi",200,-4,4,200,-4,4);
282  TH2D *tphi_gphi_EtaCut_CircularCut = new TH2D("tphi_gphi_EtaCut_CircularCut","tphi vs gphi",200,-4,4,200,-4,4);
283 
284  TH2D *ttheta_gtheta_EtaCut = new TH2D("ttheta_gtheta_EtaCut","ttheta vs gtheta",200,ttheta_gtheta_min,ttheta_gtheta_max,200,ttheta_gtheta_min,ttheta_gtheta_max);
285  TH2D *ttheta_gtheta_EtaCut_CircularCut = new TH2D("ttheta_gtheta_EtaCut_CircularCut","ttheta vs gtheta",200,ttheta_gtheta_min,ttheta_gtheta_max,200,ttheta_gtheta_min,ttheta_gtheta_max);
286 
287  TH1D *counts_towerCounts_EtaCut = new TH1D("counts_towerCounts_EtaCut","n_towers",200,-1,towerCounts_max);
288  TH1D *counts_towerCounts_EtaCut_CircularCut = new TH1D("counts_towerCounts_EtaCut_CircularCut","n_towers",200,-1,towerCounts_CircularCut_max);
289 
290  TH1D *te_aggregate_EtaCut_CircularCut = new TH1D("te_aggregate_EtaCut_CircularCut","",200,0, 30);
291  TH1D *hist_geta = new TH1D("hist_geta","",200,-4,4);
292 
293  TH2D *dphi_dtheta_EtaCut = new TH2D("dphi_dtheta_EtaCut","dphi vs dtheta",200,-1,1,200,-1,1);
294  auto *mean_te_geta_EtaCut = new TProfile("mean_te_geta_EtaCut","te vs geta",200,-4,4,0,5);
295 
296  Double_t theta_min = 0;
297  if (detector == "EEMC"){
298  theta_min = 2;
299  }
300 
301  auto *mean_te_gtheta_EtaCut = new TProfile("mean_te_gtheta_EtaCut","te vs gtheta",200,theta_min,M_PI,0,5);
302 
303  T1->SetBranchAddress("DST#EvalTTree_" + detector,&evaltree1);
304 
305  for(int i=0; i<T1->GetEntries(); i++){
306 
307  T1->GetEntry(i);
308 
309  if(debug==1){
310  std::cout<<"\n\n\n------------------------------------------\nParticle: "<<i<<"\n\n";
311  std::cout<<"Initial Parameters "<<"\n";
312  }
313 
314  Double_t geta1 = evaltree1->get_geta();
315  if(debug==1){
316  std::cout<<"geta: "<<geta1<<"\n";
317  }
318 
319  if(geta1>=eta_min && geta1<=eta_max){
320 
321  if(debug==1){
322  cout<<"\ngeta cut applied (1.3, 3.3)"<<"\n\n";
323  }
324 
325  Double_t ge = evaltree1->get_ge();
326  if(debug==1){
327  std::cout<<"ge: "<<ge<<"\n";
328  }
329 
330  Double_t gphi = evaltree1->get_gphi();
331  if(debug==1){
332  std::cout<<"gphi: "<<gphi<<"\n";
333  }
334 
335  Double_t gtheta = evaltree1->get_gtheta();
336  if(debug==1){
337  std::cout<<"gtheta: "<<gtheta<<"\n";
338  }
339 
340  total_ge += ge;
341  if(debug==1){
342  std::cout<<"total_ge till now = "<<total_ge<<"\n";
343  }
344 
345  int twr_count = 0;
346  int twr_count_CircularCut = 0;
347  Double_t te_aggregate = 0;
348  Double_t te_aggregate_CircularCut = 0;
349 
350  for (int j=0; j<evaltree1->get_ntowers(); j++){
351 
352  if(debug==1){
353  std::cout<<"\n"<<detector<<" Tower: "<<j<<"\n";
354  }
355 
356  EvalTower *twr1 = evaltree1->get_tower(j);
357  if (twr1){
358  twr_count += 1;
359  if(debug==1){
360  cout<<"non-empty "<<detector<<" tower\n";
361  }
362  Double_t tphi = twr1->get_tphi();
363  if(debug==1){
364  std::cout<<"tphi: "<<tphi<<"\n";
365  }
366  Double_t ttheta = twr1->get_ttheta();
367  if(debug==1){
368  std::cout<<"ttheta: "<<ttheta<<"\n";
369  }
370  Double_t dphi = tphi - gphi;
371  if(debug==1){
372  std::cout<<"tphi-gphi: "<<dphi<<"\n";
373  }
374  Double_t dtheta = ttheta - gtheta;
375  if(debug==1){
376  std::cout<<"ttheta-gtheta: "<<dtheta<<"\n";
377  }
378  Double_t te = twr1->get_te();
379  if(debug==1){
380  std::cout<<"te: "<<te<<"\n";
381  }
382 
383  if(te > energyCut){
384  te_aggregate += te;
385 
386  dphi_dtheta_EtaCut->Fill(dtheta, dphi);
387  tphi_gphi_EtaCut->Fill(gphi, tphi, te);
388  ttheta_gtheta_EtaCut->Fill(gtheta, ttheta, te);
389 
390  if(debug==1){
391  std::cout<<"pow(dphi/y_radius,2)+pow(dtheta/x_radius,2): "<<pow(dphi/y_radius,2)+pow(dtheta/x_radius,2)<<"\n";
392  }
393 
394  if (pow(dphi/y_radius,2)+pow(dtheta/x_radius,2)<=1){
395  if(debug==1){
396  cout<<"Tower included after circular cut\n";
397  }
398  twr_count_CircularCut += 1;
399  te_aggregate_CircularCut += te;
400  tphi_gphi_EtaCut_CircularCut->Fill(gphi, tphi, te);
401  ttheta_gtheta_EtaCut_CircularCut->Fill(gtheta, ttheta, te);
402  }
403 
404  if(debug==1){
405  cout<<"te_aggregate till now = "<<te_aggregate<<"\n";
406  cout<<"te_aggregate_CircularCut till now = "<<te_aggregate_CircularCut<<"\n";
407  std::cout<<"te += "<<twr1->get_te()<<"\n";
408  }
409  }
410  }
411  }
412 
413  total_te += te_aggregate;
414  total_te_CircularCut += te_aggregate_CircularCut;
415  if(debug==1){
416  cout<<"total_te till now = "<<total_te<<"\n";
417  cout<<"total_te_CircularCut till now = "<<total_te_CircularCut<<"\n\n";
418  }
419 
420  if(te_aggregate_CircularCut > energyCutAggregate){
421 
422  te_aggregate_EtaCut_CircularCut->Fill(te_aggregate_CircularCut);
423  te_aggregate_by_ge->Fill(ge, te_aggregate_CircularCut);
424  te_aggregate_by_geta->Fill(geta1, te_aggregate_CircularCut);
425  hist_geta->Fill(geta1, te_aggregate_CircularCut);
426  te_minus_ge_by_ge_ge_EtaCut->Fill(ge, (te_aggregate-ge)/ge);
427  if(debug==1){
428  cout<<"(ge, (te_aggregate-ge)/ge): ("<<ge<<", "<<(te_aggregate-ge)/ge<<")\n";
429  }
430  te_minus_ge_by_ge_ge_EtaCut_CircularCut->Fill(ge, (te_aggregate_CircularCut-ge)/ge);
431  if(debug==1){
432  cout<<"(ge, (te_aggregate_CircularCut-ge)/ge): ("<<ge<<", "<<(te_aggregate_CircularCut-ge)/ge<<")\n";
433  }
434 
435  te_by_ge_ge_EtaCut->Fill(ge, te_aggregate/ge);
436  if(debug==1){
437  cout<<"(ge, te_aggregate/ge): ("<<ge<<", "<<te_aggregate/ge<<")\n";
438  }
439 
440  te_by_ge_ge_EtaCut_CircularCut->Fill(ge, te_aggregate_CircularCut/ge);
441  mean_te_by_ge_ge_EtaCut_CircularCut->Fill(ge, te_aggregate_CircularCut/ge);
442 
443  if(debug==1){
444  cout<<"(ge, te_aggregate_CircularCut/ge): ("<<ge<<", "<<te_aggregate_CircularCut/ge<<")\n";
445  }
446 
447  counts_towerCounts_EtaCut->Fill(twr_count);
448  mean_te_geta_EtaCut->Fill(geta1, te_aggregate_CircularCut);
449  mean_te_gtheta_EtaCut->Fill(gtheta, te_aggregate_CircularCut);
450  if(debug==1){
451  cout<<"(twr_count): ("<<twr_count<<")\n";
452  }
453 
454  counts_towerCounts_EtaCut_CircularCut->Fill(twr_count_CircularCut);
455  if(debug==1){
456  cout<<"(twr_count_CircularCut): ("<<twr_count_CircularCut<<")\n";
457  }
458  }
459  }
460  }
461 
462  if(debug==1){
463  cout<<"\neta cut if ends"<<"\n";
464  }
465 
466 
467  TString arr[nSlicesx];
468  for(int sno = 0; sno < nSlicesx; sno++){
469  arr[sno] = TString::Itoa(sno + 1, 10); // Used for naming of slices
470  }
471 
472  if(debug==1){
473  std::cout<<"\nGenerating sigma_e vs ge plots\n\n";
474  }
475 
476 
477  Double_t recalibrationArr[nSlicesx];
478 
479  recalibrationArr[0] = recalibration_firstSlice;
480  cout << "Recalibration factor for slice " << 1 << " is: " << recalibrationArr[0] << endl;
481 
482  for(int binIter = 2; binIter <= nSlicesx; binIter++){
483  recalibrationArr[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut->GetBinContent(binIter);
484  cout << "Recalibration factor for slice " << binIter << " is: " << recalibrationArr[binIter-1] << endl;
485  }
486 
487  for(int i=0; i<T1->GetEntries(); i++){
488 
489  T1->GetEntry(i);
490 
491  Double_t geta1 = evaltree1->get_geta();
492  Double_t gphi = evaltree1->get_gphi();
493  Double_t gtheta = evaltree1->get_gtheta();
494  Double_t ge = evaltree1->get_ge();
495  Double_t te_aggregate_CircularCut = 0;
496 
497  if(geta1>=eta_min && geta1<=eta_max){
498 
499 
500  for (int j=0; j<evaltree1->get_ntowers(); j++){
501 
502  EvalTower *twr1 = evaltree1->get_tower(j);
503  if (twr1){
504 
505  Double_t tphi = twr1->get_tphi();
506  Double_t ttheta = twr1->get_ttheta();
507  Double_t dphi = tphi - gphi;
508  Double_t dtheta = ttheta - gtheta;
509  Double_t te = twr1->get_te();
510 
511  if(te > energyCut){
512 
513  if(debug==1){
514  std::cout<<"pow(dphi/y_radius,2)+pow(dtheta/x_radius,2): "<<pow(dphi/y_radius,2)+pow(dtheta/x_radius,2)<<"\n";
515  }
516 
517  if (pow(dphi/y_radius,2)+pow(dtheta/x_radius,2)<=1){
518  te_aggregate_CircularCut += te;
519  }
520  }
521  }
522  }
523 
524  if(te_aggregate_CircularCut > energyCutAggregate){
525  int recalibration_factor = ceil((ge/30.0)*(Double_t)nSlicesx)- 1;
526  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Fill(ge, ((te_aggregate_CircularCut/recalibrationArr[recalibration_factor])-ge)/ge);
527  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->Fill(ge, ((te_aggregate_CircularCut/recalibrationArr[recalibration_factor])-ge)/ge);
528  if(debug==1){
529  cout<<"(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge: ("<<ge<<", "<<((te_aggregate_CircularCut/recalibrationArr[recalibration_factor])-ge)/ge<<")\n";
530  }
531  }
532  }
533  }
534 
535 
536  // Initialising fit functions
537  TF1 *fit = new TF1("fit",fitting_function);
538  TF1 *fit1 = new TF1("fit1",fitting_function,fit_min,fit_max);
539  TEllipse *el1 = new TEllipse(0,0,x_radius,y_radius);
540  TF1 *fExp = new TF1("fExp",eRes,0,30);
541  TF1 *fTrue = new TF1("fTrue","[0] + [1]/sqrt(x) + [2]/x",0,30);
542 
543  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->FitSlicesY(0, 0, -1, 0, "QN");
544  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2 = (TH2D*)gDirectory->Get("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2");
545  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1 = (TH2D*)gDirectory->Get("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1");
546  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2 = (TH2D*)gDirectory->Get("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2");
547 
548 
549  // Generating individual slices
550 
551  TH1D* slices[nSlicesx];
552 
553  for(int sno = 0; sno < nSlicesx; sno++){
554  int plusOne = sno+1;
555  TString sname = "slice " + arr[sno];
556  slices[sno] = te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->ProjectionY(sname, plusOne, plusOne);
557  }
558 
559  if(debug==1){
560  std::cout<<"\nHistogram Formatting\n\n";
561  }
562 
563  counts_towerCounts_EtaCut->GetXaxis()->SetTitle("Number of Towers");
564  counts_towerCounts_EtaCut->GetXaxis()->SetLabelSize(0.045);
565  counts_towerCounts_EtaCut->GetXaxis()->SetTitleSize(0.045);
566  counts_towerCounts_EtaCut->GetYaxis()->SetTitle("Counts");
567  counts_towerCounts_EtaCut->GetYaxis()->SetLabelSize(0.045);
568  counts_towerCounts_EtaCut->GetYaxis()->SetTitleSize(0.045);
569 
570  counts_towerCounts_EtaCut_CircularCut->GetXaxis()->SetTitle("Number of Towers");
571  counts_towerCounts_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.045);
572  counts_towerCounts_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.045);
573  counts_towerCounts_EtaCut_CircularCut->GetYaxis()->SetTitle("Counts");
574  counts_towerCounts_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.045);
575  counts_towerCounts_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.045);
576 
577  mean_te_geta_EtaCut->GetXaxis()->SetTitle("geta");
578  mean_te_geta_EtaCut->GetXaxis()->SetLabelSize(0.045);
579  mean_te_geta_EtaCut->GetXaxis()->SetTitleSize(0.045);
580  mean_te_geta_EtaCut->GetYaxis()->SetTitle("te_{agg}");
581  mean_te_geta_EtaCut->GetYaxis()->SetLabelSize(0.045);
582  mean_te_geta_EtaCut->GetYaxis()->SetTitleSize(0.045);
583 
584  mean_te_gtheta_EtaCut->GetXaxis()->SetTitle("gtheta");
585  mean_te_gtheta_EtaCut->GetXaxis()->SetLabelSize(0.045);
586  mean_te_gtheta_EtaCut->GetXaxis()->SetTitleSize(0.045);
587  mean_te_gtheta_EtaCut->GetYaxis()->SetTitle("te_{agg}");
588  mean_te_gtheta_EtaCut->GetYaxis()->SetLabelSize(0.045);
589  mean_te_gtheta_EtaCut->GetYaxis()->SetTitleSize(0.045);
590 
591  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
592  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetLabelSize(0.05);
593  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetTitleSize(0.05);
594  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
595  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetLabelSize(0.05);
596  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetTitleSize(0.05);
597 
598  te_aggregate_by_ge->GetXaxis()->SetTitle("Generated Energy (GeV)");
599  te_aggregate_by_ge->GetXaxis()->SetLabelSize(0.05);
600  te_aggregate_by_ge->GetXaxis()->SetTitleSize(0.05);
601  te_aggregate_by_ge->GetYaxis()->SetTitle("te_{agg} (GeV)");
602  te_aggregate_by_ge->GetYaxis()->SetLabelSize(0.05);
603  te_aggregate_by_ge->GetYaxis()->SetTitleSize(0.05);
604 
605  te_aggregate_by_geta->GetXaxis()->SetTitle("Generated #eta");
606  te_aggregate_by_geta->GetXaxis()->SetLabelSize(0.05);
607  te_aggregate_by_geta->GetXaxis()->SetTitleSize(0.05);
608  te_aggregate_by_geta->GetYaxis()->SetTitle("te_{agg} (GeV)");
609  te_aggregate_by_geta->GetYaxis()->SetLabelSize(0.05);
610  te_aggregate_by_geta->GetYaxis()->SetTitleSize(0.05);
611 
612  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
613  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
614  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
615  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
616  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
617  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
618 
619  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetTitle("Generated Energy (GeV)");
620  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetLabelSize(0.05);
621  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetTitleSize(0.05);
622  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
623  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetLabelSize(0.05);
624  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetTitleSize(0.05);
625 
626  te_by_ge_ge_EtaCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
627  te_by_ge_ge_EtaCut->GetXaxis()->SetLabelSize(0.05);
628  te_by_ge_ge_EtaCut->GetXaxis()->SetTitleSize(0.05);
629  te_by_ge_ge_EtaCut->GetYaxis()->SetTitle("te_{agg}/ge");
630  te_by_ge_ge_EtaCut->GetYaxis()->SetLabelSize(0.05);
631  te_by_ge_ge_EtaCut->GetYaxis()->SetTitleSize(0.05);
632 
633  te_aggregate_EtaCut_CircularCut->GetXaxis()->SetTitle("te_{agg} (GeV)");
634  te_aggregate_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
635  te_aggregate_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
636  te_aggregate_EtaCut_CircularCut->GetYaxis()->SetTitle("Counts");
637  te_aggregate_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
638  te_aggregate_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
639 
640  hist_geta->GetXaxis()->SetTitle("Generated #eta");
641  hist_geta->GetXaxis()->SetLabelSize(0.05);
642  hist_geta->GetXaxis()->SetTitleSize(0.05);
643  hist_geta->GetYaxis()->SetTitle("Counts");
644  hist_geta->GetYaxis()->SetLabelSize(0.05);
645  hist_geta->GetYaxis()->SetTitleSize(0.05);
646 
647  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
648  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
649  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
650  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("te_{agg}/ge");
651  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
652  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
653 
654  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
655  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
656  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
657  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("Mean of te_{agg}/ge");
658  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
659  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
660 
661  tphi_gphi_EtaCut->GetXaxis()->SetTitle("Generated #Phi");
662  tphi_gphi_EtaCut->GetXaxis()->SetLabelSize(0.05);
663  tphi_gphi_EtaCut->GetXaxis()->SetTitleSize(0.05);
664  tphi_gphi_EtaCut->GetYaxis()->SetTitle("Tower #Phi");
665  tphi_gphi_EtaCut->GetYaxis()->SetLabelSize(0.05);
666  tphi_gphi_EtaCut->GetYaxis()->SetTitleSize(0.05);
667 
668  tphi_gphi_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated #Phi");
669  tphi_gphi_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
670  tphi_gphi_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
671  tphi_gphi_EtaCut_CircularCut->GetYaxis()->SetTitle("Tower #Phi");
672  tphi_gphi_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
673  tphi_gphi_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
674 
675  ttheta_gtheta_EtaCut->GetXaxis()->SetTitle("Generated #theta");
676  ttheta_gtheta_EtaCut->GetXaxis()->SetLabelSize(0.05);
677  ttheta_gtheta_EtaCut->GetXaxis()->SetTitleSize(0.05);
678  ttheta_gtheta_EtaCut->GetYaxis()->SetTitle("Tower #theta");
679  ttheta_gtheta_EtaCut->GetYaxis()->SetLabelSize(0.05);
680  ttheta_gtheta_EtaCut->GetYaxis()->SetTitleSize(0.05);
681 
682  ttheta_gtheta_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated #theta");
683  ttheta_gtheta_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
684  ttheta_gtheta_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
685  ttheta_gtheta_EtaCut_CircularCut->GetYaxis()->SetTitle("Tower #theta");
686  ttheta_gtheta_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
687  ttheta_gtheta_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
688 
689  dphi_dtheta_EtaCut->GetXaxis()->SetTitle("ttheta-gtheta");
690  dphi_dtheta_EtaCut->GetXaxis()->SetLabelSize(0.05);
691  dphi_dtheta_EtaCut->GetXaxis()->SetTitleSize(0.05);
692  dphi_dtheta_EtaCut->GetYaxis()->SetTitle("tphi-gphi");
693  dphi_dtheta_EtaCut->GetYaxis()->SetLabelSize(0.05);
694  dphi_dtheta_EtaCut->GetYaxis()->SetTitleSize(0.05);
695 
696  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetTitle("Generated Energy (GeV)");
697  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetLabelSize(0.05);
698  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetTitleSize(0.05);
699  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetTitle("#sigma_{e_{agg}}");
700  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetLabelSize(0.05);
701  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetTitleSize(0.05);
702  //te_minus_ge_by_ge_ge_EtaCut_2->SetTitle("#sigma_{e_{agg}} vs true_e" + cut_text);
703 
704  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetTitle("Generated Energy (GeV)");
705  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetLabelSize(0.05);
706  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetTitleSize(0.05);
707  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetTitle("Reduced_#chi^{2}_{e_{agg}}");
708  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetLabelSize(0.05);
709  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetTitleSize(0.04);
710  //te_minus_ge_by_ge_ge_EtaCut_chi2->SetTitle("#chi^{2}_{e_{agg}} vs true_e" + cut_text);
711 
712  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetTitle("Generated Energy (GeV)");
713  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetLabelSize(0.05);
714  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetTitleSize(0.05);
715  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetTitle("Mean_{e_{agg}}");
716  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetLabelSize(0.05);
717  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetTitleSize(0.05);
718  //te_minus_ge_by_ge_ge_EtaCut_1->SetTitle("mean_{e_{agg}} vs true_e" + cut_text);
719 
720  for(int sno = 0; sno < nSlicesx; sno++){
721  slices[sno]->GetXaxis()->SetTitle("#Delta e^{agg}/ ge");
722  slices[sno]->GetXaxis()->SetLabelSize(0.05);
723  slices[sno]->GetXaxis()->SetTitleSize(0.05);
724  slices[sno]->GetYaxis()->SetTitle("Counts");
725  slices[sno]->GetYaxis()->SetLabelSize(0.05);
726  slices[sno]->GetYaxis()->SetTitleSize(0.05);
727  }
728 
729  if(debug==1){
730  std::cout<<"\nWriting Histograms to File\n";
731  }
732 
733  TFile *f = new TFile("Energy_verification_EtaCut_CircularCut_" + detector + ".root","RECREATE");
734 
735  f->GetList()->Add(mean_te_geta_EtaCut);
736  f->GetList()->Add(mean_te_gtheta_EtaCut);
737  f->GetList()->Add(counts_towerCounts_EtaCut);
738  f->GetList()->Add(counts_towerCounts_EtaCut_CircularCut);
739  f->GetList()->Add(te_by_ge_ge_EtaCut);
740  f->GetList()->Add(te_by_ge_ge_EtaCut_CircularCut);
741  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut);
742  f->GetList()->Add(tphi_gphi_EtaCut);
743  f->GetList()->Add(tphi_gphi_EtaCut_CircularCut);
744  f->GetList()->Add(ttheta_gtheta_EtaCut);
745  f->GetList()->Add(ttheta_gtheta_EtaCut_CircularCut);
746  f->GetList()->Add(dphi_dtheta_EtaCut);
747  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut);
748  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut);
749  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated);
750  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2);
751  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1);
752  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2);
753  f->GetList()->Add(te_aggregate_EtaCut_CircularCut);
754  f->GetList()->Add(hist_geta);
755  f->GetList()->Add(te_aggregate_by_geta);
756  f->GetList()->Add(te_aggregate_by_ge);
757 
758  for(int sno = 0; sno < nSlicesx; sno++){
759  f->GetList()->Add(slices[sno]);
760  }
761 
762  f->Write();
763 
764  gStyle -> SetOptStat(11);
765  gStyle -> SetOptFit(112);
766 
767  if(print==1){
768  if(debug==1){
769  std::cout<<"\nSaving Histograms as .png\n";
770  }
771 
772  gStyle -> SetOptStat(0);
773 
774  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->SetAxisRange(mean_min,mean_max,"Y");
775  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->Draw();
776  c->Print(detector + "_meanE_ge_EtaCut_CircularCut.png");
777 
778  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->SetAxisRange(chi2_min,chi2_max,"Y");
779  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->Draw();
780  c->Print(detector + "_chi2E_ge_EtaCut_CircularCut.png");
781 
782  gStyle -> SetOptStat(0);
783  gStyle -> SetOptFit(0);
784 
785  fExp->SetLineColor(4); //38
786  fTrue->SetLineColor(2); //46
787 
788  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerStyle(kFullCircle);
789  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerColor(46); //30
790  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerSize(0.75);
791  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetAxisRange(sigma_min,sigma_max,"Y");
792 
793  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->Fit("fTrue", "M+");
794  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->Draw("same");
795  fExp->Draw("same");
796 
797  TLegend* legend = new TLegend(1.75,1.75);
798  legend->SetHeader("Legend", "C");
799  legend->AddEntry(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2, "#sigma_{e_{agg}} vs Generated Energy", "flep");
800  legend->AddEntry((TObject*)0,"","");
801  legend->AddEntry(fTrue, "p_{0} + p_{1}/#sqrt{ge} + p_{2}/ge (Fitted)", "l");
802  legend->AddEntry((TObject*)0,"","");
803  legend->AddEntry(fExp, eRes1, "l");
804  legend->AddEntry((TObject*)0,"(Requirement)","" );
805  legend->SetTextSize(0.033);
806  legend->Draw();
807 
808  std::cout<<"reduced_chi2 of fit: "<<fTrue->GetChisquare()/fTrue->GetNDF()<<"\n";
809 
810  c->Print(detector + "_sigmaE_ge_EtaCut_CircularCut.png");
811 
812  gStyle -> SetOptStat(11);
813  gStyle -> SetOptFit(112);
814  for(int sno = 0; sno < nSlicesx; sno++){ // sno = 1
815  TString plusOne = (TString)(sno + 1);
816  TString nameF = detector + "_sigmaE_slice" + arr[sno] + "_EtaCut_CircularCut.png";
817  slices[sno] -> Fit("fit", "M+");
818  slices[sno] -> Draw("hist same");
819  c->Print(nameF);
820  }
821 
822  gStyle -> SetOptStat(1);
823 
824  counts_towerCounts_EtaCut->Draw("colz");
825  c->Print(detector + "_counts_towerCounts_EtaCut.png");
826  counts_towerCounts_EtaCut_CircularCut->Draw("colz");
827  c->Print(detector + "_counts_towerCounts_EtaCut_CircularCut.png");
828 
829  te_by_ge_ge_EtaCut->Draw("colz");
830  c->Print(detector + "_te_by_ge_ge_EtaCut.png");
831  te_by_ge_ge_EtaCut_CircularCut->Draw("colz");
832  c->Print(detector + "_te_by_ge_ge_EtaCut_CircularCut.png");
833 
834  gStyle -> SetOptStat(0);
835 
836  mean_te_by_ge_ge_EtaCut_CircularCut->Draw();
837  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut.png");
838 
839  gStyle -> SetOptStat(1);
840 
841  tphi_gphi_EtaCut->Draw("colz");
842  c->Print(detector + "_tphi_gphi_EtaCut.png");
843  tphi_gphi_EtaCut_CircularCut->Draw("colz");
844  c->Print(detector + "_tphi_gphi_EtaCut_CircularCut.png");
845 
846  ttheta_gtheta_EtaCut->Draw("colz");
847  c->Print(detector + "_ttheta_gtheta_EtaCut.png");
848  ttheta_gtheta_EtaCut_CircularCut->Draw("colz");
849  c->Print(detector + "_ttheta_gtheta_EtaCut_CircularCut.png");
850 
851  te_minus_ge_by_ge_ge_EtaCut->Draw("colz");
852  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut.png");
853  te_minus_ge_by_ge_ge_EtaCut_CircularCut->Draw("colz");
854  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut_CircularCut.png");
855  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Draw("colz");
856  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated.png");
857 
858  c->SetCanvasSize(700,700);
859  dphi_dtheta_EtaCut->Draw("colz");
860  el1->SetFillColorAlpha(0, 0);
861  el1->SetLineColorAlpha(2, 1);
862  el1->SetLineWidth(3);
863  el1->Draw("same");
864  c->Print(detector + "_dphi_dtheta_EtaCut_CircularCut.png");
865 
866  c->SetCanvasSize(500,700);
867  }
868 
869  if (mips == 1) {
870  gStyle -> SetOptStat(11);
871  gStyle -> SetOptFit(112);
872 
873  te_aggregate_EtaCut_CircularCut->Draw();
874  c->Print(detector + "_te_aggregate_EtaCut_CircularCut.png");
875 
876  hist_geta->Draw("colz");
877  c->Print(detector + "_hist_geta.png");
878 
879  gStyle -> SetOptStat(0);
880  gStyle -> SetOptFit(0);
881 
882  te_aggregate_by_ge->Draw("colz");
883  c->Print(detector + "_te_aggregate_by_ge.png");
884 
885  te_aggregate_by_geta->Draw("colz");
886  c->Print(detector + "_te_aggregate_by_geta.png");
887 
888  TF1 *polFit = new TF1("polFit","pol4",eta_min,eta_max);
889  TF1 *polFit2 = new TF1("polFit2","pol4",2*atan(exp(-eta_min)),2*atan(exp(-eta_max)));
890  if (detector == "FHCAL"){
891  polFit = new TF1("polFit","pol7",eta_min+0.2,eta_max-0.55);
892  polFit2 = new TF1("polFit2","pol6",2*atan(exp(-eta_min-0.2)),2*atan(exp(-eta_max+0.55)));
893  }
894  if (detector == "EEMC"){
895  polFit = new TF1("polFit","pol3",eta_min+0.15,eta_max-0.05);
896  polFit2 = new TF1("polFit2","pol3",2*atan(exp(-eta_min-0.15)),2*atan(exp(-eta_max+0.05)));
897  }
898  if (detector == "FEMC"){
899  polFit = new TF1("polFit","pol4",eta_min+0.1,eta_max-0.3);
900  polFit2 = new TF1("polFit2","pol4",2*atan(exp(-eta_min-0.1)),2*atan(exp(-eta_max+0.3)));
901  }
902 
903  TF1 *HcalFit = new TF1("HcalFit","pol2",-0.96,-0.65);
904  TF1 *HcalFit2 = new TF1("HcalFit2","pol2",-0.65,0.7);
905  TF1 *HcalFit3 = new TF1("HcalFit3","pol2",0.69,0.92);
906  TF1 *HcalFit4 = new TF1("HcalFit4","pol2",2*atan(exp(0.96)),2*atan(exp(0.65)));
907  TF1 *HcalFit5 = new TF1("HcalFit5","pol2",2*atan(exp(0.65)),2*atan(exp(-0.7)));
908  TF1 *HcalFit6 = new TF1("HcalFit6","pol2",2*atan(exp(-0.68)),2*atan(exp(-0.93)));
909 
910  if (detector == "HCALOUT"){
911  mean_te_geta_EtaCut->Fit("HcalFit", "ER+");
912  mean_te_geta_EtaCut->Fit("HcalFit2", "ER+");
913  mean_te_geta_EtaCut->Fit("HcalFit3", "ER+");
914  mean_te_geta_EtaCut->Draw("colz");
915  c->Print(detector + "_mean_te_geta_EtaCut.png");
916  std::cout<<"reduced_chi2 of eta fits from left to right: "<<HcalFit->GetChisquare()/HcalFit->GetNDF()<<", "<<HcalFit2->GetChisquare()/HcalFit2->GetNDF()<<", "<<HcalFit3->GetChisquare()/HcalFit3->GetNDF()<<"\n";
917 
918  mean_te_gtheta_EtaCut->Fit("HcalFit4", "ER+");
919  mean_te_gtheta_EtaCut->Fit("HcalFit5", "ER+");
920  mean_te_gtheta_EtaCut->Fit("HcalFit6", "ER+");
921  mean_te_gtheta_EtaCut->Draw("colz");
922  c->Print(detector + "_mean_te_gtheta_EtaCut.png");
923  std::cout<<"reduced_chi2 of theta fits from left to right: "<<HcalFit4->GetChisquare()/HcalFit4->GetNDF()<<", "<<HcalFit5->GetChisquare()/HcalFit5->GetNDF()<<", "<<HcalFit6->GetChisquare()/HcalFit6->GetNDF()<<"\n";
924  }
925 
926  else{
927  mean_te_geta_EtaCut->Fit("polFit", "ER+");
928  mean_te_geta_EtaCut->Draw("colz");
929  c->Print(detector + "_mean_te_geta_EtaCut.png");
930  std::cout<<"reduced_chi2 of eta fit: "<<polFit->GetChisquare()/polFit->GetNDF()<<"\n";
931 
932  mean_te_gtheta_EtaCut->Fit("polFit2", "ER+");
933  mean_te_gtheta_EtaCut->Draw("colz");
934  c->Print(detector + "_mean_te_gtheta_EtaCut.png");
935  std::cout<<"reduced_chi2 of theta fit: "<<polFit2->GetChisquare()/polFit2->GetNDF()<<"\n";
936  }
937  }
938 
939  c->Close();
940  std::cout << "The total_te is: " << total_te << endl;
941  std::cout << "The total_te_CircularCut is: " << total_te_CircularCut << endl;
942  std::cout << "The total_ge is: " << total_ge << endl;
943  std::cout<<"\n\nDone\n----------------------------------------------------------------------\n\n";
944 
945 }