nue_fnex_vs_caf_one_to_one.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 plots for the nue analysis for the following cases for all three cvn pid bins
4 //1. ND data
5 //
6 //
7 
8 #include <Header.h>
9 #include <HistogramAttr.h>
10 
12  gStyle->SetOptStat(" ");
13 
14  //TVectorT<double> par;
15 
16  TCanvas *c1;
17 
18  TPad* pad1;
19  TPad* pad2;
20 
21  TLegend* legend;
22  TLegend* legend_ratio;
23 
24  TLine* one;;
25 
26  double nd_pot = 0.;
27  double fd_pot = 0.;
28  double dcp = 0.;
29  double th13 = 0.;
30  double th23 = 0.;
31  double dmsq32 = 0.;
32  double ymax = 0.;
33  double ratio = 0.;
34  double ratio_error = 0.;
35 
36  int first_intr = 0;
37  int last_intr = 0;
38  int intr = 0; // User: this value will decide which interaction to consider.
39  //intr = 0 for ND fnex ee and caf nue
40  //inrt = 1 for ND fnex mm and caf numu
41  //intr = 2 for ND fnex nc and caf nc
42  //
43  //intr = 3 for FD fnex me and caf nue_app
44  //intr = 4 for FD fnex ee and caf nue_surv
45  //intr = 5 for FD fnex mm and caf numu_surv
46  //intr = 6 for FD fnex mt and caf tau_cc_all
47  //intr = 7 for FD fnex nc and caf nc
48 
49  TString det_caf;
50  TString datamc;
51 
52  TString PID[] = {"Low", "Mid", "High"};
53  //int ipid = 2;
54  //User:
55  //ipid = 0 for Low PID
56  //ipid = 1 for Mid PID
57  //ipid = 2 for High PID
58 
59  TString fnex_intr[] = {"ee", "mm", "nc", "me", "ee", "mm", "mt", "nc"};
60  TString caf_intr[] = {"nue", "numu", "nc", "nue_app", "nue_surv", "numu_surv", "tau_cc_all", "nc"};
61 
62  //TString det = "ND"; //User set ND or FD
63  TString det[] = {"ND", "FD"};
64  //int idet = 1;
65  //User:
66  //idet = 0 for ND
67  //idet = 1 for FD
68  //loop over detector type
69 
70  TFile* f1;
71  TFile* f2;
72 
73  TH1F* hfnex;
74  TH1D* hcaf;
75  TH1D* h2;
76  TH1F* hist_ratio;
77 
78  for(int idet = 0; idet < sizeof(det)/sizeof(det[0]); idet++){
79  if(det[idet]=="ND")
80  det_caf = "nd";
81  else if(det[idet]=="FD")
82  det_caf = "fd";
83  datamc = "data"; //User set data or mc
84 
85  if(det[idet]=="ND"){
86  first_intr = 0;
87  last_intr = 3;
88  }
89  else if(det[idet]=="FD"){
90  first_intr = 3;
91  last_intr = sizeof(caf_intr)/sizeof(caf_intr[0]);
92  }
93  for(int ipid = 0; ipid < sizeof(PID)/sizeof(PID[0]); ipid++){
94  //loop over all three PIDs
95  //for(int intr = first_intr; intr < last_intr; intr++){ //loop over all interactions
96  //fnex
97  f1 = new TFile("/nova/ana/users/prabhjot/FNEX/nue_analysis/fnex_plotpoint_hist.root", "read");
98 
99  if(datamc=="data") hfnex = (TH1F*)f1->Get("plot/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/CorrectedStacks/"+det[idet]+"DataCorrected");
100 
101  else if(datamc=="mc") hfnex = (TH1F*)f1->Get("plot/MultiExperiment/NOvA_NuESel_"+PID[ipid]+"PID_NuMuSel/CorrectedStacks/"+det[idet]+datamc+"_"+fnex_intr[intr]+"Corrected");
102  hfnex->Rebin(2); //rebin fnex because caf bin is 0.5GeV/bin and fnex bin is 0.25GeV/bin
103  //caf
104  f2 = new TFile("/nova/ana/users/prabhjot/FNEX/nue_analysis/caf/histograms_forFNEX_CVN_prop_stats.root", "read");
105 
106  if(datamc=="data")
107  h2 = (TH1D*)f2->Get(det_caf+"_"+datamc);
108  else if(datamc=="mc" && det[idet]=="ND")
109  h2 = (TH1D*)f2->Get(det_caf+"_"+datamc+"/"+caf_intr[intr]);
110  else if(datamc=="mc" && det[idet]=="FD")
111  h2 = (TH1D*)f2->Get(det_caf+"_extrap_oscil_BF/"+caf_intr[intr]);
112 
113  hcaf = CorrectRange(h2, PID[ipid]);
114  //scale fnex by ratio of caf/fnex
115  hfnex->Scale(hcaf->Integral()/hfnex->Integral());
116 
117  //POT and best fit information.
118  //std::cout << " " << std::endl;
119  //std::cout << "--------------------- POT and BF parameters (START) ----------------" << std::endl;
120  TVectorT<double> pot = (TVectorT<double>)f2->Get("pot_nd_fd");
121  //pot->Print();
122  nd_pot = pot[0];
123  fd_pot = pot[1];
124  // cout << "ND POT = " << nd_pot << endl;
125  // cout << "FD POT = " << fd_pot << endl;
126 
127  TVectorT<double> par = (TVectorT<double>)f2->Get("best_dcp_th13_th23_dmsq32");
128  //par->Print();
129  dcp = par[0];
130  th13 = par[1];
131  th23 = par[2];
132  dmsq32 = par[3];
133  // cout << "dCP = " << dcp << endl;
134  // cout << "th31 = " << th13 << endl;
135  // cout << "th23 = " << th23 << endl;
136  // cout << "dmsq32 = " << dmsq32 << endl;
137  //std::cout << "--------------------- POT and BF parameters (END) ---------------" << std::endl;
138  //std::cout << " " << std::endl;
139 
140  //new TCanvas;
141  c1 = new TCanvas("c1"," ",800,600);
142  // gPad->SetLogy();
143  //c1->SetGridy();
144  pad1 = new TPad("pad1","",0.05,0.20,0.95,1);
145  pad2 = new TPad("pad2","",0.05,0.05,0.95,0.25);
146  pad1->Draw();
147  pad2->Draw();
148  pad1->cd();
149  gPad->SetLeftMargin(0.15);
150 
151  //set y maximum range
152  if(hfnex->GetBinContent(hfnex->GetMaximumBin()) > hcaf->GetBinContent(hcaf->GetMaximumBin()) )
153  ymax = 1.1*hfnex->GetBinContent(hfnex->GetMaximumBin());
154  else if(hfnex->GetBinContent(hfnex->GetMaximumBin()) <= hcaf->GetBinContent(hcaf->GetMaximumBin()) )
155  ymax = 1.1*hcaf->GetBinContent(hcaf->GetMaximumBin());
156 
157  hfnex->Draw("HIST same");
158  hcaf->Draw("HIST same");
159  hfnex->GetXaxis()->SetRangeUser(0, 5);
160  hfnex->GetYaxis()->SetRangeUser(0, ymax);
161  hfnex->GetYaxis()->SetTitle("Events");
162  hfnex->GetYaxis()->CenterTitle();
163  hfnex->GetYaxis()->SetTitleSize(0.06);
164  hfnex->GetYaxis()->SetTitleOffset(1.);
165  hfnex->GetYaxis()->SetLabelSize(0.05);
166  hfnex->GetYaxis()->SetDecimals();
167  hfnex->GetXaxis()->SetDecimals();
168 
169  hcaf->SetLineColor(kRed);
170  hcaf->SetFillColor(0);
171  hcaf->SetLineWidth(2);
172  hcaf->SetLineStyle(2);
173 
174  hfnex->SetLineColor(kBlue);
175  hfnex->SetFillColor(0);
176  hfnex->SetLineWidth(2);
177 
178  legend = new TLegend(0.61, 0.59, 0.86, 0.87);
179  legend->SetBorderSize(0);
180  legend->SetFillStyle(0);
181  legend->AddEntry((TObject*)0, PID[ipid]+ " "+"CVN"," ");
182  if(datamc=="mc")
183  legend->AddEntry((TObject*)0, det[idet]+" "+datamc+ " "+caf_intr[intr]," ");
184  else if(datamc=="data")
185  legend->AddEntry((TObject*)0, det[idet]+" "+datamc," ");
186  legend->AddEntry(hfnex, "FNEX", "l");
187  legend->AddEntry(hcaf, "CAFAna", "l");
188  legend->Draw();
189  //cout << "bins fnex = " << hfnex->GetNbinsX() << endl;
190  //cout << "bins caf = " << hcaf->GetNbinsX() << endl;
191 
192  //ratio plot
193  hist_ratio = new TH1F("hist_ratio", "", hcaf->GetNbinsX(), hcaf->GetXaxis()->GetXmin(), hcaf->GetXaxis()->GetXmax());
194  hist_ratio->Sumw2();
195  legend->AddEntry(hist_ratio, "FNEX/CAF", "lep");
196  ratio = 0.;
197  ratio_error = 0.;
198  for(int i = 0; i < hcaf->GetNbinsX()+1; i++){
199  if(hcaf->GetBinContent(i)==0) ratio = 0.;
200  else ratio = hfnex->GetBinContent(i)/hcaf->GetBinContent(i);
201 
202  if(hcaf->GetBinContent(i)==0) ratio_error = 0.;
203  else ratio_error = error_division(hfnex->GetBinContent(i), hcaf->GetBinContent(i));
204 
205  /* cout << "bin = " << i << endl;
206  cout << "fnex content = " << hfnex->GetBinContent(i) << endl;
207  cout << "hcaf content = " << hcaf->GetBinContent(i) << endl;
208  cout << "ratio (fnex/caf) = " << ratio << endl;
209  cout << "ratio error (fnex/caf) = " << ratio_error << endl;
210  cout << " " << endl;
211  */
212  hist_ratio->SetBinContent(i, ratio);
213  hist_ratio->SetBinError(i, ratio_error);
214  }//end of bin by bin ratio
215 
216  pad2->cd();
217  gPad->SetBottomMargin(0.30);
218  gPad->SetLeftMargin(0.15);
219  hist_ratio->Draw("P");
220  hist_ratio->SetLineColor(kBlack);
221  hist_ratio->SetMarkerColor(kBlack);
222  hist_ratio->SetMarkerStyle(20);
223  hist_ratio->GetYaxis()->SetTitle("Ratio");
224  hist_ratio->GetYaxis()->CenterTitle();
225  hist_ratio->GetYaxis()->SetTitleSize(0.23);
226  hist_ratio->GetYaxis()->SetTitleOffset(0.20);
227  hist_ratio->GetYaxis()->SetLabelSize(0.17);
228  hist_ratio->GetYaxis()->SetDecimals();
229  hist_ratio->GetYaxis()->SetRangeUser(0, 3);
230  hist_ratio->GetXaxis()->SetTitle("E^{reco}_{#nu}");
231  hist_ratio->GetXaxis()->CenterTitle();
232  hist_ratio->GetXaxis()->SetTitleSize(0.23);
233  hist_ratio->GetXaxis()->SetTitleOffset(0.45);
234  hist_ratio->GetXaxis()->SetLabelSize(0.17);
235 
236  legend_ratio = new TLegend(0.65, 0.45, 0.93, 0.97);
237  legend_ratio->SetBorderSize(0);
238  legend_ratio->SetFillStyle(0);
239  legend_ratio->Draw();
240  one = new TLine(0, 1, 5, 1);
241  one->SetLineWidth(2);
242  one->SetLineStyle(7);
243  one->Draw();
244 
245  c1->cd();
246  if(datamc=="mc")
247  c1->SaveAs("Images/fnexvscaf/"+PID[ipid]+"_PID_"+det[idet]+"_"+datamc+"_"+caf_intr[intr]+".pdf");
248  else if(datamc=="data")
249  c1->SaveAs("Images/fnexvscaf/"+PID[ipid]+"_PID_"+det[idet]+"_"+datamc+".pdf");
250 
251  // }//end of loop over all interactions.
252  }//end of loop over all PIDs
253  }//end of loop over detector type
254 }//end of main program
255 //--------------------------------------------------------------//
256 TH1D* CorrectRange(TH1D* h, TString PID){
257  int firstbin = 0;
258  if(PID=="Low") firstbin = 0;
259  if(PID=="Mid") firstbin = 10;
260  if(PID=="High") firstbin = 20;
261  TH1D* hlow = new TH1D("hlow", "", 10, 0, 5);
262  for(int i = 0; i < h->GetNbinsX()++1; i++){
263  hlow->SetBinContent(i, h->GetBinContent(i+firstbin));
264  }
265  return hlow;
266 }//end of CorrectRange
267 //---------------------------------------------------------------//
268 double error_division(double num, double deno){
269 
270  return (num/deno)*std::sqrt(std::pow(std::sqrt(num)/num, 2) + std::pow(std::sqrt(deno)/deno, 2));
271 
272 }//end of error_division(double num, double deno)
double th23
Definition: runWimpSim.h:98
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
constexpr T pow(T x)
Definition: pow.h:75
Int_t par
Definition: SimpleIterate.C:24
Float_t f2
bool datamc
Definition: fnexvscaf.C:20
nue_fnex_vs_caf_one_to_one()
Double_t ymax
Definition: plot.C:25
#define pot
Float_t f1
TH1F * h2
Definition: plot.C:45
double dmsq32
int num
Definition: f2_nu.C:119
TH1D * CorrectRange(TH1D *h, TString PID)
c1
Definition: demo5.py:24
auto one()
Definition: PMNS.cxx:49
TPad * pad2
Definition: analysis.C:13
double error_division(double num, double deno)
double th13
Definition: runWimpSim.h:98
TPad * pad1
Definition: analysis.C:13