modularextrap_demo_numu.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "3FlavorAna/Extrap/NumuExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "TFile.h"
#include "OscLib/OscCalcPMNSOpt.h"
#include "CAFAna/Core/Spectrum.h"
#include "TH1.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "TCanvas.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Vars/FitVars.h"

using namespace ana;

void modularextrap_demo_numu()
{

  // POT for fake FD data
  const double pot(1e20);

  // Specify which files to load
  std::string haddDir("/nova/ana/users/tamsett/data/sam_definitions/");
  SpectrumLoader loaderNDData(haddDir + "hadd_prod_decaf_S15-05-22a_"
    "nd_numi_numu_contain_goodruns.root");
  SpectrumLoader loaderNDMC(haddDir + "hadd_prod_decaf_S15-05-22a_"
    "nd_genie_fhc_nonswap_numu_contain_s00-04_without_zero_entry_files.root");
  SpectrumLoader loaderFDNonSwap(haddDir + "hadd_prod_decaf_S15-05-22a_"
    "fd_genie_fhc_nonswap_numu_contain_matchedruns.root");
  SpectrumLoader loaderFDSwap(haddDir + "hadd_prod_decaf_S15-05-22a_"
    "fd_genie_fhc_fluxswap_numu_contain_matchedruns.root");
  NullLoader loaderFDTauSwap;
  NullLoader loaderFDCosmic;

  // Variables (from Vars/Vars.h) and binning for nue and ND numu
  const HistAxis axis( "Quasielastic Energy Estimator (GeV)",
    Binning::Simple(30, 0.5, 3.5), kQEE );

  // Analysis cuts (from Cuts/Cuts.h), example only
  const Cut cutFD( kNumuQEFD );
  const Cut cutND( kNumuQEND );

  // Setup a nue and a numu decomposition
  NumuDecomp decomp( loaderNDMC, loaderNDData, axis, cutND );

  // Setup the extrapolation
  NumuExtrap extrap(
    loaderNDMC, loaderFDSwap, loaderFDNonSwap, loaderFDTauSwap,
    decomp, axis, cutFD, cutND );

  // Make a prediction based on the extrapolation
  PredictionExtrap pred( &extrap );
  // Make a prediction based on the MC for fake FD data
  PredictionNoExtrap predFDfake(
    loaderFDNonSwap, loaderFDSwap,
    axis.label, axis.bins, axis.var, cutFD );

  // Load from files
  loaderNDData.Go();
  loaderNDMC.Go();
  loaderFDNonSwap.Go();
  loaderFDSwap.Go();
  loaderFDTauSwap.Go();
  loaderFDCosmic.Go();

  // Open output file
  TFile file("modularextrap_demo_numu.root","RECREATE");
  // Save things in output file
  extrap.SaveTo(& file, "numu_extrap") ;
  extrap.SavePlotsNumu( file.mkdir("numu_extrap_plots"), pot );
  pred.SaveTo(& file, "numu_pred") ;
  predFDfake.SaveTo(& file, "numu_pred_noextrap") ;

  // Set up oscillation calculator with default values
  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);
  calc.SetTh13(asin(sqrt(.10))/2);
  calc.SetTh23(TMath::Pi()/4);
  calc.SetdCP(0);

  // Make fake FD data
  Spectrum fakeFDMC( predFDfake.Predict(&calc) );
  Spectrum fakeFDData( fakeFDMC.FakeData(pot) );

  // An experiment is just a prediction plus data
  SingleSampleExperiment expt( &pred, fakeFDData );

  TCanvas c;

  // Plot a chi-squared surface with best fit and contours
  FrequentistSurface surf( &expt, &calc,
    &kFitSinSq2Theta23, 30, .8, 1.,
    &kFitDmSq32Scaled, 30, 1.5, 3.5 );
  surf.Draw();
  surf.DrawBestFit(kRed);
  surf.DrawContour( Gaussian68Percent2D(surf), 7, kRed );
  surf.DrawContour( Gaussian90Percent2D(surf), 1, kRed );

  c.Write("surface");

  file.Close();
  std::cout << "output saved to modularextrap_demo_numu.root" << std::endl;

}