demo5.C
Go to the documentation of this file.
1 // Use that prediction to make a contour for a really dumb numu analysis
2 
4 
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Fit/Fit.h"
8 #include "CAFAna/Vars/FitVars.h"
9 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Analysis/Plots.h"
13 #include "CAFAna/Core/Spectrum.h"
16 
18 
19 #include "TCanvas.h"
20 
21 
22 using namespace ana;
23 
24 
25 void demo5()
26 {
27  const std::string fname = "prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_numu2020";
28  const std::string fnameSwap = "prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1_numu2020";
29 
30  SpectrumLoader loader(fname);
31  SpectrumLoader loaderSwap(fnameSwap);
32 
33  const Binning bins = Binning::Simple(100, 0, 1000);
34 
35  // All interaction types
36  PredictionNoExtrap pred(loader, loaderSwap, "Number of hits in slice", bins, kNHit, kNoCut);
37  // Another way of defining a PredNoExtrap. Note there are different `ProdXLoaders` depending on Ana year.
38  /*
39  Header: #include "CAFAna/Analysis/Prod4Loaders.h"
40 
41  Prod4NomLoaders loader( <ECAFType::>, <Loaders::kFHC/Loaders::kRHC>, <std::string Period> );
42 
43  PredictionNoExtrap pred( loader, "Number of hits per slice", bins, kNHit, kNoCut);
44  */
45  // Finding Vars:
46  // All CAFAna Vars should be in `CAFAna/Vars`. The most general should be in `CAFAna/Vars/Vars.h`.
47  // A quick search using grep/find/<etc> in that directory should hopefully yield the Var you want.
48 
49 
50  // Do it!
51  loader.Go();
52  loaderSwap.Go();
53 
54  // How to scale histograms
55  const double pot = 18e20;
56 
57  // Make a calculator. This is the fastest variant
59 
60  calc.SetL(810);
61  calc.SetRho(2.75);
62  calc.SetDmsq21(7.6e-5);
63  calc.SetDmsq32(2.35e-3);
64  calc.SetTh12(asin(sqrt(.87))/2);
65  calc.SetTh13(asin(sqrt(.10))/2);
66  calc.SetTh23(TMath::Pi()/4);
67  calc.SetdCP(0);
68 
69  new TCanvas;
70  // What does it predict at these oscillation parameters?
71  Spectrum spred = pred.Predict(&calc);
72  // Add some poisson fluctuations, for fun
73  Spectrum mock = spred.MockData(pot);
74 
75  // See all the components
76  DataMCComparisonComponents(mock, &pred, &calc);
77  // http://nusoft.fnal.gov/nova/novasoft/doxygen/html/namespaceana.html#a19a474dbc249aab1b2c3fb5ae723c22e
78  // Personally, I (Karl W), tend to expand the contents of that function into my macro, but it's a useful function.
79  // Also note that you can specify different;
80  // Flavours - http://nusoft.fnal.gov/nova/novasoft/doxygen/html/namespaceana_1_1Flavors.html#aa60f64652048d17d131d9da6e38df26eaf6d7acd93b8a01cbaf2d7a86a7858124
81  // Currents - http://nusoft.fnal.gov/nova/novasoft/doxygen/html/namespaceana_1_1Current.html#a0284685bc4ae90572e7c9ac760ec6873af9f37ae1299b9b253572ac1de34dc4ad
82  // Signs - http://nusoft.fnal.gov/nova/novasoft/doxygen/html/namespaceana.html#a46dc35ff760761f9c65207f8663d3487a022f571a6e5da9f570c50e5df4ef9849
83 
84  // What did it look like unoscillated?
85  pred.PredictUnoscillated().ToTH1(pot)->Draw("hist same");
86 
87  // Not too complicated to make a contour too
88  new TCanvas;
89  SingleSampleExperiment expt(&pred, mock);
90  FrequentistSurface surf(&expt, &calc,
91  &kFitSinSqTheta23, 30, .35, .65,
92  &kFitDmSq32Scaled, 30, 2.2, 2.8);
93  surf.Draw();
94  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
95  surf.DrawBestFit(kRed);
96 }
enum BeamMode kRed
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh13(const T &th13) override
void SetL(double L) override
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:148
string fnameSwap
Definition: demo4.py:6
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
virtual Spectrum PredictUnoscillated() const
Definition: IPrediction.cxx:33
spred
Definition: demo4.py:26
osc::OscCalcDumb calc
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Log-likelihood scan across two parameters.
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
expt
Definition: demo5.py:34
void demo5()
Definition: demo5.C:25
Spectrum Predict(osc::IOscCalc *calc) const override
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
const Var kNHit
Definition: Vars.cxx:71
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
#define pot
virtual void Go() override
Load all the registered spectra.
void SetTh23(const T &th23) override
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
loader
Definition: demo0.py:10
void SetDmsq21(const T &dmsq21) override
void SetDmsq32(const T &dmsq32) override
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH1 * DataMCComparisonComponents(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc)
Plot MC broken down into flavour components, overlayed with data.
Definition: Plots.cxx:114
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SetdCP(const T &dCP) override
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
void SetTh12(const T &th12) override
Prediction that just uses FD MC, with no extrapolation.
mock
Definition: demo4.py:28
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
void SetRho(double rho) override
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
enum BeamMode string