plot_NDvsFD_REW.C
Go to the documentation of this file.
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/Cut.h"
5 #include "CAFAna/Core/HistAxis.h"
7 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Cuts/NumuCuts.h"
13 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Vars/NueVars.h"
21 #include "CAFAna/Vars/Vars.h"
22 
25 #include "CAFAna/Analysis/Style.h"
26 
27 //#include "CAFAna/nue/SecondAna/draw_plots_util.h"
28 
29 using namespace ana;
30 
31 #include "TFile.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TCanvas.h"
35 #include "TGaxis.h"
36 #include "TLatex.h"
37 #include "TLegend.h"
38 #include "TLine.h"
39 #include "TSystem.h"
40 #include <iostream>
41 #include <iomanip>
42 #include <fstream>
43 
44 //make a pretty legend
45 void Legend(double x0, double y0, double x1, double y1)
46 {
47  // Doesn't handle log-y well
48  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
49  TLegend* leg = new TLegend(x0, y0, x1, y1);
50  leg->SetFillStyle(0);
51 
52  TH1* dummy = new TH1F("", "", 1, 0, 1);
53  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
54  dummyFill->SetLineColor(kTotalMCColor);
55  dummyFill->SetFillColor(kTotalMCErrorBandColor);
56 
57  dummy->SetMarkerStyle(kFullCircle);
58 
59  dummy->SetLineColor(kGray+2);
60  dummy->SetMarkerColor(kGray+2);
61  dummy->SetLineColor(kTotalMCColor);
62  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} Total MC", "l");
63 
64  dummy->SetLineColor(kNueSignalColor);
65  leg->AddEntry(dummy->Clone(), "FD signal #nu_{e}", "l");
66 
67  leg->Draw();
68 }
69 
70 void plot_NDvsFD_REW(const std::string weight_name, const std::string fname_ND, const std::string fname_FD)
71 {
72  const double kPOT = 8.09e20;
73 
74  struct HistDef
75  {
77  };
78 
79  //which variables do you actually want to look at?
80  const int kNumVars = 5;
81  const HistDef defs[kNumVars] = {
82  "trueQ2",
83  "recoQ2",
84  "trueW2",
85  "PtP",
86  "CosNumi"
87  };
88 
89  //which cut tiers do you want to look at?
90  const int kNumSels = 1;
91  const std::string selNamesND[kNumSels] = {
92  "numuND"
93  };
94 
95  const std::string selNamesFD[kNumSels] = {
96  "AllSamples"
97  };
98 
99  IPrediction* predsND[kNumSels][kNumVars];
100 
103  ResetOscCalcToDefault(&mycalc);
104 
105  mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
106  mycalc.SetdCP(3*TMath::Pi()/2);
107  mycalc.SetTh23(asin(sqrt(0.404)));
108 
109  //pick up spectra/prediction objects from the input root file
110  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
111  const char* selName = selNamesND[selIdx].c_str();
112  const char* wName = weight_name.c_str();
113 
114  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
115  const char* varName = defs[varIdx].name.c_str();
116 
117  std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
118 
119  predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s_%s", selName, wName, varName).Data()).release();
120 
121  if(selIdx == 0 && varIdx == 0)
122  std::cout << "ND MC POT: "
123  << predsND[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
124 
125  } // end for varIdx
126  } // end for selIdx
127 
128 
129  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
130  TDirectory * dpred_no[kNumSels][kNumVars];
131  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
132  const std::string selName = selNamesFD[selIdx];
133  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
134  const std::string varName = defs[varIdx].name;
135  if(varName == "recoQ2") {
136  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
137  } else if(varName == "CosNumi"){
138  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
139  } else {
140  dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
141  }
142  }
143  }
144 
145  TH1* hFullNDHists[kNumVars-1];
146  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
147  const std::string selName = selNamesFD[selIdx];
148  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
149 
150  const std::string varName = defs[varIdx].name;
151 
152  //play with the binning choice for given variables:
153  int rebinFactor=1;
154 
155  //create a canvas
156  TCanvas *c2 =new TCanvas("c2","c2",500,400);
157 
158 
159  TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
160  padAbs -> SetBottomMargin(0);
161  padAbs->SetFillStyle(0);
162  padAbs -> SetLeftMargin(0.15);
163 
164  padAbs -> Draw();
165  padAbs -> cd();
166 
167  double legendxlow=0.55;
168  double legendxhigh=0.95;
169 
170  if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.2; legendxhigh=0.5;}
171 
172  //Pull down the variable we want
173 
174  IPrediction* predND = predsND[selIdx][varIdx];
175  auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);
176 
177  //Get components of the prediction and store them in histograms
178 
179  TH1* hFD = predFD->PredictComponent(&mycalc,
181  Current::kCC,
182  Sign::kBoth).ToTH1(kPOT);
183 
184  TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);
185 
186  TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
187 
188  hND->SetLineColor(kTotalMCColor);
189 
190  hFD->SetLineColor(kNueSignalColor);
191  hND->Rebin(rebinFactor);
192 
193  std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
194 
195  hND->GetXaxis()->CenterTitle();
196  hND->GetYaxis()->CenterTitle();
197  hND->GetXaxis()->SetDecimals();
198  hND->GetYaxis()->SetDecimals();
199  hND->GetYaxis()->SetTitleOffset(1.15);
200  hND->GetXaxis()->SetLabelSize(.0);
201  hND->GetYaxis()->SetLabelSize(.05);
202  hND->SetMaximum(1.1*TMath::Max(hND->GetMaximum(), hFD->GetMaximum()));
203  hND->GetYaxis()->SetTitle(TString::Format("Events"));
204 
205  if(varName=="recoQ2" || varName=="trueQ2"){
206  hND->GetXaxis()->SetRangeUser(0,2);
207  }
208 
209  if(varName=="recoE"){
210  hND->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
211  }
212  hND->GetXaxis()->SetNdivisions(405,kFALSE);
213  hND->SetMinimum(.001);
214 
215  hND ->DrawCopy("hist");
216 
217  hND->Draw("hist same");
218  hFD->Draw("hist same");
219 
220  //Fix up the axis
221  padAbs->RedrawAxis();
222  padAbs->Update();
223 
224  //Add the legend and preliminary
225  Legend(legendxlow, .70, legendxhigh, .85);
226  std::string wTitle = weight_name+" reweighted";
227 
228  TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
229  selTitle->SetNDC();
230  selTitle->SetTextSize(2/30.);
231  selTitle->SetTextAlign(32);
232  selTitle->Draw();
233 
234  c2->cd();
235  TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
236  padRatio -> SetTopMargin(0);
237  padRatio -> SetBottomMargin(0.25);
238  padRatio -> SetFillStyle(0);
239  padRatio -> SetLeftMargin(0.15);
240 
241  padRatio -> Draw();
242  padRatio -> cd();
243 
244  TH1* hRatio = (TH1*)hFD -> Clone("hRatio");
245 
246  double b=5;
247  if(varName=="trueW2") {b = 2.0;}
248  if(varName=="trueQ2") {b = 2.0;}
249  if(varName=="recoQ2") {b = 2.0;}
250  if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
251 
252  TLine* lOne = new TLine(0,1.0,b,1.0);
253  lOne -> SetLineStyle(2);
254 
255  hRatio -> Divide(hND);
256 
257  hRatio->GetXaxis()->CenterTitle();
258  hRatio->GetYaxis()->CenterTitle();
259  hRatio->GetYaxis()->SetLabelSize(4.0/30.);
260  hRatio->GetXaxis()->SetLabelSize(4.0/30.);
261  hRatio->GetXaxis()->SetLabelOffset(.04);
262  hRatio->GetXaxis()->SetTickSize(.06);
263  hRatio->GetYaxis()->SetTitleSize(4.5/30.);
264  hRatio->GetXaxis()->SetTitleSize(4.5/30.);
265  hRatio->GetYaxis()->SetRangeUser(0.5,1.5);
266  hRatio->SetMarkerSize(.1);
267  hRatio->GetXaxis()->SetDecimals();
268  hRatio->GetYaxis()->SetDecimals();
269  hRatio->GetYaxis()->SetTitleOffset(0.4);
270  hRatio->GetYaxis()->SetTitle("Sig. FD/#nu_{#mu} ND");
271 
272  if(varName=="recoE"){
273  hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
274  }
275  if(varName=="recoQ2"){
276  hRatio->GetXaxis()->SetTitle("ND reco Q^{2} vs. FD true Q^{2} (GeV^{2})");
277  }
278 
279  if(varName=="recoQ2" || varName=="trueQ2"){
280  hRatio->GetXaxis()->SetRangeUser(0,2);
281  }
282 
283  hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
284  hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
285 
286  hRatio-> Draw("hist");
287  lOne -> Draw("same");
288 
289  c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".png").c_str());
290 
291  } // end for varIdx
292  } // end for selIdx
293 }
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.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
void plot_NDvsFD_REW(const std::string weight_name, const std::string fname_ND, const std::string fname_FD)
Float_t x1[n_points_granero]
Definition: compare.C:5
hmean SetLineStyle(2)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void ResetOscCalcToDefault(osc::IOscCalculatorAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
void SetTh23(const T &th23) override
const Color_t kTotalMCErrorBandColor
Definition: Style.h:18
ntuple SetFillStyle(1001)
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
void SetDmsq32(const T &dmsq32) override
const Color_t kNueSignalColor
Definition: Style.h:19
c2
Definition: demo5.py:33
void SetdCP(const T &dCP) override
Charged-current interactions.
Definition: IPrediction.h:39
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
void Legend(double x0, double y0, double x1, double y1)
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
virtual T GetDmsq32() const
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)
T asin(T number)
Definition: d0nt_math.hpp:60