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().

{
  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

(

)

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().

{
  //  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
}