nue_fnex_vs_caf_noextrapolation.C
Go to the documentation of this file.
1 //Author: Prabhjot Singh (prabhjot@fnal.gov).
2 //Dated: Sep 20, 2016 15:03.
3 //This program is a template macro to plot and compare FNEX vs. CAFAna NOT extrapolated oscillated spectrum 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  ofstream summary;
27  summary .open("Images/fnexvscaf/summary.txt");
28 
29  TFile* outfile = new TFile("Images/fnexvscaf/nue_fnex_vs_caf_histograms_noextrapolation.root", "RECREATE");
30 
31  double nd_pot = 0.;
32  double fd_pot = 0.;
33  double dcp = 0.;
34  double th13 = 0.;
35  double th23 = 0.;
36  double dmsq32 = 0.;
37  double ymax = 0.;
38  double ratio = 0.;
39  double ratio_error = 0.;
40  double totalevents_fnex = 0.;
41  double totalevents_caf = 0.;
42 
43  int first_intr = 0;
44  int last_intr = 0;
45 
46  TString det_caf;
47  TString fit_or_plot = "plot"; // "plot" or "fit/Fit"
48  TString corr_uncorr = "Uncorrected"; // "Corrected" or "Uncorrected"
49 
50  TString datamc[] = { "mc", "data"};
51  TString PID[] = {"Low", "Mid", "High"};
52 
53  TString fnex_intr[] = {"ee", "mm" , "nc", "me" , "ee" , "mm" , "mt" , "nc", "em" , "mebar" , "eebar" , "mmbar" , "cosmic"};
54  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"};
55 
56  TString det[] = {"ND", "FD"};
57 
58  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"); //fnex file
59 
60  // TFile* file_caf = new TFile("/nova/ana/users/prabhjot/FNEX/nue_analysis/caf/histograms_forFNEX_CVN_prop_stats.root", "read"); //cafana file
61  TFile* file_caf = new TFile("/nova/ana/users/prabhjot/SRT/development/results/histograms_forFNEX_CVN_prop_stats.root", "read");
62 
63  TH1F* hfnex;
64  TH1D* hcaf;
65  TH1D* h2;
66  TH1F* hist_ratio;
67  TH1F* hist_fd_totaldata_fnex = new TH1F("hist_fd_totaldata_fnex", " ", 60, 0 , 30);
68  TH1F* hist_fd_totaldata_caf = new TH1F("hist_fd_totaldata_caf", " ", 10, 0 , 5 );
69  TH1F* hist_fd_totalmc_fnex = new TH1F("hist_fd_totalmc_fnex", " ", 60, 0 , 30);
70  TH1F* hist_fd_totalmc_caf = new TH1F("hist_fd_totalmc_caf", " ", 10, 0 , 5 );
71  TH1F* hist_nd_totaldata_fnex = new TH1F("hist_nd_totaldata_fnex", " ", 60, 0 , 30);
72  TH1F* hist_nd_totaldata_caf = new TH1F("hist_nd_totaldata_caf", " ", 10, 0 , 5 );
73  TH1F* hist_nd_totalmc_fnex = new TH1F("hist_nd_totalmc_fnex", " ", 60, 0 , 30);
74  TH1F* hist_nd_totalmc_caf = new TH1F("hist_nd_totalmc_caf", " ", 10, 0 , 5 );
75 
76  //POT and best fit information.
77  std::cout << " " << std::endl;
78  std::cout << "--------------------- POT and BF parameters (START) ----------------" << std::endl;
79  TVectorT<double> pot = (TVectorT<double>)file_caf->Get("pot_nd_fd");
80  // //pot->Print();
81  nd_pot = pot[0];
82  fd_pot = pot[1];
83  cout << "ND POT = " << nd_pot << endl;
84  cout << "FD POT = " << fd_pot << endl;
85 
86  TVectorT<double> par = (TVectorT<double>)file_caf->Get("best_dcp_th13_th23_dmsq32");
87  dcp = par[0];
88  th13 = par[1];
89  th23 = par[2];
90  dmsq32 = par[3];
91  cout << "dCP = " << dcp << endl;
92  cout << "th31 = " << th13 << endl;
93  cout << "th23 = " << th23 << endl;
94  cout << "dmsq32 = " << dmsq32 << endl;
95  std::cout << "--------------------- POT and BF parameters (END) ---------------" << std::endl;
96  std::cout << " " << std::endl;
97 
98  //loop over ND and FD detectors.
99  //idet = 0 for ND
100  //idet = 1 for FD
101  for(int idet = 0; idet < sizeof(det)/sizeof(det[0]); idet++){
102  if(det[idet]=="ND"){
103  det_caf = "nd";
104  first_intr = 0;
105  last_intr = 3;
106  }
107  else if(det[idet]=="FD"){
108  det_caf = "fd";
109  first_intr = 3;
110  last_intr = sizeof(caf_intr)/sizeof(caf_intr[0]);
111  }
112  //loop over all three PIDs
113  //ipid = 0 for Low PID
114  //ipid = 1 for Mid PID
115  //ipid = 2 for High PID
116  for(int ipid = 0; ipid < sizeof(PID)/sizeof(PID[0]); ipid++){
117  TH1F* h_fnex_pid = new TH1F("h_fnex_pid", " ", 60, 0., 30.);
118  TH1F* h_caf_pid = new TH1F("h_caf_pid", " ", 10, 0., 5.);
119 
120  TH1F* h_fnex_pid_background = new TH1F("h_fnex_pid_background", " ", 60, 0., 30.);
121  TH1F* h_caf_pid_background = new TH1F("h_caf_pid_background", " ", 10, 0., 5.); //loop over data or mc plots
122  for(int idatamc = 0; idatamc < sizeof(datamc)/sizeof(datamc[0]); idatamc++){
123 
124  //loop over all interactions
125  //intr = 0 for ND fnex ee and caf nue
126  //inrt = 1 for ND fnex mm and caf numu
127  //intr = 2 for ND fnex nc and caf nc
128  //
129  //intr = 3 for FD fnex me and caf nue_app
130  //intr = 4 for FD fnex ee and caf nue_surv
131  //intr = 5 for FD fnex mm and caf numu_surv
132  //intr = 6 for FD fnex mt and caf tau_cc_all
133  //intr = 7 for FD fnex nc and caf nc
134  for(int intr = first_intr; intr < last_intr; intr++){
135  //fnex Uncorrected Histograms
136  if(datamc[idatamc]=="data") hfnex = (TH1F*)file_fnex->Get(fit_or_plot+"/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/"+corr_uncorr+"Stacks/"+det[idet]+"Data"+corr_uncorr);
137 
138  else if(datamc[idatamc]=="mc" && fnex_intr[intr] != "cosmic")
139  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);
140 
141  else if(datamc[idatamc]=="mc" && fnex_intr[intr] == "cosmic")
142  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);
143 
144  if(!hfnex){
145  if(datamc[idatamc]=="data")
146  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;
147 
148  else if(datamc[idatamc]=="mc" && fnex_intr[intr] != "cosmic")
149  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;
150 
151  else if(datamc[idatamc]=="mc" && fnex_intr[intr] == "cosmic")
152  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;
153  std::cout << "---------------------->>>>>>>>>>>> ABORTING >>>>>>>>>>>>>>>>--------------------------" << std::endl;
154  abort();
155  }//condition to check if fnex histogram exists or not
156 
157  if(det[idet] == "FD" && datamc[idatamc]=="mc" && fnex_intr[intr] == "mt"){
158  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));
159  }
160 
161  hfnex->Rebin(2); //rebin fnex because caf bin is 0.5GeV/bin and fnex bin is 0.25GeV/bin
162 
163  //caf
164  if(datamc[idatamc]=="data")
165  h2 = (TH1D*)file_caf->Get(det_caf+"_"+datamc[idatamc]);
166  else if(datamc[idatamc]=="mc" && det[idet]=="ND")
167  h2 = (TH1D*)file_caf->Get(det_caf+"_"+datamc[idatamc]+"/"+caf_intr[intr]);
168  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] != "cosmic")
169  h2 = (TH1D*)file_caf->Get(det_caf+"_extrap_oscil_BF/"+caf_intr[intr]);
170  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] == "cosmic")
171  h2 = (TH1D*)file_caf->Get(det_caf+"_"+caf_intr[intr]);
172 
173  if(!h2){
174  if(datamc[idatamc]=="data")
175  std::cout << "CAF histogram: " << det_caf+"_"+datamc[idatamc] << " does not exists." << std::endl;
176  else if(datamc[idatamc]=="mc" && det[idet]=="ND")
177  std::cout << "CAF histogram: " << det_caf+"_"+datamc[idatamc]+"/"+caf_intr[intr] << " does not exists." << std::endl;
178  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] != "cosmic")
179  std::cout << "CAF histogram: " << det_caf+"_extrap_oscil_BF/"+caf_intr[intr] << " does not exists." << std::endl;
180  else if(datamc[idatamc]=="mc" && det[idet]=="FD" && fnex_intr[intr] == "cosmic")
181  std::cout << "CAF histogram: " << det_caf+"_"+caf_intr[intr] << " does not exists." << std::endl;
182  std::cout << "---------------------->>>>>>>>>>>> ABORTING >>>>>>>>>>>>>>>>--------------------------" << std::endl;
183  abort();
184  }//condition to check if CAF histogram exists or not
185 
186  hcaf = CorrectRange(h2, PID[ipid]);
187 
188  //scale fnex by ratio of caf/fnex
189  // if(datamc[idatamc]=="mc" || det[idet]=="ND")
190  // hfnex->Scale(hcaf->Integral()/hfnex->Integral());
191 
192  //Write integrals of different histograms in a summary file
193  summary << "//---------------------------------------------------------------//" << std::endl;
194  summary << "PID : " << PID[ipid] << " PID" << std::endl;
195  summary << "Detector : " << det[idet] << std::endl;
196  summary << "Data or MC : " << datamc[idatamc] << std::endl;
197  if(datamc[idatamc]=="mc")
198  summary << "Interaction : " << fnex_intr[intr] << std::endl;
199  summary << "CAF events : " << hcaf->Integral() << std::endl;
200  summary << "FNEX Corrected events : " << hfnex->Integral() << std::endl;
201  summary << "//---------------------------------------------------------------//" << std::endl;
202 
203  if(datamc[idatamc]=="mc")
204  CompareHistos(
205  hcaf,
206  hfnex,
207  "E^{reco}_{#nu} (GeV)",
208  "Events",
209  632,
210  600,
211  2,
212  1,
213  0.64,
214  0.53,
215  0.89,
216  0.88,
217  det[idet]+" "+datamc[idatamc]+ " "+caf_intr[intr],
218  "CAF",
219  "FNEX",
220  PID[ipid]+"_PID_"+det[idet]+"_"+datamc[idatamc]+"_"+caf_intr[intr]+"_noextrapolation"
221  );
222 
223  else if(datamc[idatamc]=="data")
224  CompareHistos( hcaf,
225  hfnex,
226  "E^{reco}_{#nu} (GeV)",
227  "Events",
228  632,
229  600,
230  2,
231  1,
232  0.64,
233  0.53,
234  0.89,
235  0.88,
236  det[idet]+" "+datamc[idatamc],
237  "CAF",
238  "FNEX",
239  PID[ipid]+"_PID_"+det[idet]+"_"+datamc[idatamc]+"_noextrapolation"
240  );
241 
242  if(datamc[idatamc]=="data" && det[idet]=="FD"){
243  hist_fd_totaldata_fnex->Add(hfnex);
244  hist_fd_totaldata_caf->Add(hcaf);
245  }
246  if(datamc[idatamc]=="mc" && det[idet]=="FD"){
247  hist_fd_totalmc_fnex->Add(hfnex);
248  hist_fd_totalmc_caf->Add(hcaf);
249 
250  h_fnex_pid->Add(hfnex);
251  h_caf_pid->Add(hcaf);
252 
253  if( caf_intr[intr] != "nue_app"){
254  h_fnex_pid_background->Add(hfnex);
255  h_caf_pid_background->Add(hcaf);
256  }
257  }
258  if(datamc[idatamc]=="data" && det[idet]=="ND"){
259  hist_nd_totaldata_fnex->Add(hfnex);
260  hist_nd_totaldata_caf->Add(hcaf);
261  }
262  if(datamc[idatamc]=="mc" && det[idet]=="ND"){
263  hist_nd_totalmc_fnex->Add(hfnex);
264  hist_nd_totalmc_caf->Add(hcaf);
265 
266  h_fnex_pid->Add(hfnex);
267  h_caf_pid->Add(hcaf);
268  }
269 
270  delete hfnex;
271  delete hcaf;
272  delete h2;
273 
274  //loop only once if its a data case
275  if(datamc[idatamc]=="data") break;
276  }//end of loop over all interactions.
277  if(datamc[idatamc]=="mc") CompareHistos(h_caf_pid, h_fnex_pid, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, det[idet]+" "+datamc[idatamc]+" "+PID[ipid], "CAF", "FNEX", det[idet]+"_"+datamc[idatamc]+"_"+PID[ipid]);
278 
279  if(datamc[idatamc]=="mc" && det[idet]=="FD")
280  CompareHistos(h_caf_pid_background, h_fnex_pid_background, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, det[idet]+" "+datamc[idatamc]+" "+PID[ipid]+ " Background", "CAF", "FNEX", det[idet]+"_"+datamc[idatamc]+"_"+PID[ipid]+"_background_noextrapolation");
281 
282  if(h_caf_pid) delete h_caf_pid;
283  if(h_fnex_pid) delete h_fnex_pid;
284  if(h_caf_pid_background) delete h_caf_pid_background;
285  if(h_fnex_pid_background) delete h_fnex_pid_background;
286  }//end of loop over data or mc
287  }//end of loop over all PIDs
288  }//end of loop over detector type
289  CompareHistos(hist_fd_totaldata_caf, hist_fd_totaldata_fnex, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "FD data total", "CAF", "FNEX", "FD_data_total_noextrapolation");
290  CompareHistos(hist_fd_totalmc_caf, hist_fd_totalmc_fnex, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "FD MC total", "CAF", "FNEX", "FD_mc_total_noextrapolation");
291 
292  CompareHistos(hist_nd_totaldata_caf, hist_nd_totaldata_fnex, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "ND data total", "CAF", "FNEX", "ND_data_total_noextrapolation");
293  CompareHistos(hist_nd_totalmc_caf, hist_nd_totalmc_fnex, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "ND mc total", "CAF", "FNEX", "ND_mc_total_noextrapolation");
294 
295  CompareHistos(hist_fd_totalmc_caf, hist_fd_totaldata_caf, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "FD CAF total", "MC", "data", "CAF_FD_databymc_total_noextrapolation");
296  CompareHistos(hist_fd_totalmc_fnex, hist_fd_totaldata_fnex, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "FD FNEX total", "MC", "data", "FNEX_FD_databymc_total_noextrapolation");
297 
298  CompareHistos(hist_nd_totalmc_caf, hist_nd_totaldata_caf, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "ND CAF total", "MC", "data", "CAF_ND_databymc_total_noextrapolation");
299  CompareHistos(hist_nd_totalmc_fnex, hist_nd_totaldata_fnex, "E^{reco}_{#nu} (GeV)", "Events", 632, 600, 2, 1, 0.64, 0.53, 0.89, 0.88, "ND FNEX total", "MC", "data", "FNEX_ND_databymc_total_noextrapolation");
300 
301  summary .close();
302 
303  outfile->Close();
304 }//end of main program
305 //--------------------------------------------------------------//
306 TH1D* CorrectRange(TH1D* h, TString PID){
307  int firstbin = 0;
308  if(PID=="Low") firstbin = 0;
309  if(PID=="Mid") firstbin = 10;
310  if(PID=="High") firstbin = 20;
311  TH1D* hlow = new TH1D("hlow", "", 10, 0, 5);
312  for(unsigned int i = 1; i < h->GetNbinsX()+1; i++){
313  hlow->SetBinContent(i, h->GetBinContent(i+firstbin));
314  }
315  return hlow;
316 }//end of CorrectRange
317 //---------------------------------------------------------------//
318 
double th23
Definition: runWimpSim.h:98
TH1 * ratio(TH1 *h1, TH1 *h2)
Int_t par
Definition: SimpleIterate.C:24
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)
Double_t ymax
Definition: plot.C:25
#define pot
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
ofstream summary
void nue_fnex_vs_caf_noextrapolation()
double dmsq32
double th13
Definition: runWimpSim.h:98
FILE * outfile
Definition: dump_event.C:13
TH1D * CorrectRange(TH1D *h, TString PID)