nue_fnex_vs_caf_events.C
Go to the documentation of this file.
1 //Author: Prabhjot Singh (prabhjot@fnal.gov).
2 //Dated: Feb 03 2017.
3 //This program is a template macro find events for FNEX vs. CAFAna predictions or fits for the nue analysis for the following cases for all three cvn pid bins
4 //1. ND data
5 //2. ND MC intrinsic beam νe
6 //3. ND MC νμ CC
7 //4. ND MC neutral currents, NC
8 //
9 //5. FD data
10 //6. FD MC nue appearance, νμ → νe (SIGNAL)
11 //7. FD MC nue survival , νe → νe (BACKGROUND)
12 //8. FD MC numu survival , νμ → νμ (BACKGROUND)
13 //9. FD MC NC → NC (BACKGROUND)
14 //10. FD MC tau appearance (BACKGROUND)
15 //
16 //
17 //Set fit_or_plot = "plot" for Prediction Comparisons and
18 // = "fit" for Fir comparison
19 
20 #include <Header.h>
21 #include <HistogramAttr.h>
22 
24  gStyle->SetOptStat(" ");
25 
26  //fnex file
27  TFile* file_fnex = new TFile("/nova/ana/users/prabhjot/FNEX/nue_analysis/realdata/plotpoint_job_caf_bestfitpoint_with_cosmic_nh/fnex_plotpoint_hist.root", "read");
28 
29  //caf file
30  TFile* file_caf = new TFile("/nova/ana/users/prabhjot/SRT/development/results/histograms_forFNEX_CVN_prop_stats.root", "read");
31 
32  double nd_pot = 0.;
33  double fd_pot = 0.;
34  double dcp = 0.;
35  double th13 = 0.;
36  double th23 = 0.;
37  double dmsq32 = 0.;
38  double ymax = 0.;
39 
40  float caf_data_total = 0.;
41  float caf_mc_total = 0.;
42  float caf_signal_total = 0.;
43  float caf_background_total = 0.;
44  float fnex_data_total = 0.;
45  float fnex_mc_total = 0.;
46  float fnex_signal_total = 0.;
47  float fnex_background_total = 0.;
48 
49  //[3] = "low", "mid", "high"
50  float caf_data_pid[3] = {0.};
51  float caf_mc_pid[3] = {0.};
52  float caf_signal_pid[3] = {0.};
53  float caf_background_pid[3] = {0.};
54  float fnex_data_pid[3] = {0.};
55  float fnex_mc_pid[3] = {0.};
56  float fnex_signal_pid[3] = {0.};
57  float fnex_background_pid[3]= {0.};
58 
59 
60  //POT and best fit information.
61  std::cout << " " << std::endl;
62  std::cout << "--------------------- POT and BF parameters (START) ----------------" << std::endl;
63  TVectorT<double> pot = (TVectorT<double>)file_caf->Get("pot_nd_fd");
64  // //pot->Print();
65  nd_pot = pot[0];
66  fd_pot = pot[1];
67  cout << "ND POT = " << nd_pot << endl;
68  cout << "FD POT = " << fd_pot << endl;
69 
70  TVectorT<double> par = (TVectorT<double>)file_caf->Get("best_dcp_th13_th23_dmsq32");
71  dcp = par[0];
72  th13 = par[1];
73  th23 = par[2];
74  dmsq32 = par[3];
75  cout << "dCP (pi rad) = " << dcp/TMath::Pi() << endl;
76  cout << "th31 = " << th13 << endl;
77  cout << "th23 = " << th23 << " sinsqth23 = " << sin(th23)*sin(th23) << endl;
78  cout << "dmsq32 = " << dmsq32 << endl;
79  std::cout << "--------------------- POT and BF parameters (END) ---------------" << std::endl;
80  std::cout << " " << std::endl;
81 
82  ofstream summary;
83 // summary .open("Images/fnexvscaf/summary.txt");
84  summary .open("Images/fnexvscaf/summary_"+TString::Format("dCP_%2.2f",dcp)+TString::Format("_sinsqth23_%2.2f",sin(th23)*sin(th23))+".txt");
85 
86  ofstream summary_datamc_pids;
87 // summary_datamc_pids .open("Images/fnexvscaf/summary_datamc_pids.txt");
88  summary_datamc_pids .open("Images/fnexvscaf/summary_"+TString::Format("dCP_%2.2f",dcp)+TString::Format("_sinsqth23_%2.2f",sin(th23)*sin(th23))+"_datamc_pids.txt");
89 
90  ofstream summary_datamc_total;
91 // summary_datamc_total .open("Images/fnexvscaf/summary_datamc_total.txt");
92  summary_datamc_total .open("Images/fnexvscaf/summary_"+TString::Format("dCP_%2.2f",dcp)+TString::Format("_sinsqth23_%2.2f",sin(th23)*sin(th23))+"_datamc_total.txt");
93 
94  int first_intr = 0;
95  int last_intr = 0;
96 
97  TString det_caf;
98  TString fit_or_plot = "plot"; // "plot" or "fit/Fit"
99  TString corr_uncorr = "Corrected"; // "Corrected" or "Uncorrected"
100 
101  TString datamc[] = { "mc", "data"};
102  TString PID[] = {"Low", "Mid", "High"};
103 
104  TString fnex_intr[] = {"ee", "mm" , "nc", "me" , "ee" , "mm" , "mt" , "nc", "em" , "mebar" , "eebar" , "mmbar" , "cosmic"};
105  TString caf_intr[] = {"nue", "numu", "nc", "nue_app", "nue_surv", "numu_surv", "tau_cc_all", "nc", "numu_app", "anue_app", "anue_surv", "anumu_surv", "cosmic"};
106 
107  TString det[] = {"ND", "FD"};
108 
109  TH1F* hfnex;
110  TH1D* hcaf;
111  TH1D* h2;
112 
113  //POT and best fit information.
114  summary << " " << std::endl;
115  summary << "--------------------- POT and BF parameters (START) ----------------" << std::endl;
116  summary << "ND POT = " << nd_pot << std::endl;
117  summary << "FD POT = " << fd_pot << std::endl;
118  summary << "dCP (pi rad) = " << dcp/TMath::Pi() << std::endl;
119  summary << "th31 = " << th13 << std::endl;
120  summary << "th23 = " << th23 << " sinsqth23 = " << sin(th23)*sin(th23) << endl;
121  summary << "dmsq32 = " << dmsq32 << std::endl;
122  summary << "--------------------- POT and BF parameters (END) ---------------" << std::endl;
123 
124  //POT and best fit information.
125  summary_datamc_total << " " << std::endl;
126  summary_datamc_total << "--------------------- POT and BF parameters (START) ----------------" << std::endl;
127  summary_datamc_total << "ND POT = " << nd_pot << std::endl;
128  summary_datamc_total << "FD POT = " << fd_pot << std::endl;
129  summary_datamc_total << "dCP (pi rad) = " << dcp/TMath::Pi() << std::endl;
130  summary_datamc_total << "th31 = " << th13 << std::endl;
131  summary_datamc_total << "th23 = " << th23 << " sinsqth23 = " << sin(th23)*sin(th23) << endl;
132  summary_datamc_total << "dmsq32 = " << dmsq32 << std::endl;
133  summary_datamc_total << "--------------------- POT and BF parameters (END) ---------------" << std::endl;
134 
135  //POT and best fit information.
136  summary_datamc_pids << " " << std::endl;
137  summary_datamc_pids << "--------------------- POT and BF parameters (START) ----------------" << std::endl;
138  summary_datamc_pids << "ND POT = " << nd_pot << std::endl;
139  summary_datamc_pids << "FD POT = " << fd_pot << std::endl;
140  summary_datamc_pids << "dCP (pi rad) = " << dcp/TMath::Pi() << std::endl;
141  summary_datamc_pids << "th31 = " << th13 << std::endl;
142  summary_datamc_pids << "th23 = " << th23 << " sinsqth23 = " << sin(th23)*sin(th23) << endl;
143  summary_datamc_pids << "dmsq32 = " << dmsq32 << std::endl;
144  summary_datamc_pids << "--------------------- POT and BF parameters (END) ---------------" << std::endl;
145 
146  summary_datamc_total << " " << std::endl;
147  summary_datamc_total << "//---------------------------------------------------------------//" << std::endl;
148  summary_datamc_total
149  << std::setw(20)
150  << "CAF"
151  << std::setw(10)
152  << "FNEX"
153  << std::setw(20)
154  << "Change %"
155  << std::endl;
156 
157  summary_datamc_pids << " " << std::endl;
158  summary_datamc_pids << "//---------------------------------------------------------------//" << std::endl;
159  summary_datamc_pids
160  << std::setw(20)
161  << "CAF"
162  << std::setw(10)
163  << "FNEX"
164  << std::setw(20)
165  << "Change %"
166  << std::endl;
167 
168  //loop over ND and FD detectors.
169  //idet = 0 for ND
170  //idet = 1 for FD
171  for(int idet = 0; idet < sizeof(det)/sizeof(det[0]); idet++){
172  if(det[idet]=="ND"){
173  det_caf = "nd";
174  first_intr = 0;
175  last_intr = 3;
176  }
177  else if(det[idet]=="FD"){
178  det_caf = "fd";
179  first_intr = 3;
180  last_intr = sizeof(caf_intr)/sizeof(caf_intr[0]);
181  }
182 
183  //loop over data or MC
184  for(int idatamc = 0; idatamc < sizeof(datamc)/sizeof(datamc[0]); idatamc++){
185 
186  //loop over all interactions
187  //intr = 0 for ND fnex ee and caf nue
188  //inrt = 1 for ND fnex mm and caf numu
189  //intr = 2 for ND fnex nc and caf nc
190  //
191  //intr = 3 for FD fnex me and caf nue_app
192  //intr = 4 for FD fnex ee and caf nue_surv
193  //intr = 5 for FD fnex mm and caf numu_surv
194  //intr = 6 for FD fnex mt and caf tau_cc_all
195  //intr = 7 for FD fnex nc and caf nc
196  for(int intr = first_intr; intr < last_intr; intr++){
197 
198  float caf_events_interactionwise = 0.;
199  float fnex_events_interactionwise = 0.;
200  float change = 0.;
201  //Write integrals of different histograms in a summary file
202  summary << " " << std::endl;
203  summary << "//---------------------------------------------------------------//" << std::endl;
204  if(datamc[idatamc]=="mc")
205  summary
206  << "Detector : "
207  << det[idet]
208  << " Data/MC : "
209  << datamc[idatamc]
210  << " Interaction : "
211  << fnex_intr[intr]
212  << std::endl;
213  if(datamc[idatamc]=="data")
214  summary
215  << "Detector : "
216  << det[idet]
217  << " Data/MC : "
218  << datamc[idatamc]
219  << std::endl;
220  summary
221  << std::setw(20)
222  << "CAF"
223  << std::setw(10)
224  << "FNEX"
225  << std::setw(20)
226  << "Change %"
227  << std::endl;
228 
229  for(int ipid = 0; ipid < sizeof(PID)/sizeof(PID[0]); ipid++){
230 
231  TString sign = " ";
232 
233  //fnex Corrected Histograms
234  if(datamc[idatamc]=="data")
235  hfnex = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"Data"+corr_uncorr);
236 
237  else if(datamc[idatamc]=="mc" && fnex_intr[intr] != "cosmic")
238  hfnex = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+datamc[idatamc]+"_"+fnex_intr[intr]+corr_uncorr);
239 
240  else if(datamc[idatamc]=="mc" && fnex_intr[intr] == "cosmic")
241  hfnex = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+fnex_intr[intr]+corr_uncorr);
242 
243  if(!hfnex){
244  if(datamc[idatamc]=="data")
245  std::cout << "FNEX histogram: " << "fit/Fit/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"Data"+corr_uncorr << " does not exists." << std::endl;
246 
247  else if(datamc[idatamc]=="mc" && fnex_intr[intr] != "cosmic")
248  std::cout << "FNEX histogram: " << "fit/Fit/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+datamc[idatamc]+"_"+fnex_intr[intr]+corr_uncorr << " does not exists." << std::endl;
249 
250  else if(datamc[idatamc]=="mc" && fnex_intr[intr] == "cosmic")
251  std::cout << "FNEX histogram: " << "fit/Fit/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+fnex_intr[intr]+corr_uncorr << " does not exists." << std::endl;
252  std::cout << "---------------------->>>>>>>>>>>> ABORTING >>>>>>>>>>>>>>>>--------------------------" << std::endl;
253  abort();
254  }//condition to check if fnex histogram exists or not
255 
256  if(det[idet] == "FD" && datamc[idatamc]=="mc" && fnex_intr[intr] == "mt"){
257  hfnex->Add((TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+datamc[idatamc]+"_mtbar"+corr_uncorr));
258  }
259 
260  hfnex->Rebin(2); //rebin fnex because caf bin is 0.5GeV/bin and fnex bin is 0.25GeV/bin
261 
262  //caf
263  if(datamc[idatamc]=="data")
264  h2 = (TH1D*)file_caf->Get(det_caf+"_"+datamc[idatamc]);
265  else if(datamc[idatamc]=="mc" && det[idet]=="ND")
266  h2 = (TH1D*)file_caf->Get(det_caf+"_"+datamc[idatamc]+"/"+caf_intr[intr]);
267  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] != "cosmic")
268  h2 = (TH1D*)file_caf->Get(det_caf+"_extrap_oscil_BF/"+caf_intr[intr]);
269  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] == "cosmic")
270  h2 = (TH1D*)file_caf->Get(det_caf+"_"+caf_intr[intr]);
271 
272  if(!h2){
273  if(datamc[idatamc]=="data")
274  std::cout << "CAF histogram: " << det_caf+"_"+datamc[idatamc] << " does not exists." << std::endl;
275  else if(datamc[idatamc]=="mc" && det[idet]=="ND")
276  std::cout << "CAF histogram: " << det_caf+"_"+datamc[idatamc]+"/"+caf_intr[intr] << " does not exists." << std::endl;
277  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] != "cosmic")
278  std::cout << "CAF histogram: " << det_caf+"_extrap_oscil_BF/"+caf_intr[intr] << " does not exists." << std::endl;
279  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] == "cosmic")
280  std::cout << "CAF histogram: " << det_caf+"_"+caf_intr[intr] << " does not exists." << std::endl;
281  std::cout << "---------------------->>>>>>>>>>>> ABORTING >>>>>>>>>>>>>>>>--------------------------" << std::endl;
282  abort();
283  }//condition to check if CAF histogram exists or not
284 
285  hcaf = CorrectRange(h2, PID[ipid]);
286 
287  //set y maximum range
288  if(hfnex->GetBinContent(hfnex->GetMaximumBin()) > hcaf->GetBinContent(hcaf->GetMaximumBin()) )
289  ymax = 1.1*hfnex->GetBinContent(hfnex->GetMaximumBin());
290  else if(hfnex->GetBinContent(hfnex->GetMaximumBin()) <= hcaf->GetBinContent(hcaf->GetMaximumBin()) )
291  ymax = 1.1*hcaf->GetBinContent(hcaf->GetMaximumBin());
292 
293  change = std::abs(hcaf->Integral() - hfnex->Integral() )/hcaf->Integral() *100.;
294 
295  if(hcaf->Integral() > hfnex->Integral() )
296  sign = "-";
297  else if(hcaf->Integral() < hfnex->Integral() )
298  sign = "+";
299 
300  if(PID[ipid]=="Low")
301  summary
302  << std::setw(10)
303  << "LOW : "
304  << std::setw(10)
305  << std::fixed << std::setprecision(2) << hcaf->Integral()
306  << std::setw(10)
307  << std::fixed << std::setprecision(2) << hfnex->Integral()
308  << std::setw(10)
309  << sign
310  << " "
311  << std::fixed << std::setprecision(1) << change
312  << std::endl;
313 
314  else if(PID[ipid]=="Mid")
315  summary
316  << std::setw(10)
317  << "MID : "
318  << std::setw(10)
319  << std::fixed << std::setprecision(2) << hcaf->Integral()
320  << std::setw(10)
321  << std::fixed << std::setprecision(2) << hfnex->Integral()
322  << std::setw(10)
323  << sign
324  << " "
325  << std::fixed << std::setprecision(1) << change
326  << std::endl;
327 
328  else if(PID[ipid]=="High")
329  summary
330  << std::setw(10)
331  << "High : "
332  << std::setw(10)
333  << std::fixed << std::setprecision(2) << hcaf->Integral()
334  << std::setw(10)
335  << std::fixed << std::setprecision(2) << hfnex->Integral()
336  << std::setw(10)
337  << sign
338  << " "
339  << std::fixed << std::setprecision(1) << change
340  << std::endl;
341 
342  else
343  summary
344  << "Unknown PID"
345  << std::endl;
346 
347  caf_events_interactionwise += hcaf ->Integral();
348  fnex_events_interactionwise += hfnex ->Integral();
349 
350  //total data
351  if(det[idet] == "FD" && datamc[idatamc]=="data"){
352  caf_data_total += hcaf ->Integral();
353  fnex_data_total += hfnex ->Integral();
354 
355  if(PID[ipid]=="Low"){
356  caf_data_pid[0] += hcaf ->Integral();
357  fnex_data_pid[0] += hfnex ->Integral();
358  }
359  else if (PID[ipid]=="Mid"){
360  caf_data_pid[1] += hcaf ->Integral();
361  fnex_data_pid[1] += hfnex ->Integral();
362  }
363  else if (PID[ipid]=="High"){
364  caf_data_pid[2] += hcaf ->Integral();
365  fnex_data_pid[2] += hfnex ->Integral();
366  }
367  }
368  //total MC
369  if(det[idet] == "FD" && datamc[idatamc]=="mc" && fnex_intr[intr] != "mebar"){
370  caf_mc_total += hcaf ->Integral();
371  fnex_mc_total += hfnex ->Integral();
372  if(PID[ipid]=="Low"){
373  caf_mc_pid[0] += hcaf ->Integral();
374  fnex_mc_pid[0] += hfnex ->Integral();
375  }
376  else if (PID[ipid]=="Mid"){
377  caf_mc_pid[1] += hcaf ->Integral();
378  fnex_mc_pid[1] += hfnex ->Integral();
379  }
380  else if (PID[ipid]=="High"){
381  caf_mc_pid[2] += hcaf ->Integral();
382  fnex_mc_pid[2] += hfnex ->Integral();
383  }
384 }
385  //Signal
386  //with mebar
387 // if(det[idet] == "FD" && datamc[idatamc]=="mc" && (fnex_intr[intr] == "me" || fnex_intr[intr] == "mebar")){
388  //without mebar
389  if(det[idet] == "FD" && datamc[idatamc]=="mc" && fnex_intr[intr] == "me"){
390  caf_signal_total += hcaf ->Integral();
391  fnex_signal_total += hfnex ->Integral();
392  if(PID[ipid]=="Low"){
393  caf_signal_pid[0] += hcaf ->Integral();
394  fnex_signal_pid[0] += hfnex ->Integral();
395  }
396  else if (PID[ipid]=="Mid"){
397  caf_signal_pid[1] += hcaf ->Integral();
398  fnex_signal_pid[1] += hfnex ->Integral();
399  }
400  else if (PID[ipid]=="High"){
401  caf_signal_pid[2] += hcaf ->Integral();
402  fnex_signal_pid[2] += hfnex ->Integral();
403  }
404 }
405  //backgrounds
406  if(det[idet] == "FD" && datamc[idatamc]=="mc" && fnex_intr[intr] != "me" && fnex_intr[intr] != "mebar"){
407  caf_background_total += hcaf ->Integral();
408  fnex_background_total += hfnex ->Integral();
409  if(PID[ipid]=="Low"){
410  caf_background_pid[0] += hcaf ->Integral();
411  fnex_background_pid[0] += hfnex ->Integral();
412  }
413  else if (PID[ipid]=="Mid"){
414  caf_background_pid[1] += hcaf ->Integral();
415  fnex_background_pid[1] += hfnex ->Integral();
416  }
417  else if (PID[ipid]=="High"){
418  caf_background_pid[2] += hcaf ->Integral();
419  fnex_background_pid[2] += hfnex ->Integral();
420  }
421 }
422 
423  delete hfnex;
424  delete hcaf;
425  delete h2;
426 
427  }//end of loop over all PIDs
428 
429  if(caf_events_interactionwise > fnex_events_interactionwise)
430  sign = "-";
431  if(caf_events_interactionwise < fnex_events_interactionwise)
432  sign = "+";
433 
434  summary
435  << std::setw(10)
436  << "ALL :"
437  << std::setw(10)
438  << std::fixed << std::setprecision(2) << caf_events_interactionwise
439  << std::setw(10)
440  << std::fixed << std::setprecision(2) << fnex_events_interactionwise
441  << std::setw(10)
442  << sign
443  << " "
444  << std::fixed << std::setprecision(1) << std::abs(caf_events_interactionwise - fnex_events_interactionwise)/caf_events_interactionwise * 100.
445  << std::endl;
446 
447  //loop only once if its a data case
448  if(datamc[idatamc]=="data") break;
449 
450  }//end of loop over all interactions.
451  }//end of loop over data or mc
452  }//end of loop over detector type
453 
454  if(caf_data_total > fnex_data_total) sign = "-";
455  else if(caf_data_total < fnex_data_total) sign = "+";
456  summary_datamc_total
457  << "Data : "
458  << std::setw(10) << std::fixed << std::setprecision(2) << caf_data_total
459  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_data_total
460  << std::setw(10)
461  << sign
462  << " "
463  << std::fixed << std::setprecision(1) << std::abs(caf_data_total - fnex_data_total)/caf_data_total*100.
464  << std::endl;
465 
466  if(caf_mc_total > fnex_mc_total) sign = "-";
467  else if(caf_mc_total < fnex_mc_total) sign = "+";
468  summary_datamc_total
469  << "Prediction : "
470  << std::setw(10) << std::fixed << std::setprecision(2) << caf_mc_total
471  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_mc_total
472  << std::setw(10)
473  << sign
474  << " "
475  << std::fixed << std::setprecision(1) << std::abs(caf_mc_total - fnex_mc_total)/caf_mc_total*100.
476  << std::endl;
477 
478  if(caf_signal_total > fnex_signal_total) sign = "-";
479  else if(caf_signal_total < fnex_signal_total) sign = "+";
480  summary_datamc_total
481  << "Signal : "
482  << std::setw(10) << std::fixed << std::setprecision(2) << caf_signal_total
483  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_signal_total
484  << std::setw(10)
485  << sign
486  << " "
487  << std::fixed << std::setprecision(1) << std::abs(caf_signal_total - fnex_signal_total)/caf_signal_total*100.
488  << std::endl;
489 
490  if(caf_background_total > fnex_background_total) sign = "-";
491  else if(caf_background_total < fnex_background_total) sign = "+";
492  summary_datamc_total
493  << "Background : "
494  << std::setw(10) << std::fixed << std::setprecision(2) << caf_background_total
495  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_background_total
496  << std::setw(10)
497  << sign
498  << " "
499  << std::fixed << std::setprecision(1) << std::abs(caf_background_total - fnex_background_total)/caf_background_total*100.
500  << std::endl;
501 
502  //cout prints for all pids
503  for(int i_pid = 0; i_pid < 3; ++i_pid){
504  summary_datamc_pids << " " << std::endl;
505  if(i_pid==0) summary_datamc_pids << "LOW" << std::endl;
506  if(i_pid==1) summary_datamc_pids << "MID" << std::endl;
507  if(i_pid==2) summary_datamc_pids << "HIGH" << std::endl;
508 
509  if(caf_data_pid[i_pid] > fnex_data_pid[i_pid]) sign = "-";
510  else if(caf_data_pid[i_pid] < fnex_data_pid[i_pid]) sign = "+";
511  summary_datamc_pids
512  << "Data : "
513  << std::setw(10) << std::fixed << std::setprecision(2) << caf_data_pid[i_pid]
514  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_data_pid[i_pid]
515  << std::setw(10)
516  << sign
517  << " "
518  << std::fixed << std::setprecision(1) << std::abs(caf_data_pid[i_pid] - fnex_data_pid[i_pid])/caf_data_pid[i_pid]*100.
519  << std::endl;
520 
521  if(caf_mc_pid[i_pid] > fnex_mc_pid[i_pid]) sign = "-";
522  else if(caf_mc_pid[i_pid] < fnex_mc_pid[i_pid]) sign = "+";
523  summary_datamc_pids
524  << "Prediction : "
525  << std::setw(10) << std::fixed << std::setprecision(2) << caf_mc_pid[i_pid]
526  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_mc_pid[i_pid]
527  << std::setw(10)
528  << sign
529  << " "
530  << std::fixed << std::setprecision(1) << std::abs(caf_mc_pid[i_pid] - fnex_mc_pid[i_pid])/caf_mc_pid[i_pid]*100.
531  << std::endl;
532 
533  if(caf_signal_pid[i_pid] > fnex_signal_pid[i_pid]) sign = "-";
534  else if(caf_signal_pid[i_pid] < fnex_signal_pid[i_pid]) sign = "+";
535  summary_datamc_pids
536  << "Signal : "
537  << std::setw(10) << std::fixed << std::setprecision(2) << caf_signal_pid[i_pid]
538  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_signal_pid[i_pid]
539  << std::setw(10)
540  << sign
541  << " "
542  << std::fixed << std::setprecision(1) << std::abs(caf_signal_pid[i_pid] - fnex_signal_pid[i_pid])/caf_signal_pid[i_pid]*100.
543  << std::endl;
544 
545  if(caf_background_pid[i_pid] > fnex_background_pid[i_pid]) sign = "-";
546  else if(caf_background_pid[i_pid] < fnex_background_pid[i_pid]) sign = "+";
547  summary_datamc_pids
548  << "Background : "
549  << std::setw(10) << std::fixed << std::setprecision(2) << caf_background_pid[i_pid]
550  << std::setw(10) << std::fixed << std::setprecision(2) << fnex_background_pid[i_pid]
551  << std::setw(10)
552  << sign
553  << " "
554  << std::fixed << std::setprecision(1) << std::abs(caf_background_pid[i_pid] - fnex_background_pid[i_pid])/caf_background_pid[i_pid]*100.
555  << std::endl;
556  }//end of loop over all pids
557 
558  summary .close();
559  summary_datamc_pids .close();
560  summary_datamc_total .close();
561 
562 }//end of main program
563 //--------------------------------------------------------------//
564 TH1D* CorrectRange(TH1D* h, TString PID){
565  int firstbin = 0;
566  if(PID=="Low") firstbin = 0;
567  if(PID=="Mid") firstbin = 10;
568  if(PID=="High") firstbin = 20;
569  TH1D* hlow = new TH1D("hlow", "", 10, 0, 5);
570  for(unsigned int i = 1; i < h->GetNbinsX()+1; i++){
571  hlow->SetBinContent(i, h->GetBinContent(i+firstbin));
572  }
573  return hlow;
574 }//end of CorrectRange
575 //---------------------------------------------------------------//
double th23
Definition: runWimpSim.h:98
Int_t par
Definition: SimpleIterate.C:24
float abs(float number)
Definition: d0nt_math.hpp:39
bool datamc
Definition: fnexvscaf.C:20
Double_t ymax
Definition: plot.C:25
void nue_fnex_vs_caf_events()
#define pot
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
ofstream summary
T sin(T number)
Definition: d0nt_math.hpp:132
double dmsq32
TH1D * CorrectRange(TH1D *h, TString PID)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
double th13
Definition: runWimpSim.h:98
def sign(x)
Definition: canMan.py:197