caf_vs_fnex_data_mc.C
Go to the documentation of this file.
1 //Author Prabhjot Singh
2 //Dated: June 23, 2017
3 //Comapres foloowing for both nue only and numu-nue joint analysis
4 //Need to pass in correct CAF and FNEX input files for nue only or numu-nue anaysis
5  //1. FNEX ND Data vs. CAF ND Data,
6  //2. FNEX ND MC vs. CAF ND MC
7  //3. FNEX ND Data vs. FNEX ND MC
8  //4. CAF ND Data vs. CAF ND MC
9  //5. FNEX ND Data/MC vs. CAF ND Data/MC and their double ratios
10  //
11  //
12  //6. FNEX FD Data vs. CAF FD Data,
13  //7. FNEX FD MC vs. CAF FD MC
14  //8. FNEX FD Data vs. FNEX FD MC
15  //9. CAF FD Data vs. CAF FD MC
16  //10. FNEX FD Data/MC vs. CAF FD Data/MC and their double ratios
17 
18 
19 #include <Header.h>
20 #include <HistogramAttr.h>
21 
22 
24  gStyle->SetOptStat(" ");
25 
26  //FNEX file
27  //1. nue only stats only
28  //1a. FNEX_prediction_at_CAF_BFP
29  //std::string f_fnex = "/nova/ana/users/prabhjot/FNEX/nue_analysis/FNEX_prediction_at_CAF_BFP/Stats/revision_25698/fnex_plotpoint_hist.root";
30  //1b. FNEX BFP
31  std::string f_fnex = "/nova/ana/users/prabhjot/FNEX/nue_analysis/BestFitPoint/StatOnly/RealData/UO/Tolerance_1/fnex_fitpoint_hist.root";
32 
33  //2. nue only + systs
34  //2a. FNEX_prediction_at_CAF_BFP
35  //std::string f_fnex = " ";
36  //2b. FNEX BFP
37  //std::string f_fnex = " ";
38 
39  //3. joint fit stats only
40  //3a. FNEX_prediction_at_CAF_BFP
41  //std::string f_fnex = "/nova/ana/users/prabhjot/FNEX/numu_nue_joint_analysis/FNEX_Prediction_at_CAF_BFP/StatOnly/nominal_cap/revision_25698/fnex_plotpoint_hist.root";
42  //3b. FNEX BFP
43  //std::string f_fnex = "/nova/ana/users/prabhjot/FNEX/numu_nue_joint_analysis/BestFitPoint/StatOnly/RealData/NH/LO/from_interactive/Tolerance_1/Brian_fix_revision_25698/fnex_fitpoint_hist.root";
44 
45  //4. joint fit + syst
46  //4a. FNEX_prediction_at_CAF_BFP
47  //std::string f_fnex = " ";
48  //4b. FNEX BFP
49  //std::string f_fnex = " ";
50 
51  TFile* file_fnex = new TFile(f_fnex.c_str());
52 
53  //CAF file
54  //nue only stats
55  std::string f_caf = "/nova/ana/users/prabhjot/FNEX/CAF_Files_by_Erika/root_files/histograms_forFNEX_CVN_nueOnly_stats.root";
56 
57  //nue only systs
58  //std::string f_caf = "";
59 
60  //joint fit stats
61  //std::string f_caf = "";
62 
63  //joint fit + systs
64  //std::string f_caf = "";
65 
66  TFile* file_caf = new TFile(f_caf.c_str());
67 
68  TString det_caf; //CAF Detector
69  TString savefile; //Output pdf or png file
70 
71  TString fit_or_plot;
72  TString analysis;
73  TString syst;
74 
75  TString corr_uncorr = "Corrected"; // "Corrected" or "Uncorrected"
76  TString selection = "nue"; // nue or numu Selection type
77 
78  ofstream summary; //Summary file with all the information about the number of events
79 
80  if(f_fnex.find("plotpoint") != std::string::npos) fit_or_plot = "plot" ; else fit_or_plot = "fit/Fit";
81  if(f_caf .find("nueOnly") != std::string::npos) analysis = "nue" ; else analysis = "joint";
82  if(f_caf .find("stats") != std::string::npos) syst = "stats" ; else syst = "systs";
83 
84  summary .open("/nova/app/users/prabhjot/mrb_build/srcs/novasoft/FNEX/macros/Images/fnexvscaf/summary_"+analysis+"_"+syst+"_fnex_vs_caf_predictions.txt"); //Name of the summary files
85 
86  TString datamc[] = {"MC", "Data"}; //MC or Data type FNEX spectrum
87  TString PID[] = {"Low", "Mid", "High"}; //LOW, MID or HIGH PID bins
88 
89  TString det[] = {"ND", "FD"};
90 
91  TH1F* hfnex; //FNEX histograms
92  TH1D* h2; //Three binned CAF histograms
93  TH1D* hcaf;
94 
95  TH1F* hfnex_data;
96  TH1F* hfnex_mc;
97  TH1F* hfnex_datamc_ratio;
98  TH1D* h2_data;
99  TH1D* h2_mc;
100  TH1D* h2_datamc_ratio;
101  TH1D* hcaf_data;
102  TH1D* hcaf_mc;
103  TH1D* hcaf_datamc_ratio;
104 
105  //loop over ND and FD detectors.
106  for(int idet = 0; idet < sizeof(det)/sizeof(det[0]); idet++){
107 
108  if(det[idet]=="ND") det_caf = "nd";
109  else if(det[idet]=="FD") det_caf = "fd";
110 
111  //Canvas for Data/MC ratio and double ratios
112  TCanvas* can_datamc_ratios = new TCanvas("can_"+det_caf+"_datamcratios", " ", 1200, 600);
113  if(analysis=="nue") can_datamc_ratios ->Divide(3, 1);
114  else can_datamc_ratios ->Divide(4, 1);
115 
116  //Canvas for CAF based Data and MC and their ratios
117  TCanvas* can_datamc_caf = new TCanvas("can_caf_"+det_caf+"_datamc", " ", 1200, 600);
118  if(analysis=="nue") can_datamc_caf ->Divide(3, 1);
119  else can_datamc_caf ->Divide(4, 1);
120 
121  //Canvas for FNEX based Data and MC and their ratios
122  TCanvas* can_datamc_fnex = new TCanvas("can_fnex_"+det_caf+"_datamc", " ", 1200, 600);
123  if(analysis=="nue") can_datamc_fnex ->Divide(3, 1);
124  else can_datamc_fnex ->Divide(4, 1);
125 
126  //Canvas for FNEX Data vs. CAF Data and their ratios
127  TCanvas* can_data = new TCanvas("can_fnex_"+det_caf+"_data", " ", 1200, 600);
128  if(analysis=="nue") can_data ->Divide(3, 1);
129  else can_data ->Divide(4, 1);
130 
131  //Canvas for FNEX MC vs. CAF MC and their ratios
132  TCanvas* can_mc = new TCanvas("can_fnex_"+det_caf+"_mc", " ", 1200, 600);
133  if(analysis=="nue") can_mc ->Divide(3, 1);
134  else can_mc ->Divide(4, 1);
135 
136  savefile=analysis+"_"+syst+"_fnex_vs_caf_predictions_"+det[idet];
137 
138  if(analysis=="joint"){
139  //numu selected plots
140  //FNEX histograms
141  hfnex_data = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"Data"+corr_uncorr);
142  if(det[idet]=="ND") hfnex_data ->Rebin(2);
143  else if(det[idet]=="FD") hfnex_data ->Rebin(1);
144  hfnex_mc = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"MC"+corr_uncorr);
145  if(det[idet]=="ND") hfnex_mc ->Rebin(2);
146  else if(det[idet]=="FD") hfnex_mc ->Rebin(1);
147 
148  hfnex_datamc_ratio = (TH1F*)hfnex_data->Clone("hfnex_datamc_ratio");
149  hfnex_datamc_ratio ->Sumw2();
150  hfnex_datamc_ratio ->Divide(hfnex_mc);
151 
152  //CAF histograms
153  //ND data
154  if(det[idet]=="ND") h2_data = (TH1D*)file_caf->Get("ND_Numu/"+det_caf+"_data");
155 
156  //FD data
157  else if(det[idet]=="FD") h2_data = (TH1D*)file_caf->Get("FD_Data/"+det_caf+"_data_numu");
158 
159  //ND MC
160  if(det[idet]=="ND"){
161  h2_mc = (TH1D*)file_caf ->Get("ND_Numu/nue_comp");
162  h2_mc ->Add((TH1D*)file_caf->Get("ND_Numu/antinue_comp"));
163  h2_mc ->Add((TH1D*)file_caf->Get("ND_Numu/numu_comp"));
164  h2_mc ->Add((TH1D*)file_caf->Get("ND_Numu/antinumu_comp"));
165  h2_mc ->Add((TH1D*)file_caf->Get("ND_Numu/nc_comp"));
166  }//end of ND MC condition
167 
168  else if(det[idet]=="FD"){
169  h2_mc = (TH1D*)file_caf ->Get("FD_Numu_Extrap_Oscil_BF/nue_app");
170  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/nue_surv"));
171  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/numu_surv"));
172  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/numu_app"));
173  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/anue_app"));
174  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/anue_surv"));
175  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/anumu_surv"));
176  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/anumu_app"));
177  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/tau_cc_all"));
178  h2_mc ->Add((TH1D*)file_caf->Get("FD_Numu_Extrap_Oscil_BF/nc"));
179  h2_mc ->Add((TH1D*)file_caf->Get("FD_Data/"+det_caf+"_cosmics_numu"));
180  }//end of CAF FD MC histograms
181 
182  if( det[idet]=="FD" ){
183  h2_data ->Rebin(1);
184  h2_mc ->Rebin(1);
185  }
186  else if( det[idet]=="ND" ){
187  h2_data ->Rebin(5);
188  h2_mc ->Rebin(5);
189  }
190 
191  h2_datamc_ratio = (TH1D*)h2_data->Clone("h2_datamc_ratio");
192  h2_datamc_ratio ->Sumw2();
193  h2_datamc_ratio ->Divide(h2_mc);
194 
195  //Draw Data/MC of FNEX vs. Data/MC of CAF and their double ratios
196  can_datamc_ratios->cd(1);
197  std::vector<TString> legends;
198  legends.push_back(det[idet]);
199  legends.push_back("#nu_{#mu} Sel.");
201  h2_datamc_ratio,
202  hfnex_datamc_ratio,
203  "E^{reco}_{#nu} (GeV)",
204  "Data/MC",
205  632,
206  600,
207  2,
208  1,
209  0.55,
210  0.53,
211  0.80,
212  0.88,
213  legends,
214  "CAF",
215  "FNEX",
216  "Double ratio",
217  can_datamc_ratios,
218  "P",
219  "no_scaling"
220  );
221 
222  //CAF Data vs. MC
223  can_datamc_caf->cd(1);
224  std::vector<TString> legends;
225  legends.push_back("CAF");
226  legends.push_back(det[idet]);
227  legends.push_back("#nu_{#mu} Sel.");
229  h2_mc,
230  h2_data,
231  "E^{reco}_{#nu} (GeV)",
232  "Events",
233  632,
234  600,
235  2,
236  1,
237  0.55,
238  0.53,
239  0.80,
240  0.88,
241  legends,
242  "MC",
243  "Data",
244  "ratio",
245  can_datamc_caf,
246  "P",
247  "no_scaling"
248  );
249 
250  //FNEX Data vs. MC
251  can_datamc_fnex->cd(1);
252  std::vector<TString> legends;
253  legends.push_back("FNEX");
254  legends.push_back(det[idet]);
255  legends.push_back("#nu_{#mu} Sel.");
257  hfnex_mc,
258  hfnex_data,
259  "E^{reco}_{#nu} (GeV)",
260  "Events",
261  632,
262  600,
263  2,
264  1,
265  0.55,
266  0.53,
267  0.80,
268  0.88,
269  legends,
270  "MC",
271  "Data",
272  "ratio",
273  can_datamc_fnex,
274  "P",
275  "no_scaling"
276  );
277 
278  //FNEX Data vs. CAF Data
279  can_data ->cd(1);
280  std::vector<TString> legends;
281  legends.push_back(det[idet]+" Data");
282  legends.push_back("#nu_{#mu} Sel.");
284  hcaf_data,
285  hfnex_data,
286  "E^{reco}_{#nu} (GeV)",
287  "Events",
288  632,
289  600,
290  2,
291  1,
292  0.55,
293  0.53,
294  0.80,
295  0.88,
296  legends,
297  "CAF",
298  "FNEX",
299  "ratio",
300  can_data,
301  "HIST",
302  "yes_scaliing"
303  );
304 
305  //FNEX MC vs. CAF MC
306  can_mc ->cd(1);
307  std::vector<TString> legends;
308  legends.push_back(det[idet]+" MC");
309  legends.push_back("#nu_{#mu} Sel.");
311  hcaf_mc,
312  hfnex_mc,
313  "E^{reco}_{#nu} (GeV)",
314  "Events",
315  632,
316  600,
317  2,
318  1,
319  0.55,
320  0.53,
321  0.80,
322  0.88,
323  legends,
324  "CAF",
325  "FNEX",
326  "ratio",
327  can_mc,
328  "HIST",
329  "yes_scaliing"
330  );
331 
332  delete hfnex_data;
333  delete hfnex_mc;
334  delete hcaf_data;
335  delete hcaf_mc;
336  }//condition that it was joint fit
337 
338  //loop over three PIDs
339  for(int ipid = 0; ipid < sizeof(PID)/sizeof(PID[0]); ipid++){
340  //FNEX histograms
341  hfnex_data = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"Data"+corr_uncorr);
342  if(det[idet]=="ND" ) hfnex_data ->Rebin(2); //rebin fnex because caf bin is 0.5GeV/bin and fnex bin is 0.25GeV/bin
343  else if(det[idet]=="FD" ) hfnex_data ->Rebin(2);
344 
345  hfnex_mc = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"MC"+corr_uncorr);
346  if(det[idet]=="ND" ) hfnex_mc->Rebin(2);
347  else if(det[idet]=="FD" ) hfnex_mc->Rebin(2);
348 
349  hfnex_datamc_ratio = (TH1F*)hfnex_data->Clone("hfnex_datamc_ratio");
350  hfnex_datamc_ratio ->Sumw2();
351  hfnex_datamc_ratio ->Divide(hfnex_mc);
352 
353  //CAF histograms
354  //ND data
355  if(det[idet]=="ND" ) h2_data = (TH1D*)file_caf->Get("ND_Nue/"+det_caf+"_data");
356 
357  //FD data
358  else if(det[idet]=="FD" ) h2_data = (TH1D*)file_caf->Get("FD_Data/"+det_caf+"_data_"+selection);
359 
360  //ND MC
361  if(det[idet]=="ND"){
362  h2_mc = (TH1D*)file_caf ->Get("ND_Nue/nue_comp");
363  h2_mc ->Add((TH1D*)file_caf->Get("ND_Nue/antinue_comp"));
364  h2_mc ->Add((TH1D*)file_caf->Get("ND_Nue/numu_comp"));
365  h2_mc ->Add((TH1D*)file_caf->Get("ND_Nue/antinumu_comp"));
366  h2_mc ->Add((TH1D*)file_caf->Get("ND_Nue/nc_comp"));
367  }//end of ND MC condition
368 
369  //FD MC
370  else if(det[idet]=="FD"){
371  h2_mc = (TH1D*)file_caf ->Get("FD_Nue_Extrap_Oscil_BF/nue_app");
372  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/nue_surv"));
373  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/numu_surv"));
374  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/numu_app"));
375  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/anue_app"));
376  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/anue_surv"));
377  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/anumu_surv"));
378  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/anumu_app"));
379  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/tau_cc_all"));
380  h2_mc ->Add((TH1D*)file_caf->Get("FD_Nue_Extrap_Oscil_BF/nc"));
381  h2_mc ->Add((TH1D*)file_caf->Get("FD_Data/"+det_caf+"_cosmics_"+selection));
382  }//end of CAF FD MC histograms
383  hcaf_data = (TH1D*)CorrectRange(h2_data, PID[ipid]);
384  hcaf_mc = (TH1D*)CorrectRange(h2_mc, PID[ipid]);
385 
386  hcaf_datamc_ratio = (TH1D*)hcaf_data->Clone("hcaf_datamc_ratio");
387  hcaf_datamc_ratio ->Sumw2();
388  hcaf_datamc_ratio ->Divide(hcaf_mc);
389 
390  //Draw Data/MC of FNEX vs. Data/MC of CAF and their double ratios
391  if(analysis == "nue") can_datamc_ratios->cd(1+ipid);
392  else can_datamc_ratios->cd(2+ipid);
393  std::vector<TString> legends;
394  legends.push_back(det[idet]);
395  legends.push_back("#nu_{e} Sel.");
396  legends.push_back(PID[ipid]+ " CVN");
397 
399  hcaf_datamc_ratio,
400  hfnex_datamc_ratio,
401  "E^{reco}_{#nu} (GeV)",
402  "Data/MC",
403  632,
404  600,
405  2,
406  1,
407  0.55,
408  0.53,
409  0.80,
410  0.88,
411  legends,
412  "CAF",
413  "FNEX",
414  "Double ratio",
415  can_datamc_ratios,
416  "P",
417  "no_scaling"
418  );
419 
420  //CAF Data vs. MC
421  if(analysis == "nue") can_datamc_caf->cd(1+ipid);
422  else can_datamc_caf->cd(2+ipid);
423  std::vector<TString> legends;
424  legends.push_back("CAF");
425  legends.push_back(det[idet]);
426  legends.push_back("#nu_{e} Sel.");
427  legends.push_back(PID[ipid]+ " CVN");
429  hcaf_mc,
430  hcaf_data,
431  "E^{reco}_{#nu} (GeV)",
432  "Events",
433  632,
434  600,
435  2,
436  1,
437  0.55,
438  0.53,
439  0.80,
440  0.88,
441  legends,
442  "MC",
443  "Data",
444  "ratio",
445  can_datamc_caf,
446  "P",
447  "no_scaling"
448  );
449 
450  //FNEX Data vs. MC
451  if(analysis == "nue") can_datamc_fnex->cd(1+ipid);
452  else can_datamc_fnex->cd(2+ipid);
453  std::vector<TString> legends;
454  legends.push_back("FNEX");
455  legends.push_back(det[idet]);
456  legends.push_back("#nu_{e} Sel.");
457  legends.push_back(PID[ipid]+ " CVN");
459  hfnex_mc,
460  hfnex_data,
461  "E^{reco}_{#nu} (GeV)",
462  "Events",
463  632,
464  600,
465  2,
466  1,
467  0.55,
468  0.53,
469  0.80,
470  0.88,
471  legends,
472  "MC",
473  "Data",
474  "ratio",
475  can_datamc_fnex,
476  "P",
477  "no_scaling"
478  );
479 
480  //FNEX Data vs. CAF Data
481  if(analysis == "nue") can_data ->cd(1+ipid);
482  else can_data ->cd(2+ipid);
483  std::vector<TString> legends;
484  legends.push_back(det[idet]+" Data");
485  legends.push_back("#nu_{e} Sel.");
486  legends.push_back(PID[ipid]+ " CVN");
488  hcaf_data,
489  hfnex_data,
490  "E^{reco}_{#nu} (GeV)",
491  "Events",
492  632,
493  600,
494  2,
495  1,
496  0.55,
497  0.53,
498  0.80,
499  0.88,
500  legends,
501  "CAF",
502  "FNEX",
503  "ratio",
504  can_data,
505  "HIST",
506  "yes_scaling"
507  );
508 
509  //FNEX MC vs. CAF MC
510  if(analysis == "nue") can_mc ->cd(1+ipid);
511  else can_mc ->cd(2+ipid);
512  std::vector<TString> legends;
513  legends.push_back(det[idet]+" MC ");
514  legends.push_back("#nu_{e} Sel.");
515  legends.push_back(PID[ipid]+ " CVN");
517  hcaf_mc,
518  hfnex_mc,
519  "E^{reco}_{#nu} (GeV)",
520  "Events",
521  632,
522  600,
523  2,
524  1,
525  0.55,
526  0.53,
527  0.80,
528  0.88,
529  legends,
530  "CAF",
531  "FNEX",
532  "ratio",
533  can_mc,
534  "HIST",
535  "yes_scaling"
536  );
537 
538  delete hfnex_data;
539  //delete hfnex_mc;
540  delete hcaf_data;
541  //delete hcaf_mc;
542 
543  }//end of loop over three CVN pids
544 
545  can_datamc_ratios ->cd();
546  can_datamc_ratios ->SaveAs("/nova/app/users/prabhjot/mrb_build/srcs/novasoft/FNEX/macros/Images/fnexvscaf/"+savefile+"_datamcratios.png");
547 
548  can_datamc_caf ->cd();
549  can_datamc_caf ->SaveAs("/nova/app/users/prabhjot/mrb_build/srcs/novasoft/FNEX/macros/Images/fnexvscaf/"+savefile+"_CAF_DataMC.png");
550 
551  can_datamc_fnex ->cd();
552  can_datamc_fnex ->SaveAs("/nova/app/users/prabhjot/mrb_build/srcs/novasoft/FNEX/macros/Images/fnexvscaf/"+savefile+"_FNEX_DataMC.png");
553 
554 
555  can_data ->cd();
556  can_data ->SaveAs("/nova/app/users/prabhjot/mrb_build/srcs/novasoft/FNEX/macros/Images/fnexvscaf/"+savefile+"_Data.png");
557 
558  can_mc ->cd();
559  can_mc ->SaveAs("/nova/app/users/prabhjot/mrb_build/srcs/novasoft/FNEX/macros/Images/fnexvscaf/"+savefile+"_MC.png");
560 
561  }//end of loop over detectors
562 
563 }//end of main function
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
bool datamc
Definition: fnexvscaf.C:20
void CompareHistos(TH1 *hbase_orig, TH1 *hsec_orig, TString xlabel, TString ylabel, int cbase, int csec, int lbase, int lsec, float x1, float y1, float x2, float y2, TString legtitle, TString legbase, TString legsec, TString savefile)
TH1 * CorrectRange(TH1 *h, TString PID)
void caf_vs_fnex_data_mc()
TH1F * h2
Definition: plot.C:45
ofstream summary
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154