NueSystFuncs.h
Go to the documentation of this file.
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Cut.h"
7 #include "CAFAna/Core/HistAxis.h"
9 #include "CAFAna/Core/Ratio.h"
10 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Cuts/Cuts.h"
14 #include "CAFAna/Cuts/NumuCuts.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
19 #include "CAFAna/Vars/NueVars.h"
20 #include "CAFAna/Vars/Vars.h"
21 
22 
24 
25 #include "TStyle.h"
26 #include "TFile.h"
27 #include "TH1.h"
28 #include "TCanvas.h"
29 #include "TGraph.h"
30 #include "TGraphAsymmErrors.h"
31 #include "TGaxis.h"
32 #include "TLatex.h"
33 #include "TLegend.h"
34 #include "TLine.h"
35 #include "TPad.h"
36 #include "TDirectory.h"
37 #include "TSystem.h"
38 #include <iostream>
39 #include <fstream>
40 #include <iomanip>
41 
42 using namespace ana;
43 
46 
47 
48 //....................................................................................
49 void ComparePredictions (PredictionExtrap* predPropMC, PredictionExtrap* predPropMCPlus, PredictionExtrap* predPropMCMinus, const std::string &predName, const std::string &systName, const std::string &varName, int varIdx)
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 }
121 
122 
123 //..........................................................................................
124 void PredRatioToNom (PredictionExtrap* predPropMC, PredictionExtrap* predPropMCPlus, PredictionExtrap* predPropMCMinus, const std::string &predName, const std::string &systName, const std::string &varName, int varIdx)
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 }
191 
192 
193 //..................................................................................
194 void CompareMCCompPrediction (PredictionExtrap* predPropMC, PredictionExtrap* predPropMCPlus, PredictionExtrap* predPropMCMinus)
195 
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 }
285 
286 
287 //..........................................................................
288 void MCCompPredictionTable(PredictionExtrap* predPropMC, PredictionExtrap* predPropMCPlus, PredictionExtrap* predPropMCMinus)
289 
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);
306  c.nc = predPropMC->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).Integral(kPOT);
307  c.beamnue = predPropMC->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).Integral(kPOT);
308  c.signal = predPropMC->PredictComponent(calc, Flavors::kNuMuToNuE, 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 }
void CompareMCCompPrediction(PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus)
Definition: NueSystFuncs.h:194
Pass neutrinos through unchanged.
Spectrum Predict(osc::IOscCalculator *calc) const override
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:122
(&#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
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetdCP(const T &dCP)=0
TLegend * leg1
Definition: plot_hist.C:105
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
void PredRatioToNom(PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus, const std::string &predName, const std::string &systName, const std::string &varName, int varIdx)
Definition: NueSystFuncs.h:124
c2
Definition: demo5.py:33
Charged-current interactions.
Definition: IPrediction.h:39
General interface to any calculator that lets you set the parameters.
void MCCompPredictionTable(PredictionExtrap *predPropMC, PredictionExtrap *predPropMCPlus, PredictionExtrap *predPropMCMinus)
Definition: NueSystFuncs.h:288
const double j
Definition: BetheBloch.cxx:29
osc::IOscCalculatorAdjustable * calc
Definition: NueSystFuncs.h:45
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
virtual void SetTh13(const T &th13)=0
Take the output of an extrapolation and oscillate it as required.
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
void ComparePredictions(IPrediction *pred_no, IPrediction *pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=false)
T asin(T number)
Definition: d0nt_math.hpp:60