template_basic.C
Go to the documentation of this file.
1 // most basic template. FHC only, no cosmics file, no systematics, 2 choices of theta_23
2 
4 #include "OscLib/OscCalcPMNS.h"
5 #include "OscLib/OscCalc.h"
6 
7 #include "CAFAna/Cuts/Cuts.h"
9 #include "CAFAna/Fit/Fit.h"
10 #include "CAFAna/Vars/FitVars.h"
12 #include "CAFAna/Analysis/Plots.h"
16 #include "CAFAna/Core/Utilities.h"
18 #include "CAFAna/Vars/Vars.h"
20 #include "CAFAna/Core/Binning.h"
21 
22 
23 #include "TCanvas.h"
24 #include "TFile.h"
25 #include "TGraph.h"
26 #include "TH1.h"
27 #include "TMath.h"
28 #include "TGaxis.h"
29 #include "TMultiGraph.h"
30 #include "TLegend.h"
31 #include "TLatex.h"
32 #include "TStyle.h"
33 #include "THStack.h"
34 #include "TPaveText.h"
35 
36 #include <cmath>
37 #include <iostream>
38 #include <sstream>
39 #include <string>
40 #include <fstream>
41 #include <iomanip>
42 
43 #include "Utilities/func/MathUtil.h"
44 
45 // Use with Nova style ROOT logon script
46 
47 #define location 0 // 0 for FNAL 1 for caltech
48 #define mcset 1 // 0 = ideal MC goes with pot like 18e20, or 1 = realistic conditions mc for 1st analysis POT
49 #define pot 18e20
50 
51 using namespace ana;
52 
53 void Simulation() // for some reason doesn't always remember definition in Style script
54 {
55  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Simulation");
56  prelim->SetTextColor(kGray+1);
57  prelim->SetNDC();
58  prelim->SetTextSize(2/30.);
59  prelim->SetTextAlign(32);
60  prelim->Draw();
61 }
62 
64 {
65 
66  TGaxis::SetMaxDigits(2);
67 
68  gStyle->SetTitleOffset(.85, "y");
69  gStyle->SetTitleOffset(.84, "x");
70  gStyle->SetTitleSize(0.05,"x");
71 
73 
74  if(mcset==1) {
75 
76  if(location==1) {
77  fname = "$NOVA_DATA/mc/forFA/fdgenie.nonswap.decaf.root";
78  fnameSwap = "$NOVA_DATA/mc/forFA/fdgenie.swap.decaf.root";
79  }
80  else if(location==0) {
81  fname = "prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_numu_contain";
82  fnameSwap = "prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_numu_contain";
83  }
84  }
85 
86  else if(mcset==0) {
87 
88  if(location==1) {
89  fname = "prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_ideal_numu_contain";
90  fnameSwap = "prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_ideal_numu_contain";
91  }
92  else if(location==0) {
93  fname = "prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_ideal_numu_contain";
94  fnameSwap = "prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_ideal_numu_contain";
95  }
96  }
97 
99 
100  calc.SetL(810);
101  calc.SetRho(0); // No matter effects
102  calc.SetDmsq21(7.59e-5);
103  calc.SetDmsq32(2.4e-3);
104  calc.SetTh12(.601);
105  calc.SetTh13(.1567);
106  // calc.SetdCP(TMath::Pi()/2);
107  calc.SetdCP(0);
108  calc.SetTh23(TMath::Pi()/4); // sin^2(2theta) = 1
109 
110  const std::string filename = fname;
111  const std::string filenameSwap = fnameSwap;
112 
113  SpectrumLoader loader(filename);
114  SpectrumLoader loaderSwap(filenameSwap);
115 
116  // Load up the necessary histograms, seperated by flavour etc
117  PredictionNoExtrap predCC(loader, loaderSwap,
118  "CC Neutrino Energy (GeV)", kNumuEnergyBinning,
119  kCCE, kNumuFD);
120 
121  loader.Go();
122  loaderSwap.Go();
123 
124  Spectrum obsCC = predCC.Predict(&calc); // the prediction
125  Spectrum fakeCC = obsCC.FakeData(pot); // the data, faked
126  SingleSampleExperiment exptCC(&predCC, fakeCC); // made into an 'experiment'
127 
128  TLatex text;
129  text.SetNDC(true);
130  text.SetTextSize(0.04);
131 
132  // now, the same but with a different truth value of theta_23
133  calc.SetTh23(.67264); // sin^2(2theta) = 0.95
134 
135  Spectrum obsCC2 = predCC.Predict(&calc);
136  Spectrum fakeCC2 = obsCC2.FakeData(pot);
137  SingleSampleExperiment exptCC2(&predCC, fakeCC2);
138 
139  new TCanvas("comp1b","CC spectra .95"); // plot spectra for each sample
140  DataMCComparison(fakeCC2, obsCC2);
141  TH1* ccu2 = predCC.PredictUnoscillated().ToTH1(pot);
142  ccu2->SetLineStyle(7);
143  ccu2->Draw("hist same");
144  gPad->Print("plots/comp1b.ps");
145 
146  calc.SetTh23(1); // Seed value
147 
148  TH1F *hcc = new TH1F("cc", "cc", 10, 0, 1);
149  TH1F *hqe = new TH1F("qe", "qe", 10, 0, 1);
150  TH1F *huc = new TH1F("uc", "uc", 10, 0, 1);
151  TH1F *hall = new TH1F("all", "all", 10, 0, 1);
152 
153  hcc->SetLineWidth(2);
154  hqe->SetLineWidth(2);
155  huc->SetLineWidth(2);
156  hall->SetLineWidth(2);
157  hcc->SetLineColor(2);
158  hqe->SetLineColor(4);
159  huc->SetLineColor(kMagenta);
160  hall->SetLineColor(1);
161 
162  //The log-likelihood surface: FHC
163  new TCanvas("contoursfhc1","contoursfhc1") ;
164 
165  FrequentistSurface surfCC(&exptCC, &calc,
166  &kFitSinSq2Theta23, 30, .9, 1,
167  &kFitDmSq32Scaled, 30, 2.2, 2.7);
168  surfCC.DrawBestFit(kBlack);
169  surfCC.DrawContour(Gaussian90Percent2D(surfCC), kSolid, kBlack);
170 
171  Simulation();
172 
173  TLegend *leg = new TLegend(0.11,0.11,0.42,0.35);
174  leg->SetFillColor(10);
175  leg->SetFillStyle(0);
176  leg->SetBorderSize(0);
177  if(pot==18e20) leg->SetHeader("18e20 POT, FHC, 14 kton: 90% CL");
178  else if(pot==6e20) leg->SetHeader("6e20 POT, FHC, 14 kton: 90% CL");
179  else if(pot==2e20) leg->SetHeader("2e20 POT, FHC, 14 kton: 90% CL");
180  else if(pot==1e20) leg->SetHeader("1e20 POT, FHC, 14 kton: 90% CL");
181 
182  leg->AddEntry(hall,"Contained CC","L");
183  leg->Draw("same");
184  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
185 
186  gPad->Print("plots/basic.1.ps");
187  gPad->Print("plots/basic.1.C");
188 
189  //2
190 
191  new TCanvas("contoursfhc2","contoursfhc2") ;
192 
193  FrequentistSurface surfCC2(&exptCC2, &calc,
194  &kFitSinSq2Theta23, 30, .87, 1,
195  &kFitDmSq32Scaled, 30, 2.2, 2.7);
196  surfCC2.DrawBestFit(kRed);
197  surfCC2.DrawContour(Gaussian90Percent2D(surfCC2), kSolid, kRed);
198 
199  leg->Draw("same");
200  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
201 
202  Simulation();
203 
204  gPad->Print("plots/basic.fhc.2.ps");
205  gPad->Print("plots/basic.fhc.2.C");
206 
207  new TCanvas("contourscombine","contours90percent") ;
208  surfCC2.DrawBestFit(kRed);
209  surfCC2.DrawContour(Gaussian2Sigma2D(surfCC2), kSolid, kRed);
210  surfCC2.DrawContour(Gaussian68Percent2D(surfCC2), 7, kRed);
211 
212  surfCC.DrawBestFit(kBlack);
213  surfCC.DrawContour(Gaussian2Sigma2D(surfCC), kSolid, kBlack);
214  surfCC.DrawContour(Gaussian68Percent2D(surfCC), 7, kBlack);
215 
216  TH1F *solid = new TH1F("90% CL", "90% CL", 10, 0, 1);
217  TH1F *dashed = new TH1F("68% CL", "68% CL", 10, 0, 1);
218 
219  solid->SetLineWidth(2);
220  dashed->SetLineWidth(2);
221  dashed->SetLineStyle(7);
222  solid->SetLineColor(1);
223  dashed->SetLineColor(1);
224 
225  TLegend *leg4 = new TLegend(0.12,0.12,0.5,0.4);
226  leg4->SetFillColor(10);
227  leg4->SetFillStyle(0);
228  leg4->SetHeader("18e20 POT FHC, 14 kton");
229 
230  leg4->AddEntry(solid, "2 #sigma","L");
231  leg4->AddEntry(dashed,"1 #sigma","L");
232  leg4->SetBorderSize(0);
233  leg4->Draw("");
234  text.DrawLatex(0.3,0.92,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
235 
236  gPad->Print("plots/basic.overlay.ps");
237  gPad->Print("plots/basic.overlay.C");
238 }
void SetDmsq21(const T &dmsq21) override
Definition: OscCalcPMNS.h:35
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
#define location
void SetTh12(const T &th12) override
Definition: OscCalcPMNS.h:37
void SetTh23(const T &th23) override
Definition: OscCalcPMNS.h:39
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
Adapt the PMNS calculator to standard interface.
Definition: StanTypedefs.h:28
const FitSinSq2Theta23 kFitSinSq2Theta23
string fnameSwap
Definition: demo4.py:6
void Simulation()
void SetL(double L) override
Definition: OscCalcPMNS.h:33
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
void template_basic()
virtual Spectrum PredictUnoscillated() const
Definition: IPrediction.cxx:33
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
string filename
Definition: shutoffs.py:106
const Cut kNumuFD
Definition: NumuCuts.h:53
osc::OscCalcDumb calc
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
const Binning kNumuEnergyBinning
Definition: Binnings.cxx:13
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:341
void SetdCP(const T &dCP) override
Definition: OscCalcPMNS.h:40
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcPMNS.h:36
const Var kCCE
Definition: NumuVars.h:21
#define pot
virtual void Go() override
Load all the registered spectra.
void SetTh13(const T &th13) override
Definition: OscCalcPMNS.h:38
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
loader
Definition: demo0.py:10
TLatex * prelim
Definition: Xsec_final.C:133
void SetRho(double rho) override
Definition: OscCalcPMNS.h:34
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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.
Prediction that just uses FD MC, with no extrapolation.
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
#define mcset
Compare a single data spectrum to the MC + cosmics expectation.