Functions
test_ana.C File Reference
#include "OscLib/OscCalcPMNSOpt.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/OscCurve.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Core/Utilities.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TMath.h"
#include <cmath>
#include <iostream>
#include "Utilities/func/MathUtil.h"

Go to the source code of this file.

Functions

void cc ()
 
void test_ana ()
 

Function Documentation

void cc ( )

Definition at line 28 of file test_ana.C.

References std::asin(), ana::bins, calc, om::cout, ana::DataMCComparisonComponents(), ana::ISurface::Draw(), ana::ISurface::DrawBestFit(), ana::ISurface::DrawContour(), e, allTimeWatchdog::endl, demo5::expt, ana::Spectrum::FakeData(), ana::IFitter::Fit(), ana::Gaussian68Percent2D(), ana::Gaussian90Percent2D(), osc::_IOscCalcAdjustable< T >::GetDmsq32(), osc::_IOscCalcAdjustable< T >::GetTh23(), ana::Loaders::Go(), ana::kCaloE, ana::kCrudeCCSel, ana::Loaders::kFHC, ana::kFitDmSq32, ana::kFitSinSq2Theta23, loaders, demo4::mock, plot_validation_datamc::pred, ana::PredictionExtrap::Predict(), ana::IPrediction::PredictUnoscillated(), ana::RatioPlot(), osc::_OscCalcPMNSOpt< T >::SetdCP(), osc::_OscCalcPMNSOpt< T >::SetDmsq21(), osc::_OscCalcPMNSOpt< T >::SetDmsq32(), osc::_OscCalcPMNSOpt< T >::SetL(), osc::_OscCalcPMNSOpt< T >::SetRho(), osc::_OscCalcPMNSOpt< T >::SetTh12(), osc::_OscCalcPMNSOpt< T >::SetTh13(), osc::_OscCalcPMNSOpt< T >::SetTh23(), ana::Binning::Simple(), std::sin(), util::sqr(), std::sqrt(), demo5::surf, and ana::Spectrum::ToTH1().

Referenced by Compare_NoExtrap(), Compare_Spectra(), CompareDecompDataMC(), ComparisonPlots_Data(), ComparisonPlots_MC(), Evaluate_BDTMLP_Algs_PredNoExtrap(), Evaluate_BDTMLP_Algs_Spectra(), GetBkgd(), ana::NuWROSyst::InitializeHistograms(), make_cm_pullterm_pdf(), MakeCutFlow(), MakeCutFlowPlots(), ana::MakeCutVec(), MakeNuEnergyPlots(), MakeNueTable(), MakeNuMuTable(), NoExtrap(), NuMu2019_BasicPIDPlots_FD(), NuMu2019_BasicPIDPlots_ND(), NuMu2019_BasicPIDPlots_Spectrum(), NuMu2020_BasicPIDPlots_FD(), NuMu2020_BasicPIDPlots_ND(), NuMu2020_BasicPIDPlots_Spectrum(), Plot_Evals_BDTMLP_Algs(), plot_michel_final(), PlotMCOnly(), PlotNME(), jmshower::JMTrackMerge::produce(), jmshower::JMClusterMerge::produce(), SystematicComp(), and WriteOutEntries().

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
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 
74 
75 
76  calc.SetTh23(1); // Seed value
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;
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 }
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
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:148
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.
virtual T GetTh23() const
Definition: IOscCalc.h:89
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:40
Log-likelihood scan across two parameters.
expt
Definition: demo5.py:34
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:341
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:69
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Binning bins
void SetTh23(const T &th23) override
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
void SetDmsq21(const T &dmsq21) override
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
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
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
Prediction that just uses FD MC, with no extrapolation.
mock
Definition: demo4.py:28
Float_t e
Definition: plot.C:35
TCanvas * RatioPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption, std::vector< TH1 * > bottomHistos, std::vector< TString > bottomOption, TString pidtag, AxisType pidaxis=kNue1bin)
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
void test_ana ( )

Definition at line 99 of file test_ana.C.

References std::asin(), calc, om::cout, ana::DataMCComparisonComponents(), ana::Loaders::DisableLoader(), ana::ISurface::Draw(), ana::ISurface::DrawBestFit(), ana::ISurface::DrawContour(), e, allTimeWatchdog::endl, demo5::expt, ana::Spectrum::FakeData(), ana::IFitter::Fit(), ana::Gaussian68Percent2D(), ana::Gaussian90Percent2D(), osc::_IOscCalcAdjustable< T >::GetdCP(), osc::_IOscCalcAdjustable< T >::GetTh13(), ana::Loaders::Go(), ana::kCosmic, caf::kFARDET, ana::Loaders::kFHC, ana::kFitDeltaInPiUnits, ana::kFitSinSq2Theta13, ana::kFitSinSq2Theta23, ana::kLEM, ana::Loaders::kMC, ana::kPassesPresel, loaders, demo4::mock, plot_validation_datamc::pred, ana::PredictionExtrap::Predict(), ana::IPrediction::PredictUnoscillated(), osc::_OscCalcPMNSOpt< T >::SetdCP(), osc::_OscCalcPMNSOpt< T >::SetDmsq21(), osc::_OscCalcPMNSOpt< T >::SetDmsq32(), osc::_OscCalcPMNSOpt< T >::SetL(), osc::_OscCalcPMNSOpt< T >::SetRho(), osc::_OscCalcPMNSOpt< T >::SetTh12(), osc::_OscCalcPMNSOpt< T >::SetTh13(), osc::_OscCalcPMNSOpt< T >::SetTh23(), ana::FitSinSq2Theta13::SetValue(), ana::FitSinSq2Theta23::SetValue(), ana::Binning::Simple(), ana::SimpleFOM(), std::sin(), util::sqr(), std::sqrt(), demo5::surf, ana::OscCurve::ToTH1(), and ana::Spectrum::ToTH1().

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 
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 
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
void SetTh13(const T &th13) override
void SetL(double L) override
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.
const Color_t kMC
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:165
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:40
Log-likelihood scan across two parameters.
expt
Definition: demo5.py:34
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:341
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:69
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
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
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
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
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