EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LoopEvalFR.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LoopEvalFR.C
1 /*
2 > LoopEvalFR.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 FEMC+FHCAL (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 FEMC alone and FEMC+FHCAL 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 FEMC+FHCAL
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 FEMC to eliminate MIPs, based on the parametric equation provided in the code
11 - Output file - energy_verification_EtaCut_CircularCut_FEMC_FHCAL.root, and the .png plots if generated
12 */
13 
14 /*
15  Authors - Siddhant Rathi (me190003061@iiti.ac.in)
16  Sagar Joshi (ee190002054@iiti.ac.in)
17 
18 
19  version - 2
20 
21 */
22 
23 #include <iostream>
24 #include <stdexcept>
25 #include <eicqa_modules/EvalRootTTree.h>
26 #include <eicqa_modules/EvalHit.h>
27 #include "TMath.h"
28 #include "TStyle.h"
29 #include <unistd.h>
30 
31 R__LOAD_LIBRARY(libeicqa_modules.so)
32 
33 void LoopEvalFR(int print = 1, int debug = 0, Double_t energyCutAggregate = 0.1, Double_t energyCut = 0.0, int MIP_theta_parametrisation = 1)
34 {
35 
36  Double_t EMC_cut = 0.0;
37  //TF1 *mip_pmzn_energy_cut_ftheta = new TF1("mip_pmzn_energy_cut_ftheta", "(2.49951e-02) + (2.15109e-01)*x - (1.15673)*(x^2) + (2.75562)*(x^3) - (2.31809)*(x^4)");
38  TF1 *mip_pmzn_energy_cut_ftheta = new TF1("mip_pmzn_energy_cut_ftheta", "0");
39 
40  if(MIP_theta_parametrisation == 1 && energyCut != 0){
41  throw std::invalid_argument("We do not currently support theta-parametrized \nMIP cut on EMC simultaneously with individual tower cuts \non other detectors.:(;");
42  }
43 
44  //variable binning code
45  int arraySizeBins = 0; //size of the bins array
46  Double_t arraySizeBinsLoop = 0; //loop variable
47  int lowThresholdBins = 2; //number of bins below the threshold bin
48  int highThresholdBins = 9; //number of bins above the threshold bin
49  Double_t maxEnergyBin = 30; //highest energy
50  Double_t lowEnergyBin = 3; //threshold energy
51 
52  //calculating the number of bins required
53  while(arraySizeBinsLoop < maxEnergyBin){
54  arraySizeBins += 1;
55  if(arraySizeBinsLoop < lowEnergyBin){
56  arraySizeBinsLoop += lowEnergyBin/lowThresholdBins;
57  }
58  else{
59  arraySizeBinsLoop += (maxEnergyBin - lowEnergyBin)/highThresholdBins;
60  }
61  }
62 
63  Double_t binLimitArray[arraySizeBins + 1];
64 
65  //filing the array containing bin lower limits
66  for(int i = 0; i <= arraySizeBins; i++){
67  if(i <= lowThresholdBins){
68  binLimitArray[i] = lowEnergyBin*i/lowThresholdBins;
69  }
70  else{
71  binLimitArray[i] = lowEnergyBin + ((maxEnergyBin - lowEnergyBin)*(i - lowThresholdBins)/highThresholdBins);
72  }
73  }
74 
75  //print array
76  for(int i = 0; i < arraySizeBins; i++){
77  cout<<"element "<<i<<" of the bin array is "<<binLimitArray[i]<<" \n";
78  }
79 
80  //variable binning code end
81 
82  TString detector = "FEMC_FHCAL";
83  TFile *f1 = new TFile("merged_Eval_FHCAL.root","READ");
84  TFile *f2 = new TFile("merged_Eval_FEMC.root","READ");
85  //TFile *f1 = new TFile("../../../simrankaur/Fun4All-Calorimeter-Macros_copy/final/Pion_fhcalcorrected/Eval_FHCAL.root","READ");
86  //TFile *f2 = new TFile("../../../simrankaur/Fun4All-Calorimeter-Macros_copy/final/Pion_fhcalcorrected/Eval_FEMC.root","READ");
87  //TFile *f1 = new TFile("../../July21/pi_30_2gev/merged_Eval_FHCAL.root","READ");
88  //TFile *f2 = new TFile("../../July21/pi_30_2gev/merged_Eval_FEMC.root","READ");
89 
90  TTree* T1 = (TTree*)f1->Get("T");
91  EvalRootTTree *evaltree1 = nullptr;
92 
93  TTree* T2 = (TTree*)f2->Get("T");
94  EvalRootTTree *evaltree2 = nullptr;
95 
96  gStyle->SetCanvasPreferGL(kTRUE);
97  TCanvas *c = new TCanvas();
98  c->SetTickx();
99  c->SetTicky();
100 
101  // Modifying default plotting style
102  gStyle->SetOptTitle(0);
103  gStyle->SetOptFit(102);
104  gStyle->SetTitleXOffset(1);
105  gStyle->SetTitleYOffset(1);
106  gStyle->SetLabelSize(0.05);
107  gStyle->SetTitleXSize(0.05);
108  gStyle->SetTitleYSize(0.05);
109 
110  long double total_te = 0;
111  long double total_te_CircularCut = 0;
112  long double total_ge = 0;
113 
114  int nSlicesx = arraySizeBins; // Number of ge-axis slices taken for making sigma_e vs ge plot
115  int nSlicesy = 350;
116  Double_t eta_min, eta_max; // Eta range of detectors
117  Double_t x_radius_FHCAL = 0.15; //0.15; Length of semi-minor axis of circular (elliptical) in FHCAL
118  Double_t y_radius_FHCAL = 0.45; //0.45; Length of semi-major axis of circular (elliptical) in FHCAL
119  Double_t x_radius_FEMC = 0.13; //0.13; Length of semi-minor axis of circular (elliptical) in FEMC
120  Double_t y_radius_FEMC = 0.35; //0.35; Length of semi-major axis of circular (elliptical) in FEMC
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 = 20.746;
127  Double_t recalibration_firstSlice_FHCAL = 0.2079; // Recalibration factor of first slice (needed to be done manually because of low statistics)
128  Double_t recalibration_firstSlice_FEMC = 0.1498;
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 = 2;
132  fit_max = 15;
133  eta_min = 1.4;
134  eta_max = 3.0;
135  sigma_min = 0.0;
136  sigma_max = 1.0;
137  mean_min = -0.5;
138  mean_max = 0.5;
139  chi2_min = 0;
140  chi2_max = 2.52;
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 = " {1.4 < geta < 3.0} ";
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 
153  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_FEMC = new TH2D("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_FEMC","#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);
154  TH2D *te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_FHCAL = new TH2D("te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_FHCAL","#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);
155 
156  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
157 
158  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
159 
160  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);
161  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,1.5);
162 
163  TH2D *te_by_ge_ge_EtaCut_CircularCut_FHCAL = new TH2D("te_by_ge_ge_EtaCut_CircularCut_FHCAL","te_{agg}/ge vs ge",200,0,30,200,-1,2);
164  TH2D *te_by_ge_ge_EtaCut_CircularCut_FEMC = new TH2D("te_by_ge_ge_EtaCut_CircularCut_FEMC","te_{agg}/ge vs ge",200,0,30,200,-1,2);
165 
166  TH1D *te_aggregate_EtaCut_CircularCut_FEMC = new TH1D("te_aggregate_EtaCut_CircularCut_FEMC","",200,0,1);
167 
168  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);
169  auto *mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL","Mean_{te/ge}",nSlicesx,binLimitArray,-0.5,35);
170  auto *mean_te_by_ge_ge_EtaCut_CircularCut_FEMC = new TProfile("mean_te_by_ge_ge_EtaCut_CircularCut_FEMC","Mean_{te/ge}",nSlicesx,binLimitArray,-0.5,35);
171 
172  T1->SetBranchAddress("DST#EvalTTree_FHCAL",&evaltree1);
173  T2->SetBranchAddress("DST#EvalTTree_FEMC",&evaltree2);
174 
175  for(int i=0; i<T1->GetEntries(); i++) // We assume same no. of entries, since no cuts are applied
176  {
177  T1->GetEntry(i);
178  T2->GetEntry(i);
179 
180  if(debug==1){
181  std::cout<<"\n\n\n------------------------------------------\nParticle: "<<i<<"\n\n";
182  std::cout<<"Initial Parameters "<<"\n";
183  }
184 
185  Double_t geta1 = evaltree1->get_geta();
186 
187  if(debug==1){
188  std::cout<<"geta: "<<geta1<<"\n";
189  }
190 
191  if(geta1>=eta_min && geta1<=eta_max){
192 
193  if(debug==1){
194  cout<<"\ngeta cut applied (1.3, 3.3)"<<"\n\n";
195  }
196  Double_t ge = evaltree1->get_ge();
197  if(debug==1){
198  std::cout<<"ge: "<<ge<<"\n";
199  }
200 
201  Double_t gphi = evaltree1->get_gphi();
202  if(debug==1){
203  std::cout<<"gphi: "<<gphi<<"\n";
204  }
205 
206  Double_t gtheta = evaltree1->get_gtheta();
207  if(debug==1){
208  std::cout<<"gtheta: "<<gtheta<<"\n";
209  }
210 
211  total_ge += ge;
212  if(debug==1){
213  std::cout<<"total_ge till now = "<<total_ge<<"\n";
214  }
215 
216  Double_t te_aggregate = 0;
217  Double_t te_aggregate_CircularCut = 0;
218  Double_t te_aggregate_FHCAL_CircularCut = 0;
219 
220  for (int j=0; j<evaltree1->get_ntowers(); j++){
221 
222  if(debug==1){
223  std::cout<<"\nFHCAL Tower: "<<j<<"\n";
224  }
225 
226  EvalTower *twr1 = evaltree1->get_tower(j);
227  if (twr1){
228 
229  if(debug==1){
230  cout<<"non-empty FHCAL tower\n";
231  }
232 
233  Double_t tphi = twr1->get_tphi();
234  if(debug==1){
235  std::cout<<"tphi: "<<tphi<<"\n";
236  }
237 
238  Double_t ttheta = twr1->get_ttheta();
239  if(debug==1){
240  std::cout<<"ttheta: "<<ttheta<<"\n";
241  }
242 
243  Double_t dphi = tphi - gphi;
244  if(debug==1){
245  std::cout<<"tphi-gphi: "<<dphi<<"\n";
246  }
247 
248  Double_t dtheta = ttheta - gtheta;
249  if(debug==1){
250  std::cout<<"ttheta-gtheta: "<<dtheta<<"\n";
251  }
252 
253  Double_t te = twr1->get_te();
254  if(debug==1){
255  std::cout<<"te: "<<te<<"\n";
256  }
257 
258  if(te > energyCut){
259  te_aggregate += te;
260 
261  if(debug==1){
262  std::cout<<"FHCAL: pow(dphi/y_radius,2)+pow(dtheta/x_radius,2) = "<<pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2)<<"\n";
263  }
264 
265  if (pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2)<=1){
266  if(debug==1){
267  cout<<"FHCAL Tower included after circular cut\n";
268  }
269  te_aggregate_CircularCut += te;
270  te_aggregate_FHCAL_CircularCut += te;
271  }
272 
273  if(debug==1){
274  cout<<"te_aggregate till now = "<<te_aggregate<<"\n";
275  std::cout<<"te += "<<twr1->get_te()<<"\n";
276  cout<<"te_aggregate_CircularCut till now = "<<te_aggregate_CircularCut<<"\n";
277  }
278  }
279  }
280  }
281 
282  for (int j=0; j<evaltree2->get_ntowers(); j++){
283 
284 
285  if(debug==1){
286  std::cout<<"\nFEMC Tower: "<<j<<"\n";
287  }
288 
289  EvalTower *twr2 = evaltree2->get_tower(j);
290  if (twr2){
291 
292  if(debug==1){
293  cout<<"non-empty FEMC tower\n";
294  }
295 
296  Double_t tphi = twr2->get_tphi();
297  if(debug==1){
298  std::cout<<"tphi: "<<tphi<<"\n";
299  }
300 
301  Double_t ttheta = twr2->get_ttheta();
302  if(debug==1){
303  std::cout<<"ttheta: "<<ttheta<<"\n";
304  }
305 
306  Double_t dphi = tphi - gphi;
307  if(debug==1){
308  std::cout<<"tphi-gphi: "<<dphi<<"\n";
309  }
310 
311  Double_t dtheta = ttheta - gtheta;
312  if(debug==1){
313  std::cout<<"ttheta-gtheta: "<<dtheta<<"\n";
314  }
315 
316  Double_t te = twr2->get_te();
317  if(debug==1){
318  std::cout<<"te: "<<te<<"\n";
319  }
320 
321  if(MIP_theta_parametrisation == 1){
322  EMC_cut = mip_pmzn_energy_cut_ftheta->Eval(gtheta);
323  }
324 
325  if(te > energyCut + EMC_cut){
326  te_aggregate += te;
327 
328  if(debug==1){
329  std::cout<<"FEMC: pow(dphi/y_radius,2)+pow(dtheta/x_radius,2): "<<pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2)<<"\n";
330  }
331 
332  if (pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2)<=1){
333  if(debug==1){
334  cout<<"FEMC Tower included after circular cut\n";
335  }
336  te_aggregate_CircularCut += te;
337  }
338 
339  if(debug==1){
340  cout<<"te_aggregate till now = "<<te_aggregate<<"\n";
341  cout<<"te_aggregate_CircularCut till now = "<<te_aggregate_CircularCut<<"\n";
342  std::cout<<"te += "<<twr2->get_te()<<"\n";
343  }
344  }
345  }
346  }
347 
348  total_te += te_aggregate;
349  total_te_CircularCut += te_aggregate_CircularCut;
350 
351  if(debug==1){
352  cout<<"total_te till now = "<<total_te<<"\n";
353  cout<<"total_te_CircularCut till now = "<<total_te_CircularCut<<"\n\n";
354  }
355 
356  if(te_aggregate_CircularCut > energyCutAggregate){
357 
358  if(te_aggregate_CircularCut < 0.1) {
359 
360  cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
361 
362  }
363 
364  te_minus_ge_by_ge_ge_EtaCut->Fill(ge, (te_aggregate-ge)/ge);
365 
366  if(debug==1){
367  cout<<"(ge, (te_aggregate-ge)/ge): ("<<ge<<", "<<(te_aggregate-ge)/ge<<")\n";
368  }
369 
370  te_minus_ge_by_ge_ge_EtaCut_CircularCut->Fill(ge, (te_aggregate_CircularCut-ge)/ge);
371  if(debug==1){
372  cout<<"(ge, (te_aggregate_CircularCut-ge)/ge): ("<<ge<<", "<<(te_aggregate_CircularCut-ge)/ge<<")\n";
373  }
374 
375 
376  /*te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Fill(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge);
377  if(debug==1){
378  cout<<"(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge: ("<<ge<<", "<<((te_aggregate_CircularCut/recalibration_factor)-ge)/ge<<")\n";
379  }
380  te_minus_ge_by_ge_ge_EtaCut_Recalibrated_temp->Fill(ge, ((te_aggregate_CircularCut/recalibration_factor)-ge)/ge);*/
381 
382  te_by_ge_ge_EtaCut->Fill(ge, te_aggregate/ge);
383 
384  if(debug==1){
385  cout<<"(ge, te_aggregate/ge): ("<<ge<<", "<<te_aggregate/ge<<")\n";
386  }
387 
388  te_by_ge_ge_EtaCut_CircularCut->Fill(ge, te_aggregate_CircularCut/ge);
389 
390  te_aggregate_EtaCut_CircularCut_FEMC->Fill(te_aggregate_CircularCut - te_aggregate_FHCAL_CircularCut);
391 
392  te_by_ge_ge_EtaCut_CircularCut_FHCAL->Fill(ge, te_aggregate_FHCAL_CircularCut/ge);
393  te_by_ge_ge_EtaCut_CircularCut_FEMC->Fill(ge, (te_aggregate_CircularCut - te_aggregate_FHCAL_CircularCut)/ge);
394 
395  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->Fill(ge, te_aggregate_FHCAL_CircularCut/ge);
396  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->Fill(ge, (te_aggregate_CircularCut-te_aggregate_FHCAL_CircularCut)/ge);
397 
398 
399  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_FHCAL->Fill(ge, (te_aggregate_FHCAL_CircularCut-ge)/ge);
400  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_FEMC->Fill(ge, (te_aggregate_CircularCut - te_aggregate_FHCAL_CircularCut-ge)/ge);
401 
402 
403  if(debug==1){
404  cout<<"(ge, te_aggregate_CircularCut/ge): ("<<ge<<", "<<te_aggregate_CircularCut/ge<<")\n";
405  cout<<"(ge, te_aggregate_FHCAL_CircularCut/ge): ("<<ge<<", "<<te_aggregate_FHCAL_CircularCut/ge<<")\n";
406  cout<<"(ge, (te_aggregate_CircularCut-te_aggregate_FHCAL_CircularCut)/ge): ("<<ge<<", "<<(te_aggregate_CircularCut-te_aggregate_FHCAL_CircularCut)/ge<<")\n";
407  }
408  }
409  }
410  }
411 
412  if(debug==1){
413  cout<<"\neta cut if ends"<<"\n";
414  }
415 
416 
417 
418  TString arr[nSlicesx]; // Used for naming fitted slices used in sigma_e vs ge
419  for(int sno = 0; sno < nSlicesx; sno++){
420  arr[sno] = TString::Itoa(sno + 1, 10);
421  }
422 
423  if(debug==1){
424  std::cout<<"\nGenerating sigma_e vs ge plots\n\n";
425  }
426 
427  Double_t recalibrationArr1[nSlicesx];
428  Double_t rf_integral_FHCAL = 0;
429  Double_t weight_FHCAL = te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetMean(2);
430 
431  //recalibrationArr1[0] = recalibration_firstSlice_FHCAL;
432  //rf_integral_FHCAL += recalibrationArr1[0];
433  //cout << "Recalibration factor for slice " << 1 << " of FHCAL is: " << recalibrationArr1[0] << endl;
434 
435  for(int binIter = 1; binIter <= nSlicesx; binIter++){
436  recalibrationArr1[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetBinContent(binIter);
437  rf_integral_FHCAL += recalibrationArr1[binIter-1];
438 
439  if(debug == 1){
440  cout<<"rf_integral_FHCAL += "<<recalibrationArr1[binIter-1]<<"\n";
441  }
442 
443  cout << "Recalibration factor for slice " << binIter << " of FHCAL is: " << recalibrationArr1[binIter-1] << endl;
444  }
445 
446  if(debug == 1){
447  cout<<"rf_integral_FHCAL = "<<rf_integral_FHCAL<<"\n\n";
448  }
449 
450  Double_t recalibrationArr2[nSlicesx];
451  Double_t rf_integral_FEMC = 0;
452  Double_t weight_FEMC = te_by_ge_ge_EtaCut_CircularCut_FEMC->GetMean(2);
453  //recalibrationArr2[0] = recalibration_firstSlice_FEMC;
454  //rf_integral_FEMC += recalibrationArr2[0];
455  //cout << "Recalibration factor for slice " << 1 << " of FEMC is: " << recalibrationArr2[0] << endl;
456 
457  for(int binIter = 1; binIter <= nSlicesx; binIter++){
458  recalibrationArr2[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetBinContent(binIter);
459  rf_integral_FEMC += recalibrationArr2[binIter-1];
460 
461  if(debug == 1){
462  cout<<"rf_integral_FEMC += "<<recalibrationArr2[binIter-1]<<"\n";
463  }
464 
465  cout << "Recalibration factor for slice " << binIter << " of FEMC is: " << recalibrationArr2[binIter-1] << endl;
466  }
467 
468  if(debug == 1){
469  cout<<"rf_integral_FEMC = "<<rf_integral_FEMC<<"\n\n";
470  }
471 
472  for(int i=0; i<T1->GetEntries(); i++){
473 
474  T1->GetEntry(i);
475  T2->GetEntry(i);
476 
477  Double_t geta1 = evaltree1->get_geta();
478  Double_t gphi = evaltree1->get_gphi();
479  Double_t gtheta = evaltree1->get_gtheta();
480  Double_t ge = evaltree1->get_ge();
481  Double_t te_aggregate_CircularCut = 0;
482  Double_t te_aggregate_FHCAL_CircularCut = 0;
483  Double_t te_aggregate_CircularCut_normalised = 0;
484 
485  int recalibration_factor = 0;
486 
487 
488  if(ge > lowEnergyBin){
489  Double_t eRangeBin = (maxEnergyBin - lowEnergyBin)/highThresholdBins;
490  recalibration_factor = (lowThresholdBins - 1) + ceil(ge/eRangeBin)- 1;
491  }
492 
493  else{
494  recalibration_factor = ceil((ge/lowEnergyBin)*(Double_t)lowThresholdBins)- 1;
495  }
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_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2): "<<pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2)<<"\n";
515  }
516 
517  if (pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2)<=1){
518  te_aggregate_CircularCut += te;
519  te_aggregate_FHCAL_CircularCut += te;
520  te_aggregate_CircularCut_normalised += te*weight_FHCAL/recalibrationArr1[recalibration_factor];
521  }
522  }
523  }
524  }
525 
526  for (int j=0; j<evaltree2->get_ntowers(); j++){
527 
528  EvalTower *twr2 = evaltree2->get_tower(j);
529  if (twr2){
530 
531  Double_t tphi = twr2->get_tphi();
532  Double_t ttheta = twr2->get_ttheta();
533  Double_t dphi = tphi - gphi;
534  Double_t dtheta = ttheta - gtheta;
535  Double_t te = twr2->get_te();
536 
537  if(MIP_theta_parametrisation == 1){
538  EMC_cut = mip_pmzn_energy_cut_ftheta->Eval(gtheta);
539  }
540 
541  if(te > energyCut + EMC_cut){
542 
543  if(debug==1){
544  std::cout<<"pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2): "<<pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2)<<"\n";
545  }
546 
547  if (pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2)<=1){
548  te_aggregate_CircularCut += te;
549  te_aggregate_CircularCut_normalised += te*weight_FEMC/recalibrationArr2[recalibration_factor];
550  }
551  }
552  }
553  }
554 
555 
556  if(te_aggregate_CircularCut > energyCutAggregate){
557  mean_te_by_ge_ge_EtaCut_CircularCut->Fill(ge, te_aggregate_CircularCut_normalised/ge);
558  if(debug == 1){
559  cout<<"(ge, te_aggregate_CircularCut_normalised/ge): ("<<ge<<", "<<te_aggregate_CircularCut_normalised/ge<<")\n";
560  }
561  }
562  }
563  }
564 
565 
566  Double_t recalibrationArr[nSlicesx];
567 
568  //recalibrationArr[0] = recalibration_firstSlice;
569  //cout << "Recalibration factor for slice " << 1 << " is: " << recalibrationArr[0] << endl;
570 
571  for(int binIter = 1; binIter <= nSlicesx; binIter++){
572  recalibrationArr[binIter-1] = mean_te_by_ge_ge_EtaCut_CircularCut->GetBinContent(binIter);
573  cout << "Recalibration factor for slice " << binIter << " of is: " << recalibrationArr[binIter-1] << endl;
574  }
575 
576 
577  for(int i=0; i<T1->GetEntries(); i++){
578 
579  T1->GetEntry(i);
580  T2->GetEntry(i);
581 
582  Double_t geta1 = evaltree1->get_geta();
583  Double_t gphi = evaltree1->get_gphi();
584  Double_t gtheta = evaltree1->get_gtheta();
585  Double_t ge = evaltree1->get_ge();
586  Double_t te_aggregate_CircularCut = 0.0;
587  Double_t te_aggregate_FHCAL_CircularCut = 0.0;
588  Double_t te_aggregate_CircularCut_normalised = 0.0;
589 
590  int recalibration_factor = 0;
591 
592 
593  if(ge > lowEnergyBin){
594  Double_t eRangeBin = (maxEnergyBin - lowEnergyBin)/highThresholdBins;
595  recalibration_factor = (lowThresholdBins - 1) + ceil(ge/eRangeBin)- 1;
596  }
597 
598  else{
599  recalibration_factor = ceil((ge/lowEnergyBin)*(Double_t)lowThresholdBins)- 1;
600  }
601 
602 
603  if(geta1>=eta_min && geta1<=eta_max){
604 
605 
606  for (int j=0; j<evaltree1->get_ntowers(); j++){
607 
608  EvalTower *twr1 = evaltree1->get_tower(j);
609  if (twr1){
610 
611  Double_t tphi = twr1->get_tphi();
612  Double_t ttheta = twr1->get_ttheta();
613  Double_t dphi = tphi - gphi;
614  Double_t dtheta = ttheta - gtheta;
615  Double_t te = twr1->get_te();
616 
617  if(te > energyCut){
618 
619  if(debug==1){
620  std::cout<<"pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2): "<<pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2)<<"\n";
621  }
622 
623  if (pow(dphi/y_radius_FHCAL,2)+pow(dtheta/x_radius_FHCAL,2)<=1){
624  te_aggregate_CircularCut += te;
625  te_aggregate_FHCAL_CircularCut += te;
626  te_aggregate_CircularCut_normalised += te*weight_FHCAL/recalibrationArr1[recalibration_factor];
627  }
628  }
629  }
630  }
631 
632  for (int j=0; j<evaltree2->get_ntowers(); j++){
633 
634  EvalTower *twr2 = evaltree2->get_tower(j);
635  if (twr2){
636 
637  Double_t tphi = twr2->get_tphi();
638  Double_t ttheta = twr2->get_ttheta();
639  Double_t dphi = tphi - gphi;
640  Double_t dtheta = ttheta - gtheta;
641  Double_t te = twr2->get_te();
642 
643  if(MIP_theta_parametrisation == 1){
644  EMC_cut = mip_pmzn_energy_cut_ftheta->Eval(gtheta);
645  }
646 
647  if(te > energyCut + EMC_cut){
648 
649  if(debug==1){
650  std::cout<<"pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2): "<<pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2)<<"\n";
651  }
652 
653  if (pow(dphi/y_radius_FEMC,2)+pow(dtheta/x_radius_FEMC,2)<=1){
654  te_aggregate_CircularCut += te;
655  te_aggregate_CircularCut_normalised += te*weight_FEMC/recalibrationArr2[recalibration_factor];
656  }
657  }
658  }
659  }
660 
661 
662  if(te_aggregate_CircularCut > energyCutAggregate){
663  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Fill(ge, ((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge);
664  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->Fill(ge, ((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge);
665 
666  if(debug==1){
667  cout<<"(ge, ((te_aggregate_CircularCut_normalised/recalibration_factor)-ge)/ge: ("<<ge<<", "<<((te_aggregate_CircularCut_normalised/recalibrationArr[recalibration_factor])-ge)/ge<<")\n";
668  }
669  }
670  }
671  }
672 
673 
674  // Initialising fit functions
675  TF1 *fit = new TF1("fit", "gaus");
676  TF1 *fit1 = new TF1("fit1","gaus",fit_min,fit_max);
677  TF1 *fExp = new TF1("fExp","0.1 + 0.5/sqrt(x)",0,30);
678  TF1 *fTrue = new TF1("fTrue","[0] + [1]/sqrt(x)",0,30);
679 
680 
681  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->FitSlicesY(0, 1, -1, 0, "QN");
682  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");
683  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");
684  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");
685 
686 
687 
688  TH1D* slices[nSlicesx];
689 
690  // Generating plots for individual slices
691  for(int sno = 0; sno < nSlicesx; sno++){
692  int plusOne = sno+1;
693  TString sname = "slice " + arr[sno];
694  slices[sno] = te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp->ProjectionY(sname, plusOne, plusOne);
695  }
696 
697  if(debug==1){
698  cout<<"\nHistogram Formatting\n\n";
699  }
700 
701  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
702  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetLabelSize(0.05);
703  te_minus_ge_by_ge_ge_EtaCut->GetXaxis()->SetTitleSize(0.05);
704  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
705  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetLabelSize(0.05);
706  te_minus_ge_by_ge_ge_EtaCut->GetYaxis()->SetTitleSize(0.05);
707 
708  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
709  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
710  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
711  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
712  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
713  te_minus_ge_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
714 
715  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetTitle("Generated Energy (GeV)");
716  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetLabelSize(0.05);
717  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetXaxis()->SetTitleSize(0.05);
718  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetTitle("(te_{agg}-ge)/ge");
719  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetLabelSize(0.05);
720  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->GetYaxis()->SetTitleSize(0.05);
721 
722  te_by_ge_ge_EtaCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
723  te_by_ge_ge_EtaCut->GetXaxis()->SetLabelSize(0.05);
724  te_by_ge_ge_EtaCut->GetXaxis()->SetTitleSize(0.05);
725  te_by_ge_ge_EtaCut->GetYaxis()->SetTitle("te_{agg}/ge");
726  te_by_ge_ge_EtaCut->GetYaxis()->SetLabelSize(0.05);
727  te_by_ge_ge_EtaCut->GetYaxis()->SetTitleSize(0.05);
728 
729  te_aggregate_EtaCut_CircularCut_FEMC->GetXaxis()->SetTitle("te_{agg} (GeV)");
730  te_aggregate_EtaCut_CircularCut_FEMC->GetXaxis()->SetLabelSize(0.05);
731  te_aggregate_EtaCut_CircularCut_FEMC->GetXaxis()->SetTitleSize(0.05);
732  te_aggregate_EtaCut_CircularCut_FEMC->GetYaxis()->SetTitle("Counts");
733  te_aggregate_EtaCut_CircularCut_FEMC->GetYaxis()->SetLabelSize(0.05);
734  te_aggregate_EtaCut_CircularCut_FEMC->GetYaxis()->SetTitleSize(0.05);
735 
736  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
737  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
738  te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
739  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("te_{agg}/ge");
740  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
741  te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
742 
743  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitle("Generated Energy (GeV)");
744  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetLabelSize(0.05);
745  mean_te_by_ge_ge_EtaCut_CircularCut->GetXaxis()->SetTitleSize(0.05);
746  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitle("Mean of te_{agg}/ge");
747  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetLabelSize(0.05);
748  mean_te_by_ge_ge_EtaCut_CircularCut->GetYaxis()->SetTitleSize(0.05);
749 
750  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetXaxis()->SetTitle("Generated Energy (GeV)");
751  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetXaxis()->SetLabelSize(0.05);
752  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetXaxis()->SetTitleSize(0.05);
753  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetYaxis()->SetTitle("Mean of te_{agg}/ge (FHCAL)");
754  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetYaxis()->SetLabelSize(0.05);
755  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->GetYaxis()->SetTitleSize(0.05);
756 
757  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetXaxis()->SetTitle("Generated Energy (GeV)");
758  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetXaxis()->SetLabelSize(0.05);
759  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetXaxis()->SetTitleSize(0.05);
760  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetYaxis()->SetTitle("Mean of te_{agg}/ge (FEMC)");
761  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetYaxis()->SetLabelSize(0.05);
762  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->GetYaxis()->SetTitleSize(0.05);
763 
764  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetTitle("Generated Energy (GeV)");
765  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetLabelSize(0.05);
766  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetXaxis()->SetTitleSize(0.05);
767  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetTitle("#sigma_{e_{agg}}");
768  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetLabelSize(0.05);
769  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->GetYaxis()->SetTitleSize(0.05);
770  //te_minus_ge_by_ge_ge_EtaCut_2->SetTitle("#sigma_{e_{agg}} vs true_e" + cut_text);
771 
772  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetTitle("Generated Energy (GeV)");
773  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetLabelSize(0.05);
774  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetXaxis()->SetTitleSize(0.05);
775  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetTitle("Reduced_#chi^{2}_{e_{agg}}");
776  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetLabelSize(0.05);
777  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->GetYaxis()->SetTitleSize(0.04);
778  //te_minus_ge_by_ge_ge_EtaCut_chi2->SetTitle("#chi^{2}_{e_{agg}} vs true_e" + cut_text);
779 
780  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetTitle("Generated Energy (GeV)");
781  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetLabelSize(0.05);
782  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetXaxis()->SetTitleSize(0.05);
783  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetTitle("Mean_{e_{agg}}");
784  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetLabelSize(0.05);
785  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->GetYaxis()->SetTitleSize(0.05);
786  //te_minus_ge_by_ge_ge_EtaCut_1->SetTitle("mean_{e_{agg}} vs true_e" + cut_text);
787 
788  for(int sno = 0; sno < nSlicesx; sno++){
789  slices[sno]->GetXaxis()->SetTitle("#Delta e^{agg}/ ge");
790  slices[sno]->GetXaxis()->SetLabelSize(0.05);
791  slices[sno]->GetXaxis()->SetTitleSize(0.05);
792  slices[sno]->GetYaxis()->SetTitle("Counts");
793  slices[sno]->GetYaxis()->SetLabelSize(0.05);
794  slices[sno]->GetYaxis()->SetTitleSize(0.05);
795  }
796 
797  if(debug==1){
798  std::cout<<"\nWrite Histograms to File\n";
799  }
800 
801  TFile *f = new TFile("energy_verification_EtaCut_CircularCut_" + detector + ".root","RECREATE");
802 
803  f->GetList()->Add(te_by_ge_ge_EtaCut);
804  f->GetList()->Add(te_by_ge_ge_EtaCut_CircularCut);
805  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut);
806  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut);
807  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut);
808  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL);
809  f->GetList()->Add(mean_te_by_ge_ge_EtaCut_CircularCut_FEMC);
810  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated);
811  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2);
812  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1);
813  f->GetList()->Add(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2);
814  f->GetList()->Add(te_aggregate_EtaCut_CircularCut_FEMC);
815 
816  for(int sno = 0; sno < nSlicesx; sno++){
817  f->GetList()->Add(slices[sno]);
818  }
819  f->Write();
820 
821  gStyle -> SetOptStat(11);
822  gStyle -> SetOptFit(112);
823 
824  /* int sno = 0;
825  TString plusOne = (TString)(sno + 1);
826  TString nameF = detector + "_sigmaE_slice" + arr[sno] + "_EtaCut_CircularCut.png";
827  slices[sno] -> Fit("fit1", "R+");
828  slices[sno] -> Draw("hist same");
829  c->Print(nameF);
830 
831  double_t mean = fit1->GetParameter(1);
832  double_t mean_error = fit1->GetParError(1);
833  double_t sigma = fit1->GetParameter(2);
834  double_t sigma_error = fit1->GetParError(2);
835  double_t chi2 = (fit1->GetChisquare())/(fit1->GetNDF());
836 
837  TLine *Mean = new TLine(0,mean,30.0/nSlicesx,mean);
838  TLine *Sigma = new TLine(0,sigma,30.0/nSlicesx,sigma);
839  TLine *Chi2 = new TLine(0,chi2,30.0/nSlicesx,chi2);
840  TLine *Periphery_Chi2 = new TLine(30.0/nSlicesx,0,30.0/nSlicesx,chi2);
841  TLine *Sigma_error = new TLine(30.0/(2.0*nSlicesx),sigma-(sigma_error),30.0/(2.0*nSlicesx),sigma+(sigma_error));
842  TLine *Mean_error = new TLine(30.0/(2.0*nSlicesx),mean-(mean_error),30.0/(2.0*nSlicesx),mean+(mean_error));
843  TLine *Point = new TLine(30.0/(2.0*nSlicesx),sigma,30.0/(2.0*nSlicesx),sigma);*/
844 
845  if(print==1){
846  if(debug==1){
847  std::cout<<"\nSaving Histograms as .png\n";
848  }
849 
850  gStyle -> SetOptStat(0);
851 
852  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->SetAxisRange(mean_min,mean_max,"Y");
853  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_1->Draw();
854  /*Mean->SetLineColor(1);
855  Mean->SetLineWidth(1);
856  Mean_error->SetLineColor(1);
857  Mean_error->SetLineWidth(1);
858  Mean->Draw("same");
859  Mean_error->Draw("same");*/
860  c->Print(detector + "_meanE_ge_EtaCut_CircularCut.png");
861 
862  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->SetAxisRange(chi2_min,chi2_max,"Y");
863  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_chi2->Draw();
864  /*Chi2->SetLineColor(1);
865  Chi2->SetLineWidth(1);
866  Periphery_Chi2->SetLineColor(1);
867  Periphery_Chi2->SetLineWidth(1);
868  Chi2->Draw("same");
869  Periphery_Chi2->Draw("same");*/
870  c->Print(detector + "_chi2E_ge_EtaCut_CircularCut.png");
871 
872  gStyle -> SetOptStat(0);
873  gStyle -> SetOptFit(0);
874 
875  /* Sigma->SetLineColor(1);
876  Sigma->SetLineWidth(1);
877  Sigma_error->SetLineColor(1);
878  Sigma_error->SetLineWidth(1);
879  Point->SetLineColor(3);
880  Point->SetLineWidth(5);*/
881  fExp->SetLineColor(4); //38
882  fTrue->SetLineColor(2); //46
883 
884  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerStyle(kFullCircle);
885  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerColor(46); //30
886  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetMarkerSize(0.75);
887  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->SetAxisRange(sigma_min,sigma_max,"Y");
888 
889  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->Fit("fTrue", "M+");
890  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2->Draw("same");
891  fExp->Draw("same");
892  /*Sigma->Draw("same");
893  Sigma_error->Draw("same");
894  Point->Draw("same");*/
895 
896  TLegend* legend = new TLegend(1.75,1.75);
897  legend->SetHeader("Legend", "C");
898  //legend->AddEntry(te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated_temp_2, "#sigma_{e_{agg}} vs Generated Energy", "flep");
899  //legend->AddEntry((TObject*)0,"","");
900  legend->AddEntry(fTrue, "p_{0} + p_{1}/#sqrt{ge} (Fitted)", "l");
901  legend->AddEntry((TObject*)0,"","");
902  legend->AddEntry(fExp, "0.1 + 0.5/#sqrt{ge} (Requirement)", "l");
903  legend->SetTextSize(0.033);
904  legend->Draw();
905 
906  std::cout<<"reduced_chi2 of fit: "<<fTrue->GetChisquare()/fTrue->GetNDF()<<"\n";
907 
908  c->Print(detector + "_sigmaE_ge_EtaCut_CircularCut.png");
909  usleep(5e6);
910 
911  gStyle -> SetOptStat(11);
912  gStyle -> SetOptFit(112);
913 
914  for(int sno = 0; sno < nSlicesx; sno++){
915  TString plusOne = (TString)(sno + 1);
916  TString nameF = detector + "_sigmaE_slice" + arr[sno] + "_EtaCut_CircularCut.png";
917  slices[sno] -> Fit("fit", "M+");
918  slices[sno] -> Draw("hist same");
919  c->Print(nameF);
920  }
921 
922  te_aggregate_EtaCut_CircularCut_FEMC->Draw();
923  c->Print("te_aggregate_EtaCut_CircularCut_FEMC.png");
924  usleep(5e6);
925 
926  gStyle -> SetOptStat(1);
927 
928  te_by_ge_ge_EtaCut->Draw("colz");
929  c->Print(detector + "_te_by_ge_ge_EtaCut.png");
930  te_by_ge_ge_EtaCut_CircularCut->Draw("colz");
931  c->Print(detector + "_te_by_ge_ge_EtaCut_CircularCut.png");
932 
933  te_minus_ge_by_ge_ge_EtaCut->Draw("colz");
934  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut.png");
935  te_minus_ge_by_ge_ge_EtaCut_CircularCut->Draw("colz");
936  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut_CircularCut.png");
937  te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated->Draw("colz");
938  c->Print(detector + "_te_minus_ge_by_ge_ge_EtaCut_CircularCut_Recalibrated.png");
939 
940  gStyle -> SetOptStat(0);
941 
942  mean_te_by_ge_ge_EtaCut_CircularCut->Draw();
943  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut.png");
944 
945  mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL->Draw();
946  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut_FHCAL.png");
947 
948  mean_te_by_ge_ge_EtaCut_CircularCut_FEMC->Draw();
949  c->Print(detector + "_mean_te_by_ge_ge_EtaCut_CircularCut_FEMC.png");
950 
951  gStyle -> SetOptStat(1);
952  c->Close();
953 
954  }
955 
956  std::cout << "The total te is: " << total_te << endl;
957  std::cout << "The total te_CircularCut is: " << total_te_CircularCut << endl;
958  std::cout << "The total ge is: " << total_ge << endl;
959  std::cout<<"\n\nDone\n----------------------------------------------------------------------\n\n";
960 
961 }