fit_3flavor_withsysts.C
Go to the documentation of this file.
1 /*
2  * fit_3flavor_withsysts.C:
3  * Perform a 3-flavor oscillation fit to dm32, ssth23, dcp
4  * with all systematic uncertainty parameters
5  * using MCMC.
6  *
7  * Original author: J. Wolcott <jwolcott@fnal.gov>
8  * Date: September 2020
9  */
10 
11 
12 #include "TMath.h"
13 
14 #include "CAFAna/Core/SystShifts.h"
15 #include "CAFAna/Fit/StanConfig.h"
16 #include "CAFAna/Fit/StanFitter.h"
17 
18 #include "CAFAna/3flavor/Ana2020/joint_fit_2020_loader_tools.h"
19 
20 #include "OscLib/OscCalcAnalytic.h"
21 #include "OscLib/OscCalcDMP.h"
22 
25 
26 using namespace ana;
27 
28 namespace mcmc
29 {
30  const double MOCKDATA_TH23 = 41 * TMath::DegToRad();
31  const double MOCKDATA_DM32 = 0.0024; // normal hierarchy
32  const double MOCKDATA_DCP = (5./4.) * TMath::Pi(); // halfway between no CPV and max CPV
33 
34 }
35 
36 
37 // ==============================================================
38 
39 void fit_3flavor_withsysts(const std::string & fakeDataFile,
40  std::string fakeDataSpectraDir,
41  bool extrapolation=false,
42  bool fitSysts=true,
43  std::string dirPrefix=".",
44  std::string samplesFilename="samples_3flavor.root")
45 {
46 
47  std::cout << "Saving MCMC samples to file: " << mcmc::FullFilename(dirPrefix, samplesFilename) << std::endl;
48 
49  std::unique_ptr<osc::IOscCalcAdjustable> calc;
50  std::unique_ptr<Spectrum> fakeData;
51 
52  // build my fit variables...
53  std::vector<const IFitVar*> fitVars{&fitSsqTh23_UniformPriorSsqTh23,
56 
57  calc = std::make_unique<osc::OscCalcDMP>();
58  mcmc::ResetCalculator(*calc);
59 
60  // try the fake data
61  auto spectra = mcmc::LoadFakeData(fakeDataFile, fakeDataSpectraDir);
62 
63  // load the predictions and build a MultiExperiment out of them
64  auto preds = mcmc::LoadPreds();
65  auto expts = mcmc::BuildExperiments(spectra, preds);
67 
68  std::set<const ana::ISyst*> systs;
69  if (fitSysts)
70  systs = mcmc::SupportedSysts(expt.get());
71  std::vector<const ana::ISyst*> systVector(systs.begin(), systs.end());
72 
73  auto shifts = std::make_unique<GaussianPriorSystShifts>();
74 
76  cfg.num_warmup = 500;
77  cfg.num_samples = 1000;
78  cfg.max_depth = 15;
79  cfg.verbosity = StanConfig::Verbosity::kQuiet;
80 // cfg.verbosity = StanConfig::Verbosity::kEverything;
81  StanFitter fitter(expt.get(), fitVars, systVector);
82  fitter.SetStanConfig(cfg);
83  fitter.Fit(calc.get(), *shifts);
84 
85  std::unique_ptr<ana::MCMCSamples> warmup = fitter.ExtractSamples(ana::MemoryTupleWriter::WhichSamples::kWarmup);
86  std::unique_ptr<ana::MCMCSamples> samples = fitter.ExtractSamples(ana::MemoryTupleWriter::WhichSamples::kPostWarmup);
87 
88  mcmc::SaveToFile(mcmc::FullFilename(dirPrefix, samplesFilename),
89  fakeData.get(),
90  warmup.get(),
91  samples.get());
92 
93 
94 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string FullFilename(const std::string &dir, std::string file)
Definition: MCMC3FShared.h:53
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::unique_ptr< ana::IExperiment > BuildMultiExperiment(const std::vector< const ana::IExperiment * > &expts)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP
void SaveToFile(const std::string &fileName, const ana::Spectrum *fakeData, const ana::MCMCSamples *warmup, const ana::MCMCSamples *samples, const ana::SystShifts *trueSystPulls, const osc::IOscCalcAdjustable *trueOscParams)
Configuration parameters for the Stan MCMC fitter, bundled up together for easy storage and parameter...
Definition: StanConfig.h:6
std::vector< const ana::IExperiment * > ExptPtrs(const std::vector< std::unique_ptr< ana::IExperiment >> &exptVector)
Definition: MCMC3FShared.h:77
const double MOCKDATA_DCP
osc::OscCalcDumb calc
void SetStanConfig(const StanConfig &cfg)
Change the config used for Stan. See the StanConfig struct documentation for ideas.
Definition: StanFitter.h:224
expt
Definition: demo5.py:34
void fit_3flavor_withsysts(const std::string &fakeDataFile, std::string fakeDataSpectraDir, bool extrapolation=false, bool fitSysts=true, std::string dirPrefix=".", std::string samplesFilename="samples_3flavor.root")
std::set< const ana::ISyst * > SupportedSysts(const ana::IExperiment *expt)
Which systs does this experiment support?
const double MOCKDATA_TH23
std::vector< std::unique_ptr< ana::IExperiment > > BuildExperiments(const std::map< std::string, ana::Spectrum > &fakeData, const std::vector< ana::predictions > &preds)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior
OStream cout
Definition: OStream.cxx:6
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
const double MOCKDATA_DM32
void ResetCalculator(osc::IOscCalcAdjustable &calc)
std::vector< ana::predictions > LoadPreds(bool extrap)
Just wrap up all the args somewhere centralized for the moment.
Fitter type that bolts the Stan fitting tools onto CAFAna.
Definition: StanFitter.h:126
enum BeamMode string