plot_syst_contours.C
Go to the documentation of this file.
1 #include "OscLib/OscCalcDumb.h"
2 
4 #include "CAFAna/Fit/Fit.h"
7 
8 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Core/IFitVar.h"
12 
13 #include "CAFAna/Vars/FitVars.h"
14 
20 
21 
23 
24 #include "Utilities/func/MathUtil.h"
25 
26 #include "TCanvas.h"
27 #include "TLine.h"
28 #include "TGraph.h"
29 #include "TH2.h"
30 #include "TH1.h"
31 #include "TFile.h"
32 #include "TLegend.h"
33 #include "TRandom3.h"
34 
35 #include <iostream>
36 
37 using namespace ana;
38 
39 double pot = 9.48e20;
40 
41 
42 void SetHierOct(osc::IOscCalcAdjustable* calc, bool NH, bool UO)
43 {
44  double dmsq32 = NH ? 2.7e-3 : -2.7e-3;
45  kFitDmSq32.SetValue(calc,dmsq32);
46  double ssth23 = UO ? 0.62 : 0.40;
47  kFitSinSqTheta23.SetValue(calc,ssth23);
48 }
49 
52  std::vector<const ISyst*> systs = {})
53 {
54  double minChiSq = 1e5;
56  for (bool Hier : {true, false}){
57  for (bool Oct : {true, false}){
58  int NPts = 4;
59  for (double dCP = 0.5/NPts; dCP < 2; dCP += 2./(NPts+1)){
60  SetHierOct(calc,Hier,Oct);
62  minChiSq = std::min(minChiSq,fit.Fit(calc))->EvalMetricVal();
63  }
64  }
65  }
66  return minChiSq;
67 }
68 
69 
71 {
72  h->SetTitle(title.c_str());
73  h->GetYaxis()->SetTitle("#Delta#chi^{2}");
74  h->GetYaxis()->CenterTitle();
75  h->GetXaxis()->SetTitle("#delta_{CP}/#pi");
76  h->GetXaxis()->CenterTitle();
77 }
78 
79 void MakePlot(IPrediction *pred, Spectrum mockdata, int sIdx=-1)
80 {
81 
82  if (sIdx >= int(kNue2017Systs.size())){
83  std::cout << "This was not a valid syst index, try again." << std::endl;
84  return; // Don't assert, other calls to function may be valid
85  }
86 
88  const std::vector<const IFitVar*> slcProfVars = {
91  &kFitDmSq32
92  };
93 
94 
95  std::vector<const ISyst*> systs;
96  if (sIdx == -1)
97  for (const ISyst *syst : kNue2017Systs) systs.push_back(syst);
98  else
99  systs.push_back(kNue2017Systs[sIdx]);
100 
101  SingleSampleExperiment nue(pred,mockdata);
102 
103  // Until numu has files, just use their SA result
104  std::string fNumuName = ("/nova/ana/nu_e_ana/SecondAna/SAoutputs_Numu.root");
105  AtmConstraint *numu = new AtmConstraint(fNumuName,
106  "2dchisqmapstatonly",
107  "2dchisqmapstatonlyinv");
108  const ReactorExperiment th13Constr(0.085, 0.005);
109 
110  MultiExperiment expt({&nue,&th13Constr,numu});
111 
112 
113  // Set up the slice plots with and without systs
114  std::vector<TH1*> slcStat;
115  std::vector<TH1*> slcSyst;
116  double minChiSqStatSlc = GetMinChiSqSlc(&expt,calc);
117  double minChiSqSystSlc = GetMinChiSqSlc(&expt,calc,systs);
118  for (bool Hier : {true, false}){
119  for (bool Oct : {false, true}){
120  // Make sure that our hierarchy and octant are set right
121  SetHierOct(calc,Hier,Oct);
122  // Add to our vector of slices
123  TH1 *curStat = Profile(&expt,calc,&kFitDeltaInPiUnits,80,0,2,
124  minChiSqStatSlc,slcProfVars);
125  slcStat.push_back(curStat);
126  TH1 *curSyst = Profile(&expt,calc,&kFitDeltaInPiUnits,80,0,2,
127  minChiSqSystSlc,slcProfVars,systs);
128  slcSyst.push_back(curSyst);
129  }
130  }
131 
132  // Set up surfaces for NH and IH with and without systs
133  // Make sure we're testing the NH
134  kFitDmSq32.SetValue(calc,+2.7e-3);
135  FrequentistSurface surfSystNH(&expt,calc,
136  &kFitDeltaInPiUnits,30,0,2,
137  &kFitSinSqTheta23,30,0.28,0.72,
139  kFitDmSq32.SetValue(calc,+2.7e-3);
140  FrequentistSurface surfStatNH(&expt,calc,
141  &kFitDeltaInPiUnits,30,0,2,
142  &kFitSinSqTheta23,30,0.28,0.72,
143  {&kFitDmSq32,&kFitSinSq2Theta13});
144  // Make sure we're testing the IH
145  kFitDmSq32.SetValue(calc,-2.7e-3);
146  FrequentistSurface surfSystIH(&expt,calc,
147  &kFitDeltaInPiUnits,30,0,2,
148  &kFitSinSqTheta23,30,0.28,0.72,
149  {&kFitDmSq32,&kFitSinSq2Theta13}, systs);
150  kFitDmSq32.SetValue(calc,-2.7e-3);
151  FrequentistSurface surfStatIH(&expt,calc,
152  &kFitDeltaInPiUnits,30,0,2,
153  &kFitSinSqTheta23,30,0.28,0.72,
154  {&kFitDmSq32,&kFitSinSq2Theta13});
155  double minChiSqStatCont = std::min(surfStatNH.BestLikelihood(), surfStatIH.BestLikelihood());
156  double minChiSqSystCont = std::min(surfSystNH.BestLikelihood(), surfSystIH.BestLikelihood());
157 
158 
159  TCanvas *canSyst = new TCanvas(UniqueName().c_str(),"",1000,500);
160  canSyst->Divide(3,2);
161 
162  // List out the true params here to draw in the hist titles
163  std::vector<std::string> titles = {"NH - LO","NH - UO","IH - LO","IH - UO"};
164 
165  // Now, swicth over to plotting mode. Start with the slices
166  for (int i = 0; i < 4; i++){
167  int paneIdx = i < 2 ? i+1 : i+2;
168  canSyst->cd(paneIdx);
169  slcStat[i]->SetLineColor(kRed);
170  slcSyst[i]->SetLineColor(kBlue);
171  HandleTitles(slcStat[i],titles[i]);
172  slcStat[i]->SetMinimum(0);
173  slcStat[i]->Draw("c");
174  slcSyst[i]->Draw("c,same");
175 
176  if (i==0){
177  TLegend *legSlc = AutoPlaceLegend(0.2,0.2);
178  legSlc->AddEntry(slcStat[i],"Stat Only","l");
179  legSlc->AddEntry(slcSyst[i],"Systs","l");
180  legSlc->Draw();
181  }
182  }
183 
184 
185  canSyst->cd(3);
186  surfSystNH.DrawContour(Gaussian95Percent2D(surfSystIH),
187  1, kBlue, minChiSqSystCont);
188  surfStatNH.DrawContour(Gaussian95Percent2D(surfStatNH),
189  1, kRed, minChiSqStatCont);
190 
191  canSyst->cd(6);
192  surfSystIH.DrawContour(Gaussian95Percent2D(surfSystIH),
193  1, kBlue, minChiSqSystCont);
194  surfStatIH.DrawContour(Gaussian95Percent2D(surfStatIH),
195  1, kRed, minChiSqStatCont);
196 
197 
198  std::string idxStr = "_"+std::to_string(sIdx);
199  canSyst->SaveAs(("StatOnlySyst_NuePlots_"+idxStr+".png").c_str());
200 
201 
202 }
203 
204 
205 void plot_syst_contours(std::string fMockData="out_mockdata_for_syst_contours.root",int systIdx=-1)
206 {
207  std::string predFName = "/nova/ana/users/ecatanom/Ana2017/Predictions/v2/pred_nue_ana2017_full_v2.root";
208 
209  std::vector<const ISyst*> systs;
210  for (const ISyst *syst : kNue2017Systs) systs.push_back(syst);
211 
212 
213  const std::vector<const IFitVar*> fitvars = {
216  &kFitDmSq32,
218  };
219 
221  ENue2017ExtrapType extrap = kCombo;
222  PredictionSystNue2017 pred(extrap, calc, kNue2017Systs, predFName);
223 
224 
225  Spectrum mock = *(LoadFromFile<Spectrum>(fMockData,"mockdata").release());
226 
227 
228 
229  MakePlot(&pred,mock,systIdx);
230 
231  for (int i = 0; i <= int(systs.size()); i++)
232  MakePlot(&pred,mock,i);
233 
234 
235 
236 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double ssth23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void plot_syst_contours(std::string fMockData="out_mockdata_for_syst_contours.root", int systIdx=-1)
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const std::vector< const DummyNue2017Syst * > kNue2017Systs
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Log-likelihood scan across two parameters.
TH2 * Gaussian95Percent2D(const FrequentistSurface &s)
Up-value surface for 95% confidence in 2D in gaussian approximation.
Loads shifted spectra from files.
expt
Definition: demo5.py:34
void SetHierOct(osc::IOscCalcAdjustable *calc, bool NH, bool UO)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double dCP
double GetMinChiSqSlc(const IExperiment *expt, osc::IOscCalcAdjustable *calc, std::vector< const ISyst * > systs={})
double pot
Very simple model allowing inclusion of reactor constraints.
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double dmsq32
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Base class defining interface for experiments.
Definition: IExperiment.h:14
void HandleTitles(TH1 *h, std::string title)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)
mock
Definition: demo4.py:28
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
Compare a single data spectrum to the MC + cosmics expectation.
void MakePlot(IPrediction *pred, Spectrum mockdata, int sIdx=-1)
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17