Functions | Variables
NueSystFuncs.h File Reference
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"
#include "OscLib/IOscCalc.h"
#include "TStyle.h"
#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TPad.h"
#include "TDirectory.h"
#include "TSystem.h"
#include <iostream>
#include <fstream>
#include <iomanip>

Go to the source code of this file.

Functions

void ComparePredictions (PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus, const std::string &predName, const std::string &systName, const std::string &varName, int varIdx)
 
void PredRatioToNom (PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus, const std::string &predName, const std::string &systName, const std::string &varName, int varIdx)
 
void CompareMCCompPrediction (PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus)
 
void MCCompPredictionTable (PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus)
 

Variables

osc::NoOscillations noosc
 
osc::IOscCalcAdjustablecalc = ana::DefaultOscCalc()
 

Function Documentation

void CompareMCCompPrediction ( PredictionExtrap predPropMC,
PredictionExtrap predPropMCPlus,
PredictionExtrap predPropMCMinus 
)

Definition at line 194 of file NueSystFuncs.h.

References demo5::c1, om::cout, allTimeWatchdog::endl, MECModelEnuComparisons::i, calib::j, ana::Flavors::kAll, ana::Flavors::kAllNuMu, ana::Sign::kBoth, ana::Current::kCC, ana::Current::kNC, ana::Flavors::kNuEToNuE, ana::Flavors::kNuMuToNuE, ana::kPOT, leg1, ana::PredictionExtrap::PredictComponent(), and ana::Spectrum::ToTH1().

196 {
197  const double kPOT= 9.0000e20;//since there's no data to compare it doesn't really matter what pot you normalize to
198 
199  TH1* hNueSigMC = predPropMC->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(kPOT);
200  TH1* hBeamNueMC = predPropMC->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(kPOT);
201  TH1* hNCMC = predPropMC->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kPOT);
202  TH1* hNumuccMC = predPropMC->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kPOT);
203 
204 
205 
206  TH1* hNueSigMCPlus = predPropMCPlus->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(kPOT);
207  TH1* hBeamNueMCPlus = predPropMCPlus->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(kPOT);
208  TH1* hNCMCPlus = predPropMCPlus->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kPOT);
209  TH1* hNumuccMCPlus = predPropMCPlus->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kPOT);
210 
211 
212 
213  TH1* hNueSigMCMinus = predPropMCMinus->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).ToTH1(kPOT);
214  TH1* hBeamNueMCMinus = predPropMCMinus->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(kPOT);
215  TH1* hNCMCMinus = predPropMCMinus->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kPOT);
216  TH1* hNumuccMCMinus = predPropMCMinus->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kPOT);
217 
218  std::cout<<"Making plots"<<std::endl;
219 
220  hNueSigMC->SetLineColor(kBlack);
221  hBeamNueMC->SetLineColor(kBlack);
222  hNCMC->SetLineColor(kBlack);
223  hNumuccMC->SetLineColor(kBlack);
224 
225  hNueSigMCPlus->SetLineColor(kRed);
226  hBeamNueMCPlus->SetLineColor(kRed);
227  hNCMCPlus->SetLineColor(kRed);
228  hNumuccMCPlus->SetLineColor(kRed);
229 
230  hNueSigMCMinus->SetLineColor(kBlue);
231  hBeamNueMCMinus->SetLineColor(kBlue);
232  hNCMCMinus->SetLineColor(kBlue);
233  hNumuccMCMinus->SetLineColor(kBlue);
234 
235  std::vector <TH1*> vecMC = {hNueSigMC, hBeamNueMC, hNCMC, hNumuccMC};
236  std::vector <TH1*> vecMCPlus = {hNueSigMCPlus, hBeamNueMCPlus, hNCMCPlus, hNumuccMCPlus};
237  std::vector <TH1*> vecMCMinus = {hNueSigMCMinus,hBeamNueMCMinus,hNCMCMinus,hNumuccMCMinus};
238 
239 
240  const int kNumComps = 4;
241  const std::string mccompNames[kNumComps] = {"NueSig", "BeamNue", "NC", "NumuCC"};
242 
243  for (unsigned int i=0; i < vecMC.size(); ++i)
244  {
245  for (unsigned int j=0; j < vecMCPlus.size(); ++j)
246  {
247  for (unsigned int k=0; k < vecMCMinus.size(); ++k)
248  {
249  if (i == j && j == k)
250  {
251 
252  TCanvas *c1 = new TCanvas(mccompNames[i].c_str(), mccompNames[i].c_str());
253  TLegend* leg1 = new TLegend(0.50,0.65,0.80,0.85);
254 
255  if(i == 2)
256  (vecMCMinus[k]->GetYaxis()->SetRangeUser(0, 0.7));
257 
258  vecMCMinus[k]->GetYaxis()->SetTitle("Events / 9x10^{20} POT");
259  vecMCMinus[k]->GetYaxis()->CenterTitle();
260  vecMCMinus[k]->GetYaxis()->SetTitleSize(0.045);
261  vecMCMinus[k]->GetYaxis()->SetTitleOffset(0.99);
262  vecMCMinus[k]->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV.)");
263  vecMCMinus[k]->GetXaxis()->CenterTitle();
264  vecMCMinus[k]->GetXaxis()->SetTitleSize(0.050);
265  vecMCMinus[k]->GetXaxis()->SetTitleOffset(0.85);
266 
267  vecMCMinus[k]->Draw("hist");
268  vecMC[i]->Draw("hist same");
269  vecMCPlus[j]->Draw("hist same");
270 
271  leg1->AddEntry(vecMC[i], (mccompNames[i]+" for MC(pred.) Nominal").c_str(),"l");
272  leg1->AddEntry(vecMCPlus[j], (mccompNames[j]+" for MC(pred.) +ve Shift").c_str(),"l");
273  leg1->AddEntry(vecMCMinus[k],(mccompNames[k]+" for MC(pred.) -ve Shift").c_str(),"l");
274  leg1->SetLineColor(10);
275  leg1->SetFillStyle(0);
276  leg1->SetTextSize(0.035);
277  leg1->Draw();
278 
279  }
280  }
281  }
282  }
283 
284 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
(&#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:209
(&#39;beam &#39;)
Definition: IPrediction.h:15
TLegend * leg1
Definition: plot_hist.C:105
Charged-current interactions.
Definition: IPrediction.h:39
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
osc::IOscCalcAdjustable * calc
Definition: NueSystFuncs.h:45
c1
Definition: demo5.py:24
All neutrinos, any flavor.
Definition: IPrediction.h:26
void ComparePredictions ( PredictionExtrap predPropMC,
PredictionExtrap predPropMCPlus,
PredictionExtrap predPropMCMinus,
const std::string &  predName,
const std::string &  systName,
const std::string &  varName,
int  varIdx 
)

Definition at line 49 of file NueSystFuncs.h.

References std::asin(), demo5::c1, e, ana::kPOT, leg1, ana::PredictionExtrap::Predict(), osc::_IOscCalcAdjustable< T >::SetdCP(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh13(), osc::_IOscCalcAdjustable< T >::SetTh23(), std::sqrt(), and ana::Spectrum::ToTH1().

50 {
51  const double kPOT= 9.0000e20;//since there's no data to compare it doesn't really matter what pot you normalize to
52 
53  struct PlotDef
54  {
55  std::string xlabel;
56  std::string ylabel;
57  };
58 
59  std::vector <PlotDef> plotdef;
60 
61  plotdef.push_back ({"Shower Calorimetric Energy (GeV)", "Events / 9x10^{20} POT"});
62  plotdef.push_back ({"Hadronic Calorimetric Energy (GeV)", "Events / 9x10^{20} POT"});
63  plotdef.push_back ({"Reconstructed Neutrino Energy (GeV)", "Events / 9x10^{20} POT"});
64  plotdef.push_back ({"CVN #nu_{e}", "Events / 9x10^{20} POT"});
65  plotdef.push_back ({"Nue Energy / CVN", "Events / 9x10^{20} POT"});
66 
67 
68  calc->SetDmsq32(2.44e-3);
69  calc->SetTh13(asin(sqrt(.085))/2);
70  calc->SetTh23(asin(sqrt(0.45)));
71  calc->SetdCP(1.59*TMath::Pi());
72 
73 
74 
75  TCanvas *c1 = new TCanvas((predName+"_"+systName+"_"+varName).c_str(), (predName+"_"+systName+"_"+varName).c_str());
76  TLegend* leg1 = new TLegend(0.52,0.62,0.82,0.85);
77 
78 
79  TH1* hPredMCMinus = predPropMCMinus->Predict(calc).ToTH1(kPOT);
80  hPredMCMinus->SetMarkerStyle(2);
81  hPredMCMinus->SetMarkerSize(0.5);
82  hPredMCMinus->SetMarkerColor(kBlue);
83  hPredMCMinus->SetLineColor(kBlue);
84  hPredMCMinus->GetXaxis()->SetTitle((plotdef[varIdx].xlabel).c_str());
85  hPredMCMinus->GetXaxis()->CenterTitle();
86  hPredMCMinus->GetXaxis()->SetTitleSize(0.045);
87  hPredMCMinus->GetXaxis()->SetTitleOffset(0.95);
88  hPredMCMinus->GetYaxis()->SetTitle((plotdef[varIdx].ylabel).c_str());
89  hPredMCMinus->GetYaxis()->CenterTitle();
90  hPredMCMinus->GetYaxis()->SetTitleSize(0.045);
91  hPredMCMinus->GetYaxis()->SetTitleOffset(0.85);
92  hPredMCMinus->Draw("hist");
93 
94 
95  TH1* hPredMC = predPropMC->Predict(calc).ToTH1(kPOT);
96  hPredMC->SetMarkerStyle(2);
97  hPredMC->SetMarkerSize(0.5);
98  hPredMC->SetMarkerColor(kBlack);
99  hPredMC->SetLineColor(kBlack);
100  hPredMC->Draw("hist same");
101 
102 
103  TH1* hPredMCPlus = predPropMCPlus->Predict(calc).ToTH1(kPOT);
104  hPredMCPlus->SetMarkerStyle(2);
105  hPredMCPlus->SetMarkerSize(0.5);
106  hPredMCPlus->SetMarkerColor(kRed);
107  hPredMCPlus->SetLineColor(kRed);
108  hPredMCPlus->Draw("hist same");
109 
110 
111  leg1->AddEntry(hPredMC,"MC(pred.) Nominal","l");
112  leg1->AddEntry(hPredMCPlus,"MC(pred.) +ve Calib Shift","l");
113  leg1->AddEntry(hPredMCMinus,"MC(pred.) -ve Calib Shift","l");
114 
115  leg1->SetLineColor(10);
116  leg1->SetFillStyle(0);
117  leg1->SetTextSize(0.035);
118  leg1->Draw();
119 
120 }
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetTh13(const T &th13)=0
TLegend * leg1
Definition: plot_hist.C:105
virtual void SetDmsq32(const T &dmsq32)=0
Spectrum Predict(osc::IOscCalc *calc) const override
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual void SetTh23(const T &th23)=0
osc::IOscCalcAdjustable * calc
Definition: NueSystFuncs.h:45
c1
Definition: demo5.py:24
Float_t e
Definition: plot.C:35
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
void MCCompPredictionTable ( PredictionExtrap predPropMC,
PredictionExtrap predPropMCPlus,
PredictionExtrap predPropMCMinus 
)

Definition at line 288 of file NueSystFuncs.h.

References plot_validation_datamc::c, om::cout, allTimeWatchdog::endl, ana::Spectrum::Integral(), ana::Flavors::kAll, ana::Flavors::kAllNuMu, ana::Sign::kBoth, ana::Current::kCC, ana::Current::kNC, ana::Flavors::kNuEToNuE, ana::Flavors::kNuMuToNuE, ana::kPOT, make_root_from_grid_output::nc, and ana::PredictionExtrap::PredictComponent().

290 {
291  struct current
292  {
293  double numucc;
294  double nc;
295  double beamnue;
296  double signal;
297  };
298 
299  std::vector <current> noofevt;
300 
301  const double kPOT= 9.0000e20;
302 
303  current c;
304 
305  c.numucc = predPropMC->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).Integral(kPOT);
307  c.beamnue = predPropMC->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).Integral(kPOT);
309 
310  noofevt.push_back(c);
311 
312  std::cout << "MC.numucc: " << c.numucc << std::endl;
313  std::cout << "MC.nc: " << c.nc << std::endl;
314  std::cout << "MC.beamnue: " << c.beamnue << std::endl;
315  std::cout << "MC.signal: " << c.signal << std::endl;
316 
317  c.numucc = predPropMCPlus->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).Integral(kPOT);
318  c.nc = predPropMCPlus->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).Integral(kPOT);
319  c.beamnue = predPropMCPlus->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).Integral(kPOT);
320  c.signal = predPropMCPlus->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).Integral(kPOT);
321 
322  noofevt.push_back(c);
323 
324  std::cout << "MCPlus.numucc: " << c.numucc << std::endl;
325  std::cout << "MCPlus.nc: " << c.nc << std::endl;
326  std::cout << "MCPlus.beamnue: " << c.beamnue << std::endl;
327  std::cout << "MCPlus.signal: " << c.signal << std::endl;
328 
329  c.numucc = predPropMCMinus->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).Integral(kPOT);
330  c.nc = predPropMCMinus->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).Integral(kPOT);
331  c.beamnue = predPropMCMinus->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).Integral(kPOT);
332  c.signal = predPropMCMinus->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth).Integral(kPOT);
333 
334  noofevt.push_back(c);
335 
336  std::cout << "MCMinus.numucc: " << c.numucc << std::endl;
337  std::cout << "MCMinus.nc: " << c.nc << std::endl;
338  std::cout << "MCMinus.beamnue: " << c.beamnue << std::endl;
339  std::cout << "MCMinus.signal: " << c.signal << std::endl;
340 
341 
342 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
(&#39; appearance&#39;)
Definition: IPrediction.h:18
(&#39;beam &#39;)
Definition: IPrediction.h:15
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:333
Charged-current interactions.
Definition: IPrediction.h:39
OStream cout
Definition: OStream.cxx:6
osc::IOscCalcAdjustable * calc
Definition: NueSystFuncs.h:45
All neutrinos, any flavor.
Definition: IPrediction.h:26
void PredRatioToNom ( PredictionExtrap predPropMC,
PredictionExtrap predPropMCPlus,
PredictionExtrap predPropMCMinus,
const std::string &  predName,
const std::string &  systName,
const std::string &  varName,
int  varIdx 
)

Definition at line 124 of file NueSystFuncs.h.

References std::asin(), demo5::c2, e, make_syst_table_plots::line, ana::PredictionExtrap::Predict(), osc::_IOscCalcAdjustable< T >::SetdCP(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh13(), osc::_IOscCalcAdjustable< T >::SetTh23(), std::sqrt(), and ana::Ratio::ToTH1().

Referenced by plot_nue_filesyst_pred(), and plot_nue_xsec_pred().

125 {
126  struct PlotDef
127  {
128  std::string xlabel;
129  std::string ylabel;
130  };
131 
132  std::vector <PlotDef> plotdef;
133 
134  plotdef.push_back ({"Shower Calorimetric Energy (GeV)", "Ratio"});
135  plotdef.push_back ({"Hadronic Calorimetric Energy (GeV)", "Ratio"});
136  plotdef.push_back ({"Reconstructed Neutrino Energy (GeV)", "Ratio"});
137  plotdef.push_back ({"CVN #nu_{e}", "Ratio"});
138  plotdef.push_back ({"Nue Energy / CVN", "Ratio"});
139 
140  calc->SetDmsq32(2.44e-3);
141  calc->SetTh13(asin(sqrt(.085))/2);
142  calc->SetTh23(asin(sqrt(0.45)));
143  calc->SetdCP(1.59*TMath::Pi());
144 
145 
146  Ratio* ratioPredPlus = new Ratio (predPropMCPlus->Predict(calc), predPropMC->Predict(calc));
147  Ratio* ratioPredMinus = new Ratio (predPropMCMinus->Predict(calc), predPropMC->Predict(calc));
148 
149  TH1* hRatioPredPlus = ratioPredPlus->ToTH1();
150  TH1* hRatioPredMinus = ratioPredMinus->ToTH1();
151 
152 
153  TCanvas *c2 = new TCanvas((predName+"_"+systName+"RATIO"+varName).c_str(), (predName+"_"+systName+"RATIO"+varName).c_str());
154  TLegend* leg2 = new TLegend(0.45,0.70,0.70,0.85);
155 
156 
157  hRatioPredMinus->GetYaxis()->SetRangeUser(0.5, 1.5);
158  hRatioPredMinus->SetLineColor(kBlue);
159  hRatioPredMinus->GetXaxis()->SetTitle((plotdef[varIdx].xlabel).c_str());
160  hRatioPredMinus->GetXaxis()->CenterTitle();
161  hRatioPredMinus->GetXaxis()->SetTitleSize(0.045);
162  hRatioPredMinus->GetXaxis()->SetTitleOffset(0.95);
163  hRatioPredMinus->GetYaxis()->SetTitle((plotdef[varIdx].ylabel).c_str());
164  hRatioPredMinus->GetYaxis()->CenterTitle();
165  hRatioPredMinus->GetYaxis()->SetTitleSize(0.050);
166  hRatioPredMinus->GetYaxis()->SetTitleOffset(0.85);
167 
168  hRatioPredPlus->SetLineColor(kRed);
169 
170 
171  hRatioPredMinus->Draw();
172  hRatioPredPlus->Draw("same");
173 
174 
175  gPad->Update();
176  TLine* line = new TLine(gPad->GetUxmin(),1,gPad->GetUxmax(),1);
177  line->SetLineColor(kBlack);
178  line->SetLineStyle(2);
179  line->SetLineWidth(2);
180  line->Draw();
181  line->Draw();
182 
183  leg2->AddEntry(hRatioPredMinus,"Ratio (MC -ve Calib Shift / MC Nom.)","lpe");
184  leg2->AddEntry(hRatioPredPlus, "Ratio (MC +ve Calib Shift / MC Nom.)","lpe");
185  leg2->SetLineColor(10);
186  leg2->SetFillStyle(0);
187  leg2->SetTextSize(0.035);
188  leg2->Draw();
189 
190 }
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:68
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetTh13(const T &th13)=0
virtual void SetDmsq32(const T &dmsq32)=0
c2
Definition: demo5.py:33
Spectrum Predict(osc::IOscCalc *calc) const override
Represent the ratio between two spectra.
Definition: Ratio.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual void SetTh23(const T &th23)=0
osc::IOscCalcAdjustable * calc
Definition: NueSystFuncs.h:45
Float_t e
Definition: plot.C:35
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60

Variable Documentation

Definition at line 45 of file NueSystFuncs.h.