test_ana.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Vars/FitVars.h"
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/OscCurve.h"
13 #include "CAFAna/Core/Utilities.h"
14 
15 #include "TCanvas.h"
16 #include "TFile.h"
17 #include "TGraph.h"
18 #include "TH1.h"
19 #include "TMath.h"
20 
21 #include <cmath>
22 #include <iostream>
23 
24 #include "Utilities/func/MathUtil.h"
25 
26 using namespace ana;
27 
28 void cc()
29 {
30  Loaders loaders("$NOVA_DATA/mc/S13-06-18/hadd/", Loaders::kFHC);
31 
32  // Define our binning
33  Binning bins = Binning::Simple(40, 0, 4);
34 
35  // Load up the necessary histograms, seperated by flavour etc
36  PredictionNoExtrap pred(loaders,
37  "Calorimetric energy", bins,
39 
40  loaders.Go();
41 
42  // Make a calculator. Any of the variants is fine
44 
45  calc.SetL(810);
46  calc.SetRho(0); // No matter effects
47  calc.SetDmsq21(7.6e-5);
48  calc.SetDmsq32(2.35e-3);
49  calc.SetTh12(asin(sqrt(.87))/2); // PDG
50  calc.SetTh13(asin(sqrt(.10))/2);
51  calc.SetTh23(TMath::Pi()/4);
52  calc.SetdCP(0);
53 
54  // What we expect to observe at those oscillation parameters
55  Spectrum obs = pred.Predict(&calc);
56 
57  // Make it into mock data
58  // Spectrum mock = obs.MockData(18e20);
59  Spectrum mock = obs.FakeData(18e20);
60 
61  Spectrum unosc = pred.PredictUnoscillated();
62 
63  new TCanvas;
64  // DataMCComparison(mock, obs);
65  DataMCComparisonComponents(mock, &pred, &calc);
66  TH1* uoh = unosc.ToTH1(18e20);
67  uoh->SetLineColor(kBlue);
68  uoh->Draw("same");
69 
70  new TCanvas;
71  RatioPlot(mock, unosc, obs);
72 
73  SingleSampleExperiment expt(&pred, mock);
74 
75 
76  calc.SetTh23(1); // Seed value
77  MinuitFitter fit(&expt, {&kFitSinSq2Theta23, &kFitDmSq32});
78  const double chi = fit.Fit(&calc)->EvalMetricVal();
79  std::cout << "Best fit, dmsq = " << calc.GetDmsq32() << " theta = " << calc.GetTh23() << " (ss2th23 = " << util::sqr(sin(2*calc.GetTh23())) << ") LL = " << chi << std::endl;
80 
81  // The log-likelihood surface
82  new TCanvas;
83  FrequentistSurface surf(&expt, &calc,
84  &kFitSinSq2Theta23, 20, .9, 1,
85  &kFitDmSq32, 20, 2e-3, 3e-3);
86  surf.Draw();
87  surf.DrawBestFit(kRed);
88 
89  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
90  surf.DrawContour(Gaussian90Percent2D(surf), kSolid, kRed);
91 
92  new TCanvas;
93  surf.DrawBestFit(kRed);
94 
95  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
96  surf.DrawContour(Gaussian90Percent2D(surf), kSolid, kRed);
97 }
98 
99 void test_ana()
100 {
101  // cc();
102  // return;
103 
104  Loaders loaders("$NOVA_DATA/mc/S13-06-18/hadd/", Loaders::kFHC);
106 
107  // Make a calculator. Any of the variants is fine
109 
110  calc.SetL(810);
111  calc.SetRho(2.75);
112  calc.SetDmsq21(7.6e-5);
113  calc.SetDmsq32(2.35e-3);
114  calc.SetTh12(asin(sqrt(.87))/2); // PDG
115  calc.SetTh13(asin(sqrt(.10))/2);
116  calc.SetTh23(TMath::Pi()/4);
117  calc.SetdCP(0);
118 
119  // Make sure we know what we're doing for sure
120  kFitSinSq2Theta13.SetValue(&calc, 0.1);
121  kFitSinSq2Theta23.SetValue(&calc, 1);
122 
123  // We can make oscillation curves
124  OscCurve curve(&calc, 14, 12);
125  new TCanvas;
126  curve.ToTH1()->Draw();
127 
128  OscCurve curve2(&calc, 14, 14);
129  new TCanvas;
130  curve2.ToTH1()->Draw();
131 
132 
133 
134  Binning bins2 = Binning::Simple(30, -.25, 1.25);
135 
136  PredictionNoExtrap pred(loaders,
137  "LEM PID", bins2,
138  // "LEM PID", 14, -.2, 1.2,
140  // kJmID, kPassesPresel);
141 
142  loaders.Go();
143 
144  // What we expect to observe at those oscillation parameters
145  Spectrum obs = pred.Predict(&calc);
146 
147  // Make it into mock data
148  // Spectrum mock = obs.MockData(18e20);
149  Spectrum mock = obs.FakeData(18e20);
150 
151  new TCanvas;
152  // DataMCComparison(mock, obs);
153  DataMCComparisonComponents(mock, &pred, &calc);
154 
155  // Overlay prediction at no oscillations
156  TH1* hExp = pred.PredictUnoscillated().ToTH1(18e20);
157  hExp->SetLineColor(kRed);
158  hExp->Draw("same");
159  hExp->SetLineStyle(7);
160  calc.SetTh13(0);
161  // And at no theta13
162  TH1* hExp2 = pred.Predict(&calc).ToTH1(18e20);
163  hExp2->SetLineColor(kGreen+2);
164  hExp2->Draw("same");
165 
166  // std::cout << "LL to no osc = " << LogLikelihood(pred.PredictUnoscillated(), mock) << std::endl;
167  // std::cout << "LL to input = " << LogLikelihood(obs, mock) << std::endl;
168  std::cout << "FOM = " << SimpleFOM(mock, pred.Predict(&calc)) << std::endl;
169 
170 
171  SingleSampleExperiment expt(&pred, mock);
172 
174  const double chi = fit.Fit(&calc)->EvalMetricVal();
175  std::cout << "Best fit, th13 = " << calc.GetTh13() << " (ss2th13 = " << util::sqr(sin(2*calc.GetTh13())) << ") delta = " << calc.GetdCP() << " LL = " << chi << std::endl;
176 
177  // The log-likelihood surface
178  new TCanvas;
179  FrequentistSurface surf(&expt, &calc,
180  &kFitSinSq2Theta13, 50, 0, .4,
181  &kFitDeltaInPiUnits, 50, 0, 2);
182  surf.Draw();
183  surf.DrawBestFit(kRed);
184 
185  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
186  surf.DrawContour(Gaussian90Percent2D(surf), 1, kRed); // solid
187 
188  // Inverted hierarchy
189  calc.SetDmsq32(-2.35e-3);
190 
191  FrequentistSurface surfInv(&expt, &calc,
192  &kFitSinSq2Theta13, 50, 0, .4,
193  &kFitDeltaInPiUnits, 50, 0, 2);
194  surfInv.DrawBestFit(kBlue);
195 
196  surfInv.DrawContour(Gaussian68Percent2D(surf), 7, kBlue); // dashed
197  surfInv.DrawContour(Gaussian90Percent2D(surf), 1, kBlue); // solid
198 }
const Var kLEM
PID
Definition: Vars.cxx:26
Far Detector at Ash River.
Definition: SREnums.h:11
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
const Cut kCrudeCCSel([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks==0) return false;return sr->trk.kalman.tracks[0].len > 200;})
Extremely crude CC selection, only for use in some old demos.
Definition: Cuts.h:17
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:209
void test_ana()
Definition: test_ana.C:99
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
const FitSinSq2Theta23 kFitSinSq2Theta23
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
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.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:165
virtual Spectrum PredictUnoscillated() const
Definition: IPrediction.cxx:56
void DisableLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Definition: Loaders.cxx:65
virtual T GetTh23() const
Definition: IOscCalc.h:89
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:48
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
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
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:64
TH1D * ToTH1(bool title=false) const
Definition: OscCurve.cxx:59
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Cut kPassesPresel(PassesPreselFunc)
Does this event pass the nue preselection?
Definition: Cuts.h:14
void SetTh23(const T &th23) override
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
void SetDmsq21(const T &dmsq21) override
virtual T GetdCP() const
Definition: IOscCalc.h:90
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
virtual T GetTh13() const
Definition: IOscCalc.h:88
const Binning bins
Definition: NumuCC_CPiBin.h:8
Transition probability for any one channel as a function of energy.
Definition: OscCurve.h:11
T sin(T number)
Definition: d0nt_math.hpp:132
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
void cc()
Definition: test_ana.C:28
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
void SetdCP(const T &dCP) override
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void SetTh12(const T &th12) override
double SimpleFOM(const Spectrum &obs, const Spectrum &unosc, double pot)
Figure-of-merit with no systematics, for binned data.
Definition: Fit.cxx:23
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Prediction that just uses FD MC, with no extrapolation.
mock
Definition: demo4.py:28
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
TCanvas * RatioPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption, std::vector< TH1 * > bottomHistos, std::vector< TString > bottomOption, TString pidtag, bool pidaxis=false)
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
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17