plot_ND_numu_NOMvsREW.C
Go to the documentation of this file.
2 
4 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Cuts/NueCutsFirstAna.h"
8 #include "CAFAna/Cuts/NueCutsSecondAna.h"
9 
10 #include "CAFAna/Vars/Vars.h"
12 #include "CAFAna/Analysis/Style.h"
13 using namespace ana;
14 
15 #include "OscLib/IOscCalc.h"
16 
17 #include "TCanvas.h"
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TLatex.h"
21 #include "TLine.h"
22 #include "TLegend.h"
23 #include "TPad.h"
24 #include "TStyle.h"
25 
26 #include <iostream>
27 #include <fstream>
28 
29 //make a pretty legend
30 void Legend(double x0, double y0, double x1, double y1)
31 {
32  // Doesn't handle log-y well
33  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
34  TLegend* leg = new TLegend(x0, y0, x1, y1);
35  leg->SetFillStyle(0);
36 
37  TH1* dummy = new TH1F("", "", 1, 0, 1);
38  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
39  dummyFill->SetLineColor(kTotalMCColor);
40  dummyFill->SetFillColor(kTotalMCErrorBandColor);
41 
42  dummy->SetMarkerStyle(kFullCircle);
43 
44  leg->AddEntry(dummy->Clone(), "ND data NOM.", "p");
45  dummy->SetMarkerStyle(24);
46  leg->AddEntry(dummy->Clone(), "ND data REW.", "p");
47  dummy->SetLineColor(kGray+2);
48  dummy->SetMarkerColor(kGray+2);
49  // leg->AddEntry(dummy->Clone(), "ND data stagger", "lep");
50  dummy->SetLineColor(kTotalMCColor);
51  leg->AddEntry(dummy->Clone(), "Total MC", "l");
52  //leg->AddEntry(dummyFill->Clone(), "Flux Uncert.", "fl");
53  dummy->SetLineColor(kNCBackgroundColor);
54  leg->AddEntry(dummy->Clone(), "NC", "l");
55  dummy->SetLineColor(kBeamNueBackgroundColor);
56  leg->AddEntry(dummy->Clone(), "Beam #nu_{e} CC", "l");
57  dummy->SetLineColor(kNumuBackgroundColor);
58  leg->AddEntry(dummy->Clone(), "#nu_{#mu} CC", "l");
59  dummy->SetLineColor(1);
60  leg->AddEntry(dummy->Clone(), "NOMINAL", "l");
61  dummy->SetLineStyle(3);
62  leg->AddEntry(dummy->Clone(), "REWEIGHTED", "l");
63 
64  leg->Draw();
65 }
66 
67 void plot_ND_numu_NOMvsREW(const std::string weight_name, const std::string fname_nom = "datamc_ND_numu_kinematics.root", const std::string fname = "datamc_ND_numu_kinematics_REW.root")
68 {
69  struct HistDef
70  {
72  };
73 
74  //which variables do you actually want to look at?
75  const int kNumVars =4;
76  const HistDef defs[kNumVars] = {
77  "recoE",
78  "recoQ2",
79  "PtP",
80  "CosNumi"
81  };
82 
83  //which cut tiers do you want to look at?
84  const int kNumSels = 1;
85  const std::string selNames[kNumSels] = {
86  "numuND"
87  };
88 
89  Spectrum* spects_nom[kNumSels][kNumVars];
90  IPrediction* preds_nom[kNumSels][kNumVars];
91 
92  Spectrum* spects[kNumSels][kNumVars];
93  IPrediction* preds[kNumSels][kNumVars];
94 
96 
97  //pick up spectra/prediction objects from the input NOMINAL root file
98  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
99  const char* selName = selNames[selIdx].c_str();
100  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
101  const char* varName = defs[varIdx].name.c_str();
102 
103  std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
104  spects_nom[selIdx][varIdx] = LoadFromFile<Spectrum>(fname_nom, TString::Format("%s/spect_%s", selName, varName).Data()).release();
105 
106  preds_nom[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_nom, TString::Format("%s/pred_%s", selName, varName).Data()).release();
107 
108  if(selIdx == 0 && varIdx == 0)
109  std::cout << "data, MC POT: "
110  << spects_nom[selIdx][varIdx]->POT() << " "
111  << preds_nom[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
112 
113  } // end for varIdx
114  } // end for selIdx
115 
116  //pick up spectra/prediction objects from the input REWEIGHTED root file
117  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
118  const char* selName = selNames[selIdx].c_str();
119  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
120  const char* varName = defs[varIdx].name.c_str();
121  const char* wName = weight_name.c_str();
122 
123 
124  std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
125  spects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()).release();
126 
127  preds[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/pred_%s_%s", selName, wName, varName).Data()).release();
128 
129  if(selIdx == 0 && varIdx == 0)
130  std::cout << "data, MC POT: "
131  << spects[selIdx][varIdx]->POT() << " "
132  << preds[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
133 
134  } // end for varIdx
135  } // end for selIdx
136 
137  std::ofstream tfile;
138 
139  tfile.open("plots/TableOut-numuND-REW.tex", ios::out);
140 
141  tfile << "\\begin{tabular}{llllll} \n";
142  tfile << "\\multicolumn{6}{c}{numuND} \\\\ \n";
143  tfile << "Sample & \\%Data/MC & Data & MC & NC & $\\nu_\\mu$ CC \\\\ \n";
144  tfile << "\\hline \n";
145 
146  const double kPOTFD = 8.09e20;
147 
148  //loop over the different selector/cut levels to plot each one
149  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
150  const std::string selName = selNames[selIdx];
151  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
152 
153  const std::string varName = defs[varIdx].name;
154 
155  //play with the binning choice for given variables:
156  int rebinFactor=1;
157 
158  if(varName=="recoE"){
159  rebinFactor=2;
160  } else if(varName=="recoQ2"){
161  rebinFactor=10;
162  } else {
163  rebinFactor=20;
164  }
165 
166  //create a canvas
167  TCanvas *c2 =new TCanvas("c2","c2",500,400);
168 
169  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
170  padAbs -> SetBottomMargin(0);
171  padAbs->SetFillStyle(0);
172  padAbs -> SetLeftMargin(0.15);
173 
174  padAbs -> Draw();
175  padAbs -> cd();
176 
177  double legendxlow=0.55;
178  double legendxhigh=0.9;
179 
180  //Pull down the variable we want
181  Spectrum* spect_nom = spects_nom[selIdx][varIdx];
182  IPrediction* pred_nom = preds_nom[selIdx][varIdx];
183 
184  Spectrum* spect = spects[selIdx][varIdx];
185  IPrediction* pred = preds[selIdx][varIdx];
186 
187  //Get components of the prediction and store them in histograms
188  //NOMINAL
189  TH1* hMC_nom = pred_nom->Predict(&noosc).ToTH1(kPOTFD);
190  hMC_nom->SetLineColor(kTotalMCColor);
191  hMC_nom->Rebin(rebinFactor);
192 
193  TH1* hNC_nom = pred_nom->PredictComponent(&noosc,
195  Current::kNC,
196  Sign::kBoth).ToTH1(kPOTFD);
197  hNC_nom->SetLineColor(kNCBackgroundColor);
198  hNC_nom->Rebin(rebinFactor);
199 
200  TH1* hCC_nom = pred_nom->PredictComponent(&noosc,
202  Current::kCC,
203  Sign::kBoth).ToTH1(kPOTFD);
204  hCC_nom->SetLineColor(kNumuBackgroundColor);
205  hCC_nom->Rebin(rebinFactor);
206 
207  TH1* hBeam_nom = pred_nom->PredictComponent(&noosc,
209  Current::kCC,
210  Sign::kBoth).ToTH1(kPOTFD);
211  hBeam_nom->SetLineColor(kBeamNueBackgroundColor);
212  hBeam_nom->Rebin(rebinFactor);
213 
214  TH1* hData_nom = spect_nom->ToTH1(kPOTFD);
215  hData_nom->SetMarkerStyle(kFullCircle);
216  hData_nom->Rebin(rebinFactor);
217 
218  //REWEIGHTED
219  TH1* hMC = pred->Predict(&noosc).ToTH1(kPOTFD);
220  hMC->SetLineColor(kTotalMCColor);
221  hMC->SetLineStyle(3);
222  hMC->Rebin(rebinFactor);
223 
224  TH1* hNC = pred->PredictComponent(&noosc,
226  Current::kNC,
227  Sign::kBoth).ToTH1(kPOTFD);
228  hNC->SetLineColor(kNCBackgroundColor);
229  hNC->SetLineStyle(3);
230  hNC->Rebin(rebinFactor);
231 
232  TH1* hCC = pred->PredictComponent(&noosc,
234  Current::kCC,
235  Sign::kBoth).ToTH1(kPOTFD);
236  hCC->SetLineColor(kNumuBackgroundColor);
237  hCC->SetLineStyle(3);
238  hCC->Rebin(rebinFactor);
239 
240  TH1* hBeam = pred->PredictComponent(&noosc,
242  Current::kCC,
243  Sign::kBoth).ToTH1(kPOTFD);
244  hBeam->SetLineColor(kBeamNueBackgroundColor);
245  hBeam->SetLineStyle(3);
246  hBeam->Rebin(rebinFactor);
247 
248  TH1* hData = spect->ToTH1(kPOTFD);
249  hData->SetMarkerStyle(24);
250  hData->Rebin(rebinFactor);
251 
252  double intData = round(100.0*hData->Integral())/100.0;
253  double intMC = round(100.0*hMC->Integral())/100.0;
254  double intNC = round(100.0*hNC->Integral())/100.0;
255  double intCC = round(100.0*hCC->Integral())/100.0;
256  double intBeam = round(100.0*hBeam->Integral())/100.0;
257 
258  string CompName = selName.c_str();
259  CompName.erase(0,7);
260 
261  tfile.precision(4);
262  tfile << std::fixed;
263  tfile << weight_name << " REW." << " & " << (intData/intMC) << " & ";
264  tfile.precision(0);
265  tfile << intData << " & ";
266  tfile.precision(0);
267  tfile << intMC << " & ";
268  tfile.precision(0);
269  tfile << intNC << " & ";
270  tfile.precision(0);
271  tfile << intCC << " \\\\ \n";
272  tfile << "\\hline \n";
273  tfile << "\\end{tabular}";
274  tfile.close();
275 
276  hMC->GetXaxis()->CenterTitle();
277  hMC->GetYaxis()->CenterTitle();
278  hMC->GetXaxis()->SetDecimals();
279  hMC->GetYaxis()->SetDecimals();
280  hMC->GetYaxis()->SetTitleOffset(1.15);
281  hMC->GetXaxis()->SetLabelSize(.0);
282  hMC->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", kPOTFD/1E20));
283 
284  if(varName=="recoE" || varIdx == 0){
285  hMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
286  hData->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
287  }
288 
289  if(varName=="recoQ2" || varName=="trueQ2"){
290  hMC->GetXaxis()->SetRangeUser(0,2);
291  }
292 
293  //Draw the error band, MC, Data, MC components, and then Data one last time
294 
295  if(varName=="cvne_zoom"||varName=="shw_len"){
296  legendxlow=0.2;
297  legendxhigh=0.55;
298  }
299 
300  hMC->GetXaxis()->SetNdivisions(405,kFALSE);
301 
302  hMC ->SetMinimum(.001);
303  hMC ->SetMaximum(1.1*TMath::Max(TMath::Max(hData->GetMaximum(),hMC->GetMaximum()),TMath::Max(hData_nom->GetMaximum(),hMC_nom->GetMaximum())));
304 
305  hMC ->DrawCopy("hist");
306  hData_nom->Draw("ep same");
307 
308  hNC_nom->Draw("hist same");
309  hCC_nom->Draw("hist same");
310  hBeam_nom->Draw("hist same");
311  hMC_nom->Draw("hist same");
312  hData_nom->Draw("ep same");
313 
314  hData->Draw("ep same");
315 
316  hNC->Draw("hist same");
317  hCC->Draw("hist same");
318  hBeam->Draw("hist same");
319  hMC->Draw("hist same");
320  hData->Draw("ep same");
321 
322  //Fix up the axis
323  padAbs->RedrawAxis();
324  padAbs->Update();
325 
326  //Add the legend and preliminary
327  Legend(legendxlow, .35, legendxhigh, .85);
328 
329  std::string wTitle = "numuND "+weight_name+" reweighted";
330 
331  TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
332  selTitle->SetNDC();
333  selTitle->SetTextSize(2/30.);
334  selTitle->SetTextAlign(32);
335  selTitle->Draw();
336 
337  c2->cd();
338  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
339  padRatio -> SetTopMargin(0);
340  padRatio -> SetBottomMargin(0.25);
341  padRatio -> SetFillStyle(0);
342  padRatio -> SetLeftMargin(0.15);
343 
344  padRatio -> Draw();
345  padRatio -> cd();
346 
347  TH1* hRatioMC = (TH1*)hData -> Clone("hRatioMC");
348  TH1* hRatioMC_nom = (TH1*)hData_nom -> Clone("hRatioMC_nom");
349  double b = 5.0;
350 
351  if(varName=="recoQ2"){b=2.0;}
352  if(varName=="CosNumi" || varName=="PtP"){b=1.0;}
353 
354  TLine* lOne = new TLine(0,1.0,b,1.0);
355  lOne -> SetLineStyle(2);
356 
357  hRatioMC -> Divide(hMC);
358  hRatioMC_nom -> Divide(hMC_nom);
359 
360  hRatioMC -> SetLineStyle(3);
361 
362  //TH1* hRatio = (TH1*)hRatioMC -> Clone("hRatio");
363 
364  hRatioMC->GetXaxis()->CenterTitle();
365  hRatioMC->GetYaxis()->CenterTitle();
366  hRatioMC->GetYaxis()->SetLabelSize(3.7/30.);
367  hRatioMC->GetXaxis()->SetLabelSize(3.7/30.);
368  hRatioMC->GetXaxis()->SetLabelOffset(.04);
369  hRatioMC->GetXaxis()->SetTickSize(.06);
370  hRatioMC->GetYaxis()->SetTitleSize(4.5/30.);
371  hRatioMC->GetXaxis()->SetTitleSize(4.5/30.);
372  hRatioMC->GetYaxis()->SetRangeUser(0.5,1.5);
373  hRatioMC->SetMarkerSize(.1);
374  hRatioMC_nom->SetMarkerSize(.1);
375  hRatioMC->GetXaxis()->SetDecimals();
376  hRatioMC->GetYaxis()->SetDecimals();
377  hRatioMC->GetYaxis()->SetTitleOffset(0.4);
378  hRatioMC->GetYaxis()->SetTitle("Data/MC");
379 
380  if(varName=="recoQ2" || varName=="trueQ2"){
381  hRatioMC->GetXaxis()->SetRangeUser(0,2);
382  }
383 
384  hRatioMC->GetYaxis()->SetNdivisions(502,kFALSE);
385  hRatioMC->GetXaxis()->SetNdivisions(405,kFALSE);
386 
387  hRatioMC -> DrawCopy("ep");
388  hRatioMC_nom -> Draw("ep same");
389  hRatioMC -> Draw("ep same");
390  lOne -> Draw("same");
391 
392  c2->Print(("plots/"+selName+"_"+varName+"_REW.png").c_str());
393  } // end for varIdx
394  } // end for selIdx
395 }
tree Draw("slc.nhit")
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
ratio_hxv Divide(hxv, goal_hxv)
const int kNumVars
Definition: vars.h:14
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Float_t y1[n_points_granero]
Definition: compare.C:5
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:209
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
hmean SetLineStyle(2)
std::string name
Definition: NuePlotLists.h:12
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
ntuple SetFillStyle(1001)
const Color_t kNumuBackgroundColor
Definition: Style.h:30
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
c2
Definition: demo5.py:33
void Legend(double x0, double y0, double x1, double y1)
Charged-current interactions.
Definition: IPrediction.h:39
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
double POT() const
Definition: Spectrum.h:231
OStream cout
Definition: OStream.cxx:6
void plot_ND_numu_NOMvsREW(const std::string weight_name, const std::string fname_nom="datamc_ND_numu_kinematics.root", const std::string fname="datamc_ND_numu_kinematics_REW.root")
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const hit & b
Definition: hits.cxx:21
Neutral-current interactions.
Definition: IPrediction.h:40
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:16
All neutrinos, any flavor.
Definition: IPrediction.h:26
const std::string selNames[kNumSels]
Definition: vars.h:46
const Color_t kNCBackgroundColor
Definition: Style.h:22
c cd(1)