template_GENIE_systs.C
Go to the documentation of this file.
1 // FHC only, GENIE systematics
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"
17 #include "CAFAna/Core/Utilities.h"
18 #include "CAFAna/Vars/Vars.h"
20 #include "CAFAna/Core/Loaders.h"
22 #include "CAFAna/Core/Binning.h"
23 #include "CAFAna/Systs/Systs.h"
24 
25 #include "TCanvas.h"
26 #include "TFile.h"
27 #include "TGraph.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TMath.h"
31 #include "TGaxis.h"
32 #include "TMultiGraph.h"
33 #include "TLegend.h"
34 #include "TLatex.h"
35 #include "TStyle.h"
36 #include "THStack.h"
37 #include "TPaveText.h"
38 
39 #include <cmath>
40 #include <iostream>
41 #include <sstream>
42 #include <string>
43 #include <fstream>
44 #include <iomanip>
45 
46 #include "Utilities/func/MathUtil.h"
47 
48 // Use with Nova style ROOT logon script
49 
50 #define pot 18e20
51 
52 using namespace ana;
53 
54 void Simulation() // for some reason doesn't always remember definition in Style script
55 {
56  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Simulation");
57  prelim->SetTextColor(kGray+1);
58  prelim->SetNDC();
59  prelim->SetTextSize(2/30.);
60  prelim->SetTextAlign(32);
61  prelim->Draw();
62 }
63 
65 {
66 
67  TGaxis::SetMaxDigits(2);
68 
69  gStyle->SetTitleOffset(.85, "y");
70  gStyle->SetTitleOffset(.84, "x");
71  gStyle->SetTitleSize(0.05,"x"); // these are needed to make the X axis label fit, as it has super and subscripts and tends to cramp or get cut off
72 
73  TLatex text; // I'll write some text on plots later
74  text.SetNDC(true);
75  text.SetTextSize(0.04);
76 
77  TH1F *hblack = new TH1F("black", "black", 10, 0, 1); // for adding to legends later
78  TH1F *hred = new TH1F("red", "red", 10, 0, 1);
79  TH1F *hblue = new TH1F("blue", "blue", 10, 0, 1);
80  TH1F *hmagenta = new TH1F("magenta", "magenta", 10, 0, 1);
81  TH1F *hgreen = new TH1F("green", "green", 10, 0, 1);
82 
83  hblack->SetLineWidth(2); // don't count on NOvA style script being loaded
84  hred->SetLineWidth(2);
85  hblue->SetLineWidth(2);
86  hmagenta->SetLineWidth(2);
87  hgreen->SetLineWidth(2);
88  hblack->SetLineColor(1);
89  hred->SetLineColor(2);
90  hblue->SetLineColor(4);
91  hmagenta->SetLineColor(kMagenta);
92  hgreen->SetLineColor(kGreen);
93 
94  // Caltech
95  // const std::string fname = "$NOVA_DATA/mc/forFA/fd_genie_fhc_nonswap.caf.root";
96  // const std::string fnameSwap = "$NOVA_DATA/mc/forFA/fd_genie_fhc_swap.caf.root";
97  // const std::string fnameCosmic = "/nfs/raid5/kirk/mc/fd_cry/fd_cry.153.caf.root";
98 
99  // FNAL
100  const std::string fname = "prod_caf_FA14-10-28_fd_genie_fhc_nonswap";
101  const std::string fnameSwap = "prod_caf_FA14-10-28_fd_genie_fhc_fluxswap";
102  const std::string fnameCosmic = "prod_caf_FA14-10-28_fd_cry_all"; // this may take a long time to load, it is a lot of files. copy/hadd a subset for faster processing
103 
104  osc::OscCalcPMNS calc; // our standard calculator we'll use
105 
106  calc.SetL(810);
107  calc.SetRho(0); // No matter effects
108  calc.SetDmsq21(7.59e-5);
109  calc.SetDmsq32(2.4e-3);
110  calc.SetTh12(.601);
111  calc.SetTh13(.1567);
112  calc.SetdCP(0);
113  calc.SetTh23(TMath::Pi()/4); // sin^2(2theta) = 1
114 
115  osc::OscCalcPMNS calc2; // this is just used once in the prediction, don't touch it
116 
117  calc2.SetL(810);
118  calc2.SetRho(0); // No matter effects
119  calc2.SetDmsq21(7.59e-5);
120  calc2.SetDmsq32(2.4e-3);
121  calc2.SetTh12(.601);
122  calc2.SetTh13(.1567);
123  calc2.SetdCP(0);
124  calc2.SetTh23(TMath::Pi()/4); // sin^2(2theta) = 1
125 
130 
131  GenieRwgtTableSyst kMaCCQE(rwgt::fReweightMaCCQE); // defining these here makes it a lot simpler every time we use them
132  GenieRwgtTableSyst kMaCCRES(rwgt::fReweightMaCCRES);
133 
134  // Load up the necessary histograms, seperated by flavour etc. First vector is list of systematics so it will know how to use them later.
135  NoExtrapGenerator predGenNonQE(
136  HistAxis("Non-QE Neutrino Energy (GeV)", kNumuEnergyBinning,kCCE),
138  PredictionInterp predCC({&kMaCCQE,&kMaCCRES}, &calc2, predGenNonQE, loaders);
139 
140  NoExtrapGenerator predGenQE(
141  HistAxis("Non-QE Neutrino Energy (GeV)", kNumuEnergyBinning,kQEE),
143  PredictionInterp predQE({&kMaCCQE, &kMaCCRES}, &calc2, predGenQE, loaders);
144 
145  loaders.Go();
146 
147  Spectrum obsCC = predCC.Predict(&calc); // the prediction
148  Spectrum obsQE = predQE.Predict(&calc);
149  Spectrum fakeCC = obsCC.FakeData(pot); // the data, faked
150  Spectrum fakeQE = obsQE.FakeData(pot);
151  SingleSampleExperiment exptCC(&predCC, fakeCC); // made into an 'experiment'
152  SingleSampleExperiment exptQE(&predQE, fakeQE);
153 
154  std::vector<const IExperiment*> exptsAll;
155  exptsAll.push_back(&exptCC);
156  exptsAll.push_back(&exptQE);
157  MultiExperiment exptAll(exptsAll);
158 
159  calc.SetTh23(1); // Seed value
160 
161  // let's look at some basic spectra and see if the shifts do anything
162 
163  new TCanvas("QEspectra","QEspectra") ;
164 
165  TH1* hNom = predQE.Predict(&calc).ToTH1(pot);
166  TH1* hQEMa = predQE.PredictSyst(&calc, SystShifts(&kMaCCQE, +2)).ToTH1(pot);
167  TH1* hRESMa = predQE.PredictSyst(&calc, SystShifts(&kMaCCRES, +2)).ToTH1(pot);
168 
169  hNom->SetLineWidth(2);
170  hNom->SetLineColor(4);
171  hQEMa->SetLineWidth(2);
172  hQEMa->SetLineColor(2);
173  hRESMa->SetLineWidth(2);
174  hRESMa->SetLineColor(kGreen);
175 
176  hNom->GetXaxis()->SetTitle("Reconstructed QE Energy (GeV)");
177  hNom->GetYaxis()->SetTitle("Events / 18e20 POT");
178 
179  hNom->Draw();
180  hQEMa->Draw("same");
181  hRESMa->Draw("same");
182 
183  TLegend *legspec = new TLegend(0.6,0.5,0.89,0.79);
184  legspec->SetFillColor(10);
185  legspec->SetFillStyle(0);
186  legspec->SetBorderSize(0);
187  legspec->SetHeader("18e20 POT, FHC, 14 kton, QE");
188 
189  legspec->AddEntry(hNom, "No systematics","L");
190  legspec->AddEntry(hQEMa,"CC QE Ma +2 sigma","L");
191  legspec->AddEntry(hRESMa,"CC RES Ma +2 sigma","L");
192  legspec->Draw("same");
193  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
194 
195  gPad->Print("plots/spec.qema.ps");
196  gPad->Print("plots/spec.qema.C");
197 
198  new TCanvas("CCspectra","CCspectra") ;
199 
200  hNom = predCC.Predict(&calc).ToTH1(pot);
201  hQEMa = predCC.PredictSyst(&calc, SystShifts(&kMaCCQE, +2)).ToTH1(pot);
202  hQEMa->SetLineColor(2);
203  hRESMa = predCC.PredictSyst(&calc, SystShifts(&kMaCCRES, +2)).ToTH1(pot);
204  hRESMa->SetLineColor(kGreen);
205 
206  hNom->GetXaxis()->SetTitle("Reconstructed CC Energy (GeV)");
207 
208  hNom->Draw();
209  hQEMa->Draw("same");
210  hRESMa->Draw("same");
211 
212  legspec->SetHeader("18e20 POT, FHC, 14 kton, CC");
213 
214  legspec->Draw("same");
215  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
216 
217  gPad->Print("plots/spec.ccma.ps");
218  gPad->Print("plots/spec.ccma.C");
219 
220  //The log-likelihood surface: FHC, no systs
221 
222  new TCanvas("contoursfhc1","contoursfhc1") ;
223 
224  FrequentistSurface surfCC(&exptCC, &calc,
225  &kFitSinSq2Theta23, 30, .9, 1,
226  &kFitDmSq32Scaled, 30, 2, 2.8);
227  surfCC.DrawBestFit(kRed);
228  surfCC.DrawContour(Gaussian90Percent2D(surfCC), kSolid, kRed);
229 
230  FrequentistSurface surfQE(&exptQE, &calc,
231  &kFitSinSq2Theta23, 30, .93, 1,
232  &kFitDmSq32Scaled, 30, 2.2, 2.7);
233  surfQE.DrawBestFit(kBlue);
234  surfQE.DrawContour(Gaussian90Percent2D(surfQE), kSolid, kBlue);
235 
236  FrequentistSurface surfAll(&exptAll, &calc,
237  &kFitSinSq2Theta23, 30, .96, 1,
238  &kFitDmSq32Scaled, 30, 2.2, 2.6);
239  surfAll.DrawBestFit(kBlack);
240  surfAll.DrawContour(Gaussian90Percent2D(surfAll), kSolid, kBlack);
241 
242  Simulation();
243 
244  TLegend *leg = new TLegend(0.11,0.11,0.5,0.45);
245  leg->SetFillColor(10);
246  leg->SetFillStyle(0);
247  leg->SetBorderSize(0);
248  leg->SetHeader("18e20 POT, FHC, 14 kton: 90% CL, no Systs");
249 
250  leg->AddEntry(hred, "Contained non-QE CC","L");
251  leg->AddEntry(hblue,"Contained QE CC","L");
252  leg->AddEntry(hblack,"Combined","L");
253  leg->Draw("same");
254  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
255 
256  gPad->Print("plots/withcosmics.1.ps");
257  gPad->Print("plots/withcosmics.1.C");
258 
259  //The log-likelihood surface: FHC, MaQE
260  new TCanvas("contoursfhcMaQE","contoursfhcMaQE") ;
261 
262  FrequentistSurface surfCCqe(&exptCC, &calc,
263  &kFitSinSq2Theta23, 30, .9, 1,
264  &kFitDmSq32Scaled, 30, 2, 2.8,{},{&kMaCCQE});
265  surfCCqe.DrawBestFit(kRed);
266  surfCCqe.DrawContour(Gaussian90Percent2D(surfCCqe), kSolid, kRed);
267 
268  FrequentistSurface surfQEqe(&exptQE, &calc,
269  &kFitSinSq2Theta23, 30, .92, 1,
270  &kFitDmSq32Scaled, 30, 2.2, 2.7,{},{&kMaCCQE});
271  surfQEqe.DrawBestFit(kBlue);
272  surfQEqe.DrawContour(Gaussian90Percent2D(surfQEqe), kSolid, kBlue);
273 
274  FrequentistSurface surfAllqe(&exptAll, &calc,
275  &kFitSinSq2Theta23, 30, .96, 1,
276  &kFitDmSq32Scaled, 30, 2.2, 2.6,{},{&kMaCCQE});
277  surfAllqe.DrawBestFit(kBlack);
278  surfAllqe.DrawContour(Gaussian90Percent2D(surfAllqe), kSolid, kBlack);
279 
280  Simulation();
281 
282  leg->SetHeader("18e20 POT, FHC, 14 kton: 90% CL, MaQE sys");
283 
284  leg->Draw("same");
285  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
286 
287  gPad->Print("plots/maqe.ps");
288  gPad->Print("plots/maqe.C");
289 
290  //The log-likelihood surface: FHC, MaRes
291  new TCanvas("contoursfhcMaRES","contoursfhcMaRES") ;
292 
293  FrequentistSurface surfCCres(&exptCC, &calc,
294  &kFitSinSq2Theta23, 30, .9, 1,
295  &kFitDmSq32Scaled, 30, 2, 2.8,{},{&kMaCCRES});
296  surfCCres.DrawBestFit(kRed);
297  surfCCres.DrawContour(Gaussian90Percent2D(surfCCres), kSolid, kRed);
298 
299  FrequentistSurface surfQEres(&exptQE, &calc,
300  &kFitSinSq2Theta23, 30, .92, 1,
301  &kFitDmSq32Scaled, 30, 2.2, 2.7,{},{&kMaCCRES});
302  surfQEres.DrawBestFit(kBlue);
303  surfQEres.DrawContour(Gaussian90Percent2D(surfQEqe), kSolid, kBlue);
304 
305  FrequentistSurface surfAllres(&exptAll, &calc,
306  &kFitSinSq2Theta23, 30, .96, 1,
307  &kFitDmSq32Scaled, 30, 2.2, 2.6,{},{&kMaCCRES});
308  surfAllres.DrawBestFit(kBlack);
309  surfAllres.DrawContour(Gaussian90Percent2D(surfAllqe), kSolid, kBlack);
310 
311  Simulation();
312 
313  leg->SetHeader("18e20 POT, FHC, 14 kton: 90% CL, MaRes sys");
314 
315  leg->Draw("same");
316  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
317 
318  gPad->Print("plots/mares.ps");
319  gPad->Print("plots/mares.C");
320 
321  //The log-likelihood surface: FHC, MaQE and MaRes
322  new TCanvas("contoursfhcMaBoth","contoursfhcMaBoth") ;
323 
324  FrequentistSurface surfCCboth(&exptCC, &calc,
325  &kFitSinSq2Theta23, 30, .9, 1,
326  &kFitDmSq32Scaled, 30, 2, 2.8,{},{&kMaCCQE,&kMaCCRES});
327  surfCCboth.DrawBestFit(kRed);
328  surfCCboth.DrawContour(Gaussian90Percent2D(surfCCboth), kSolid, kRed);
329 
330  FrequentistSurface surfQEboth(&exptQE, &calc,
331  &kFitSinSq2Theta23, 30, .91, 1,
332  &kFitDmSq32Scaled, 30, 2.2, 2.7,{},{&kMaCCQE,&kMaCCRES});
333  surfQEboth.DrawBestFit(kBlue);
334  surfQEboth.DrawContour(Gaussian90Percent2D(surfQEboth), kSolid, kBlue);
335 
336  FrequentistSurface surfAllboth(&exptAll, &calc,
337  &kFitSinSq2Theta23, 30, .96, 1,
338  &kFitDmSq32Scaled, 30, 2.2, 2.6,{},{&kMaCCQE,&kMaCCRES});
339  surfAllboth.DrawBestFit(kBlack);
340  surfAllboth.DrawContour(Gaussian90Percent2D(surfAllboth), kSolid, kBlack);
341 
342  Simulation();
343 
344  leg->SetHeader("18e20 POT, FHC, 14 kton: 90% CL, MaQE and MaRes systs");
345 
346  leg->Draw("same");
347  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");
348 
349  gPad->Print("plots/maqeres.ps");
350  gPad->Print("plots/maqeres.C");
351 
352 }
Implements systematic errors by interpolation between shifted templates.
enum BeamMode kRed
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Far Detector at Ash River.
Definition: SREnums.h:11
void SetDmsq21(const T &dmsq21) override
Definition: OscCalcPMNS.h:35
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh12(const T &th12) override
Definition: OscCalcPMNS.h:37
void SetTh23(const T &th23) override
Definition: OscCalcPMNS.h:39
Adapt the PMNS calculator to standard interface.
Definition: StanTypedefs.h:28
const FitSinSq2Theta23 kFitSinSq2Theta23
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
string fnameSwap
Definition: demo4.py:6
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
void SetL(double L) override
Definition: OscCalcPMNS.h:33
#define pot
const Cut kNumuCosmicRej([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.535 && sr->slc.nhit< 400);})
Definition: NumuCuts.h:32
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
osc::OscCalcDumb calc
Generates FD-only predictions (no extrapolation)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void template_GENIE_systs()
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 FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
const Cut kNumuQESel([](const caf::SRProxy *sr){std::cout<< "WARNING! Attempting to access kNumuQESel which relies on qepid (which no longer exists.) Aborting..."<< std::endl;abort(); return false;})
Definition: NumuCuts.h:38
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
static const unsigned char kMaCCQE
Definition: VarVals.h:70
const Var kCCE
Definition: NumuVars.h:21
void SetTh13(const T &th13) override
Definition: OscCalcPMNS.h:38
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
TLatex * prelim
Definition: Xsec_final.C:133
void SetRho(double rho) override
Definition: OscCalcPMNS.h:34
Combine multiple component experiments.
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void Simulation()
const Var kQEE
Energy estimator for quasielastic CC events.
Definition: NumuVars.cxx:28
Float_t e
Definition: plot.C:35
enum BeamMode kGreen
enum BeamMode kBlue
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
Compare a single data spectrum to the MC + cosmics expectation.
enum BeamMode string