plot_datamc_ND_numu.C
Go to the documentation of this file.
2 
4 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Spectrum.h"
9 
10 #include "CAFAna/Vars/Vars.h"
11 #include "CAFAna/Vars/NueVars.h"
12 #include "CAFAna/Analysis/Style.h"
13 using namespace ana;
14 
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", "lep");
45  dummy->SetLineColor(kGray+2);
46  dummy->SetMarkerColor(kGray+2);
47  // leg->AddEntry(dummy->Clone(), "ND data stagger", "lep");
48  dummy->SetLineColor(kTotalMCColor);
49  leg->AddEntry(dummy->Clone(), "Total MC", "l");
50  //leg->AddEntry(dummyFill->Clone(), "Flux Uncert.", "fl");
51  dummy->SetLineColor(kNCBackgroundColor);
52  leg->AddEntry(dummy->Clone(), "NC", "l");
53  dummy->SetLineColor(kBeamNueBackgroundColor);
54  leg->AddEntry(dummy->Clone(), "Beam #nu_{e} CC", "l");
55  dummy->SetLineColor(kNumuBackgroundColor);
56  leg->AddEntry(dummy->Clone(), "#nu_{#mu} CC", "l");
57 
58  leg->Draw();
59 }
60 
61 void plot_datamc_ND_numu(std::string fname = "datamc_ND_numu_kinematics.root")
62 {
63  struct HistDef
64  {
66  };
67 
68  //which variables do you actually want to look at?
69  const int kNumVars =4;
70  const HistDef defs[kNumVars] = {
71  "recoE",
72  "recoQ2",
73  "PtP",
74  "CosNumi"
75  };
76 
77  //which cut tiers do you want to look at?
78  const int kNumSels = 1;
79  const std::string selNames[kNumSels] = {
80  "numuND"
81  };
82 
84  IPrediction* preds[kNumSels][kNumVars];
85 
87 
88  //pick up spectra/prediction objects from the input root file
89  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
90  const char* selName = selNames[selIdx].c_str();
91  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
92  const char* varName = defs[varIdx].name.c_str();
93 
94  std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
95  spects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s", selName, varName).Data()).release();
96 
97  preds[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/pred_%s", selName, varName).Data()).release();
98 
99  if(selIdx == 0 && varIdx == 0)
100  std::cout << "data, MC POT: "
101  << spects[selIdx][varIdx]->POT() << " "
102  << preds[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
103 
104  } // end for varIdx
105  } // end for selIdx
106 
107  std::ofstream tfile;
108 
109  tfile.open("plots/TableOut-numuND.tex", ios::out);
110 
111  tfile << "\\begin{tabular}{llllll} \n";
112  tfile << "\\multicolumn{6}{c}{numuND} \\\\ \n";
113  tfile << "Sample & \\%Data/MC & Data & MC & NC & $\\nu_\\mu$ CC \\\\ \n";
114  tfile << "\\hline \n";
115 
116  //loop over the different selector/cut levels to plot each one
117  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
118  const std::string selName = selNames[selIdx];
119  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
120 
121  const std::string varName = defs[varIdx].name;
122 
123  //play with the binning choice for given variables:
124  int rebinFactor=1;
125 
126  if(varName=="recoE"){
127  rebinFactor=2;
128  } else if(varName=="recoQ2"){
129  rebinFactor=10;
130  } else {
131  rebinFactor=20;
132  }
133 
134  //create a canvas
135  TCanvas *c2 =new TCanvas("c2","c2",500,400);
136 
137 
138  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
139  padAbs -> SetBottomMargin(0);
140  padAbs->SetFillStyle(0);
141  padAbs -> SetLeftMargin(0.15);
142 
143  padAbs -> Draw();
144  padAbs -> cd();
145 
146  double legendxlow=0.6;
147  double legendxhigh=0.95;
148  const double kPOTFD = 8.09e20;
149 
150  //Pull down the variable we want
151  Spectrum* spect = spects[selIdx][varIdx];
152 
153  IPrediction* pred = preds[selIdx][varIdx];
154 
155  //Get components of the prediction and store them in histograms
156  TH1* hMC = pred->Predict(&noosc).ToTH1(kPOTFD);
157  hMC->SetLineColor(kTotalMCColor);
158  hMC->Rebin(rebinFactor);
159 
160  TH1* hNC = pred->PredictComponent(&noosc,
162  Current::kNC,
163  Sign::kBoth).ToTH1(kPOTFD);
164  hNC->SetLineColor(kNCBackgroundColor);
165  hNC->Rebin(rebinFactor);
166 
167  TH1* hCC = pred->PredictComponent(&noosc,
169  Current::kCC,
170  Sign::kBoth).ToTH1(kPOTFD);
171  hCC->SetLineColor(kNumuBackgroundColor);
172  hCC->Rebin(rebinFactor);
173 
174  TH1* hBeam = pred->PredictComponent(&noosc,
176  Current::kCC,
177  Sign::kBoth).ToTH1(kPOTFD);
178  hBeam->SetLineColor(kBeamNueBackgroundColor);
179  hBeam->Rebin(rebinFactor);
180 
181  TH1* hData = spect->ToTH1(kPOTFD);
182  hData->SetMarkerStyle(kFullCircle);
183  hData->Rebin(rebinFactor);
184 
185  double intData = round(100.0*hData->Integral())/100.0;
186  double intMC = round(100.0*hMC->Integral())/100.0;
187  double intNC = round(100.0*hNC->Integral())/100.0;
188  double intCC = round(100.0*hCC->Integral())/100.0;
189  double intBeam = round(100.0*hBeam->Integral())/100.0;
190 
191  tfile.precision(4);
192  tfile << std::fixed;
193  tfile << "NOMINAL" << " & " << (intData/intMC) << " & ";
194  tfile.precision(0);
195  tfile << intData << " & ";
196  tfile.precision(0);
197  tfile << intMC << " & ";
198  tfile.precision(0);
199  tfile << intNC << " & ";
200  tfile.precision(0);
201  tfile << intCC << " \\\\ \n";
202  tfile << "\\hline \n";
203  tfile << "\\end{tabular}";
204  tfile.close();
205 
206  hMC->GetXaxis()->CenterTitle();
207  hMC->GetYaxis()->CenterTitle();
208  hMC->GetXaxis()->SetDecimals();
209  hMC->GetYaxis()->SetDecimals();
210  hMC->GetYaxis()->SetTitleOffset(1.15);
211  hMC->GetXaxis()->SetLabelSize(.0);
212  hMC->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", kPOTFD/1E20));
213 
214  if(varName=="recoE"){
215  hMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
216  }
217 
218  if(varName=="recoQ2" || varName=="trueQ2"){
219  hMC->GetXaxis()->SetRangeUser(0,2);
220  }
221 
222  //Draw the error band, MC, Data, MC components, and then Data one last time
223 
224  if(varName=="cvne_zoom"||varName=="shw_len"){
225  legendxlow=0.2;
226  legendxhigh=0.55;
227  }
228 
229  hMC->GetXaxis()->SetNdivisions(405,kFALSE);
230 
231  hMC ->SetMinimum(.001);
232  hMC ->SetMaximum(1.1*TMath::Max(hData->GetMaximum(),hMC->GetMaximum()));
233 
234  hMC ->DrawCopy("hist");
235  hData->Draw("ep same");
236 
237  hNC->Draw("hist same");
238  hCC->Draw("hist same");
239  hBeam->Draw("hist same");
240  hMC->Draw("hist same");
241  hData->Draw("ep same");
242 
243  //Fix up the axis
244  padAbs->RedrawAxis();
245  padAbs->Update();
246 
247  //Add the legend and preliminary
248  Legend(legendxlow, .5, legendxhigh, .85);
249 
250  TLatex* selTitle = new TLatex(.935, .95, "numuND");
251  selTitle->SetNDC();
252  selTitle->SetTextSize(2/30.);
253  selTitle->SetTextAlign(32);
254  selTitle->Draw();
255 
256  c2->cd();
257  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
258  padRatio -> SetTopMargin(0);
259  padRatio -> SetBottomMargin(0.25);
260  padRatio -> SetFillStyle(0);
261  padRatio -> SetLeftMargin(0.15);
262 
263  padRatio -> Draw();
264  padRatio -> cd();
265 
266  TH1* hRatioMC = (TH1*)hData -> Clone("hRatioMC");
267  double b = 5.0;
268 
269  if(varName=="recoQ2"){b=2.0;}
270  if(varName=="CosNumi" || varName=="PtP"){b=1.0;}
271 
272  TLine* lOne = new TLine(0,1.0,b,1.0);
273  lOne -> SetLineStyle(2);
274 
275  hRatioMC -> Divide(hMC);
276 
277  TH1* hRatio = (TH1*)hRatioMC -> Clone("hRatio");
278 
279  hRatioMC->GetXaxis()->CenterTitle();
280  hRatioMC->GetYaxis()->CenterTitle();
281  hRatioMC->GetYaxis()->SetLabelSize(3.7/30.);
282  hRatioMC->GetXaxis()->SetLabelSize(3.7/30.);
283  hRatioMC->GetXaxis()->SetLabelOffset(.04);
284  hRatioMC->GetXaxis()->SetTickSize(.06);
285  hRatioMC->GetYaxis()->SetTitleSize(4.5/30.);
286  hRatioMC->GetXaxis()->SetTitleSize(4.5/30.);
287  hRatioMC->GetYaxis()->SetRangeUser(0.5,1.5);
288  hRatioMC->SetMarkerSize(.1);
289  hRatioMC->GetXaxis()->SetDecimals();
290  hRatioMC->GetYaxis()->SetDecimals();
291  hRatioMC->GetYaxis()->SetTitleOffset(0.4);
292  hRatioMC->GetYaxis()->SetTitle("Data/MC");
293 
294  if(varName=="recoQ2" || varName=="trueQ2"){
295  hRatioMC->GetXaxis()->SetRangeUser(0,2);
296  }
297 
298  if(varName=="recoE"){
299  hRatioMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
300  }
301 
302  hRatioMC->GetYaxis()->SetNdivisions(502,kFALSE);
303  hRatioMC->GetXaxis()->SetNdivisions(405,kFALSE);
304 
305  hRatioMC -> Draw("ep");
306  lOne -> Draw("same");
307 
308  c2->Print(("plots/"+selName+"_"+varName+".png").c_str());
309 
310  } // end for varIdx
311  } // end for selIdx
312 }
Pass neutrinos through unchanged.
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
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
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
hmean SetLineStyle(2)
const Color_t kTotalMCErrorBandColor
Definition: Style.h:18
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
ntuple SetFillStyle(1001)
const Color_t kNumuBackgroundColor
Definition: Style.h:31
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
c2
Definition: demo5.py:33
Charged-current interactions.
Definition: IPrediction.h:39
void plot_datamc_ND_numu(std::string fname="datamc_ND_numu_kinematics.root")
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: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
const hit & b
Definition: hits.cxx:21
recTree Draw("energy.numu.trkccE>>precosmics","fdcuts&&preshutdown")
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 Legend(double x0, double y0, double x1, double y1)
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
void spects(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC)
Definition: sensitivity.C:100
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)