plot_datamc_ND_numu_REW.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_REW(const std::string weight_name, const std::string fname = "datamc_ND_numu_kinematics_REW.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  const char* wName = weight_name.c_str();
94 
95 
96  std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
97  spects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()).release();
98 
99  preds[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/pred_%s_%s", selName, wName, varName).Data()).release();
100 
101  if(selIdx == 0 && varIdx == 0)
102  std::cout << "data, MC POT: "
103  << spects[selIdx][varIdx]->POT() << " "
104  << preds[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
105 
106  } // end for varIdx
107  } // end for selIdx
108 
109  std::ofstream tfile;
110 
111  tfile.open("plots/TableOut-numuND-REW.tex", ios::out);
112 
113  tfile << "\\begin{tabular}{llllll} \n";
114  tfile << "\\multicolumn{6}{c}{numuND} \\\\ \n";
115  tfile << "Sample & \\%Data/MC & Data & MC & NC & $\\nu_\\mu$ CC \\\\ \n";
116  tfile << "\\hline \n";
117 
118  //loop over the different selector/cut levels to plot each one
119  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
120  const std::string selName = selNames[selIdx];
121  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
122 
123  const std::string varName = defs[varIdx].name;
124 
125  //play with the binning choice for given variables:
126  int rebinFactor=1;
127 
128  if(varName=="recoE"){
129  rebinFactor=2;
130  } else if(varName=="recoQ2"){
131  rebinFactor=10;
132  } else {
133  rebinFactor=20;
134  }
135 
136  //create a canvas
137  TCanvas *c2 =new TCanvas("c2","c2",500,400);
138 
139  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
140  padAbs -> SetBottomMargin(0);
141  padAbs->SetFillStyle(0);
142  padAbs -> SetLeftMargin(0.15);
143 
144  padAbs -> Draw();
145  padAbs -> cd();
146 
147  double legendxlow=0.6;
148  double legendxhigh=0.95;
149 
150  //Pull down the variable we want
151  Spectrum* spect = spects[selIdx][varIdx];
152  IPrediction* pred = preds[selIdx][varIdx];
153  const double kPOTFD = 8.09e20;
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  string CompName = selName.c_str();
192  CompName.erase(0,7);
193 
194  tfile.precision(4);
195  tfile << std::fixed;
196  tfile << weight_name << " REW." << " & " << (intData/intMC) << " & ";
197  tfile.precision(0);
198  tfile << intData << " & ";
199  tfile.precision(0);
200  tfile << intMC << " & ";
201  tfile.precision(0);
202  tfile << intNC << " & ";
203  tfile.precision(0);
204  tfile << intCC << " \\\\ \n";
205  tfile << "\\hline \n";
206  tfile << "\\end{tabular}";
207  tfile.close();
208 
209  hMC->GetXaxis()->CenterTitle();
210  hMC->GetYaxis()->CenterTitle();
211  hMC->GetXaxis()->SetDecimals();
212  hMC->GetYaxis()->SetDecimals();
213  hMC->GetYaxis()->SetTitleOffset(1.15);
214  hMC->GetXaxis()->SetLabelSize(.0);
215  hMC->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", kPOTFD/1E20));
216 
217  if(varName=="recoE" || varIdx == 0){
218  hMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
219  hData->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
220  }
221 
222  if(varName=="recoQ2" || varName=="trueQ2"){
223  hMC->GetXaxis()->SetRangeUser(0,2);
224  }
225 
226  //Draw the error band, MC, Data, MC components, and then Data one last time
227 
228  if(varName=="cvne_zoom"||varName=="shw_len"){
229  legendxlow=0.2;
230  legendxhigh=0.55;
231  }
232 
233  hMC->GetXaxis()->SetNdivisions(405,kFALSE);
234 
235  hMC ->SetMinimum(.001);
236  hMC ->SetMaximum(1.1*TMath::Max(hData->GetMaximum(),hMC->GetMaximum()));
237 
238  hMC ->DrawCopy("hist");
239  hData->Draw("ep same");
240 
241  hNC->Draw("hist same");
242  hCC->Draw("hist same");
243  hBeam->Draw("hist same");
244  hMC->Draw("hist same");
245  hData->Draw("ep same");
246 
247  //Fix up the axis
248  padAbs->RedrawAxis();
249  padAbs->Update();
250 
251  //Add the legend and preliminary
252  Legend(legendxlow, .5, legendxhigh, .85);
253 
254  std::string wTitle = "numuND "+weight_name+" reweighted";
255 
256  TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
257  selTitle->SetNDC();
258  selTitle->SetTextSize(2/30.);
259  selTitle->SetTextAlign(32);
260  selTitle->Draw();
261 
262  c2->cd();
263  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
264  padRatio -> SetTopMargin(0);
265  padRatio -> SetBottomMargin(0.25);
266  padRatio -> SetFillStyle(0);
267  padRatio -> SetLeftMargin(0.15);
268 
269  padRatio -> Draw();
270  padRatio -> cd();
271 
272  TH1* hRatioMC = (TH1*)hData -> Clone("hRatioMC");
273  double b = 5.0;
274 
275  if(varName=="recoQ2"){b=2.0;}
276  if(varName=="CosNumi" || varName=="PtP"){b=1.0;}
277 
278  TLine* lOne = new TLine(0,1.0,b,1.0);
279  lOne -> SetLineStyle(2);
280 
281  hRatioMC -> Divide(hMC);
282 
283  TH1* hRatio = (TH1*)hRatioMC -> Clone("hRatio");
284 
285  hRatioMC->GetXaxis()->CenterTitle();
286  hRatioMC->GetYaxis()->CenterTitle();
287  hRatioMC->GetYaxis()->SetLabelSize(3.7/30.);
288  hRatioMC->GetXaxis()->SetLabelSize(3.7/30.);
289  hRatioMC->GetXaxis()->SetLabelOffset(.04);
290  hRatioMC->GetXaxis()->SetTickSize(.06);
291  hRatioMC->GetYaxis()->SetTitleSize(4.5/30.);
292  hRatioMC->GetXaxis()->SetTitleSize(4.5/30.);
293  hRatioMC->GetYaxis()->SetRangeUser(0.5,1.5);
294  hRatioMC->SetMarkerSize(.1);
295  hRatioMC->GetXaxis()->SetDecimals();
296  hRatioMC->GetYaxis()->SetDecimals();
297  hRatioMC->GetYaxis()->SetTitleOffset(0.4);
298  hRatioMC->GetYaxis()->SetTitle("Data/MC");
299 
300  if(varName=="recoQ2" || varName=="trueQ2"){
301  hRatioMC->GetXaxis()->SetRangeUser(0,2);
302  }
303 
304  hRatioMC->GetYaxis()->SetNdivisions(502,kFALSE);
305  hRatioMC->GetXaxis()->SetNdivisions(405,kFALSE);
306 
307  hRatioMC -> Draw("ep");
308  lOne -> Draw("same");
309 
310  c2->Print(("plots/"+selName+"_"+varName+"_REW.png").c_str());
311  } // end for varIdx
312  } // end for selIdx
313 }
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)
void plot_datamc_ND_numu_REW(const std::string weight_name, const std::string fname="datamc_ND_numu_kinematics_REW.root")
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
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
void Legend(double x0, double y0, double x1, double y1)
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 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)