test_stanfit_statsonly.C

Go to the documentation of this file.

/*
 * test_stanfit_statsonly.C:
 *    Test the use of the MCMC tool Stan for fitting in CAFAna.
 *    This version tests an actual neutrino oscillation fit,
 *    albeit a statistics-only one.
 *
 *    Original author:  J. Wolcott <jwolcott@fnal.gov>
 *               Date:  November 2018
 */

#include "TCanvas.h"
#include "TGraph.h"
#include "TMarker.h"

#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Prod5Loaders.h"
#include "CAFAna/Core/FitVarWithPrior.h"
#include "CAFAna/Core/SystShifts.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/Bayesian1DMarginal.h"
#include "CAFAna/Fit/BayesianSurface.h"
#include "CAFAna/Fit/Priors.h"
#include "CAFAna/Fit/StanConfig.h"
#include "CAFAna/Fit/StanFitter.h"
#include "CAFAna/Fit/MCMCSamples.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcAnalytic.h"
#include "OscLib/OscCalcDMP.h"

#include "Utilities/func/MathUtil.h"

namespace test
{
  double POT = 10e20;
  const std::string SAVED_PRED_FILE = "/nova/ana/users/jwolcott/scratch/pred.root";
  const std::string SAVED_SAMPLES_FILE = "/nova/ana/users/jwolcott/scratch/mcmcsamples.root";

  double MOCKDATA_TH23 = TMath::Pi()/4; // 0.72;     // 0.72 radians --> 41 degrees
  double MOCKDATA_DM32 = 0.0025;

  void ResetCalculator(osc::IOscCalcAdjustable & calc)
  {
    calc.SetL(810);
    calc.SetRho(2.75);
    calc.SetDmsq21(7.6e-5);
    calc.SetDmsq32(2.35e-3);
    calc.SetTh12(asin(sqrt(.87))/2);
    calc.SetTh13(asin(sqrt(.10))/2);
//    calc.SetTh23(TMath::Pi()/4);
    calc.SetTh23(TMath::DegToRad() * 43);
    calc.SetdCP(0);
  }
}

void test_stanfit_statsonly(bool loadPredFromFile=true,
                            bool savePredToFile=false,
                            bool loadSamplesFromFile=true,
                            bool saveSamplesToFile=false)
{
  auto TFText = [](bool val) { return val ? "TRUE" : "FALSE"; };
  std::cout << std::endl;
  std::cout << "Loading predictions from file: " << TFText(loadPredFromFile) << std::endl;
  std::cout << "Storing predictions to file: " << TFText(savePredToFile) << std::endl;
  std::cout << "Loading MCMC samples from file: " << TFText(loadSamplesFromFile) << std::endl;
  std::cout << "Saving MCMC samples to file: " << TFText(saveSamplesToFile) << std::endl;
  std::cout << std::endl;

  // let's try a nice "easy" problem: numu disappearance.
  // start with pretty basic NOvA parameters
  // (incl. max mixing)
//  osc::OscCalcPMNSOpt calc;
//  osc::OscCalcAnalytic calc;
  osc::OscCalcDMP calc;
  test::ResetCalculator(calc);

  std::unique_ptr<ana::PredictionExtrap> pred;
  if (!loadPredFromFile)
  {
    ana::Prod5NomLoaders loaders(ana::kNumuConcat, ana::Loaders::kFHC);
    pred = std::make_unique<ana::PredictionNoExtrap>(
      loaders,
      ana::kNumuCCOptimisedAxis,
      ana::kNumu2020FD,
      ana::kNoShift,
      ana::kXSecCVWgt2020 * ana::kPPFXFluxCVWgt);
    loaders.Go();
  }
  else
  {
    TFile inf(test::SAVED_PRED_FILE.c_str());
    pred = ana::PredictionExtrap::LoadFrom(dynamic_cast<TDirectory*>(inf.Get("pred")));
  }

  if (savePredToFile && !loadPredFromFile)
  {
    TFile outf(test::SAVED_PRED_FILE.c_str(), "recreate");
    pred->SaveTo(&outf, "pred");
  }

  // pick a few test values for some mock data...
  calc.SetTh23(test::MOCKDATA_TH23);
  calc.SetDmsq32(test::MOCKDATA_DM32); //
  ana::Spectrum spec_pred = pred->Predict(&calc);
//  ana::Spectrum fakeData = spec_pred.MockData(test::POT);  // mock data, of course, means that best fit may not reproduce the seed parameters
  ana::Spectrum fakeData = spec_pred.FakeData(test::POT);
  ana::SingleSampleExperiment expt(pred.get(), fakeData);

  // now put the calc back to normal
  test::ResetCalculator(calc);

  // build my fit variables...
  ana::ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23(&ana::kFitSinSqTheta23, ana::PriorUniformInFitVar, "FlatSSTh23");
  ana::ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior(&ana::kFitDmSq32Scaled, ana::PriorUniformInFitVar, "FlatDmSq32");

  std::unique_ptr<ana::MCMCSamples> samples;
  if (!loadSamplesFromFile)
  {
    ana::StanFitter fitter(&expt, {&fitSsqTh23_UniformPriorSsqTh23, &fitDmSq32Scaled_UniformPrior}, {});
    ana::StanConfig config;
    config.num_warmup =  500;
    config.num_samples = 2000;
//    config.verbosity = ana::StanConfig::Verbosity::kEverything;
    config.verbosity = ana::StanConfig::Verbosity::kQuiet;

    fitter.SetStanConfig(config);
    ana::SystShifts systs;
    fitter.Fit(&calc, systs);
    samples = fitter.ExtractSamples();

    if (saveSamplesToFile)
    {
      TFile outf(test::SAVED_SAMPLES_FILE.c_str(), "recreate");
      samples->SaveTo(&outf, "samples");
    }
  }
  else
  {
    TFile inf(test::SAVED_SAMPLES_FILE.c_str());
    samples = ana::MCMCSamples::LoadFrom(&inf, "samples");
  }

  TCanvas c;
  ana::BayesianSurface surf = samples->MarginalizeTo(&fitSsqTh23_UniformPriorSsqTh23, 30, .3, .7,
                                                     &fitDmSq32Scaled_UniformPrior, 30, 2.2, 2.8);
  surf.Draw();
  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian1Sigma), 7, kGreen+2); // dashed
  surf.DrawBestFit(kGray);
  TMarker marker(util::sqr(sin(test::MOCKDATA_TH23)), test::MOCKDATA_DM32*1e3, kFullStar);
  marker.SetMarkerColor(kMagenta);
  marker.SetMarkerSize(3);
  marker.Draw();
  c.SaveAs("/nova/ana/users/jwolcott/scratch/test_stanfit_statsonly_surface_contour.png");

  auto marg_dm2 = samples->MarginalizeTo(&fitDmSq32Scaled_UniformPrior);
  auto bin_marg_dm2 = ana::Binning::Simple(100, 2, 3);
  auto h_marg_dm2 = marg_dm2.ToTH1(bin_marg_dm2);
  h_marg_dm2.Draw();
  c.SaveAs("/nova/ana/users/jwolcott/scratch/test_stanfit_statsonly_dm2_marg.png");
  std::cout << "1sigma credible interval(s) for dm2:" << std::endl;
  for (const auto & range : marg_dm2.QuantileRanges(ana::Quantile::kGaussian1Sigma, bin_marg_dm2))
    std::cout << "  (" << range.first << ", " << range.second << ")" << std::endl;

  c.Clear();
  fitSsqTh23_UniformPriorSsqTh23.SetValue(&calc, surf.GetBestFitX());
  fitDmSq32Scaled_UniformPrior.SetValue(&calc, surf.GetBestFitY());
  ana::DataMCComparison(fakeData, pred.get(), &calc, ana::kNoShift, ana::kBinDensity);
  c.SaveAs("/nova/ana/users/jwolcott/scratch/test_stanfit_statsonly_bestfitpred.png");

//  tree->Scan("*");

}