test_ana.C

Go to the documentation of this file.

#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"

using namespace ana;

void cc()
{
  Loaders loaders("$NOVA_DATA/mc/S13-06-18/hadd/", Loaders::kFHC);

  // Define our binning
  Binning bins = Binning::Simple(40, 0, 4);

  // Load up the necessary histograms, seperated by flavour etc
  PredictionNoExtrap pred(loaders,
                          "Calorimetric energy", bins,
                          kCaloE, kCrudeCCSel);

  loaders.Go();

  // Make a calculator. Any of the variants is fine
  osc::OscCalcPMNSOpt calc;

  calc.SetL(810);
  calc.SetRho(0); // No matter effects
  calc.SetDmsq21(7.6e-5);
  calc.SetDmsq32(2.35e-3);
  calc.SetTh12(asin(sqrt(.87))/2); // PDG
  calc.SetTh13(asin(sqrt(.10))/2);
  calc.SetTh23(TMath::Pi()/4);
  calc.SetdCP(0);

  // What we expect to observe at those oscillation parameters
  Spectrum obs = pred.Predict(&calc);

  // Make it into mock data
  //  Spectrum mock = obs.MockData(18e20);
  Spectrum mock = obs.FakeData(18e20);

  Spectrum unosc = pred.PredictUnoscillated();

  new TCanvas;
  //  DataMCComparison(mock, obs);
  DataMCComparisonComponents(mock, &pred, &calc);
  TH1* uoh = unosc.ToTH1(18e20);
  uoh->SetLineColor(kBlue);
  uoh->Draw("same");

  new TCanvas;
  RatioPlot(mock, unosc, obs);

  SingleSampleExperiment expt(&pred, mock);


  calc.SetTh23(1); // Seed value
  MinuitFitter fit(&expt, {&kFitSinSq2Theta23, &kFitDmSq32});
  const double chi = fit.Fit(&calc)->EvalMetricVal();
  std::cout << "Best fit, dmsq = " << calc.GetDmsq32() << " theta = " << calc.GetTh23() << " (ss2th23 = " << util::sqr(sin(2*calc.GetTh23())) << ") LL = " << chi << std::endl;

  // The log-likelihood surface
  new TCanvas;
  FrequentistSurface surf(&expt, &calc,
               &kFitSinSq2Theta23, 20, .9, 1,
               &kFitDmSq32, 20, 2e-3, 3e-3);
  surf.Draw();
  surf.DrawBestFit(kRed);

  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
  surf.DrawContour(Gaussian90Percent2D(surf), kSolid, kRed);

  new TCanvas;
  surf.DrawBestFit(kRed);

  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
  surf.DrawContour(Gaussian90Percent2D(surf), kSolid, kRed);
}

void test_ana()
{
  //  cc();
  //  return;

  Loaders loaders("$NOVA_DATA/mc/S13-06-18/hadd/", Loaders::kFHC);
  loaders.DisableLoader(caf::kFARDET, Loaders::kMC, ana::kCosmic);

  // Make a calculator. Any of the variants is fine
  osc::OscCalcPMNSOpt calc;

  calc.SetL(810);
  calc.SetRho(2.75);
  calc.SetDmsq21(7.6e-5);
  calc.SetDmsq32(2.35e-3);
  calc.SetTh12(asin(sqrt(.87))/2); // PDG
  calc.SetTh13(asin(sqrt(.10))/2);
  calc.SetTh23(TMath::Pi()/4);
  calc.SetdCP(0);

  // Make sure we know what we're doing for sure
  kFitSinSq2Theta13.SetValue(&calc, 0.1);
  kFitSinSq2Theta23.SetValue(&calc, 1);

  // We can make oscillation curves
  OscCurve curve(&calc, 14, 12);
  new TCanvas;
  curve.ToTH1()->Draw();

  OscCurve curve2(&calc, 14, 14);
  new TCanvas;
  curve2.ToTH1()->Draw();



  Binning bins2 = Binning::Simple(30, -.25, 1.25);

  PredictionNoExtrap pred(loaders,
                          "LEM PID", bins2,
                          //                          "LEM PID", 14, -.2, 1.2,
                          kLEM, kPassesPresel);
  //                      kJmID, kPassesPresel);

  loaders.Go();

  // What we expect to observe at those oscillation parameters
  Spectrum obs = pred.Predict(&calc);

  // Make it into mock data
  //  Spectrum mock = obs.MockData(18e20);
  Spectrum mock = obs.FakeData(18e20);

  new TCanvas;
  //  DataMCComparison(mock, obs);
  DataMCComparisonComponents(mock, &pred, &calc);

  // Overlay prediction at no oscillations
  TH1* hExp = pred.PredictUnoscillated().ToTH1(18e20);
  hExp->SetLineColor(kRed);
  hExp->Draw("same");
  hExp->SetLineStyle(7);
  calc.SetTh13(0);
  // And at no theta13
  TH1* hExp2 = pred.Predict(&calc).ToTH1(18e20);
  hExp2->SetLineColor(kGreen+2);
  hExp2->Draw("same");

  //  std::cout << "LL to no osc = " << LogLikelihood(pred.PredictUnoscillated(), mock) << std::endl;
  //  std::cout << "LL to input = " << LogLikelihood(obs, mock) << std::endl;
  std::cout << "FOM = " << SimpleFOM(mock, pred.Predict(&calc)) << std::endl;


  SingleSampleExperiment expt(&pred, mock);

  MinuitFitter fit(&expt, {&kFitSinSq2Theta13, &kFitDeltaInPiUnits});
  const double chi = fit.Fit(&calc)->EvalMetricVal();
  std::cout << "Best fit, th13 = " << calc.GetTh13() << " (ss2th13 = " << util::sqr(sin(2*calc.GetTh13())) << ") delta = " << calc.GetdCP() << " LL = " << chi << std::endl;

  // The log-likelihood surface
  new TCanvas;
  FrequentistSurface surf(&expt, &calc,
               &kFitSinSq2Theta13, 50, 0, .4,
               &kFitDeltaInPiUnits, 50, 0, 2);
  surf.Draw();
  surf.DrawBestFit(kRed);

  surf.DrawContour(Gaussian68Percent2D(surf), 7, kRed); // dashed
  surf.DrawContour(Gaussian90Percent2D(surf), 1, kRed); // solid

  // Inverted hierarchy
  calc.SetDmsq32(-2.35e-3);

  FrequentistSurface surfInv(&expt, &calc,
                  &kFitSinSq2Theta13, 50, 0, .4,
                  &kFitDeltaInPiUnits, 50, 0, 2);
  surfInv.DrawBestFit(kBlue);

  surfInv.DrawContour(Gaussian68Percent2D(surf), 7, kBlue); // dashed
  surfInv.DrawContour(Gaussian90Percent2D(surf), 1, kBlue); // solid
}