test_genie_systs.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/HistAxis.h"
4 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Vars/Vars.h"
8 
9 #include "CAFAna/Systs/Systs.h"
10 
11 #include "TCanvas.h"
12 #include "TGraph.h"
13 #include "TH1.h"
14 
15 #include <cmath>
16 #include <iostream>
17 
18 #include "Utilities/func/MathUtil.h"
19 
20 using namespace ana;
21 
23 {
24  // NB: this is only about 1% of the available stats
25  SpectrumLoader loader("$NOVA_PROD/mc/S14-05-12/fardet/genie/prod_pid_caf_loosened_presel_neutrino2014nightmare_besttoavoid/fardet_genie_fhc_nonswap_none_3000_r0100000?_s01*.caf.root");
26  SpectrumLoader loaderSwap("$NOVA_PROD/mc/S14-05-12/fardet/genie/prod_pid_caf_loosened_presel_neutrino2014nightmare_besttoavoid/fardet_genie_fhc_swap_none_3000_r0100000?_s01*.caf.root");
27 
28  // Make a calculator. Any of the variants is fine
30 
31  // Standard oscillation parameters
32  calc.SetL(810);
33  calc.SetRho(2.75);
34  calc.SetDmsq21(7.6e-5);
35  calc.SetDmsq32(2.35e-3);
36  calc.SetTh12(asin(sqrt(.87))/2);
37  calc.SetTh13(asin(sqrt(.10))/2);
38  calc.SetTh23(TMath::Pi()/4);
39  calc.SetdCP(0);
40 
41  HistAxis axis("Calorimetric Energy (GeV)", 20, 0, 5, kCaloE);
42 
43  // Nominal prediction
44  PredictionNoExtrap predNom(loader, loaderSwap, axis.label, axis.bins, axis.var, kLEM > 0.8, kNoShift);
45 
46  // M_A QE +1 sigma
47  PredictionNoExtrap predUp(loader, loaderSwap, axis.label, axis.bins, axis.var, kLEM > 0.8, SystShifts(GetGenieKnobSyst(rwgt::fReweightMaCCQE), +1));
48 
49  // M_A QE -1 sigma
50  PredictionNoExtrap predDn(loader, loaderSwap, axis.label, axis.bins, axis.var, kLEM > 0.8, SystShifts(GetGenieKnobSyst(rwgt::fReweightMaCCQE), -1));
51 
52 
53  // Load the files
54  loader.Go();
55  loaderSwap.Go();
56 
57  // Draw the plot
58 
59  TH1* hNom = predNom.Predict(&calc).ToTH1(18e20);
60  hNom->Draw("hist");
61 
62  TH1* hUp = predUp.Predict(&calc).ToTH1(18e20);
63  hUp->SetLineColor(kRed);
64  hUp->Draw("hist same");
65 
66  TH1* hDn = predDn.Predict(&calc).ToTH1(18e20);
67  hDn->SetLineColor(kBlue);
68  hDn->Draw("hist same");
69 
70  new TCanvas;
71  TH1* rUp = (TH1*)hUp->Clone();
72  rUp->Divide(hNom);
73  rUp->Draw("hist ][");
74  rUp->GetYaxis()->SetTitle("Ratio to nominal");
75  rUp->GetYaxis()->SetRangeUser(.8, 1.2);
76 
77  TH1* rDn = (TH1*)hDn->Clone();
78  rDn->Divide(hNom);
79  rDn->Draw("hist same ][");
80 
81  TGraph* one = new TGraph;
82  one->SetPoint(0, 0, 1);
83  one->SetPoint(1, 100, 1);
84  one->SetLineWidth(2);
85  one->SetLineStyle(7);
86  one->Draw("same");
87 }
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
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:176
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
osc::OscCalcDumb calc
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Spectrum Predict(osc::IOscCalc *calc) const override
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:32
virtual void Go() override
Load all the registered spectra.
void SetTh23(const T &th23) override
loader
Definition: demo0.py:10
void SetDmsq21(const T &dmsq21) override
const SystShifts kNoShift
Definition: SystShifts.cxx:22
void SetDmsq32(const T &dmsq32) override
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SetdCP(const T &dCP) override
void SetTh12(const T &th12) override
Prediction that just uses FD MC, with no extrapolation.
auto one()
Definition: PMNS.cxx:57
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
enum BeamMode kBlue
void test_genie_systs()
void SetRho(double rho) override
T asin(T number)
Definition: d0nt_math.hpp:60