test_fc_nue.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/FC/FCPoint.h"
6 #include "CAFAna/FC/FCSurface.h"
8 #include "CAFAna/Vars/FitVars.h"
9 #include "CAFAna/Core/Loaders.h"
13 #include "CAFAna/Analysis/Plots.h"
14 #include "CAFAna/Core/Progress.h"
18 #include "CAFAna/Systs/Systs.h"
19 #include "CAFAna/Core/Utilities.h"
20 #include "CAFAna/Vars/Vars.h"
21 
22 using namespace ana;
23 
24 #include "TCanvas.h"
25 #include "TFile.h"
26 #include "TH2.h"
27 #include "TLegend.h"
28 #include "TRandom3.h"
29 
31 {
32  // Make a calculator. Any of the variants is fine
34 
35  calc->SetL(810);
36  calc->SetRho(2.75);
37  calc->SetDmsq21(7.6e-5);
38  calc->SetDmsq32(2.35e-3);
39  calc->SetTh12(asin(sqrt(.87))/2); // PDG
40  calc->SetTh13(asin(sqrt(.10))/2);
41  calc->SetTh23(TMath::Pi()/4);
42  calc->SetdCP(0);
43 
44  return calc;
45 }
46 
47 void load_files()
48 {
49  Loaders loaders("$NOVA_PROD/mc/S14-05-12/fardet/genie/prod_pid_caf_loosened_presel_neutrino2014nightmare_besttoavoid/", Loaders::kFHC);
50  // Don't bother with cosmics for now...
52 
53  const Binning bins = Binning::Simple(1, -1, +2);
54 
55  // PredictionNoExtrap pred(loaders, "", bins, kLEM, kLEM > 0.8);
56 
57  NoExtrapGenerator predGen(HistAxis("", bins, kLEM), kLEM > 0.8);
58  PredictionInterp pred({&kNormSyst, &kNCScaleSyst/*, kReweightMaCCQE*/},
59  default_calc(),
60  predGen,
61  loaders);
62 
63  loaders.Go();
64 
65  TFile* fout = new TFile("first_ana_fc_pred.root", "RECREATE");
66  pred.SaveTo(fout);
67 }
68 
69 void fill_col(int nTrials, bool genSysts, bool fitSysts)
70 {
71  TFile* fin = new TFile("~/dev_fc/CAFAna/first_ana_fc_pred.root");
73 
75 
76  kFitSinSq2Theta23.SetValue(calc, 1);
77 
78  FCCollection fccol;
79 
80  std::map<int, double> bestChis;
81 
82  TRandom3 rnd(0); // zero seeds randomly
83 
84  Progress prog("Performing mock experiments");
85  for(int iTrial = 0; iTrial < nTrials; ++iTrial){
86  const double trueX = rnd.Uniform(0, 0.4);
87  const double trueY = rnd.Uniform(0, 2);
88  kFitSinSq2Theta13.SetValue(calc, trueX);
89  kFitDeltaInPiUnits.SetValue(calc, trueY);
90 
91  SystShifts shifts;
92  if(genSysts){
93  shifts.SetShift(&kNormSyst, rnd.Gaus(0, 1));
94  shifts.SetShift(&kNCScaleSyst, rnd.Gaus(0, 1));
95  }
96 
97  // What we expect to observe at those oscillation parameters
98  Spectrum obs = pred->PredictSyst(calc, shifts);
99 
100  // Poisson fluctuations
101  Spectrum mock = obs.MockData(1e20);
102 
103  std::unique_ptr<TH1> hobs(obs.ToTH1(1e20));
104  std::unique_ptr<TH1> hmock(mock.ToTH1(1e20));
105 
106  SingleSampleExperiment expt(pred, mock);
107 
108  double trueChi;
109  if(fitSysts){
110  MinuitFitter fit(&expt, {}, {&kNormSyst, &kNCScaleSyst});
111  kFitSinSq2Theta13.SetValue(calc, trueX);
112  kFitDeltaInPiUnits.SetValue(calc, trueY);
113  trueChi = fit.Fit(calc, IFitter::kQuiet);
114  }
115  else{
116  trueChi = LogLikelihood(hobs.get(), hmock.get());
117  }
118 
119  // TODO, this is a dumb way to work this count out
120  TH1* h = mock.ToTH1(mock.POT());
121  const int nMock = h->Integral();
122  delete h;
123 
124  if(bestChis.find(nMock) == bestChis.end()){
125  std::vector<const ISyst*> systlist = {&kNormSyst, &kNCScaleSyst};
126  if(!fitSysts) systlist = {};
127 
128  MinuitFitter fit(&expt, {&kFitSinSq2Theta13, &kFitDeltaInPiUnits}, systlist);
129 
130  // TODO: this number doesn't seem to be stable for the N=0 case!
131  bestChis.insert(std::make_pair(nMock, fit.Fit(calc, IFitter::kQuiet)));
132  }
133 
134  const double bestChi = bestChis[nMock];
135 
136  const double bestX = 0;//kFitSinSq2Theta13.GetValue(&calc);
137  const double bestY = 0;//kFitDeltaInPiUnits.GetValue(&calc);
138 
139  FCPoint pt(trueX, trueY, trueChi,
140  bestX, bestY, bestChi,
141  0); // Don't know seed or rng state
142 
143  fccol.AddPoint(pt);
144 
145  prog.SetProgress((iTrial+1.)/nTrials);
146  }
147 
148  fccol.SaveToFile("fccol.root");
149 
150  std::cout << "#events\tbest chisq" << std::endl;
151  for(auto it: bestChis) std::cout << it.first << "\t" << it.second << std::endl;
152 
153  delete pred;
154 }
155 
156 void plots(bool genSysts, bool fitSysts)
157 {
158 #if 0
159  if(fitSysts){
160  FCCollection* fccolsyst = FCCollection::LoadFromFile("~/pbs/fc_Aug14/out/nue_syst_fit/*[12345]0/fccol.root").release();
161 
162  FCSurface* surffcsyst = new FCSurface(40, 0, .4,
163  40, 0, 2);
164  surffcsyst->Add(*fccolsyst);
165  delete fccolsyst;
166 
167  surffcsyst->SaveToFile("surffc_nue_syst_fit.root");
168  return;
169  }
170 
171  if(genSysts){
172  FCCollection* fccolsyst = FCCollection::LoadFromFile("~/pbs/fc_Aug14/out/nue_syst_nofit/*[12345]0/fccol.root").release();
173 
174  FCSurface* surffcsyst = new FCSurface(40, 0, .4,
175  40, 0, 2);
176  surffcsyst->Add(*fccolsyst);
177  delete fccolsyst;
178 
179  surffcsyst->SaveToFile("surffc_nue_syst_nofit.root");
180  return;
181  }
182 
183  // *** WARNING!!! THE FULL STATS RUNS POSITRON01 OUT OF MEMORY AND HANGS IT ***
184  FCCollection* fccol = FCCollection::LoadFromFile("~/pbs/fc_Aug14/out/nue/*[12345]0/fccol.root").release();
185 
186  FCCollection* fccolsyst;
187  if(fitSysts){
188  // fccolsyst = FCCollection::LoadFromFile("~/pbs/fc_Aug14/out/nue_syst_fit/*/fccol.root").release();
189  }
190  else{
191  // fccolsyst = FCCollection::LoadFromFile("~/pbs/fc_Aug14/out/nue_syst_nofit/*/fccol.root").release();
192  }
193 
194 
195  // FCCollection* fccol = FCCollection::LoadFromFile("fccol.1e6.root");
196 
197  FCSurface* surffc = new FCSurface(40, 0, .4,
198  40, 0, 2);
199  surffc->Add(*fccol);
200  delete fccol;
201 
202  surffc->SaveToFile("surffc_nue.root");
203  return;
204 
205 #else
206  FCSurface* surffc = FCSurface::LoadFromFile("surffc_nue.root").release();
207 #endif
208 
209  FCSurface* surffc_nofit = 0;
210  if(genSysts) surffc_nofit = FCSurface::LoadFromFile("surffc_nue_syst_nofit.root").release();
211  FCSurface* surffc_fit = 0;
212  if(fitSysts) surffc_fit = FCSurface::LoadFromFile("surffc_nue_syst_fit.root").release();
213 
214  FCSurface surffcsyst(40, 0, .4,
215  40, 0, 2);
216  // surffcsyst.Add(*fccolsyst);
217 
218  /*
219  for(int i = 10; i < 990; ++i){
220  TH2* h = surffc->SurfaceForSignificance(.001*i);
221  h->SetTitle(TString::Format("%g%%", i/10.));
222  h->Draw("colz");
223  h->SetMinimum(0);
224  gPad->Update();
225  gPad->Print(TString::Format("plots/fc_anim%03d.gif", i));
226  }
227  return;
228  */
229 
230  /*
231  surffc->GetBin(120, 80)->Draw();
232  gPad->Print("plots/bin.eps");
233 
234  new TCanvas;
235  surffc->SurfaceForSignificance(.68)->Draw("colz");
236  gPad->Print("plots/fc68.eps");
237 
238  new TCanvas;
239  surffcsyst.SurfaceForSignificance(.68)->Draw("colz");
240  gPad->Print("plots/fc68syst.eps");
241 
242  new TCanvas;
243  surffc->SurfaceForSignificance(.9)->Draw("colz");
244  gPad->Print("plots/fc90.eps");
245 
246  new TCanvas;
247  surffcsyst.SurfaceForSignificance(.9)->Draw("colz");
248  gPad->Print("plots/fc90syst.eps");
249  */
250 
251  new TCanvas;
252  TH2* fc68 = surffc->SurfaceForSignificance(.6827);
253  fc68->Draw("colz");
254  fc68->SetTitle("Critical value for 68% C.L.;sin^{2}2#theta_{13};#delta / #pi");
255  fc68->GetXaxis()->CenterTitle();
256  fc68->GetYaxis()->CenterTitle();
257  gPad->Print("plots/nue_fc68.eps");
258 
259  new TCanvas;
260  TH2* fc90 = surffc->SurfaceForSignificance(.9);
261  fc90->Draw("colz");
262  fc90->SetTitle("Critical value for 90% C.L.;sin^{2}2#theta_{13};#delta / #pi");
263  fc90->GetXaxis()->CenterTitle();
264  fc90->GetYaxis()->CenterTitle();
265  gPad->Print("plots/nue_fc90.eps");
266 
267  TH2 *fc68_nofit = 0, *fc90_nofit = 0;
268  if(genSysts){
269  new TCanvas;
270  fc68_nofit = surffc_nofit->SurfaceForSignificance(.6827);
271  fc68_nofit->Draw("colz");
272  fc68_nofit->SetTitle("Critical value for 68% C.L.;sin^{2}2#theta_{13};#delta / #pi");
273  fc68_nofit->GetXaxis()->CenterTitle();
274  fc68_nofit->GetYaxis()->CenterTitle();
275  gPad->Print("plots/nue_fc68_nofit.eps");
276 
277  new TCanvas;
278  fc90_nofit = surffc_nofit->SurfaceForSignificance(.9);
279  fc90_nofit->Draw("colz");
280  fc90_nofit->SetTitle("Critical value for 90% C.L.;sin^{2}2#theta_{13};#delta / #pi");
281  fc90_nofit->GetXaxis()->CenterTitle();
282  fc90_nofit->GetYaxis()->CenterTitle();
283  gPad->Print("plots/nue_fc90_nofit.eps");
284  }
285 
286  TH2 *fc68_fit = 0, *fc90_fit = 0;
287  if(fitSysts){
288  new TCanvas;
289  fc68_fit = surffc_fit->SurfaceForSignificance(.6827);
290  fc68_fit->Draw("colz");
291  fc68_fit->SetTitle("Critical value for 68% C.L.;sin^{2}2#theta_{13};#delta / #pi");
292  fc68_fit->GetXaxis()->CenterTitle();
293  fc68_fit->GetYaxis()->CenterTitle();
294  gPad->Print("plots/nue_fc68_fit.eps");
295 
296  new TCanvas;
297  fc90_fit = surffc_fit->SurfaceForSignificance(.9);
298  fc90_fit->Draw("colz");
299  fc90_fit->SetTitle("Critical value for 90% C.L.;sin^{2}2#theta_{13};#delta / #pi");
300  fc90_fit->GetXaxis()->CenterTitle();
301  fc90_fit->GetYaxis()->CenterTitle();
302  gPad->Print("plots/nue_fc90_fit.eps");
303  }
304 
305  new TCanvas;
306 
307  TFile* fin = new TFile("first_ana_fc_pred.root");
309 
311  Spectrum obs = pred->Predict(calc);
312 
313  SingleSampleExperiment expt(pred, obs.FakeData(1e20));
314 
315  FrequentistSurface surfStat(&expt, calc,
316  &kFitSinSq2Theta13, 40, 0, .4,
317  &kFitDeltaInPiUnits, 40, 0, 2);
318 
319  FrequentistSurface* surfSyst = &surfStat;
320  if(fitSysts){
321  surfSyst = new FrequentistSurface(&expt, calc,
322  &kFitSinSq2Theta13, 40, 0, .4,
323  &kFitDeltaInPiUnits, 40, 0, 2,
324  {}, {&kNormSyst, &kNCScaleSyst});
325  }
326 
327  surfStat.SetTitle("68% and 90% C.L.s");
328  if(!genSysts) surfStat.DrawContour(Gaussian68Percent1D(surfStat), 7, kRed);
329  surfStat.DrawContour(fc68, 7, kBlue);
330  if(genSysts) surfStat.DrawContour(fc68_nofit, 7, kGreen+2);
331  if(fitSysts) surfSyst->DrawContour(fc68_fit, 7, kMagenta);
332  // surfSyst->DrawContour(surffcsyst.SurfaceForSignificance(.6827), kSolid, kBlack);
333 
334  // gPad->Print("plots/cont68.eps");
335 
336  // new TCanvas;
337 
338  if(!genSysts) surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid, kRed);
339  surfStat.DrawContour(fc90, kSolid, kBlue);
340  if(genSysts) surfStat.DrawContour(fc90_nofit, kSolid, kGreen+2);
341  if(fitSysts) surfSyst->DrawContour(fc90_fit, kSolid, kMagenta);
342  // surfSyst->DrawContour(surffcsyst.SurfaceForSignificance(.9), kSolid, kBlack);
343  // gPad->Print("plots/cont90.eps");
344 
345  TLegend* leg = new TLegend(.7, .7, .9, .9);
346  leg->SetFillStyle(0);
347  TH1* dummy = new TH1F("", "", 1, 0, 1);
348  dummy->SetLineColor(kRed);
349  if(!genSysts) leg->AddEntry(dummy->Clone(), "Gaus", "l");
350  dummy->SetLineColor(kBlue);
351  leg->AddEntry(dummy->Clone(), "FC stat", "l");
352  dummy->SetLineColor(kGreen+2);
353  if(genSysts) leg->AddEntry(dummy->Clone(), "FC syst", "l");
354  dummy->SetLineColor(kMagenta);
355  if(fitSysts) leg->AddEntry(dummy->Clone(), "FC syst fit", "l");
356  leg->Draw();
357 
358  if(fitSysts){
359  gPad->Print("plots/nue_both_fit.eps");
360  return;
361  }
362 
363  if(genSysts){
364  gPad->Print("plots/nue_both_nofit.eps");
365  }
366  else{
367  gPad->Print("plots/nue_both.eps");
368  }
369 }
370 
371 void test_fc_nue(int nexpts = 1000000, bool genSysts = false, bool fitSysts = false)
372 {
373  // load_files();
374  // fill_col(nexpts, genSysts, fitSysts);
375  plots(genSysts, fitSysts);
376 }
TString fin
Definition: Style.C:24
TH2 * Gaussian90Percent1D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 1D in gaussian approximation.
const Var kLEM
PID
Definition: Vars.cxx:24
virtual void SetL(double L)=0
Implements systematic errors by interpolation between shifted templates.
Far Detector at Ash River.
Definition: SREnums.h:11
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
set< int >::iterator it
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:553
const FitSinSq2Theta23 kFitSinSq2Theta23
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:67
virtual void SetRho(double rho)=0
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:17
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetTitle(const char *str)
Definition: ISurface.cxx:192
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
void DisableLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Definition: Loaders.cxx:65
virtual void SetdCP(const T &dCP)=0
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:335
stan::math::var LogLikelihood(const std::vector< stan::math::var > &exp, const TH1 *obs)
Variant that handles the prediction in the form of Stan vars.
Definition: StanUtils.cxx:11
void fill_col(int nTrials, bool genSysts, bool fitSysts)
Definition: test_fc_nue.C:69
virtual void SetDmsq21(const T &dmsq21)=0
Generates FD-only predictions (no extrapolation)
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:98
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:170
static std::unique_ptr< FCSurface > LoadFromFile(const std::string &fname)
Definition: FCSurface.cxx:170
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
Log-likelihood scan across two parameters.
const NormSyst kNormSyst
Definition: Systs.cxx:13
void test_fc_nue(int nexpts=1000000, bool genSysts=false, bool fitSysts=false)
Definition: test_fc_nue.C:371
expt
Definition: demo5.py:34
General interface to any calculator that lets you set the parameters.
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:819
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:45
void SaveToFile(const std::string &fname) const
Definition: FCSurface.cxx:158
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
double POT() const
Definition: Spectrum.h:263
OStream cout
Definition: OStream.cxx:6
void load_files()
Definition: test_fc_nue.C:47
const Binning bins
Definition: NumuCC_CPiBin.h:8
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:46
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:787
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
osc::IOscCalculatorAdjustable * default_calc()
Definition: test_fc_nue.C:30
std::vector< Loaders * > loaders
Definition: syst_header.h:385
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:95
A simple ascii-art progress bar.
Definition: Progress.h:9
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
virtual void SetTh13(const T &th13)=0
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
void SaveToFile(const std::string &fname) const
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:103
mock
Definition: demo4.py:28
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
Float_t e
Definition: plot.C:35
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
const NCScaleSyst kNCScaleSyst
Definition: Systs.cxx:15
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
_OscCalculatorPMNSOpt< double > OscCalculatorPMNSOpt
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
void plots(bool genSysts, bool fitSysts)
Definition: test_fc_nue.C:156
TH2 * Gaussian68Percent1D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 1D in gaussian approximation.
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
h
Definition: demo3.py:41
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
virtual void SetTh12(const T &th12)=0
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:15