syst_plot_test.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
3 #include "CAFAna/Core/Loaders.h"
9 #include "CAFAna/Systs/Systs.h"
10 #include "CAFAna/Vars/Vars.h"
11 
12 #include "OscLib/OscCalcPMNSOpt.h"
13 
14 #include "TCanvas.h"
15 #include "TFile.h"
16 #include "TH2.h"
18 
19 using namespace ana;
20 
22 {
23  calc->SetL(810);
24  calc->SetRho(2.75);
25  calc->SetDmsq21(7.59e-5);
26  calc->SetDmsq32(2.40e-3);
27  calc->SetTh12(asin(sqrt(.87))/2); // PDG
28  calc->SetTh13(asin(sqrt(.095))/2);
29  calc->SetTh23(asin(sqrt(.99))/2);
30  calc->SetdCP(0);
31 }
32 
33 void GenFile()
34 {
35  // This is the oscillation parameters that PredictionInterp will expand
36  // around
37  osc::OscCalcPMNSOpt calcOrigin;
38  resetCalc(&calcOrigin);
39 
40  Loaders loaders("$NOVA_DATA/mc/S13-06-18/hadd/", Loaders::kFHC);
41  // This is just a test macro. These would require cosmic rejection cuts.
43 
44  // Pass the list of systematics we want to be able to handle
45  NoExtrapGenerator predGen( HistAxis("Simple Energy (GeV)",
46  Binning::Simple(20, 0, 5),
47  kCCE),
49  PredictionInterp interp({&kEnergyScaleSyst, &kNormSyst, &kNCScaleSyst},
50  &calcOrigin,
51  predGen,
52  loaders);
53 
54  loaders.Go();
55 
56  TFile fout("interp.root", "RECREATE");
57  interp.SaveTo(&fout);
58 }
59 
61 {
62  // Can comment this call after the first time, will read prediction out of
63  // the file.
64  GenFile();
65 
66 
68  resetCalc(&calc);
69 
70  // This is the interpolator that does most of the work
71  PredictionInterp* interp = LoadFrom<PredictionInterp>(new TFile("interp.root")).release();
72 
73  // Include NC background just to make the plots look cooler
74  TH1* nc = interp->PredictComponent(&calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(18e20);
75  nc->SetLineColor(kBlue);
76  nc->SetFillColor(kBlue-9);
77 
78  // Likewise some mock data
79  TH1* data = interp->Predict(&calc).MockData(18e20).ToTH1(18e20);
80  data->SetMarkerStyle(kFullCircle);
81 
82  // CC energy scale syst
83  PlotWithSystErrorBand(interp, {&kEnergyScaleSyst}, &calc, 18e20);
84  nc->Draw("hist ][ same");
85  data->Draw("ep same");
86 
87  // Normalization syst
88  new TCanvas;
89  PlotWithSystErrorBand(interp, {&kNormSyst}, &calc, 18e20);
90  nc->Draw("hist ][ same");
91  data->Draw("ep same");
92 
93  // NC background scale syst
94  new TCanvas;
95  PlotWithSystErrorBand(interp, {&kNCScaleSyst}, &calc, 18e20);
96  nc->Draw("hist ][ same");
97  data->Draw("ep same");
98 
99  // All systematics
100  new TCanvas;
101  PlotWithSystErrorBand(interp, {&kEnergyScaleSyst, &kNormSyst, &kNCScaleSyst}, &calc, 18e20);
102  nc->Draw("hist ][ same");
103  data->Draw("ep same");
104 }
virtual void SetL(double L)=0
Implements systematic errors by interpolation between shifted templates.
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void resetCalc(osc::IOscCalcAdjustable *calc)
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:149
virtual void SetDmsq21(const T &dmsq21)=0
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype, double alpha)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Cut kNumuCC
Definition: NumuCuts.h:70
virtual void SetTh13(const T &th13)=0
void DisableLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Definition: Loaders.cxx:65
void syst_plot_test()
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
osc::OscCalcDumb calc
Generates FD-only predictions (no extrapolation)
void GenFile()
virtual void SetDmsq32(const T &dmsq32)=0
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:301
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const XML_Char const XML_Char * data
Definition: expat.h:268
const NormSyst kNormSyst
Definition: Systs.cxx:14
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
const Cut kNumuQESel([](const caf::SRProxy *sr){std::cout<< "WARNING! Attempting to access kNumuQESel which relies on qepid (which no longer exists.) Aborting..."<< std::endl;abort(); return false;})
Definition: NumuCuts.h:38
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:32
const Var kCCE
Definition: NumuVars.h:21
Spectrum Predict(osc::IOscCalc *calc) const override
virtual void SetRho(double rho)=0
std::vector< Loaders * > loaders
Definition: syst_header.h:386
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
enum BeamMode nc
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
enum BeamMode kBlue
const NCScaleSyst kNCScaleSyst
Definition: Systs.cxx:16
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60