plot_3NDvsFD_RHC.C
Go to the documentation of this file.
1 
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/Cut.h"
6 #include "CAFAna/Core/HistAxis.h"
8 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Cuts/NumuCuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
21 #include "CAFAna/Vars/NueVars.h"
22 #include "CAFAna/Vars/Vars.h"
23 
26 #include "CAFAna/Analysis/Style.h"
27 
28 //#include "CAFAna/nue/SecondAna/draw_plots_util.h"
29 
30 using namespace ana;
31 
32 #include "TFile.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TCanvas.h"
36 #include "TGaxis.h"
37 #include "TLatex.h"
38 #include "TLegend.h"
39 #include "TLine.h"
40 #include "TSystem.h"
41 #include "TStyle.h"
42 #include <iostream>
43 #include <iomanip>
44 #include <fstream>
45 
46 //make a pretty legend
47 void Legend(double x0, double y0, double x1, double y1)
48 {
49  // Doesn't handle log-y well
50  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
51  TLegend* leg = new TLegend(x0, y0, x1, y1);
52  leg->SetFillStyle(0);
53  leg->SetHeader("#splitline{#bf{Antineutrino Beam}}{CORE sample: Q^{2} reweighted}");
54  TH1* dummy = new TH1F("", "", 1, 0, 1);
55  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
56  dummyFill->SetLineColor(kTotalMCColor);
57  dummyFill->SetFillColor(kTotalMCErrorBandColor);
58 
59  dummy->SetMarkerStyle(kFullCircle);
60 
61  dummy->SetLineColor(kGray+2);
62  dummy->SetMarkerColor(kGray+2);
63  dummy->SetLineColor(kTotalMCColor);
64  leg->AddEntry(dummy->Clone(), "ND #bar{#nu}_{#mu} MC reweighted", "l");
65  dummy->SetLineColor(1);
66  leg->AddEntry(dummy->Clone(), "ND #bar{#nu}_{#mu} MC nominal", "l");
67 
68  dummy->SetLineColor(kNueSignalColor);
69  leg->AddEntry(dummy->Clone(), "FD MC #bar{#nu}_{e} + WS", "l");
70 
71  leg->Draw();
72 }
73 
74 void plot_3NDvsFD_RHC(const std::string weight_name, const std::string fname_NDnon, const std::string fname_ND, const std::string fname_FD)
75 {
76  gStyle->SetPadTickX(1);
77  gStyle->SetPadTickY(0);
78 
79  TH1::AddDirectory(kFALSE);
80 
81  struct HistDef
82  {
84  };
85 
86  //which variables do you actually want to look at?
87  const int kNumVars = 5;
88  const HistDef defs[kNumVars] = {
89  "trueQ2",
90  "recoQ2",
91  "trueW2",
92  "PtP",
93  "CosNumi"
94  };
95 
96  //which cut tiers do you want to look at?
97  const int kNumSels = 1;
98  const std::string selNamesND[kNumSels] = {
99  "numuND"
100  };
101 
102  const std::string selNamesFD[kNumSels] = {
103  "CORE"
104  };
105 
106  IPrediction* predsND[kNumSels][kNumVars];
107  IPrediction* predsNDnon[kNumSels][kNumVars];
108 
110  osc::OscCalculatorDumb mycalc;
111  //osc::OscCalculatorPMNSOpt mycalc;
112  //ResetOscCalcToDefault(&mycalc);
113 
114  //mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
115  //mycalc.SetdCP(3*TMath::Pi()/2);
116  //mycalc.SetTh23(asin(sqrt(0.404)));
117 
118  //pick up spectra/prediction objects from the input root file
119  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
120  const char* selName = selNamesND[selIdx].c_str();
121  const char* selNameFD = selNamesFD[selIdx].c_str();
122  const char* wName = weight_name.c_str();
123 
124  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
125  const char* varName = defs[varIdx].name.c_str();
126 
127  std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
128 
129  predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s_%s", selNameFD, wName, varName).Data()).release();
130  predsNDnon[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_NDnon, TString::Format("%s/pred_noweight_%s", selNameFD, varName).Data()).release();
131 
132  if(selIdx == 0 && varIdx == 0)
133  std::cout << "ND MC POT: "
134  << predsNDnon[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
135 
136  } // end for varIdx
137  } // end for selIdx
138 
139 
140  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
141  TDirectory * dpred_no[kNumSels][kNumVars];
142  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
143  const std::string selName = selNamesFD[selIdx];
144  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
145  const std::string varName = defs[varIdx].name;
146  if(varName == "recoQ2") {
147  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
148  } else if(varName == "CosNumi"){
149  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
150  } else {
151  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
152  }
153  }
154  }
155 
156  TH1* hFullNDHists[kNumVars-1];
157  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
158  const std::string selName = selNamesFD[selIdx];
159  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
160 
161  const std::string varName = defs[varIdx].name;
162 
163  //play with the binning choice for given variables:
164  int rebinFactor=1;
165 
166  if(varName=="PtP") rebinFactor = 5;
167 
168  //create a canvas
169  TCanvas *c2 =new TCanvas("c2","c2",500,400);
170 
171  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
172  padAbs -> SetBottomMargin(0);
173  padAbs->SetFillStyle(0);
174  padAbs->SetFillColor(0);
175  padAbs -> SetLeftMargin(0.15);
176 
177  padAbs -> Draw();
178  padAbs -> cd();
179 
180  double legendxlow=0.45;
181  double legendxhigh=0.93;
182 
183  if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.17; legendxhigh=0.6;}
184  //if(varName=="CosNumi"){legendxlow=0.17; legendxhigh=0.6;}
185 
186  //Pull down the variable we want
187  IPrediction* predND = predsND[selIdx][varIdx];
188  IPrediction* predNDnon = predsNDnon[selIdx][varIdx];
189 
190  const double kPOT = predND->Predict(&noosc).POT();
191 
192  auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);
193 
194  //Get components of the prediction and store them in histograms
195 
196  TH1* hFD = predFD->PredictComponent(&mycalc,
198  Current::kCC,
199  Sign::kBoth).ToTH1(kPOT);
200 
201  TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);
202  TH1* hAuxnon = predNDnon->Predict(&noosc).ToTH1(kPOT);
203 
204  TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
205  TH1* hNDnon = predNDnon->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAuxnon->Integral());
206 
207  hND->SetLineColor(kTotalMCColor);
208  hNDnon->SetLineColor(1);
209 
210  hFD->SetLineColor(kNueSignalColor);
211  hFD->Rebin(rebinFactor);
212  hND->Rebin(rebinFactor);
213  hNDnon->Rebin(rebinFactor);
214 
215  std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
216 
217  hND->GetXaxis()->CenterTitle();
218  hND->GetYaxis()->CenterTitle();
219  hND->GetXaxis()->SetDecimals();
220  hND->GetYaxis()->SetDecimals();
221  hND->GetYaxis()->SetTitleOffset(.8);
222  hND->GetXaxis()->SetLabelSize(.0);
223  hND->GetYaxis()->SetLabelSize(.06);
224  hND->GetYaxis()->SetTitleSize(.06);
225  hND->SetMaximum(1.1*TMath::Max(hND->GetMaximum(), TMath::Max(hNDnon->GetMaximum(), hFD->GetMaximum())));
226  if(varName=="trueW2"){
227  hND->SetMaximum(1.2*hND->GetMaximum());
228  }
229 
230  hND->GetYaxis()->SetTitle(TString::Format("Events area normalized"));
231 
232  if(varName=="recoQ2" || varName=="trueQ2"){
233  hND->GetXaxis()->SetRangeUser(0,2);
234  }
235 
236  if(varName=="recoE"){
237  hND->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
238  }
239 
240  if(varName=="CosNumi") hND->GetXaxis()->SetRangeUser(0.5,1.0);
241 
242  hND->GetXaxis()->SetNdivisions(405,kFALSE);
243  hND->SetMinimum(.001);
244 
245  hND ->DrawCopy("h");
246 
247  hNDnon->Draw("h same");
248 
249  hND->Draw("h same");
250  hFD->Draw("hist same");
251 
252  //Fix up the axis
253  padAbs->RedrawAxis();
254  padAbs->Update();
255 
256  //Add the legend and preliminary
257  Legend(legendxlow, .47, legendxhigh, .87);
258  std::string wTitle = "FD RHC CORE sample: Q^{2} reweighted";
259 
260  if(weight_name == "Cos") wTitle = "FD RHC CORE sample: Cos reweighted";
261  if(weight_name == "PtP") wTitle = "FD RHC CORE sample: p_{T}/p reweighted";
262  wTitle = "NOvA Simulation";
263 
264  TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
265  selTitle->SetNDC();
266  selTitle->SetTextSize(2/30.);
267  selTitle->SetTextAlign(32);
268  selTitle->SetTextColor(kGray+1);
269  selTitle->Draw();
270 
271  c2->cd();
272  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.0, 1, 0.33);
273  padRatio -> SetTopMargin(0.11);
274  padRatio -> SetBottomMargin(0.33);
275  padRatio -> SetFillStyle(0);
276  padRatio -> SetLeftMargin(0.15);
277 
278  padRatio -> Draw();
279  padRatio -> cd();
280 
281  TH1* hRatio = (TH1*)hND -> Clone("hRatio");
282  TH1* hRationon = (TH1*)hNDnon -> Clone("hRationon");
283 
284  double b=5;
285  double a=0;
286  if(varName=="trueW2") {b = 2.0;}
287  if(varName=="trueQ2") {b = 2.0;}
288  if(varName=="recoQ2") {b = 2.0;}
289  if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
290  if(varName=="CosNumi") a=0.5;
291 
292  TLine* lOne = new TLine(a,1.0,b,1.0);
293  lOne -> SetLineStyle(2);
294 
295  hRatio -> Divide(hFD);
296  hRationon -> Divide(hFD);
297 
298  hRatio->GetXaxis()->CenterTitle();
299  hRatio->GetYaxis()->CenterTitle();
300  hRatio->GetYaxis()->SetLabelSize(.12);
301  hRatio->GetXaxis()->SetLabelSize(.12);
302  hRatio->GetXaxis()->SetTitleOffset(1.);
303  hRatio->GetXaxis()->SetLabelOffset(0.03);
304  hRatio->GetXaxis()->SetTickSize(.1);
305  hRatio->GetYaxis()->SetTitleSize(.12);
306  hRatio->GetXaxis()->SetTitleSize(.12);
307  hRatio->GetYaxis()->SetRangeUser(0,2);
308  if(varName=="CosNumi"){
309  hRatio->GetYaxis()->SetRangeUser(0.5,2.5);
310  }
311  hRatio->SetMarkerSize(.1);
312  hRatio->GetXaxis()->SetDecimals();
313  hRatio->GetYaxis()->SetDecimals();
314  hRatio->GetYaxis()->SetTitleOffset(0.4);
315  hRatio->GetYaxis()->SetTitle("Ratio to FD #bar{#nu}_{e}");
316 
317  if(varName=="recoE"){
318  hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
319  }
320  if(varName=="recoQ2"){
321  hRatio->GetXaxis()->SetTitle("ND Reconstructed Q^2 or FD True Q^2 (GeV^2)");
322  }
323  if(varName=="CosNumi"){
324  hRatio->GetXaxis()->SetTitle("cos #theta");
325  }
326  if(varName=="PtP"){
327  hRatio->GetXaxis()->SetTitle("p_{T}/p");
328  }
329 
330  if(varName=="recoQ2" || varName=="trueQ2"){
331  hRatio->GetXaxis()->SetRangeUser(0,2);
332  }
333  if(varName=="trueW2"){
334  hRatio->GetXaxis()->SetTitle("True W^{2} (GeV^{2})");
335  }
336 
337  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
338  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
339 
340  c2->Update();
341  hRatio-> DrawCopy("hist");
342  lOne -> Draw("same");
343  hRationon->Draw("hist same");
344  hRatio->Draw("hist same");
345  padRatio->RedrawAxis();
346  padRatio->Update();
347 
348  c2->Update();
349  c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".png").c_str());
350  c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".pdf").c_str());
351 
352  } // end for varIdx
353  } // end for selIdx
354 }
355 
Pass neutrinos through unchanged.
void Legend(double x0, double y0, double x1, double y1)
ratio_hxv Divide(hxv, goal_hxv)
const int kNumVars
Definition: vars.h:14
Oscillation analysis framework, runs over CAF files outside of ART.
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
Float_t x1[n_points_granero]
Definition: compare.C:5
hmean SetLineStyle(2)
const Color_t kTotalMCErrorBandColor
Definition: Style.h:18
ntuple SetFillStyle(1001)
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
const Color_t kNueSignalColor
Definition: Style.h:19
c2
Definition: demo5.py:33
Charged-current interactions.
Definition: IPrediction.h:39
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
const double a
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
void plot_3NDvsFD_RHC(const std::string weight_name, const std::string fname_NDnon, const std::string fname_ND, const std::string fname_FD)
double POT() const
Definition: Spectrum.h:263
std::string name
Definition: NuePlotLists.h:12
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
const hit & b
Definition: hits.cxx:21
recTree Draw("energy.numu.trkccE>>precosmics","fdcuts&&preshutdown")
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:17
c cd(1)