modularextrap_demo_numu.C
Go to the documentation of this file.
3 #include "CAFAna/Core/HistAxis.h"
4 #include "CAFAna/Vars/Vars.h"
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Cut.h"
7 #include "CAFAna/Cuts/Cuts.h"
13 #include "TFile.h"
14 #include "OscLib/OscCalcPMNSOpt.h"
15 #include "CAFAna/Core/Spectrum.h"
16 #include "TH1.h"
18 #include "TCanvas.h"
20 #include "CAFAna/Vars/FitVars.h"
21 
22 using namespace ana;
23 
25 {
26 
27  // POT for fake FD data
28  const double pot(1e20);
29 
30  // Specify which files to load
31  std::string haddDir("/nova/ana/users/tamsett/data/sam_definitions/");
32  SpectrumLoader loaderNDData(haddDir + "hadd_prod_decaf_S15-05-22a_"
33  "nd_numi_numu_contain_goodruns.root");
34  SpectrumLoader loaderNDMC(haddDir + "hadd_prod_decaf_S15-05-22a_"
35  "nd_genie_fhc_nonswap_numu_contain_s00-04_without_zero_entry_files.root");
36  SpectrumLoader loaderFDNonSwap(haddDir + "hadd_prod_decaf_S15-05-22a_"
37  "fd_genie_fhc_nonswap_numu_contain_matchedruns.root");
38  SpectrumLoader loaderFDSwap(haddDir + "hadd_prod_decaf_S15-05-22a_"
39  "fd_genie_fhc_fluxswap_numu_contain_matchedruns.root");
40  NullLoader loaderFDTauSwap;
41  NullLoader loaderFDCosmic;
42 
43  // Variables (from Vars/Vars.h) and binning for nue and ND numu
44  const HistAxis axis( "Quasielastic Energy Estimator (GeV)",
45  Binning::Simple(30, 0.5, 3.5), kQEE );
46 
47  // Analysis cuts (from Cuts/Cuts.h), example only
48  const Cut cutFD( kNumuQEFD );
49  const Cut cutND( kNumuQEND );
50 
51  // Setup a nue and a numu decomposition
52  NumuDecomp decomp( loaderNDMC, loaderNDData, axis, cutND );
53 
54  // Setup the extrapolation
55  NumuExtrap extrap(
56  loaderNDMC, loaderFDSwap, loaderFDNonSwap, loaderFDTauSwap,
57  decomp, axis, cutFD, cutND );
58 
59  // Make a prediction based on the extrapolation
60  PredictionExtrap pred( &extrap );
61  // Make a prediction based on the MC for fake FD data
62  PredictionNoExtrap predFDfake(
63  loaderFDNonSwap, loaderFDSwap,
64  axis.label, axis.bins, axis.var, cutFD );
65 
66  // Load from files
67  loaderNDData.Go();
68  loaderNDMC.Go();
69  loaderFDNonSwap.Go();
70  loaderFDSwap.Go();
71  loaderFDTauSwap.Go();
72  loaderFDCosmic.Go();
73 
74  // Open output file
75  TFile file("modularextrap_demo_numu.root","RECREATE");
76  // Save things in output file
77  extrap.SaveTo(& file, "numu_extrap") ;
78  extrap.SavePlotsNumu( file.mkdir("numu_extrap_plots"), pot );
79  pred.SaveTo(& file, "numu_pred") ;
80  predFDfake.SaveTo(& file, "numu_pred_noextrap") ;
81 
82  // Set up oscillation calculator with default values
84  calc.SetL(810);
85  calc.SetRho(2.75);
86  calc.SetDmsq21(7.6e-5);
87  calc.SetDmsq32(2.35e-3);
88  calc.SetTh12(asin(sqrt(.87))/2);
89  calc.SetTh13(asin(sqrt(.10))/2);
90  calc.SetTh23(TMath::Pi()/4);
91  calc.SetdCP(0);
92 
93  // Make fake FD data
94  Spectrum fakeFDMC( predFDfake.Predict(&calc) );
95  Spectrum fakeFDData( fakeFDMC.FakeData(pot) );
96 
97  // An experiment is just a prediction plus data
98  SingleSampleExperiment expt( &pred, fakeFDData );
99 
100  TCanvas c;
101 
102  // Plot a chi-squared surface with best fit and contours
103  FrequentistSurface surf( &expt, &calc,
104  &kFitSinSq2Theta23, 30, .8, 1.,
105  &kFitDmSq32Scaled, 30, 1.5, 3.5 );
106  surf.Draw();
107  surf.DrawBestFit(kRed);
108  surf.DrawContour( Gaussian68Percent2D(surf), 7, kRed );
109  surf.DrawContour( Gaussian90Percent2D(surf), 1, kRed );
110 
111  c.Write("surface");
112 
113  file.Close();
114  std::cout << "output saved to modularextrap_demo_numu.root" << std::endl;
115 
116 }
void modularextrap_demo_numu()
enum BeamMode kRed
void SaveTo(TDirectory *dir, const std::string &name) const override
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh13(const T &th13) override
void SetL(double L) override
const FitSinSq2Theta23 kFitSinSq2Theta23
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
osc::OscCalcDumb calc
virtual void Go() override
Load all the registered spectra.
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Cut kNumuQEND
Definition: NumuCuts.h:66
Log-likelihood scan across two parameters.
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
expt
Definition: demo5.py:34
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
void SavePlotsNumu(TDirectory *dir, double potFD) const
Spectrum Predict(osc::IOscCalc *calc) const override
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
#define pot
virtual void Go() override
Load all the registered spectra.
void SetTh23(const T &th23) override
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
void SetDmsq21(const T &dmsq21) override
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
void SetdCP(const T &dCP) override
void SetTh12(const T &th12) override
const Var kQEE
Energy estimator for quasielastic CC events.
Definition: NumuVars.cxx:28
TFile * file
Definition: cellShifts.C:17
Dummy loader that doesn&#39;t load any files.
const Cut kNumuQEFD
Definition: NumuCuts.h:65
Prediction that just uses FD MC, with no extrapolation.
Take the output of an extrapolation and oscillate it as required.
Float_t e
Definition: plot.C:35
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
void SetRho(double rho) override
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
enum BeamMode string