EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LoopEvalHR.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LoopEvalHR.C
1 /*
2 > LoopEvalHR.C(int print = 1, int debug = 0, Double_t energyCutAggregate = 0.1, Double_t energyCut = 0.0, int MIP_theta_parametrisation = 1)
3 - Creates the analysis plots for the combined forward calorimeters CEMC+HCALIN+HCALOUT (default parameters for pions)
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 CEMC alone and CEMC+HCALIN+HCALOUT as well, and polar angle based energy cuts to eliminate MIPs
5 - Arguments
6  # print - saves plots as .png files if not 0
7  # debug - prints debugging messages if not 0; Pipe the output to a file while using this
8  # energyCutAggregate - specify the value for the tower energy cut on CEMC+HCALIN+HCALOUT
9  # energyCut - specify the value for the tower energy cut on individual towers
10  # MIP_theta_parametrization - applies a polar angle dependent energy cut on individual towers of CEMC to eliminate MIPs, based on the parametric equation provided in the code
11  - Output file - energy_verification_EtaCut_CircularCut_CEMC_HCALIN_HCALOUT.root, and the .png plots if generated
12 */
13 
14 /*
15  authors - Sagar Joshi (ee190002054@iiti.ac.in)
16  Siddhant Rathi (me190003061@iiti.ac.in)
17 
18  version - 2
19 
20 */
21 
22 #include <iostream>
23 #include <stdexcept>
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 LoopEvalHR(int print = 1, int debug = 0, Double_t energyCutAggregate = 0.1, Double_t energyCut = 0.0, int MIP_theta_parametrisation = 1){
32 
33  Double_t EMC_cut = 0.0;
34  TF1 *mip_pmzn_energy_cut_ftheta = new TF1("mip_pmzn_energy_cut_ftheta", "(9.46093e-01) - 1.62771*x + 1.37776*(x^2) - (5.4996e-01)*(x^3) + (8.82673e-02)*(x^4)");
35 
36 
37  if(MIP_theta_parametrisation == 1 && energyCut != 0){
38  throw std::invalid_argument("We do not currently support theta-parametrized \nMIP cut on EMC simultaneously with individual tower cuts \non other detectors.:(;");
39  }
40 
41 
42  //variable binning code
43  int arraySizeBins = 0; //size of the bins array
44  Double_t arraySizeBinsLoop = 0; //loop variable
45  int lowThresholdBins = 2; //number of bins below the threshold bin
46  int highThresholdBins = 9; //number of bins above the threshold bin
47  Double_t maxEnergyBin = 30; //highest energy
48  Double_t lowEnergyBin = 3; //threshold energy
49 
50  //calculating the number of bins required
51  while(arraySizeBinsLoop < maxEnergyBin){
52  arraySizeBins += 1;
53  if(arraySizeBinsLoop < lowEnergyBin){
54  arraySizeBinsLoop += lowEnergyBin/lowThresholdBins;
55  }
56  else{
57  arraySizeBinsLoop += (maxEnergyBin - lowEnergyBin)/highThresholdBins;
58  }
59  }
60 
61  Double_t binLimitArray[arraySizeBins + 1];
62 
63  //filing the array containing bin lower limits
64  for(int i = 0; i <= arraySizeBins; i++){
65  if(i <= lowThresholdBins){
66  binLimitArray[i] = lowEnergyBin*i/lowThresholdBins;
67  }
68  else{
69  binLimitArray[i] = lowEnergyBin + ((maxEnergyBin - lowEnergyBin)*(i - lowThresholdBins)/highThresholdBins);
70  }
71  }
72 
73  //print array
74  for(int i = 0; i < arraySizeBins; i++){
75  cout<<"element "<<i<<" of the bin array is "<<binLimitArray[i]<<" \n";
76  }
77 
78  //variable binning code end
79 
80  TString detector = "CEMC_HCALIN_HCALOUT";
81  TFile *f1 = new TFile("merged_Eval_HCALIN.root","READ");
82  TFile *f2 = new TFile("merged_Eval_HCALOUT.root","READ");
83  TFile *f3 = new TFile("merged_Eval_CEMC.root","READ");
84 
85  TTree* T1 = (TTree*)f1->Get("T");
86  EvalRootTTree *evaltree1 = nullptr;
87 
88  TTree* T2 = (TTree*)f2->Get("T");
89  EvalRootTTree *evaltree2 = nullptr;
90 
91  TTree* T3 = (TTree*)f3->Get("T");
92  EvalRootTTree *evaltree3 = nullptr;
93 
94  gStyle->SetCanvasPreferGL(kTRUE);
95  TCanvas *c = new TCanvas();
96  c->SetTickx();
97  c->SetTicky();
98 
99  // Modifying default plotting style
100  gStyle->SetOptTitle(0);
101  gStyle->SetOptFit(102);
102  gStyle->SetTitleXOffset(1);
103  gStyle->SetTitleYOffset(1);
104  gStyle->SetLabelSize(0.05);
105  gStyle->SetTitleXSize(0.05);
106  gStyle->SetTitleYSize(0.05);
107 
108  long double total_te = 0;
109  long double total_te_CircularCut = 0;
110  long double total_ge = 0;
111 
112  int nSlicesx = arraySizeBins; // Number of ge-axis slices taken for making sigma_e vs ge plot
113  int nSlicesy = 350;
114  Double_t eta_min, eta_max; // Eta range of detectors
115  Double_t x_radius_HCALIN = 0.15; // Length of semi-minor axis of circular (elliptical) in HCALIN
116  Double_t y_radius_HCALIN = 0.25; // Length of semi-major axis of circular (elliptical) in HCALIN
117  Double_t x_radius_HCALOUT = 0.2; // Length of semi-minor axis of circular (elliptical) in HCALOUT
118  Double_t y_radius_HCALOUT = 0.3; // Length of semi-major axis of circular (elliptical) in HCALOUT
119  Double_t x_radius_CEMC = 0.1; // Length of semi-minor axis of circular (elliptical) in CEMC
120  Double_t y_radius_CEMC = 0.2; // Length of semi-minor axis of circular (elliptical) in CEMC
121  Double_t fit_min, fit_max; // Fit range of the first slice of [(te-ge)/ge vs ge] plot
122  Double_t sigma_min, sigma_max; // Range of Y-axis in sigma_e vs ge plot
123  Double_t mean_min, mean_max; // Range of Y-axis in mean_e vs ge plot
124  Double_t chi2_min, chi2_max; // Range of Y-axis in chi2_e vs ge plot
125  Double_t recalibration_factor; // Number divided from entries of [(te-ge)/ge vs ge] plot for recalibration
126  Double_t recalibration_firstSlice = 3.8147;
127  Double_t recalibration_firstSlice_FHCAL = 0.0260; // Recalibration factor of first slice (needed to be done manually because of low statistics)
128  Double_t recalibration_firstSlice_FEMC = 0.0300;
129  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
130 
131  fit_min = -0.8;
132  fit_max = 1.0;
133  eta_min = -0.96;
134  eta_max = 0.92;
135  sigma_min = 0;
136  sigma_max = 1.5;
137  mean_min = -0.7;
138  mean_max = 0.3;
139  chi2_min = 0;
140  chi2_max = 2.23;
141  // recalibration_factor = 0.7088;
142  te_minus_ge_by_ge_ge_min = -0.99;
143  te_minus_ge_by_ge_ge_max = 1.0;
144 
145 
146  TString cut_text = " {-0.96 < geta < 0.92} ";
147 
148 
149  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,-2,1);
150  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);;
151  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.5);
152  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,binLimitArray,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
153 
154  TH2D *te_minus_ge_by_ge_ge_EtaCut_temp = new TH2D("te_minus_ge_by_ge_ge_EtaCut_temp","#frac{#Delta e_{agg}}{truth e} vs truth e",nSlicesx,binLimitArray,500,-0.99,1.3); // histogram from which mean vs ge, sigma vs ge, and reduced_chi2 vs ge plots are derived
155 
156  TH2D *te_by_ge_ge_EtaCut = new TH2D("te_by_ge_ge_EtaCut","te_{agg}/ge vs ge",200,0,30,200,-0.5,1.5);
157  TH2D *te_by_ge_ge_EtaCut_CircularCut = new TH2D("te_by_ge_ge_EtaCut_CircularCut","te_{agg}/ge vs ge",200,0,30,200,-1,2);
158  TH2D *te_by_ge_ge_EtaCut_CircularCut_HCALIN = new TH2D("te_by_ge_ge_EtaCut_CircularCut_HCALIN","te_{agg}/ge vs ge",200,0,30,200,-1,2);
159  TH2D *te_by_ge_ge_EtaCut_CircularCut_HCALOUT = new TH2D("te_by_ge_ge_EtaCut_CircularCut_HCALOUT","te_{agg}/ge vs ge",200,0,30,200,-1,2);
160  TH2D *te_by_ge_ge_EtaCut_CircularCut_CEMC = new TH2D("te_by_ge_ge_EtaCut_CircularCut_CEMC","te_{agg}/ge vs ge",200,0,30,200,-1,2);
161 
162  auto *mean_te_by_ge_ge_EtaCut_CircularCut = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut","Mean_{te/ge}",nSlicesx,binLimitArray,-0.5,35);
163  auto *mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN","Mean_{te/ge}",nSlicesx,binLimitArray,-0.5,35);
164  auto *mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT","Mean_{te/ge}",nSlicesx,binLimitArray,-0.5,35);
165  auto *mean_te_by_ge_ge_EtaCut_CircularCut_CEMC = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut_CEMC","Mean_{te/ge}",nSlicesx,binLimitArray,-0.5,35);
166 
167  T1->SetBranchAddress("DST#EvalTTree_HCALIN",&evaltree1);
168  T2->SetBranchAddress("DST#EvalTTree_HCALOUT",&evaltree2);
169  T3->SetBranchAddress("DST#EvalTTree_CEMC",&evaltree3);
170 
171  //T1->GetEntries()
172  for(int i=0; i<T1->GetEntries(); i++) // We assume same no. of entries, since no cuts are applied
173  {
174  T1->GetEntry(i);
175  T2->GetEntry(i);
176  T3->GetEntry(i);
177 
178  if(debug==1){
179  std::cout<<"\n\n\n------------------------------------------\nParticle: "<<i<<"\n\n";
180  std::cout<<"Initial Parameters "<<"\n";
181  }
182 
183  Double_t geta1 = evaltree1->get_geta();
184 
185  if(debug==1){
186  std::cout<<"geta: "<<geta1<<"\n";
187  }
188 
189  if(geta1>=eta_min && geta1<=eta_max){
190 
191  if(debug==1){
192  cout<<"\ngeta cut applied (1.3, 3.3)"<<"\n\n";
193  }
194  Double_t ge = evaltree1->get_ge();
195  if(debug==1){
196  std::cout<<"ge: "<<ge<<"\n";
197  }
198 
199  Double_t gphi = evaltree1->get_gphi();
200  if(debug==1){
201  std::cout<<"gphi: "<<gphi<<"\n";
202  }
203 
204  Double_t gtheta = evaltree1->get_gtheta();
205  if(debug==1){
206  std::cout<<"gtheta: "<<gtheta<<"\n";
207  }
208 
209  total_ge += ge;
210  if(debug==1){
211  std::cout<<"total_ge till now = "<<total_ge<<"\n";
212  }
213 
214  Double_t te_aggregate = 0;
215  Double_t te_aggregate_CircularCut = 0;
216  Double_t te_aggregate_HCALIN_CircularCut = 0;
217  Double_t te_aggregate_HCALOUT_CircularCut = 0;
218 
219  for (int j=0; j<evaltree1->get_ntowers(); j++){
220 
221  if(debug==1){
222  std::cout<<"\nHCALIN Tower: "<<j<<"\n";
223  }
224 
225  EvalTower *twr1 = evaltree1->get_tower(j);
226  if (twr1){
227 
228  if(debug==1){
229  cout<<"non-empty HCALIN tower\n";
230  }
231 
232  Double_t tphi = twr1->get_tphi();
233  if(debug==1){
234  std::cout<<"tphi: "<<tphi<<"\n";
235  }
236 
237  Double_t ttheta = twr1->get_ttheta();
238  if(debug==1){
239  std::cout<<"ttheta: "<<ttheta<<"\n";
240  }
241 
242  Double_t dphi = tphi - gphi;
243  if(debug==1){
244  std::cout<<"tphi-gphi: "<<dphi<<"\n";
245  }
246 
247  Double_t dtheta = ttheta - gtheta;
248  if(debug==1){
249  std::cout<<"ttheta-gtheta: "<<dtheta<<"\n";
250  }
251 
252  Double_t te = twr1->get_te();
253  if(debug==1){
254  std::cout<<"te: "<<te<<"\n";
255  }
256 
257  if(te > energyCut){
258  te_aggregate += te;
259 
260  if(debug==1){
261  std::cout<<"HCALIN: pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALOUT,2) = "<<pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2)<<"\n";
262  }
263 
264  if (pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2)<=1){
265  if(debug==1){
266  cout<<"HCALIN Tower included after circular cut\n";
267  }
268 
269  te_aggregate_CircularCut += te;
270  te_aggregate_HCALIN_CircularCut += te;
271 
272  }
273 
274  if(debug==1){
275  cout<<"te_aggregate till now = "<<te_aggregate<<"\n";
276  std::cout<<"te += "<<twr1->get_te()<<"\n";
277  cout<<"te_aggregate_CircularCut till now = "<<te_aggregate_CircularCut<<"\n";
278  }
279  }
280  }
281  }
282 
283  for (int j=0; j<evaltree2->get_ntowers(); j++){
284 
285 
286  if(debug==1){
287  std::cout<<"\nHCALOUT Tower: "<<j<<"\n";
288  }
289 
290  EvalTower *twr2 = evaltree2->get_tower(j);
291  if (twr2){
292 
293  if(debug==1){
294  cout<<"non-empty HCALOUT tower\n";
295  }
296 
297  Double_t tphi = twr2->get_tphi();
298  if(debug==1){
299  std::cout<<"tphi: "<<tphi<<"\n";
300  }
301 
302  Double_t ttheta = twr2->get_ttheta();
303  if(debug==1){
304  std::cout<<"ttheta: "<<ttheta<<"\n";
305  }
306 
307  Double_t dphi = tphi - gphi;
308  if(debug==1){
309  std::cout<<"tphi-gphi: "<<dphi<<"\n";
310  }
311 
312  Double_t dtheta = ttheta - gtheta;
313  if(debug==1){
314  std::cout<<"ttheta-gtheta: "<<dtheta<<"\n";
315  }
316 
317  Double_t te = twr2->get_te();
318  if(debug==1){
319  std::cout<<"te: "<<te<<"\n";
320  }
321 
322  if(te > energyCut){
323  te_aggregate += te;
324 
325  if(debug==1){
326  std::cout<<"HCALOUT: pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2): "<<pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2)<<"\n";
327  }
328 
329  if (pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2)<=1){
330  if(debug==1){
331  cout<<"HCALOUT Tower included after circular cut\n";
332  }
333 
334  te_aggregate_HCALOUT_CircularCut += te;
335  te_aggregate_CircularCut += te;
336  }
337 
338  if(debug==1){
339  cout<<"te_aggregate till now = "<<te_aggregate<<"\n";
340  cout<<"te_aggregate_CircularCut till now = "<<te_aggregate_CircularCut<<"\n";
341  std::cout<<"te += "<<twr2->get_te()<<"\n";
342  }
343  }
344  }
345  }
346 
347 
348  for (int j=0; j<evaltree3->get_ntowers(); j++){
349 
350 
351  if(debug==1){
352  std::cout<<"\nCEMC Tower: "<<j<<"\n";
353  }
354 
355  EvalTower *twr3 = evaltree3->get_tower(j);
356  if (twr3){
357 
358  if(debug==1){
359  cout<<"non-empty CEMC tower\n";
360  }
361 
362  Double_t tphi = twr3->get_tphi();
363  if(debug==1){
364  std::cout<<"tphi: "<<tphi<<"\n";
365  }
366 
367  Double_t ttheta = twr3->get_ttheta();
368  if(debug==1){
369  std::cout<<"ttheta: "<<ttheta<<"\n";
370  }
371 
372  Double_t dphi = tphi - gphi;
373  if(debug==1){
374  std::cout<<"tphi-gphi: "<<dphi<<"\n";
375  }
376 
377  Double_t dtheta = ttheta - gtheta;
378  if(debug==1){
379  std::cout<<"ttheta-gtheta: "<<dtheta<<"\n";
380  }
381 
382  Double_t te = twr3->get_te();
383  if(debug==1){
384  std::cout<<"te: "<<te<<"\n";
385  }
386 
387  if(MIP_theta_parametrisation == 1){
388  EMC_cut = mip_pmzn_energy_cut_ftheta->Eval(gtheta);
389  }
390 
391  if(te > energyCut + EMC_cut){
392  te_aggregate += te;
393 
394  if(debug==1){
395  std::cout<<"CEMC: pow(dphi/y_radius,2)+pow(dtheta/x_radius,2): "<<pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2)<<"\n";
396  }
397 
398  if (pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2)<=1){
399  if(debug==1){
400  cout<<"CEMC Tower included after circular cut\n";
401  }
402  te_aggregate_CircularCut += te;
403  }
404 
405  if(debug==1){
406  cout<<"te_aggregate till now = "<<te_aggregate<<"\n";
407  cout<<"te_aggregate_CircularCut till now = "<<te_aggregate_CircularCut<<"\n";
408  std::cout<<"te += "<<twr3->get_te()<<"\n";
409  }
410  }
411  }
412  }
413 
414  total_te += te_aggregate;
415  total_te_CircularCut += te_aggregate_CircularCut;
416 
417  if(debug==1){
418  cout<<"total_te till now = "<<total_te<<"\n";
419  cout<<"total_te_CircularCut till now = "<<total_te_CircularCut<<"\n\n";
420  }
421 
422  if(te_aggregate_CircularCut > energyCutAggregate){
423 
424  te_minus_ge_by_ge_ge_EtaCut->Fill(ge, (te_aggregate-ge)/ge);
425 
426  if(debug==1){
427  cout<<"(ge, (te_aggregate-ge)/ge): ("<<ge<<", "<<(te_aggregate-ge)/ge<<")\n";
428  }
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 
436  /*te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Fill(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge);
437  r early on, it took me a while to realize what I need is best done by a weighted sum instead of trying to normalize the fit functions. if(debug==1){
438  cout<<"(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge: ("<<ge<<", "<<((te_aggregate_CircularCut/recalibration_factor)-ge)/ge<<")\n";
439  }
440  te_minus_ge_by_ge_ge_EtaCut_Recalibrated_temp->Fill(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge);*/
441 
442  te_by_ge_ge_EtaCut->Fill(ge, te_aggregate/ge);
443 
444  if(debug==1){
445  cout<<"(ge, te_aggregate/ge): ("<<ge<<", "<<te_aggregate/ge<<")\n";
446  }
447 
448  te_by_ge_ge_EtaCut_CircularCut->Fill(ge, te_aggregate_CircularCut/ge);
449  te_by_ge_ge_EtaCut_CircularCut_HCALIN->Fill(ge, te_aggregate_HCALIN_CircularCut/ge);
450  te_by_ge_ge_EtaCut_CircularCut_HCALOUT->Fill(ge, te_aggregate_HCALOUT_CircularCut/ge);
451  te_by_ge_ge_EtaCut_CircularCut_CEMC->Fill(ge, (te_aggregate_CircularCut-te_aggregate_HCALIN_CircularCut-te_aggregate_HCALOUT_CircularCut)/ge);
452 
453  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->Fill(ge, te_aggregate_HCALIN_CircularCut/ge);
454  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->Fill(ge, te_aggregate_HCALOUT_CircularCut/ge);
455  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->Fill(ge, (te_aggregate_CircularCut-te_aggregate_HCALIN_CircularCut-te_aggregate_HCALOUT_CircularCut)/ge);
456 
457 
458  if(debug==1){
459  cout<<"\n(ge, te_aggregate_CircularCut/ge): ("<<ge<<", "<<te_aggregate_CircularCut/ge<<")\n";
460  cout<<"(ge, te_aggregate_HCALIN_CircularCut/ge): ("<<ge<<", "<<te_aggregate_HCALIN_CircularCut/ge<<")\n";
461  cout<<"(ge, te_aggregate_HCALOUT_CircularCut/ge): ("<<ge<<", "<<te_aggregate_HCALOUT_CircularCut/ge<<")\n";
462  cout<<"(ge, (te_aggregate_CircularCut-te_aggregate_HCALIN_CircularCut - te_aggregate_HCALOUT_CircularCut)/ge): ("<<ge<<", "<<(te_aggregate_CircularCut-te_aggregate_HCALIN_CircularCut - te_aggregate_HCALOUT_CircularCut)/ge<<")\n\n";
463  }
464  }
465  }
466  }
467 
468  if(debug==1){
469  cout<<"\neta cut if ends"<<"\n";
470  }
471 
472 
473 
474  TString arr[nSlicesx]; // Used for naming fitted slices used in sigma_e vs ge
475  for(int sno = 0; sno < nSlicesx; sno++){
476  arr[sno] = TString::Itoa(sno + 1, 10);
477  }
478 
479  if(debug==1){
480  std::cout<<"\nGenerating sigma_e vs ge plots\n\n";
481  }
482 
483  Double_t recalibrationArr1[nSlicesx];
484  Double_t rf_integral_HCALIN = 0;
485  Double_t weight_HCALIN = te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetMean(2);
486 
487  if(debug==1){
488  cout << "HCALIN weight is : " << weight_HCALIN << endl;
489  }
490 
491  //recalibrationArr1[0] = recalibration_firstSlice_FHCAL;
492  //rf_integral_FHCAL += recalibrationArr1[0];
493  //cout << "Recalibration factor for slice " << 1 << " of FHCAL is: " << recalibrationArr1[0] << endl;
494 
495  for(int binIter = 1; binIter <= nSlicesx; binIter++){
496  recalibrationArr1[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetBinContent(binIter);
497  rf_integral_HCALIN += recalibrationArr1[binIter-1];
498 
499  if(debug == 1){
500  cout<<"rf_integral_HCALIN += "<<recalibrationArr1[binIter-1]<<"\n";
501  }
502 
503  cout << "Recalibration factor for slice " << binIter << " of HCALIN is: " << recalibrationArr1[binIter-1] << endl;
504  }
505 
506  if(debug == 1){
507  cout<<"rf_integral_HCALIN = "<<rf_integral_HCALIN<<"\n\n";
508  }
509 
510  Double_t recalibrationArr2[nSlicesx];
511  Double_t rf_integral_HCALOUT = 0;
512  Double_t weight_HCALOUT = te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetMean(2);
513 
514  if(debug==1){
515  cout << "HCALOUT weight is : " << weight_HCALOUT << endl;
516  }
517 
518  //recalibrationArr2[0] = recalibration_firstSlice_HCALOUT;
519  //rf_integral_HCALOUT += recalibrationArr2[0];
520  //cout << "Recalibration factor for slice " << 1 << " of HCALOUT is: " << recalibrationArr2[0] << endl;
521 
522  for(int binIter = 1; binIter <= nSlicesx; binIter++){
523  recalibrationArr2[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetBinContent(binIter);
524  rf_integral_HCALOUT += recalibrationArr2[binIter-1];
525 
526  if(debug == 1){
527  cout<<"rf_integral_HCALOUT += "<<recalibrationArr2[binIter-1]<<"\n";
528  }
529 
530  cout << "Recalibration factor for slice " << binIter << " of HCALOUT is: " << recalibrationArr2[binIter-1] << endl;
531  }
532 
533  if(debug == 1){
534  cout<<"rf_integral_HCALOUT = "<<rf_integral_HCALOUT<<"\n\n";
535  }
536 
537 
538  Double_t recalibrationArr3[nSlicesx];
539  Double_t rf_integral_CEMC = 0;
540  Double_t weight_CEMC = te_by_ge_ge_EtaCut_CircularCut_CEMC->GetMean(2);
541 
542  if(debug==1){
543  cout << "CEMC weight is : " << weight_CEMC << endl;
544  }
545 
546  //recalibrationArr3[0] = recalibration_firstSlice_CEMC;
547  //rf_integral_CEMC += recalibrationArr3[0];
548  //cout << "Recalibration factor for slice " << 1 << " of CEMC is: " << recalibrationArr3[0] << endl;
549 
550  for(int binIter = 1; binIter <= nSlicesx; binIter++){
551  recalibrationArr3[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetBinContent(binIter);
552  rf_integral_CEMC += recalibrationArr3[binIter-1];
553 
554  if(debug == 1){
555  cout<<"rf_integral_CEMC += "<<recalibrationArr3[binIter-1]<<"\n";
556  }
557 
558  cout << "Recalibration factor for slice " << binIter << " of CEMC is: " << recalibrationArr3[binIter-1] << endl;
559  }
560 
561  if(debug == 1){
562  cout<<"rf_integral_CEMC = "<<rf_integral_CEMC<<"\n\n";
563  }
564 
565  for(int i=0; i<T1->GetEntries(); i++){
566 
567  T1->GetEntry(i);
568  T2->GetEntry(i);
569  T3->GetEntry(i);
570 
571  Double_t geta1 = evaltree1->get_geta();
572  Double_t gphi = evaltree1->get_gphi();
573  Double_t gtheta = evaltree1->get_gtheta();
574  Double_t ge = evaltree1->get_ge();
575  Double_t te_aggregate_CircularCut = 0;
576  Double_t te_aggregate_CircularCut_normalised = 0;
577 
578  int recalibration_factor = 0;
579 
580  if(ge > lowEnergyBin){
581  Double_t eRangeBin = (maxEnergyBin - lowEnergyBin)/highThresholdBins;
582  recalibration_factor = (lowThresholdBins - 1) + ceil(ge/eRangeBin)- 1;
583  }
584 
585  else{
586  recalibration_factor = ceil((ge/lowEnergyBin)*(Double_t)lowThresholdBins)- 1;
587  }
588 
589  if(debug==1){
590  cout << "Recalibration index in Loop 2 is : " << recalibration_factor << endl;
591  }
592 
593  if(geta1>=eta_min && geta1<=eta_max){
594 
595 
596  for (int j=0; j<evaltree1->get_ntowers(); j++){
597 
598  EvalTower *twr1 = evaltree1->get_tower(j);
599  if (twr1){
600 
601  Double_t tphi = twr1->get_tphi();
602  Double_t ttheta = twr1->get_ttheta();
603  Double_t dphi = tphi - gphi;
604  Double_t dtheta = ttheta - gtheta;
605  Double_t te = twr1->get_te();
606 
607  if(te > energyCut){
608 
609  if(debug==1){
610  std::cout<<"pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2): "<<pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2)<<"\n";
611  }
612 
613  if (pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2)<=1){
614  te_aggregate_CircularCut += te;
615  te_aggregate_CircularCut_normalised += te*weight_HCALIN/recalibrationArr1[recalibration_factor];
616  }
617  }
618  }
619  }
620 
621  for (int j=0; j<evaltree2->get_ntowers(); j++){
622 
623  EvalTower *twr2 = evaltree2->get_tower(j);
624  if (twr2){
625 
626  Double_t tphi = twr2->get_tphi();
627  Double_t ttheta = twr2->get_ttheta();
628  Double_t dphi = tphi - gphi;
629  Double_t dtheta = ttheta - gtheta;
630  Double_t te = twr2->get_te();
631 
632  if(te > energyCut){
633 
634  if(debug==1){
635  std::cout<<"pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2): "<<pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2)<<"\n";
636  }
637 
638  if (pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2)<=1){
639  te_aggregate_CircularCut += te;
640  te_aggregate_CircularCut_normalised += te*weight_HCALOUT/recalibrationArr2[recalibration_factor];
641  }
642  }
643  }
644  }
645 
646  for (int j=0; j<evaltree3->get_ntowers(); j++){
647 
648  EvalTower *twr3 = evaltree3->get_tower(j);
649  if (twr3){
650 
651  Double_t tphi = twr3->get_tphi();
652  Double_t ttheta = twr3->get_ttheta();
653  Double_t dphi = tphi - gphi;
654  Double_t dtheta = ttheta - gtheta;
655  Double_t te = twr3->get_te();
656 
657  if(MIP_theta_parametrisation == 1){
658  EMC_cut = mip_pmzn_energy_cut_ftheta->Eval(gtheta);
659  }
660 
661  if(te > energyCut + EMC_cut){
662 
663  if(debug==1){
664  std::cout<<"pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2): "<<pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2)<<"\n";
665  }
666 
667  if (pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2)<=1){
668  te_aggregate_CircularCut += te;
669  te_aggregate_CircularCut_normalised += te*weight_CEMC/recalibrationArr3[recalibration_factor];
670  }
671  }
672  }
673  }
674 
675  if(te_aggregate_CircularCut > energyCutAggregate){
676  mean_te_by_ge_ge_EtaCut_CircularCut->Fill(ge, te_aggregate_CircularCut_normalised/ge);
677 
678  if(debug==1){
679  cout << "mean_te_by_ge_ge_EtaCut_CircularCut entry " << i << " is : " << te_aggregate_CircularCut_normalised/ge << endl;
680  }
681 
682  }
683 
684  if(debug == 1){
685  cout<<"(ge, te_aggregate_CircularCut_normalised/ge): ("<<ge<<", "<<te_aggregate_CircularCut_normalised/ge<<")\n";
686  }
687  }
688  }
689 
690 
691  Double_t recalibrationArr[nSlicesx];
692 
693  // recalibrationArr[0] = recalibration_firstSlice;
694  // cout << "Recalibration factor for slice " << 1 << " is: " << recalibrationArr[0] << endl;
695 
696  for(int binIter = 1; binIter <= nSlicesx; binIter++){
697  recalibrationArr[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut->GetBinContent(binIter);
698  cout << "Recalibration factor for slice " << binIter << " is: " << recalibrationArr[binIter-1] << endl;
699  }
700 
701 
702  for(int i=0; i<T1->GetEntries(); i++){
703 
704  T1->GetEntry(i);
705  T2->GetEntry(i);
706  T3->GetEntry(i);
707 
708  Double_t geta1 = evaltree1->get_geta();
709  Double_t gphi = evaltree1->get_gphi();
710  Double_t gtheta = evaltree1->get_gtheta();
711  Double_t ge = evaltree1->get_ge();
712  Double_t te_aggregate_CircularCut = 0;
713  Double_t te_aggregate_CircularCut_normalised = 0;
714 
715  int recalibration_factor = 0;
716 
717 
718  if(ge > lowEnergyBin){
719  Double_t eRangeBin = (maxEnergyBin - lowEnergyBin)/highThresholdBins;
720  recalibration_factor = (lowThresholdBins - 1) + ceil(ge/eRangeBin)- 1;
721  }
722 
723  else{
724  recalibration_factor = ceil((ge/lowEnergyBin)*(Double_t)lowThresholdBins)- 1;
725  }
726 
727  if(debug==1){
728  cout << "Recalibration index in Loop 3 is : " << recalibration_factor << endl;
729  }
730 
731  if(geta1>=eta_min && geta1<=eta_max){
732 
733 
734  for (int j=0; j<evaltree1->get_ntowers(); j++){
735 
736  EvalTower *twr1 = evaltree1->get_tower(j);
737  if (twr1){
738 
739  Double_t tphi = twr1->get_tphi();
740  Double_t ttheta = twr1->get_ttheta();
741  Double_t dphi = tphi - gphi;
742  Double_t dtheta = ttheta - gtheta;
743  Double_t te = twr1->get_te();
744 
745  if(te > energyCut){
746 
747  if(debug==1){
748  std::cout<<"pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2): "<<pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2)<<"\n";
749  }
750 
751  if (pow(dphi/y_radius_HCALIN,2)+pow(dtheta/x_radius_HCALIN,2)<=1){
752  te_aggregate_CircularCut += te;
753  te_aggregate_CircularCut_normalised += te*weight_HCALIN/recalibrationArr1[recalibration_factor];
754  }
755  }
756  }
757  }
758 
759  for (int j=0; j<evaltree2->get_ntowers(); j++){
760 
761  EvalTower *twr2 = evaltree2->get_tower(j);
762  if (twr2){
763 
764  Double_t tphi = twr2->get_tphi();
765  Double_t ttheta = twr2->get_ttheta();
766  Double_t dphi = tphi - gphi;
767  Double_t dtheta = ttheta - gtheta;
768  Double_t te = twr2->get_te();
769 
770  if(te > energyCut){
771 
772  if(debug==1){
773  std::cout<<"pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2): "<<pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2)<<"\n";
774  }
775 
776  if (pow(dphi/y_radius_HCALOUT,2)+pow(dtheta/x_radius_HCALOUT,2)<=1){
777  te_aggregate_CircularCut += te;
778  te_aggregate_CircularCut_normalised += te*weight_HCALOUT/recalibrationArr2[recalibration_factor];
779  }
780  }
781  }
782  }
783 
784  for (int j=0; j<evaltree3->get_ntowers(); j++){
785 
786  EvalTower *twr3 = evaltree3->get_tower(j);
787  if (twr3){
788 
789  Double_t tphi = twr3->get_tphi();
790  Double_t ttheta = twr3->get_ttheta();
791  Double_t dphi = tphi - gphi;
792  Double_t dtheta = ttheta - gtheta;
793  Double_t te = twr3->get_te();
794 
795  if(MIP_theta_parametrisation == 1){
796  EMC_cut = mip_pmzn_energy_cut_ftheta->Eval(gtheta);
797  }
798 
799  if(te > energyCut + EMC_cut){
800 
801  if(debug==1){
802  std::cout<<"pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2): "<<pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2)<<"\n";
803  }
804 
805  if (pow(dphi/y_radius_CEMC,2)+pow(dtheta/x_radius_CEMC,2)<=1){
806  te_aggregate_CircularCut += te;
807  te_aggregate_CircularCut_normalised += te*weight_CEMC/recalibrationArr3[recalibration_factor];
808  }
809  }
810  }
811  }
812 
813  if(te_aggregate_CircularCut > energyCutAggregate){
814  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Fill(ge, ((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge);
815  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->Fill(ge, ((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge);
816 
817  if(debug==1){
818  cout << "te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated entry " << i << " is : " << ((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge << endl;
819  }
820 
821  }
822 
823  if(debug==1){
824  cout<<"(ge, ((te_aggregate_CircularCut_normalised/recalibration_factor)-ge)/ge: ("<<ge<<", "<<((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge<<")\n";
825  }
826  }
827  }
828 
829 
830  // Initialising fit functions
831  TF1 *fit = new TF1("fit", "gaus");
832  TF1 *fit1 = new TF1("fit1","gaus",fit_min,fit_max);
833  TF1 *fExp = new TF1("fExp","0.1 + 1.0/sqrt(x)",0,30);
834  TF1 *fTrue = new TF1("fTrue","[0] + [1]/sqrt(x)",0,30);
835 
836 
837  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->FitSlicesY(0, 1, -1, 0, "QN");
838  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");
839  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");
840  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");
841 
842 
843 
844  TH1D* slices[nSlicesx];
845 
846  // Generating plots for individual slices
847  for(int sno = 0; sno < nSlicesx; sno++){
848  int plusOne = sno+1;
849  TString sname = "slice " + arr[sno];
850  slices[sno] = te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->ProjectionY(sname, plusOne, plusOne);
851  }
852 
853 
854  if(debug==1){
855  cout<<"\nHistogram Formatting\n\n";
856  }
857 
858  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
859  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetLabelSize(0.05);
860  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetTitleSize(0.05);
861  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
862  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetLabelSize(0.05);
863  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetTitleSize(0.05);
864 
865  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
866  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
867  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
868  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
869  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
870  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
871 
872  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetTitle("Generated Energy (GeV)");
873  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetLabelSize(0.05);
874  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetTitleSize(0.05);
875  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
876  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetLabelSize(0.05);
877  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetTitleSize(0.05);
878 
879  te_by_ge_ge_EtaCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
880  te_by_ge_ge_EtaCut->GetXaxis()->SetLabelSize(0.05);
881  te_by_ge_ge_EtaCut->GetXaxis()->SetTitleSize(0.05);
882  te_by_ge_ge_EtaCut->GetYaxis()->SetTitle("te_{agg}/ge");
883  te_by_ge_ge_EtaCut->GetYaxis()->SetLabelSize(0.05);
884  te_by_ge_ge_EtaCut->GetYaxis()->SetTitleSize(0.05);
885 
886  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
887  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
888  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
889  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("te_{agg}/ge");
890  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
891  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
892 
893  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
894  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
895  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
896  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("Mean of te_{agg}/ge");
897  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
898  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
899 
900  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetXaxis()->SetTitle("Generated Energy (GeV)");
901  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetXaxis()->SetLabelSize(0.05);
902  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetXaxis()->SetTitleSize(0.05);
903  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetYaxis()->SetTitle("Mean of te_{agg}/ge (HCALIN)");
904  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetYaxis()->SetLabelSize(0.05);
905  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->GetYaxis()->SetTitleSize(0.05);
906 
907  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetXaxis()->SetTitle("Generated Energy (GeV)");
908  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetXaxis()->SetLabelSize(0.05);
909  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetXaxis()->SetTitleSize(0.05);
910  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetYaxis()->SetTitle("Mean of te_{agg}/ge (HCALOUT)");
911  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetYaxis()->SetLabelSize(0.05);
912  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->GetYaxis()->SetTitleSize(0.05);
913 
914  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetXaxis()->SetTitle("Generated Energy (GeV)");
915  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetXaxis()->SetLabelSize(0.05);
916  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetXaxis()->SetTitleSize(0.05);
917  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetYaxis()->SetTitle("Mean of te_{agg}/ge (CEMC)");
918  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetYaxis()->SetLabelSize(0.05);
919  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->GetYaxis()->SetTitleSize(0.05);
920 
921  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetTitle("Generated Energy (GeV)");
922  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetLabelSize(0.05);
923  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetTitleSize(0.05);
924  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetTitle("#sigma_{e_{agg}}");
925  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetLabelSize(0.05);
926  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetTitleSize(0.05);
927  //te_minus_ge_by_ge_ge_EtaCut_2->SetTitle("#sigma_{e_{agg}} vs true_e" + cut_text);
928 
929  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetTitle("Generated Energy (GeV)");
930  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetLabelSize(0.05);
931  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetTitleSize(0.05);
932  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetTitle("Reduced_#chi^{2}_{e_{agg}}");
933  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetLabelSize(0.05);
934  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetTitleSize(0.04);
935  //te_minus_ge_by_ge_ge_EtaCut_chi2->SetTitle("#chi^{2}_{e_{agg}} vs true_e" + cut_text);
936 
937  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetTitle("Generated Energy (GeV)");
938  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetLabelSize(0.05);
939  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetTitleSize(0.05);
940  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetTitle("Mean_{e_{agg}}");
941  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetLabelSize(0.05);
942  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetTitleSize(0.05);
943  //te_minus_ge_by_ge_ge_EtaCut_1->SetTitle("mean_{e_{agg}} vs true_e" + cut_text);
944 
945  for(int sno = 0; sno < nSlicesx; sno++){
946  slices[sno]->GetXaxis()->SetTitle("#Delta e^{agg}/ ge");
947  slices[sno]->GetXaxis()->SetLabelSize(0.05);
948  slices[sno]->GetXaxis()->SetTitleSize(0.05);
949  slices[sno]->GetYaxis()->SetTitle("Counts");
950  slices[sno]->GetYaxis()->SetLabelSize(0.05);
951  slices[sno]->GetYaxis()->SetTitleSize(0.05);
952  }
953 
954  if(debug==1){
955  std::cout<<"\nWrite Histograms to File\n";
956  }
957 
958  TFile *f = new TFile("energy_verification_EtaCut_CircularCut_" + detector + ".root","RECREATE");
959 
960  f->GetList()->Add(te_by_ge_ge_EtaCut);
961  f->GetList()->Add(te_by_ge_ge_EtaCut_CircularCut);
962  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut);
963  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut);
964  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut);
965  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN);
966  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT);
967  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut_CEMC);
968  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated);
969  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2);
970  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1);
971  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2);
972 
973  for(int sno = 0; sno < nSlicesx; sno++){
974  f->GetList()->Add(slices[sno]);
975  }
976 
977  f->Write();
978 
979  gStyle -> SetOptStat(11);
980  gStyle -> SetOptFit(112);
981 
982  /*
983  int sno = 0;
984  TString plusOne = (TString)(sno + 1);
985  TString nameF = detector + "_sigmaE_slice" + arr[sno] + "_EtaCut_CircularCut.png";
986  slices[sno] -> Fit("fit1", "R+");
987  slices[sno] -> Draw("hist same");
988  c->Print(nameF);
989  */
990 
991  /* double_t mean = fit1->GetParameter(1);
992  double_t mean_error = fit1->GetParError(1);
993  double_t sigma = fit1->GetParameter(2);
994  double_t sigma_error = fit1->GetParError(2);
995  double_t chi2 = (fit1->GetChisquare())/(fit1->GetNDF());
996 
997  TLine *Mean = new TLine(0,mean,30.0/nSlicesx,mean);
998  TLine *Sigma = new TLine(0,sigma,30.0/nSlicesx,sigma);
999  TLine *Chi2 = new TLine(0,chi2,30.0/nSlicesx,chi2);
1000  TLine *Periphery_Chi2 = new TLine(30.0/nSlicesx,0,30.0/nSlicesx,chi2);
1001  TLine *Sigma_error = new TLine(30.0/(2.0*nSlicesx),sigma-(sigma_error),30.0/(2.0*nSlicesx),sigma+(sigma_error));
1002  TLine *Mean_error = new TLine(30.0/(2.0*nSlicesx),mean-(mean_error),30.0/(2.0*nSlicesx),mean+(mean_error));
1003  TLine *Point = new TLine(30.0/(2.0*nSlicesx),sigma,30.0/(2.0*nSlicesx),sigma);*/
1004 
1005  if(print==1){
1006  if(debug==1){
1007  std::cout<<"\nSaving Histograms as .png\n";
1008  }
1009 
1010  gStyle -> SetOptStat(0);
1011 
1012  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->SetAxisRange(mean_min,mean_max,"Y");
1013  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->Draw();
1014  /*Mean->SetLineColor(1);
1015  Mean->SetLineWidth(1);
1016  Mean_error->SetLineColor(1);
1017  Mean_error->SetLineWidth(1);
1018  Mean->Draw("same");
1019  Mean_error->Draw("same");*/
1020  c->Print(detector + "_meanE_ge_EtaCut_CircularCut.png");
1021 
1022  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->SetAxisRange(chi2_min,chi2_max,"Y");
1023  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->Draw();
1024  /*Chi2->SetLineColor(1);
1025  Chi2->SetLineWidth(1);
1026  Periphery_Chi2->SetLineColor(1);
1027  Periphery_Chi2->SetLineWidth(1);
1028  Chi2->Draw("same");
1029  Periphery_Chi2->Draw("same");*/
1030  c->Print(detector + "_chi2E_ge_EtaCut_CircularCut.png");
1031 
1032  gStyle -> SetOptStat(0);
1033  gStyle -> SetOptFit(0);
1034 
1035  /*Sigma->SetLineColor(1);
1036  Sigma->SetLineWidth(1);
1037  Sigma_error->SetLineColor(1);
1038  Sigma_error->SetLineWidth(1);
1039  Point->SetLineColor(3);
1040  Point->SetLineWidth(5);*/
1041  fExp->SetLineColor(4); //38
1042  fTrue->SetLineColor(2); //46
1043 
1044  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerStyle(kFullCircle);
1045  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerColor(46); //30
1046  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerSize(0.75);
1047  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetAxisRange(sigma_min,sigma_max,"Y");
1048 
1049  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->Fit("fTrue", "M+");
1050  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->Draw("same");
1051  fExp->Draw("same");
1052  /*Sigma->Draw("same");
1053  Sigma_error->Draw("same");
1054  Point->Draw("same");*/
1055 
1056  TLegend* legend = new TLegend(1.75,1.75);
1057  legend->SetHeader("Legend", "C");
1058  //legend->AddEntry(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2, "#sigma_{e_{agg}} vs Generated Energy", "flep");
1059  //legend->AddEntry((TObject*)0,"","");
1060  legend->AddEntry(fTrue, "p_{0} + p_{1}/#sqrt{ge} (Fitted)", "l");
1061  legend->AddEntry((TObject*)0,"","");
1062  legend->AddEntry(fExp, "0.1 + 1.0/#sqrt{ge} (Requirement)", "l");
1063  legend->SetTextSize(0.033);
1064  legend->Draw();
1065 
1066  std::cout<<"reduced_chi2 of fit: "<<fTrue->GetChisquare()/fTrue->GetNDF()<<"\n";
1067 
1068  c->Print(detector + "_sigmaE_ge_EtaCut_CircularCut.png");
1069 
1070  gStyle -> SetOptStat(11);
1071  gStyle -> SetOptFit(112);
1072 
1073  for(int sno = 0; sno < nSlicesx; sno++){
1074  TString plusOne = (TString)(sno + 1);
1075  TString nameF = detector + "_sigmaE_slice" + arr[sno] + "_EtaCut_CircularCut.png";
1076  slices[sno] -> Fit("fit", "M+");
1077  slices[sno] -> Draw("hist same");
1078  c->Print(nameF);
1079  }
1080 
1081  gStyle -> SetOptStat(1);
1082 
1083  te_by_ge_ge_EtaCut->Draw("colz");
1084  c->Print(detector + "_te_by_ge_ge_EtaCut.png");
1085  te_by_ge_ge_EtaCut_CircularCut->Draw("colz");
1086  c->Print(detector + "_te_by_ge_ge_EtaCut_CircularCut.png");
1087 
1088  te_minus_ge_by_ge_ge_EtaCut->Draw("colz");
1089  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut.png");
1090  te_minus_ge_by_ge_ge_EtaCut_CircularCut->Draw("colz");
1091  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut_CircularCut.png");
1092  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Draw("colz");
1093  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated.png");
1094 
1095  gStyle -> SetOptStat(0);
1096 
1097  mean_te_by_ge_ge_EtaCut_CircularCut->Draw();
1098  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut.png");
1099 
1100  mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN->Draw();
1101  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut_HCALIN.png");
1102 
1103  mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT->Draw();
1104  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut_HCALOUT.png");
1105 
1106  mean_te_by_ge_ge_EtaCut_CircularCut_CEMC->Draw();
1107  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut_CEMC.png");
1108 
1109  gStyle -> SetOptStat(1);
1110  c->Close();
1111 
1112  }
1113 
1114  std::cout << "The total te is: " << total_te << endl;
1115  std::cout << "The total te_CircularCut is: " << total_te_CircularCut << endl;
1116  std::cout << "The total ge is: " << total_ge << endl;
1117  std::cout<<"\n\nDone\n----------------------------------------------------------------------\n\n";
1118 
1119 }