syst_test.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Vars/FitVars.h"
7 #include "CAFAna/Core/Progress.h"
11 #include "CAFAna/Systs/Systs.h"
12 
14 
15 #include "OscLib/OscCalcPMNSOpt.h"
16 
17 #include "TCanvas.h"
18 #include "TFile.h"
19 #include "TGraph.h"
20 #include "TH2.h"
21 #include "TLegend.h"
22 #include "TMarker.h"
23 #include "TPad.h"
24 
25 using namespace ana;
26 
28 {
29  calc->SetL(810);
30  calc->SetRho(2.75);
31  calc->SetDmsq21(7.59e-5);
32  calc->SetDmsq32(2.40e-3);
33  calc->SetTh12(asin(sqrt(.87))/2); // PDG
34  calc->SetTh13(asin(sqrt(.095))/2);
35  calc->SetTh23(asin(sqrt(.99))/2);
36  calc->SetdCP(0);
37 }
38 
39 void syst_test()
40 {
42  loaders.SetLoaderPath("defname: prod_caf_S15-05-22_fd_genie_fhc_nonswap with stride 100", caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap);
43  loaders.SetLoaderPath("defname: prod_caf_S15-05-22_fd_genie_fhc_fluxswap with stride 100", caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
44 
45  // This is the oscillation parameters that PredictionInterp will expand
46  // around
47  osc::OscCalcPMNSOpt calcOrigin;
48  resetCalc(&calcOrigin);
49 
50  // Pass the list of systematics we want to be able to handle
51  NoExtrapGenerator predGen(
52  HistAxis("Simple Energy (GeV)", kNumuEnergyBinning, kCCE),
53  kNumuCC && /*kNumuContainFD &&*/ !kNumuQESel );
54  PredictionInterp interp({&kEnergyScaleSyst, &kNormSyst, &kNCScaleSyst},
55  &calcOrigin,
56  predGen,
57  loaders);
58 
59  loaders.Go();
60 
61 
63 
64  resetCalc(&calc);
65  // Unshifted fake data
66  Spectrum data = interp.Predict(&calc).FakeData(18e20);
67  SingleSampleExperiment expt(&interp, data);
68 
69  FrequentistSurface surf(&expt, &calc,
70  &kFitSinSq2Theta23, 20, .94, 1,
71  &kFitDmSq32, 20, 2.2e-3, 2.6e-3);
72  surf.DrawContour(Gaussian90Percent2D(surf), kSolid, kBlack);
73 
74  // Star plots
75  const std::vector<const ISyst*> systs = {&kNormSyst, &kNCScaleSyst, &kEnergyScaleSyst};
76  for(const ISyst* s: systs){
77  Progress prog("Making star plot for "+s->ShortName());
78  TGraph* g = new TGraph;
79  for(double sigma = -1; sigma < +1.05; sigma += .2){
81  systs.SetShift(s, sigma);
82  resetCalc(&calc);
83  // Systematically shifted fake data to fit to
84  Spectrum dataSyst = interp.PredictSyst(&calc, systs).FakeData(18e20);
85 
86  SingleSampleExperiment exptSyst(&interp, dataSyst);
87  MinuitFitter fitter(&exptSyst, {&kFitSinSq2Theta23, &kFitDmSq32});
88  fitter.Fit(&calc, IFitter::kQuiet);
89 
90  g->SetPoint(g->GetN(),
92  kFitDmSq32.GetValue(&calc));
93 
94  prog.SetProgress((sigma+1)/2.);
95  }
96  if(s == &kNormSyst) g->SetLineColor(kMagenta);
97  if(s == &kNCScaleSyst) g->SetLineColor(kGreen+2);
98  if(s == &kEnergyScaleSyst) g->SetLineColor(kBlue);
99  g->SetLineWidth(3);
100  g->Draw("l");
101  }
102 
103  surf.DrawBestFit(kBlack);
104 
105  // For each systematic parameter, make a contour marginalizing over it
106  FrequentistSurface surfNorm(&expt, &calc,
107  &kFitSinSq2Theta23, 20, .94, 1,
108  &kFitDmSq32, 20, 2.2e-3, 2.6e-3,
109  {},
110  {&kNormSyst});
111  surfNorm.DrawContour(Gaussian90Percent2D(surfNorm), kSolid, kMagenta);
112 
113  FrequentistSurface surfNC(&expt, &calc,
114  &kFitSinSq2Theta23, 20, .94, 1,
115  &kFitDmSq32, 20, 2.2e-3, 2.6e-3,
116  {},
117  {&kNCScaleSyst});
118  surfNC.DrawContour(Gaussian90Percent2D(surfNC), kSolid, kGreen+2);
119 
120  FrequentistSurface surfCCE(&expt, &calc,
121  &kFitSinSq2Theta23, 20, .94, 1,
122  &kFitDmSq32, 20, 2.2e-3, 2.6e-3,
123  {},
124  {&kEnergyScaleSyst});
125  surfCCE.DrawContour(Gaussian90Percent2D(surfCCE), kSolid, kBlue);
126 
127  // And one final contour marginalizing over all three
128  FrequentistSurface surfAll(&expt, &calc,
129  &kFitSinSq2Theta23, 20, .94, 1,
130  &kFitDmSq32, 20, 2.2e-3, 2.6e-3,
131  {},
132  {&kNormSyst, &kNCScaleSyst, &kEnergyScaleSyst});
133  surfAll.DrawContour(Gaussian90Percent2D(surfAll), kSolid, kRed);
134 
135  TLegend* leg = new TLegend(.125, .125, .475, .475);
136  leg->SetFillStyle(0);
137  TH1* h = new TH1F("", "", 1, 0, 1);
138  h->SetLineColor(kBlack);
139  leg->AddEntry(h->Clone(), "Stat only", "l");
140  h->SetLineColor(kMagenta);
141  leg->AddEntry(h->Clone(), "Normalization", "l");
142  h->SetLineColor(kGreen+2);
143  leg->AddEntry(h->Clone(), "NC bkg", "l");
144  h->SetLineColor(kBlue);
145  leg->AddEntry(h->Clone(), "Energy scale", "l");
146  h->SetLineColor(kRed);
147  leg->AddEntry(h->Clone(), "All systs", "l");
148  leg->Draw();
149 
150  gPad->Print("plots/cc_systs_cont.eps");
151 
152 
153  int i = 0;
154  for(TH2* m: surfAll.GetMarginalizedHists()){
155  new TCanvas;
156  m->Draw("colz");
157  const double range = std::max(fabs(m->GetMinimum()), fabs(m->GetMaximum()));
158  m->SetMinimum(-range);
159  m->SetMaximum(+range);
160  gPad->Print(TString::Format("plots/cc_systs_marg_%d.eps", i++));
161  }
162 }
T max(const caf::Proxy< T > &a, T b)
virtual void SetL(double L)=0
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
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
virtual void SetDmsq21(const T &dmsq21)=0
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
const FitSinSq2Theta23 kFitSinSq2Theta23
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Cut kNumuCC
Definition: NumuCuts.h:70
virtual void SetTh13(const T &th13)=0
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Generates FD-only predictions (no extrapolation)
virtual void SetDmsq32(const T &dmsq32)=0
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
const XML_Char const XML_Char * data
Definition: expat.h:268
Log-likelihood scan across two parameters.
const XML_Char * s
Definition: expat.h:262
const NormSyst kNormSyst
Definition: Systs.cxx:14
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
expt
Definition: demo5.py:34
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 DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
void syst_test()
Definition: syst_test.C:39
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Var kCCE
Definition: NumuVars.h:21
void resetCalc(osc::IOscCalcAdjustable *calc)
Definition: syst_test.C:27
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:159
double sigma(TH1F *hist, double percentile)
virtual void SetRho(double rho)=0
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
virtual void SetTh23(const T &th23)=0
A simple ascii-art progress bar.
Definition: Progress.h:9
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
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
const NCScaleSyst kNCScaleSyst
Definition: Systs.cxx:16
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17