test_stanfit_statsonly.C
Go to the documentation of this file.
1 /*
2  * test_stanfit_statsonly.C:
3  * Test the use of the MCMC tool Stan for fitting in CAFAna.
4  * This version tests an actual neutrino oscillation fit,
5  * albeit a statistics-only one.
6  *
7  * Original author: J. Wolcott <jwolcott@fnal.gov>
8  * Date: November 2018
9  */
10 
11 #include "TCanvas.h"
12 #include "TGraph.h"
13 #include "TMarker.h"
14 
15 #include "CAFAna/Analysis/Plots.h"
18 #include "CAFAna/Core/SystShifts.h"
23 #include "CAFAna/Fit/Priors.h"
24 #include "CAFAna/Fit/StanConfig.h"
25 #include "CAFAna/Fit/StanFitter.h"
26 #include "CAFAna/Fit/MCMCSamples.h"
28 #include "CAFAna/Vars/FitVars.h"
33 #include "CAFAna/Vars/XsecTunes.h"
34 
35 #include "OscLib/OscCalcPMNSOpt.h"
36 #include "OscLib/OscCalcAnalytic.h"
37 #include "OscLib/OscCalcDMP.h"
38 
39 #include "Utilities/func/MathUtil.h"
40 
41 namespace test
42 {
43  double POT = 10e20;
44  const std::string SAVED_PRED_FILE = "/nova/ana/users/jwolcott/scratch/pred.root";
45  const std::string SAVED_SAMPLES_FILE = "/nova/ana/users/jwolcott/scratch/mcmcsamples.root";
46 
47  double MOCKDATA_TH23 = TMath::Pi()/4; // 0.72; // 0.72 radians --> 41 degrees
48  double MOCKDATA_DM32 = 0.0025;
49 
51  {
52  calc.SetL(810);
53  calc.SetRho(2.75);
54  calc.SetDmsq21(7.6e-5);
55  calc.SetDmsq32(2.35e-3);
56  calc.SetTh12(asin(sqrt(.87))/2);
57  calc.SetTh13(asin(sqrt(.10))/2);
58 // calc.SetTh23(TMath::Pi()/4);
59  calc.SetTh23(TMath::DegToRad() * 43);
60  calc.SetdCP(0);
61  }
62 }
63 
64 void test_stanfit_statsonly(bool loadPredFromFile=true,
65  bool savePredToFile=false,
66  bool loadSamplesFromFile=true,
67  bool saveSamplesToFile=false)
68 {
69  auto TFText = [](bool val) { return val ? "TRUE" : "FALSE"; };
71  std::cout << "Loading predictions from file: " << TFText(loadPredFromFile) << std::endl;
72  std::cout << "Storing predictions to file: " << TFText(savePredToFile) << std::endl;
73  std::cout << "Loading MCMC samples from file: " << TFText(loadSamplesFromFile) << std::endl;
74  std::cout << "Saving MCMC samples to file: " << TFText(saveSamplesToFile) << std::endl;
76 
77  // let's try a nice "easy" problem: numu disappearance.
78  // start with pretty basic NOvA parameters
79  // (incl. max mixing)
80 // osc::OscCalcPMNSOpt calc;
81 // osc::OscCalcAnalytic calc;
84 
85  std::unique_ptr<ana::PredictionExtrap> pred;
86  if (!loadPredFromFile)
87  {
89  pred = std::make_unique<ana::PredictionNoExtrap>(
90  loaders,
95  loaders.Go();
96  }
97  else
98  {
99  TFile inf(test::SAVED_PRED_FILE.c_str());
100  pred = ana::PredictionExtrap::LoadFrom(dynamic_cast<TDirectory*>(inf.Get("pred")));
101  }
102 
103  if (savePredToFile && !loadPredFromFile)
104  {
105  TFile outf(test::SAVED_PRED_FILE.c_str(), "recreate");
106  pred->SaveTo(&outf, "pred");
107  }
108 
109  // pick a few test values for some mock data...
112  ana::Spectrum spec_pred = pred->Predict(&calc);
113 // ana::Spectrum fakeData = spec_pred.MockData(test::POT); // mock data, of course, means that best fit may not reproduce the seed parameters
114  ana::Spectrum fakeData = spec_pred.FakeData(test::POT);
115  ana::SingleSampleExperiment expt(pred.get(), fakeData);
116 
117  // now put the calc back to normal
118  test::ResetCalculator(calc);
119 
120  // build my fit variables...
123 
124  std::unique_ptr<ana::MCMCSamples> samples;
125  if (!loadSamplesFromFile)
126  {
127  ana::StanFitter fitter(&expt, {&fitSsqTh23_UniformPriorSsqTh23, &fitDmSq32Scaled_UniformPrior}, {});
129  config.num_warmup = 500;
130  config.num_samples = 2000;
131 // config.verbosity = ana::StanConfig::Verbosity::kEverything;
132  config.verbosity = ana::StanConfig::Verbosity::kQuiet;
133 
134  fitter.SetStanConfig(config);
136  fitter.Fit(&calc, systs);
137  samples = fitter.ExtractSamples();
138 
139  if (saveSamplesToFile)
140  {
141  TFile outf(test::SAVED_SAMPLES_FILE.c_str(), "recreate");
142  samples->SaveTo(&outf, "samples");
143  }
144  }
145  else
146  {
147  TFile inf(test::SAVED_SAMPLES_FILE.c_str());
148  samples = ana::MCMCSamples::LoadFrom(&inf, "samples");
149  }
150 
151  TCanvas c;
152  ana::BayesianSurface surf = samples->MarginalizeTo(&fitSsqTh23_UniformPriorSsqTh23, 30, .3, .7,
153  &fitDmSq32Scaled_UniformPrior, 30, 2.2, 2.8);
154  surf.Draw();
155  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian1Sigma), 7, kGreen+2); // dashed
156  surf.DrawBestFit(kGray);
157  TMarker marker(util::sqr(sin(test::MOCKDATA_TH23)), test::MOCKDATA_DM32*1e3, kFullStar);
158  marker.SetMarkerColor(kMagenta);
159  marker.SetMarkerSize(3);
160  marker.Draw();
161  c.SaveAs("/nova/ana/users/jwolcott/scratch/test_stanfit_statsonly_surface_contour.png");
162 
163  auto marg_dm2 = samples->MarginalizeTo(&fitDmSq32Scaled_UniformPrior);
164  auto bin_marg_dm2 = ana::Binning::Simple(100, 2, 3);
165  auto h_marg_dm2 = marg_dm2.ToTH1(bin_marg_dm2);
166  h_marg_dm2.Draw();
167  c.SaveAs("/nova/ana/users/jwolcott/scratch/test_stanfit_statsonly_dm2_marg.png");
168  std::cout << "1sigma credible interval(s) for dm2:" << std::endl;
169  for (const auto & range : marg_dm2.QuantileRanges(ana::Quantile::kGaussian1Sigma, bin_marg_dm2))
170  std::cout << " (" << range.first << ", " << range.second << ")" << std::endl;
171 
172  c.Clear();
173  fitSsqTh23_UniformPriorSsqTh23.SetValue(&calc, surf.GetBestFitX());
174  fitDmSq32Scaled_UniformPrior.SetValue(&calc, surf.GetBestFitY());
175  ana::DataMCComparison(fakeData, pred.get(), &calc, ana::kNoShift, ana::kBinDensity);
176  c.SaveAs("/nova/ana/users/jwolcott/scratch/test_stanfit_statsonly_bestfitpred.png");
177 
178 // tree->Scan("*");
179 
180 }
virtual void SetL(double L)=0
TFile * inf
Definition: drawXsec.C:1
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcDMP.h:55
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
virtual void SetDmsq21(const T &dmsq21)=0
double MOCKDATA_TH23
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
TH2 * QuantileSurface(Quantile quantile) const
T sqrt(T number)
Definition: d0nt_math.hpp:156
Configuration parameters for the Stan MCMC fitter, bundled up together for easy storage and parameter...
Definition: StanConfig.h:6
Helper struct for the cache. Might not need this.
Definition: StanTypedefs.h:25
Divide bin contents by bin widths.
Definition: Utilities.h:45
virtual void SetTh13(const T &th13)=0
const std::string SAVED_SAMPLES_FILE
double GetBestFitY() const
Definition: ISurface.h:31
const std::string SAVED_PRED_FILE
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
c
Definition: test.py:13
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
expt
Definition: demo5.py:34
static std::unique_ptr< MCMCSamples > LoadFrom(TDirectory *dir, const std::string &name)
Load from file.
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod5Loaders.h:101
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
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
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior
void ResetCalculator(osc::IOscCalcAdjustable &calc)
TFile * outf
Definition: testXsec.C:51
static std::unique_ptr< PredictionExtrap > LoadFrom(TDirectory *dir, const std::string &name)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
double MOCKDATA_DM32
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual void SetRho(double rho)=0
T sin(T number)
Definition: d0nt_math.hpp:132
void SetTh23(const T &th23) override
Definition: OscCalcDMP.h:58
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
std::vector< Loaders * > loaders
Definition: syst_header.h:386
virtual void SetTh23(const T &th23)=0
stan::math::var PriorUniformInFitVar(const stan::math::var &, const osc::IOscCalcAdjustableStan *)
Prior uniform in the variable that is being fitted.
Definition: Priors.cxx:10
Version of FitVarWithPrior for use with constrained FitVar_StanSupports.
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
void test_stanfit_statsonly(bool loadPredFromFile=true, bool savePredToFile=false, bool loadSamplesFromFile=true, bool saveSamplesToFile=false)
Float_t e
Definition: plot.C:35
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
Fitter type that bolts the Stan fitting tools onto CAFAna.
Definition: StanFitter.h:126
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106
double GetBestFitX() const
Definition: ISurface.h:30
surf
Definition: demo5.py:35