modularextrap_demo_nue.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"
15 #include "TFile.h"
16 #include "OscLib/OscCalcPMNSOpt.h"
17 #include "CAFAna/Core/Spectrum.h"
18 #include "TH1.h"
20 #include "TCanvas.h"
22 #include "CAFAna/Vars/FitVars.h"
24 
25 using namespace ana;
26 
28 {
29 
30  // POT for fake FD data
31  const double pot(1e20);
32 
33  // Specify which files to load
34  std::string haddDir("/nova/ana/users/tamsett/data/sam_definitions/");
35  SpectrumLoader loaderNDData(haddDir + "hadd_prod_decaf_S15-05-22a_"
36  "nd_numi_numu_contain_goodruns.root");
37  SpectrumLoader loaderNDMC(haddDir + "hadd_prod_decaf_S15-05-22a_"
38  "nd_genie_fhc_nonswap_numu_contain_s00-04_without_zero_entry_files.root");
39  SpectrumLoader loaderFDNonSwap(haddDir + "hadd_prod_decaf_S15-05-22a_"
40  "fd_genie_fhc_nonswap_numu_contain_matchedruns.root");
41  SpectrumLoader loaderFDSwap(haddDir + "hadd_prod_decaf_S15-05-22a_"
42  "fd_genie_fhc_fluxswap_numu_contain_matchedruns.root");
43  NullLoader loaderFDTauSwap;
44  NullLoader loaderFDCosmic;
45 
46  // Variables (from Vars/Vars.h) and binning for nue and ND numu
47  const HistAxis axis(
48  "Calorimetric Energy (GeV)", Binning::Simple(12, 0.5, 3.5), kCaloE );
49  const HistAxis axisNumu(
50  "Numu CC Energy Estimator (GeV)", Binning::Simple(30, 0.5, 3.5), kCCE );
51 
52  // Analysis cuts (from Cuts/Cuts.h), example only
53  const Cut cutFD( kFDContainNue && (kLEM > .64) );
54  const Cut cutND( kNDContainNue && (kLEM > .64) );
55  const Cut cutNDNumu( kNumuND );
56 
57  // Setup a nue and a numu decomposition
58  ProportionalDecomp decomp( loaderNDMC, loaderNDData, axis, cutND );
59  NumuDecomp decompNumu( loaderNDMC, loaderNDData, axisNumu, cutNDNumu );
60 
61  // Setup the extrapolation
62  auto extrap = std::make_unique<NueExtrap>(
63  loaderNDMC, loaderFDSwap, loaderFDNonSwap, loaderFDTauSwap,
64  decomp, decompNumu, axis, axisNumu, cutFD, cutND, cutNDNumu );
65 
66  // Make a prediction based on the extrapolation
67  PredictionExtrap pred( std::move(extrap) );
68  // Make a prediction based on the MC for fake FD data
69  PredictionNoExtrap predFDfake(
70  loaderFDNonSwap, loaderFDSwap,
71  axis.label, axis.bins, axis.var, cutFD );
72 
73  // Load from files
74  loaderNDData.Go();
75  loaderNDMC.Go();
76  loaderFDNonSwap.Go();
77  loaderFDSwap.Go();
78  loaderFDTauSwap.Go();
79  loaderFDCosmic.Go();
80 
81  // Open output file
82  TFile file("modularextrap_demo_nue.root","RECREATE");
83  // Save things in output file
84  extrap.SaveTo(& file, "nue_extrap") ;
85  extrap.SavePlotsNue( file.mkdir("nue_extrap_plots"), pot );
86  pred.SaveTo(& file, "nue_pred") ;
87  predFDfake.SaveTo(& file, "nue_pred_noextrap") ;
88 
89  // Set up oscillation calculator with default values
91  calc.SetL(810);
92  calc.SetRho(2.75);
93  calc.SetDmsq21(7.6e-5);
94  calc.SetDmsq32(2.35e-3);
95  calc.SetTh12(asin(sqrt(.87))/2);
96  calc.SetTh13(asin(sqrt(.10))/2);
97  calc.SetTh23(TMath::Pi()/4);
98  calc.SetdCP(0);
99 
100  // Make fake FD data
101  Spectrum fakeFDMC( predFDfake.Predict(&calc) );
102  Spectrum fakeFDData( fakeFDMC.FakeData(pot) );
103 
104  // An experiment is just a prediction plus data
105  SingleSampleExperiment expt( &pred, fakeFDData );
106 
107  TCanvas c;
108 
109  // Plot a chi-squared surface with best fit and contours
110  FrequentistSurface surf( &expt, &calc,
111  &kFitSinSq2Theta13, 50, 0., 0.5,
112  &kFitDeltaInPiUnits, 30, 0., 2. );
113  surf.Draw();
114  surf.DrawBestFit(kRed);
115  surf.DrawContour( Gaussian68Percent2D(surf), 7, kRed );
116  surf.DrawContour( Gaussian90Percent2D(surf), 1, kRed );
117 
118  c.Write("surface");
119 
120  file.Close();
121  std::cout << "output saved to modularextrap_demo_nue.root" << std::endl;
122 
123 }
const Var kLEM
PID
Definition: Vars.cxx:26
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh13(const T &th13) override
void SetL(double L) override
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
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 modularextrap_demo_nue()
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
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
const Cut kFDContainNue([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.nshwlid<=0) return false;const caf::SRShowerLIDProxy &shwlid=sr->vtx.elastic.fuzzyk.png[0].shwlid;const caf::SRVector3DProxy &start=shwlid.start;const caf::SRVector3DProxy &stop=shwlid.stop;return(std::min(start.X(), stop.X()) >-725 && std::max(start.X(), stop.X())< +725 && std::min(start.Y(), stop.Y()) >-725 && std::max(start.Y(), stop.Y())< +725 && std::min(start.Z(), stop.Z()) > 50 && std::max(start.Z(), stop.Z())< 5910);})
Definition: Cuts.h:25
const Cut kNumuND
Definition: NumuCuts.h:55
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Var kCCE
Definition: NumuVars.h:21
#define pot
virtual void Go() override
Load all the registered spectra.
void SetTh23(const T &th23) override
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
const Cut kNDContainNue([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.nshwlid<=0) return false;const caf::SRShowerLIDProxy &shwlid=sr->vtx.elastic.fuzzyk.png[0].shwlid;const caf::SRVector3DProxy &start=shwlid.start;const caf::SRVector3DProxy &stop=shwlid.stop;return(std::min(start.X(), stop.X()) >-150 && std::max(start.X(), stop.X())< +150 && std::min(start.Y(), stop.Y()) >-150 && std::max(start.Y(), stop.Y())< +150 && std::min(start.Z(), stop.Z()) > 50 && std::max(start.Z(), stop.Z())< 1230);})
Definition: Cuts.h:22
Splits Data proportionally according to MC.
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
TFile * file
Definition: cellShifts.C:17
Dummy loader that doesn&#39;t load any files.
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Prediction that just uses FD MC, with no extrapolation.
const HistAxis axisNumu
Take the output of an extrapolation and oscillate it as required.
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
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