test_nueextrapsyst.C
Go to the documentation of this file.
1 // Macro to make nominal and shifted FD nue predictions with NueExtrapSysts
2 // Result is four plots of shifted FD nue predictions and ratios to nominal
3 
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Cuts/NueCutsFirstAna.h"
8 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Vars/Vars.h"
12 #include "CAFAna/Analysis/Style.h"
13 #include "CAFAna/Systs/Systs.h"
15 #include "CAFAna/Core/HistAxis.h"
17 #include "CAFAna/Analysis/Plots.h"
18 #include "CAFAna/Analysis/Style.h"
19 #include "CAFAna/Core/Cut.h"
21 #include "CAFAna/Cuts/NueCutsSecondAna.h"
26 #include "CAFAna/Vars/Vars.h"
29 #include "OscLib/OscCalcPMNSOpt.h"
30 #include "CAFAna/Analysis/Calcs.h"
35 #include "CAFAna/Vars/HistAxes.h"
36 
37 
38 #include "TCanvas.h"
39 #include "TH1.h"
40 #include "TFile.h"
41 #include "TLegend.h"
42 #include "TLine.h"
43 #include "TLatex.h"
44 #include "TGraph.h"
45 
46 #include <iostream>
47 #include <iomanip>
48 
49 using namespace ana;
50 
51 void test_nueextrapsyst(const std::string& outfilename = "test_extrapsyst.root")
52 {
53  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;
54 
55  //Set +- 1 sigma NueBkg systematic shift
56  SystShifts bkgPlus(&kNueExtrapSystBkg2017, 1.);
57  SystShifts bkgMinus(&kNueExtrapSystBkg2017, -1.);
58 
59  //Set +- 1 sigma NueSig systematic shift
62 
63 
64  // ND loaders
65  SpectrumLoader loaderNDData("prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_period1_nue_or_numu_or_nus_contain_v1_goodruns");
66  SpectrumLoader loaderNDMC("prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_period1_nue_or_numu_or_nus_contain_v1");
67 
68 
69  // FD loaders
70  SpectrumLoader loaderFD("prod_decaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_period1_nue_or_numu_or_nus_contain_v1");
71  SpectrumLoader loaderFDSwap("prod_decaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_period1_nue_or_numu_or_nus_contain_v1");
72  SpectrumLoader loaderFDTau("prod_decaf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1");
73 
74  HistAxis axisNue("E_{reco} (GeV)", 10, 0, 5, kNueEnergy2017);
75 
77  ResetOscCalcToDefault(&calc);
78 
79  calc.SetDmsq32(+fabs(calc.GetDmsq32()));
80  calc.SetdCP(3*TMath::Pi()/2);
81  calc.SetTh23(asin(sqrt(0.404)));
82 
83  //FD bkg nominal and shifted prediction
84  PredictionNoExtrap predNom(loaderFD, loaderFDSwap, "E_{reco} (GeV)", Binning::Simple(10, 0, 5), kNueEnergy2017, kNue2017FDAllSamples, kNoShift);
85  PredictionNoExtrap predUp(loaderFD, loaderFDSwap, "E_{reco} (GeV)", Binning::Simple(10, 0, 5), kNueEnergy2017, kNue2017FDAllSamples, bkgPlus);
86  PredictionNoExtrap predDn(loaderFD, loaderFDSwap, "E_{reco} (GeV)", Binning::Simple(10, 0, 5), kNueEnergy2017, kNue2017FDAllSamples, bkgMinus);
87 
88  //numu decomps ND
89  NumuDecomp numuDecompNom(loaderNDMC, loaderNDData, axisNumu, kNumuNDCvn, kNoShift, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
90  NumuDecomp numuDecompUp(loaderNDMC, loaderNDData, axisNumu, kNumuNDCvn, sigPlus, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
91  NumuDecomp numuDecompDn(loaderNDMC, loaderNDData, axisNumu, kNumuNDCvn, sigMinus, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
92 
93  //prop decomp
94  ProportionalDecomp propDecomp(loaderNDMC, loaderNDData, axisNue, kNue2017NDPresel && kNueSecondAnaCVNeSsb);
95 
96  //extrap
97  ModularExtrap* extrapPropNom;
98  ModularExtrap* extrapPropUp;
99  ModularExtrap* extrapPropDn;
100  extrapPropNom = new NueExtrap(loaderNDMC, loaderFDSwap, loaderFD, loaderFDTau, propDecomp, numuDecompNom, axisNue, axisNumu, kNue2017FDAllSamples, kNue2017NDPresel && kNueSecondAnaCVNeSsb, kNumuNDCvn, kNoShift);
101  extrapPropUp = new NueExtrap(loaderNDMC, loaderFDSwap, loaderFD, loaderFDTau, propDecomp, numuDecompUp, axisNue, axisNumu, kNue2017FDAllSamples, kNue2017NDPresel && kNueSecondAnaCVNeSsb, kNumuNDCvn, sigPlus);
102  extrapPropDn = new NueExtrap(loaderNDMC, loaderFDSwap, loaderFD, loaderFDTau, propDecomp, numuDecompDn, axisNue, axisNumu, kNue2017FDAllSamples, kNue2017NDPresel && kNueSecondAnaCVNeSsb, kNumuNDCvn, sigMinus);
103 
104  //predict
105  PredictionExtrap predicPropNom(extrapPropNom);
106  PredictionExtrap predicPropUp(extrapPropUp);
107  PredictionExtrap predicPropDn(extrapPropDn);
108 
109  //loaders GO!
110  loaderNDData.Go();
111  loaderNDMC.Go();
112  loaderFD.Go();
113  loaderFDSwap.Go();
114  loaderFDTau.Go();
115 
116  // Draw the plot
117  TH1* hNom = predNom.Predict(&calc).ToTH1(9e20);
118  hNom->Draw("hist");
119  TH1* hUp = predUp.Predict(&calc).ToTH1(9e20);
120  hUp->SetLineColor(kRed);
121  hUp->Draw("hist same");
122  TH1* hDn = predDn.Predict(&calc).ToTH1(9e20);
123  hDn->SetLineColor(kBlue);
124  hDn->Draw("hist same");
125 
126  new TCanvas;
127  TH1* rUp = (TH1*)hUp->Clone();
128  rUp->Divide(hNom);
129  rUp->Draw("hist ][");
130  rUp->GetYaxis()->SetTitle("Ratio to nominal");
131  rUp->GetYaxis()->SetRangeUser(.8, 1.2);
132  TH1* rDn = (TH1*)hDn->Clone();
133  rDn->Divide(hNom);
134  rDn->Draw("hist same ][");
135  TGraph* one = new TGraph;
136  one->SetPoint(0, 0, 1);
137  one->SetPoint(1, 100, 1);
138  one->SetLineWidth(2);
139  one->SetLineStyle(7);
140  one->Draw("same");
141 
142  new TCanvas;
143  TH1* hNomEx = predicPropNom.Predict(&calc).ToTH1(9e20);
144  hNomEx->Draw("hist");
145  TH1* hUpEx = predicPropUp.Predict(&calc).ToTH1(9e20);
146  hUpEx->SetLineColor(kRed);
147  hUpEx->Draw("hist same");
148  TH1* hDnEx = predicPropDn.Predict(&calc).ToTH1(9e20);
149  hDnEx->SetLineColor(kBlue);
150  hDnEx->Draw("hist same");
151 
152  new TCanvas;
153  TH1* rUpEx = (TH1*)hUpEx->Clone();
154  rUpEx->Divide(hNomEx);
155  rUpEx->Draw("hist ][");
156  rUpEx->GetYaxis()->SetTitle("Ratio to nominal");
157  rUpEx->GetYaxis()->SetRangeUser(.8, 1.2);
158  TH1* rDnEx = (TH1*)hDnEx->Clone();
159  rDnEx->Divide(hNomEx);
160  rDnEx->Draw("hist same ][");
161  TGraph* oneEx = new TGraph;
162  oneEx->SetPoint(0, 0, 1);
163  oneEx->SetPoint(1, 100, 1);
164  oneEx->SetLineWidth(2);
165  oneEx->SetLineStyle(7);
166  oneEx->Draw("same");
167 } // end
void test_nueextrapsyst(const std::string &outfilename="test_extrapsyst.root")
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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:149
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
osc::OscCalcDumb calc
string outfilename
knobs that need extra care
Uses MC for NC and CC components, assigns remainder of data to CC.
Definition: NumuDecomp.h:10
const Var kNueEnergy2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2017FDFit(kCVNemE(sr), kCVNhadE(sr));})
Definition: NueEnergy2017.h:11
const Cut kNumuNDCvn
Definition: NumuCuts.h:62
Spectrum Predict(osc::IOscCalc *calc) const override
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:32
virtual void Go() override
Load all the registered spectra.
const NueExtrapSystSignalKin2017 kNueExtrapSystSignalKin2017
void SetTh23(const T &th23) override
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
const NueExtrapSystBkg2017 kNueExtrapSystBkg2017
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const HistAxis axisNue
Definition: syst_header.h:378
void SetDmsq32(const T &dmsq32) override
const Cut kNue2017NDPresel
Definition: NueCuts2017.h:285
void NueExtrap(string beam="fhc", string cvntype="oldpresel")
Definition: NueExtrap.C:28
Splits Data proportionally according to MC.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SetdCP(const T &dCP) override
Prediction that just uses FD MC, with no extrapolation.
const HistAxis axisNumu
Take the output of an extrapolation and oscillate it as required.
auto one()
Definition: PMNS.cxx:57
const Var kXSecCVWgt2017
Definition: XsecTunes.h:36
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
T asin(T number)
Definition: d0nt_math.hpp:60
const Cut kNue2017FDAllSamples
Our FD selection including all samples, for making predictions, etc.
Definition: NueCuts2017.h:155
enum BeamMode string