sensitivity.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Cuts/Cuts.h"
5 //#include "CAFAna/ExposureAssumptions.h"
6 #include "CAFAna/Fit/Fit.h"
7 #include "CAFAna/Vars/FitVars.h"
9 #include "CAFAna/Analysis/Numu.h"
10 #include "CAFAna/Analysis/Plots.h"
11 #include "CAFAna/Core/Progress.h"
15 #include "CAFAna/Core/Utilities.h"
16 
17 #include "TCanvas.h"
18 #include "TF1.h"
19 #include "TFile.h"
20 #include "TFitResult.h"
21 #include "TGaxis.h"
22 #include "TGraphErrors.h"
23 #include "TH1.h"
24 #include "TLatex.h"
25 #include "TLegend.h"
26 #include "TMath.h"
27 
28 #include <cmath>
29 #include <iostream>
30 
31 #include "Utilities/func/MathUtil.h"
32 
33 TFile* gOut = 0;
34 
35 using namespace ana;
36 using namespace ana;
37 
39 {
40  calc->SetL(810);
41  calc->SetRho(2.75);
42  calc->SetDmsq21(7.59e-5);
43  calc->SetDmsq32(2.40e-3);
44  calc->SetTh12(asin(sqrt(.87))/2); // PDG
45  calc->SetTh13(asin(sqrt(.095))/2);
46  calc->SetTh23(TMath::Pi()/4);
47  calc->SetdCP(0);
48 }
49 
51 {
52  new TCanvas(UniqueName().c_str(), UniqueName().c_str(), 800, 800);
53 
54  TH2* axes = new TH2F("", ";P_{#mue};P_{#bar{#mu}#bar{e}}", 100, 0, 0.075, 100, 0, 0.075);
55  axes->Draw();
56 
57  TGraph* one = new TGraph;
58  one->SetPoint(0, 0, 0);
59  one->SetPoint(1, 1, 1);
60  one->SetLineWidth(2);
61  one->SetLineStyle(7);
62  one->Draw("l");
63 
64  for(double E = 1; E < 3.1; E += .5){
65  for(int hie = -1; hie <= +1; hie += 2){
66  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
67 
68  TGraph* g = new TGraph;
69 
70  double avgx = 0, avgy = 0;
71  for(double delta = 0; delta < 2.01; delta += .02){
72  calc->SetdCP(delta*M_PI);
73 
74  const double P = calc->P(14, 12, E);
75  const double Pbar = calc->P(-14, -12, E);
76  g->SetPoint(g->GetN(), P, Pbar);
77 
78  avgx += P;
79  avgy += Pbar;
80  }
81 
82  avgx /= 100;
83  avgy /= 100;
84 
85  TLatex* ltx = new TLatex(avgx, avgy, TString::Format("%.1lf GeV", E));
86  ltx->SetTextAlign(22);
87  ltx->SetTextSize(0.04);
88  ltx->SetTextColor((hie > 0) ? kBlue : kRed);
89  ltx->Draw();
90 
91  g->SetLineWidth(2);
92  g->SetLineColor((hie > 0) ? kBlue : kRed);
93  g->Draw("l");
94  }
95  }
96 
97  gPad->Print("plots/biprob.eps");
98 }
99 
101  IPrediction* predFHC, IPrediction* predRHC)
102 {
103  new TCanvas;
104  calc->SetdCP(0);
105  predFHC->Predict(calc).ToTH1(18e20)->Draw/*Normalized*/();
106  calc->SetdCP(M_PI/2);
107  TH1* cpv = predFHC->Predict(calc).ToTH1(18e20);
108  cpv->SetLineColor(kRed);
109  cpv->Draw/*Normalized*/("same");
110  calc->SetdCP(-M_PI/2);
111  TH1* cpv2 = predFHC->Predict(calc).ToTH1(18e20);
112  cpv2->SetLineColor(kBlue);
113  cpv2->Draw/*Normalized*/("same");
114 
115  new TCanvas;
116  calc->SetdCP(0);
117  predRHC->Predict(calc).ToTH1(18e20)->Draw/*Normalized*/();
118  calc->SetdCP(M_PI/2);
119  cpv = predRHC->Predict(calc).ToTH1(18e20);
120  cpv->SetLineColor(kRed);
121  cpv->Draw/*Normalized*/("same");
122  calc->SetdCP(-M_PI/2);
123  cpv2 = predRHC->Predict(calc).ToTH1(18e20);
124  cpv2->SetLineColor(kBlue);
125  cpv2->Draw/*Normalized*/("same");
126 }
127 
129 {
130  // NOvANumuToy novaCC;
131  // novaCC.SetTrueValues(.96, 2.35e-3);
132  // FrequentistSurface sn(&novaCC, calc, &kFitSinSq2Theta23, 50, .9, 1, &kFitDmSq32, 50, 2e-3, 3e-3);
133  // sn.DrawContour(Gaussian68Percent2D(sn), 7, kRed);
134  // sn.DrawContour(Gaussian90Percent2D(sn), kSolid, kRed);
135  // sn.DrawBestFit(kRed);
136 
137 
138  resetCalc(calc);
139 
140  new TCanvas;
141 
143  t2k.SetPOT(7e21);
144  FrequentistSurface st(&t2k, calc, &kFitSinSq2Theta13, 20, 0, 0.2, &kFitDeltaInPiUnits, 20, 0, 2);
145  st.DrawContour(Gaussian68Percent2D(st), 7, kRed);
146  st.DrawContour(Gaussian90Percent2D(st), kSolid, kRed);
147 
148  /*
149  T2KToyExperiment t2k_anti;
150  t2k_anti.SetPOTAnti(7e21);
151  FrequentistSurface st2(&t2k_anti, calc, &kFitSinSq2Theta13, 20, 0, 0.2, &kFitDeltaInPiUnits, 20, 0, 2);
152  st2.DrawContour(Gaussian68Percent2D(st2), 7, kBlue);
153  st2.DrawContour(Gaussian90Percent2D(st2), kSolid, kBlue);
154  */
155 
156  // http://www.hep.caltech.edu/~rbpatter/theta_13/
157  ReactorExperiment react(.095, .01);
158  FrequentistSurface sr(&react, calc, &kFitSinSq2Theta13, 20, 0, 0.2, &kFitDeltaInPiUnits, 20, 0, 2);
159  sr.DrawContour(Gaussian68Percent2D(sr), 7, kMagenta);
160  sr.DrawContour(Gaussian90Percent2D(sr), kSolid, kMagenta);
161 
162  MultiExperiment multi({&t2k, /*&t2k_anti,*/ &react});
163  FrequentistSurface sm(&multi, calc, &kFitSinSq2Theta13, 40, 0, 0.2, &kFitDeltaInPiUnits, 40, 0, 2);
164  sm.DrawContour(Gaussian68Percent2D(sm), 7, kBlack);
165  sm.DrawContour(Gaussian90Percent2D(sm), kSolid, kBlack);
166 }
167 
168 /*
169 void early(osc::IOscCalculatorAdjustable* calc, IPrediction* pred, IPrediction* predRHC, double fhcScale = 1, double rhcScale = 1, TString suffix = "")
170 {
171  ana::NominalExposureAssumptions assum;
172  assum.DrawAll();
173  gPad->Print("plots/exposure.eps");
174 
175 
176  kFitSinSq2Theta13.SetValue(calc, 0);
177  const Spectrum null = pred->Predict(calc);
178  const Spectrum nullRHC = predRHC->Predict(calc);
179  kFitSinSq2Theta13.SetValue(calc, .095);
180  const Spectrum obsNH = pred->Predict(calc);
181  const Spectrum obsNHRHC = predRHC->Predict(calc);
182  calc->SetDmsq32(-2.35e-3); // inverted hierarchy
183  const Spectrum obsIH = pred->Predict(calc);
184  const Spectrum obsIHRHC = predRHC->Predict(calc);
185 
186  const double kPOTPerMonth = 18e20/(3*12);
187 
188  // Subscript is true hierarchy
189  TGraph* gExp[2] = {new TGraph, new TGraph};
190  TGraph* gDate[2] = {new TGraph, new TGraph};
191  TGraph* gDateRHC[2] = {new TGraph, new TGraph};
192 
193  const TTimeStamp kPlotBeginDate(2013, 1, 1, 0, 0, 0);
194  const TTimeStamp kBeamSwitchDate(2014, 10, 1, 0, 0, 0);
195  const TTimeStamp kPlotEndDate(2016, 1, 1, 0, 0, 0);
196 
197  for(int inv = 0; inv < 2; ++inv){
198  for(double mon = 0; mon < 8.1; mon += .05){
199  const double pot = mon*kPOTPerMonth;
200  const double sig = pot ? sqrt(LogLikelihood(null,
201  inv ? obsIH : obsNH,
202  pot)) : 0;
203  gExp[inv]->SetPoint(gExp[inv]->GetN(), mon, sig);
204  }
205  for(int d = 0; d <= 50; ++d){
206  const TTimeStamp t = kPlotBeginDate+(d/50.)*(kBeamSwitchDate-kPlotBeginDate);
207  const double pot = fhcScale*18e20*assum.ScaleFactor(t, 700, 14, 3);
208  const double sig = pot ? sqrt(LogLikelihood(null,
209  inv ? obsIH : obsNH,
210  pot)) : 0;
211  gDate[inv]->SetPoint(gDate[inv]->GetN(), t, sig);
212  }
213  // RHC
214  for(int d = 0; d <= 50; ++d){
215  const TTimeStamp t = kBeamSwitchDate+(d/50.)*(kPlotEndDate-kBeamSwitchDate);
216  const double pot = rhcScale*18e20*(assum.ScaleFactor(t, 700, 14, 3)-
217  assum.ScaleFactor(kBeamSwitchDate, 700, 14, 3));
218  const double sig = pot ? sqrt(LogLikelihood(nullRHC,
219  inv ? obsIHRHC : obsNHRHC,
220  pot)) : 0;
221  gDateRHC[inv]->SetPoint(gDateRHC[inv]->GetN(), t, sig);
222  }
223  double x, y;
224  gDate[inv]->GetPoint(gDate[inv]->GetN()-1, x, y);
225  gDate[inv]->SetPoint(gDate[inv]->GetN(), kPlotEndDate, y);
226  }
227 
228  new TCanvas;
229  TH2* axes = new TH2F("", "NO#nuA early reach. sin^{2}2#theta_{13}=0.095, sin^{2}2#theta_{23}=1.00, #delta=0;months of 700-kW, 14-kton exposure in neutrino mode;significance of observation of #nu_{#mu}#rightarrow#nu_{e} (#sigma)", 100, 0, 8, 100, 0, 5);
230  axes->Draw();
231  axes->GetXaxis()->CenterTitle();
232  axes->GetYaxis()->CenterTitle();
233  gExp[0]->SetLineWidth(3);
234  gExp[1]->SetLineWidth(3);
235  gExp[1]->SetLineStyle(7);
236  gExp[0]->Draw("l same");
237  gExp[1]->Draw("l same");
238  gOut->cd();
239  gExp[0]->Write(TString::Format("reach_vs_pot_norm%s", suffix.Data()));
240  gExp[1]->Write(TString::Format("reach_vs_pot_inv%s", suffix.Data()));
241  gPad->SetGrid();
242 
243  for(int inv = 0; inv < 2; ++inv){
244  double x, y;
245  gExp[inv]->GetPoint(gExp[inv]->GetN()/2, x, y);
246  TLatex* ltx = new TLatex(x-.5, y+.25, inv ? "IH" : "NH");
247  ltx->Draw();
248  }
249 
250  gPad->Print("plots/reach_vs_pot"+suffix+".eps");
251 
252  new TCanvas;
253  axes = new TH2F("", "NO#nuA early reach. sin^{2}2#theta_{13}=0.095, sin^{2}2#theta_{23}=1.00, #delta=0;;significance of #nu_{e} appearance (#sigma)", 100, TTimeStamp(2013, 1, 1, 0, 0, 0), TTimeStamp(2016, 1, 1, 0, 0, 0), 100, 0, 7);
254  axes->Draw();
255  TAxis* tax = axes->GetXaxis();
256  tax->SetTimeOffset(0);
257  tax->SetTimeDisplay(1);
258  tax->SetTimeFormat("%b %Y");
259  axes->GetYaxis()->CenterTitle();
260  axes->GetYaxis()->SetNdivisions(7, 5, 0, false);
261 
262  gDate[0]->SetLineColor(kBlue);
263  gDate[0]->SetLineWidth(3);
264  gDate[0]->Draw("l same");
265  gDate[1]->SetLineColor(kBlue);
266  gDate[1]->SetLineWidth(3);
267  gDate[1]->SetLineStyle(7);
268  gDate[1]->Draw("l same");
269 
270  gOut->cd();
271  gDate[0]->Write(TString::Format("reach_vs_date_fhc_norm%s", suffix.Data()));
272  gDate[1]->Write(TString::Format("reach_vs_date_fhc_inv%s", suffix.Data()));
273 
274  gDateRHC[0]->SetLineColor(kRed);
275  gDateRHC[0]->SetLineWidth(3);
276  gDateRHC[0]->Draw("l same");
277  gDateRHC[1]->SetLineColor(kRed);
278  gDateRHC[1]->SetLineWidth(3);
279  gDateRHC[1]->SetLineStyle(7);
280  gDateRHC[1]->Draw("l same");
281 
282  gOut->cd();
283  gDateRHC[0]->Write(TString::Format("reach_vs_date_rhc_norm%s", suffix.Data()));
284  gDateRHC[1]->Write(TString::Format("reach_vs_date_rhc_inv%s", suffix.Data()));
285 
286  for(int inv = 0; inv < 2; ++inv){
287  TGraph* gComb = new TGraph;
288  for(int n = 0; n < gDateRHC[inv]->GetN(); ++n){
289  double junk, y1;
290  gDate[inv]->GetPoint(gDate[inv]->GetN()-1, junk, y1);
291  double x, y2;
292  gDateRHC[inv]->GetPoint(n, x, y2);
293  gComb->SetPoint(n, x, sqrt(y1*y1+y2*y2));
294  }
295  gComb->SetLineWidth(3);
296  gComb->SetLineColor(kBlack);
297  if(inv) gComb->SetLineStyle(7);
298  gComb->Draw("l same");
299 
300  gOut->cd();
301  gComb->Write(TString::Format("reach_vs_date_comb_%s%s", inv ? "inv" : "norm", suffix.Data()));
302  }
303 
304  gPad->SetGrid();
305 
306  gPad->Print("plots/reach_vs_date"+suffix+".eps");
307 }
308 
309 class Exposure500kW: public ExposureAssumptions
310 {
311 public:
312  Exposure500kW()
313  {
314  AddMassPt(TTimeStamp(2013, 2, 15, 0, 0, 0), 0);
315  AddMassPt(TTimeStamp(2013, 6, 1, 0, 0, 0), 2.5);
316  AddMassPt(TTimeStamp(2014, 7, 26, 0, 0, 0), 14);
317  AddMassPt(TTimeStamp(2037, 1, 1, 0, 0, 0), 14);
318 
319  AddPowerPt(TTimeStamp(2013, 6, 1, 0, 0, 0), 0);
320  AddPowerPt(TTimeStamp(2013, 6, 1, 0, 0, 0), 30);
321  AddPowerPt(TTimeStamp(2013, 7, 1, 0, 0, 0), 300);
322  AddPowerPt(TTimeStamp(2013, 12, 1, 0, 0, 0), 500);
323  AddPowerPt(TTimeStamp(2037, 1, 1, 0, 0, 0), 500);
324  }
325 };
326 
327 class Exposure18kT: public ExposureAssumptions
328 {
329 public:
330  Exposure18kT()
331  {
332  AddMassPt(TTimeStamp(2013, 2, 15, 0, 0, 0), 0);
333  AddMassPt(TTimeStamp(2013, 6, 1, 0, 0, 0), 2.5);
334  AddMassPt(TTimeStamp(2015, 1, 21, 0, 0, 0), 18);
335  AddMassPt(TTimeStamp(2037, 1, 1, 0, 0, 0), 18);
336 
337  AddPowerPt(TTimeStamp(2013, 6, 1, 0, 0, 0), 0);
338  AddPowerPt(TTimeStamp(2013, 6, 1, 0, 0, 0), 30);
339  AddPowerPt(TTimeStamp(2013, 7, 1, 0, 0, 0), 300);
340  AddPowerPt(TTimeStamp(2014, 5, 1, 0, 0, 0), 700);
341  AddPowerPt(TTimeStamp(2037, 1, 1, 0, 0, 0), 700);
342  }
343 };
344 
345 void earlyAlternates(osc::IOscCalculatorAdjustable* calc,
346  IPrediction* pred, IPrediction* predRHC)
347 {
348  ana::NominalExposureAssumptions assumNom;
349  Exposure500kW assum500;
350  Exposure18kT assum18;
351 
352  ana::ExposureAssumptions* assums[3] = {&assumNom, &assum500, &assum18};
353 
354  kFitSinSq2Theta13.SetValue(calc, 0);
355  const Spectrum null = pred->Predict(calc);
356  const Spectrum nullRHC = predRHC->Predict(calc);
357  kFitSinSq2Theta13.SetValue(calc, .095);
358  const Spectrum obs = pred->Predict(calc);
359  const Spectrum obsRHC = predRHC->Predict(calc);
360 
361  const TTimeStamp kPlotBeginDate(2013, 1, 1, 0, 0, 0);
362  const TTimeStamp kPlotEndDate(2016, 1, 1, 0, 0, 0);
363 
364  new TCanvas;
365  TH2* axes = new TH2F("", "NO#nuA early reach. sin^{2}2#theta_{13}=0.095, sin^{2}2#theta_{23}=1.00, #delta=0;;significance of #nu_{e} appearance (#sigma)", 100, kPlotBeginDate, kPlotEndDate, 100, 0, 7);
366  axes->Draw();
367  TAxis* tax = axes->GetXaxis();
368  tax->SetTimeOffset(0);
369  tax->SetTimeDisplay(1);
370  tax->SetTimeFormat("%b %Y");
371  axes->GetYaxis()->CenterTitle();
372  axes->GetYaxis()->SetNdivisions(7, 5, 0, false);
373 
374  for(int assumIdx = 0; assumIdx < 2; ++assumIdx){
375  ana::ExposureAssumptions* assum = assums[assumIdx];
376  TGraph* gDate = new TGraph;
377  TGraph* gDateRHC = new TGraph;
378 
379  TTimeStamp kBeamSwitchDate;
380  if(assumIdx == 0) kBeamSwitchDate = TTimeStamp(2014, 7, 5, 0, 0, 0);
381  if(assumIdx == 1) kBeamSwitchDate = TTimeStamp(2014, 8, 25, 0, 0, 0);
382 
383  for(int d = 0; d <= 50; ++d){
384  const TTimeStamp t = kPlotBeginDate+(d/50.)*(kBeamSwitchDate-kPlotBeginDate);
385  const double pot = 18e20*assum->ScaleFactor(t, 700, 14, 3);
386  const double sig = pot ? sqrt(LogLikelihood(null, obs, pot)) : 0;
387  gDate->SetPoint(gDate->GetN(), t, sig);
388  }
389 
390  // RHC
391  for(int d = 0; d <= 50; ++d){
392  const TTimeStamp t = kBeamSwitchDate+(d/50.)*(kPlotEndDate-kBeamSwitchDate);
393  const double pot = 18e20*(assum->ScaleFactor(t, 700, 14, 3)-
394  assum->ScaleFactor(kBeamSwitchDate, 700, 14, 3));
395  const double sig = pot ? sqrt(LogLikelihood(nullRHC, obsRHC, pot)) : 0;
396  if(sig > 3) break;
397  gDateRHC->SetPoint(gDateRHC->GetN(), t, sig);
398  }
399  double x, y;
400  gDate->GetPoint(gDate->GetN()-1, x, y);
401  gDate->SetPoint(gDate->GetN(), kPlotEndDate, y);
402 
403 
404  gDate->SetLineColor(kBlue);
405  gDate->SetLineWidth(3);
406  gDate->Draw("l same");
407  if(assumIdx == 1) gDate->SetLineStyle(7);
408  if(assumIdx == 2) gDate->SetLineStyle(2);
409 
410  gOut->cd();
411  // gDate->Write(TString::Format("reach_vs_date_fhc_norm%s", suffix.Data()));
412 
413  gDateRHC->SetLineColor(kRed);
414  gDateRHC->SetLineWidth(3);
415  gDateRHC->Draw("l same");
416  if(assumIdx == 1) gDateRHC->SetLineStyle(7);
417  if(assumIdx == 2) gDateRHC->SetLineStyle(2);
418 
419  gOut->cd();
420  // gDateRHC->Write(TString::Format("reach_vs_date_rhc_norm%s", suffix.Data()));
421 
422  TGraph* gComb = new TGraph;
423  for(int n = 0; n < gDateRHC->GetN(); ++n){
424  double junk, y1;
425  gDate->GetPoint(gDate->GetN()-1, junk, y1);
426  double x, y2;
427  gDateRHC->GetPoint(n, x, y2);
428  gComb->SetPoint(n, x, sqrt(y1*y1+y2*y2));
429  }
430  gComb->SetLineWidth(3);
431  gComb->SetLineColor(kBlack);
432  gComb->Draw("l same");
433  if(assumIdx == 1) gComb->SetLineStyle(7);
434  if(assumIdx == 2) gComb->SetLineStyle(2);
435 
436  gOut->cd();
437  // gComb->Write(TString::Format("reach_vs_date_comb_%s%s", inv ? "inv" : "norm", suffix.Data()));
438  } // end for assumIdx
439 
440  gPad->SetGrid();
441 
442  // gPad->Print("plots/reach_vs_date"+suffix+".eps");
443 }
444 */
445 
446 void hierarchy(osc::IOscCalculatorAdjustable* calc, IPrediction* predFHC, IPrediction* predRHC, double potFHC, double potRHC, TString suffix)
447 {
448  for(int t2k = false; t2k <= true; ++t2k){
449  TCanvas* c = new TCanvas;
450  TH2* axes = new TH2F("", "NO#nuA hierarchy resolution. 3+3. sin^{2}2#theta_{13}=0.095, sin^{2}2#theta_{23}=1.00, #delta=0;#delta / #pi;significance of hierarchy resolution (#sigma)", 100, 0, 2, 100, 0, 3.5);
451  axes->Draw();
452  axes->GetXaxis()->CenterTitle();
453  axes->GetYaxis()->CenterTitle();
454 
455  // Fraction plot
456  TCanvas* cfrac = new TCanvas;
457  axes = new TH2F("", ";Significance of hierarchy resolution (#sigma);Fraction of #delta", 100, 0, 3.5, 100, 0, 1);
458  axes->Draw();
459  axes->GetXaxis()->CenterTitle();
460  axes->GetYaxis()->CenterTitle();
461 
462  for(int hie = -1; hie <= +1; hie += 2){
463  TGraph* g = new TGraph;
464 
465  Progress prog("Scanning delta for hierarchy plot");
466 
467  for(double d = 0; d < 2; d += .025){
468  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
469  calc->SetdCP(M_PI*d);
470 
471  const Spectrum obsFHC = predFHC->Predict(calc).FakeData(potFHC);
472  const Spectrum obsRHC = predRHC->Predict(calc).FakeData(potRHC);
473 
474  SingleSampleExperiment exptFHC(predFHC, obsFHC);
475  SingleSampleExperiment exptRHC(predRHC, obsRHC);
476 
477  MultiExperiment expts({&exptFHC, &exptRHC});
478  if(t2k){
480  e->SetPOT(7e21); // 2020 value
481  e->SetTruthParams(calc);
482  expts.Add(e);
483  }
484 
485  // Assume the wrong hierarchy in fit
486  calc->SetDmsq32(-hie*fabs(calc->GetDmsq32()));
487 
488  MinuitFitter fit(&expts, {&kFitDeltaInPiUnits});
489 
490  double chi = fit.Fit(calc, IFitter::kQuiet);
491  // Try multiples of 90 degrees away, stabilize fit
492  const double dcp = calc->GetdCP();
493  for(int n = 1; n <= 3; ++n){
494  calc->SetdCP(dcp+n*M_PI/2);
495  chi = std::min(chi, fit.Fit(calc, IFitter::kQuiet));
496  }
497 
498  g->SetPoint(g->GetN(), d, sqrt(chi));
499 
500  prog.SetProgress(d/2);
501  }
502 
503  g->SetLineWidth(3);
504  if(hie > 0){
505  g->SetLineColor(kBlue);
506  }
507  else{
508  g->SetLineColor(kRed);
509  g->SetLineStyle(7);
510  }
511  c->cd();
512  g->Draw("l same");
513  gOut->cd();
514  g->Write(TString::Format("hie_%s%s%s", (hie < 0) ? "inv" : "norm", t2k ? "_t2k" : "", suffix.Data()));
515 
516 
517  // Fraction plot
518  std::vector<double> sigs;
519  for(int i = 0; i < g->GetN()-1; ++i){
520  double x0, y0, x1, y1;
521  g->GetPoint(i, x0, y0);
522  g->GetPoint(i+1, x1, y1);
523  // Interpolate to manufacture some more points
524  for(int n = 0; n < 10; ++n)
525  sigs.push_back(y0+n/10.*(y1-y0));
526  }
527  std::sort(sigs.begin(), sigs.end());
528 
529  TGraph* frac = new TGraph;
530  const unsigned int N = sigs.size();
531  frac->SetPoint(0, 0, 1);
532  for(unsigned int i = 0; i < N; ++i){
533  frac->SetPoint(frac->GetN(), sigs[i], 1-double(i+1)/N);
534  }
535  frac->SetPoint(frac->GetN(), sigs[N-1], 0);
536  frac->SetLineWidth(3);
537  if(hie > 0){
538  frac->SetLineColor(kBlue);
539  }
540  else{
541  frac->SetLineColor(kRed);
542  frac->SetLineStyle(7);
543  }
544  cfrac->cd();
545  frac->Draw("l same");
546  gOut->cd();
547  frac->Write(TString::Format("hiefrac_%s%s%s", (hie < 0) ? "inv" : "norm", t2k ? "_t2k" : "", suffix.Data()));
548 
549  } // end for hie
550  c->cd();
551  gPad->SetGrid();
552  gPad->Print(TString::Format("plots/hie_sens%s%s.eps", t2k ? "_t2k" : "", suffix.Data()));
553 
554  cfrac->cd();
555  gPad->SetGrid();
556  gPad->Print(TString::Format("plots/hie_sens_frac%s%s.eps", t2k ? "_t2k" : "", suffix.Data()));
557  } // end for t2k
558 }
559 
560 void cpv(osc::IOscCalculatorAdjustable* calc, IPrediction* predFHC, IPrediction* predRHC, double potFHC, double potRHC, bool knownHie, TString suffix)
561 {
562  for(int t2k = false; t2k <= true; ++t2k){
563  TCanvas* c = new TCanvas;
564  TH2* axes = new TH2F("", "NO#nuA CPV determination. 3+3. sin^{2}2#theta_{13}=0.095, sin^{2}2#theta_{23}=1.00, #delta=0;#delta / #pi;significance of CP violation (#sigma)", 100, 0, 2, 100, 0, 2.5);
565  axes->Draw();
566  axes->GetXaxis()->CenterTitle();
567  axes->GetYaxis()->CenterTitle();
568 
569  // Fraction plot
570  TCanvas* cfrac = new TCanvas;
571  axes = new TH2F("", ";Significance of CP violation (#sigma);Fraction of #delta", 100, 0, 2.5, 100, 0, 1);
572  axes->Draw();
573  axes->GetXaxis()->CenterTitle();
574  axes->GetYaxis()->CenterTitle();
575 
576  for(int hie = -1; hie <= +1; hie += 2){
577  TGraph* g = new TGraph;
578 
579  Progress prog("Scanning delta for CPV plot");
580 
581  for(double d = 0; d < 2; d += .025){
582  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
583  calc->SetdCP(M_PI*d);
584 
585  const Spectrum obsFHC = predFHC->Predict(calc).FakeData(potFHC);
586  const Spectrum obsRHC = predRHC->Predict(calc).FakeData(potRHC);
587 
588  SingleSampleExperiment exptFHC(predFHC, obsFHC);
589  SingleSampleExperiment exptRHC(predRHC, obsRHC);
590 
591  MultiExperiment expts({&exptFHC, &exptRHC});
592  if(t2k){
594  e->SetPOT(7e21); // 2020 value
595  e->SetTruthParams(calc);
596  expts.Add(e);
597  }
598 
599  calc->SetdCP(0);
600  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
601  const double chiA = expts.ChiSq(calc, {});
602  calc->SetdCP(M_PI);
603  const double chiB = expts.ChiSq(calc, {});
604  calc->SetdCP(0);
605  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
606  const double chiC = expts.ChiSq(calc, {});
607  calc->SetdCP(M_PI);
608  const double chiD = expts.ChiSq(calc, {});
609 
610  if(knownHie){
611  if(hie == +1){
612  g->SetPoint(g->GetN(), d, sqrt(std::max(0., std::min(chiA, chiB))));
613  }
614  else{
615  g->SetPoint(g->GetN(), d, sqrt(std::max(0., std::min(chiC, chiD))));
616  }
617  }
618  else{
619  g->SetPoint(g->GetN(), d, sqrt(std::max(0., std::min(std::min(chiA, chiB), std::min(chiC, chiD)))));
620  }
621 
622  prog.SetProgress(d/2);
623  }
624 
625  g->SetLineWidth(3);
626  if(hie > 0){
627  g->SetLineColor(kBlue);
628  }
629  else{
630  g->SetLineColor(kRed);
631  g->SetLineStyle(7);
632  }
633  c->cd();
634  g->Draw("l same");
635  gOut->cd();
636  g->Write(TString::Format("cpv_%s%s%s", (hie < 0) ? "inv" : "norm", t2k ? "_t2k" : "", suffix.Data()));
637 
638 
639  // Fraction plot
640  std::vector<double> sigs;
641  for(int i = 0; i < g->GetN()-1; ++i){
642  double x0, y0, x1, y1;
643  g->GetPoint(i, x0, y0);
644  g->GetPoint(i+1, x1, y1);
645  // Interpolate to manufacture some more points
646  for(int n = 0; n < 10; ++n)
647  sigs.push_back(y0+n/10.*(y1-y0));
648  }
649  std::sort(sigs.begin(), sigs.end());
650 
651  TGraph* frac = new TGraph;
652  const unsigned int N = sigs.size();
653  frac->SetPoint(0, 0, 1);
654  for(unsigned int i = 0; i < N; ++i){
655  frac->SetPoint(frac->GetN(), sigs[i], 1-double(i+1)/N);
656  }
657  frac->SetPoint(frac->GetN(), sigs[N-1], 0);
658  frac->SetLineWidth(3);
659  if(hie > 0){
660  frac->SetLineColor(kBlue);
661  }
662  else{
663  frac->SetLineColor(kRed);
664  frac->SetLineStyle(7);
665  }
666  cfrac->cd();
667  frac->Draw("l same");
668  gOut->cd();
669  frac->Write(TString::Format("cpvfrac_%s%s%s", (hie < 0) ? "inv" : "norm", t2k ? "_t2k" : "", suffix.Data()));
670  }
671  c->cd();
672  gPad->SetGrid();
673  gPad->Print(TString::Format("plots/cpv_sens%s%s.eps", t2k ? "_t2k" : "", suffix.Data()));
674 
675  cfrac->cd();
676  gPad->SetGrid();
677  gPad->Print(TString::Format("plots/cpv_sens_frac%s%s.eps", t2k ? "_t2k" : "", suffix.Data()));
678  } // end for t2k
679 }
680 
682 {
683 public:
684  HackExpt(IChiSqExperiment* e, double trueDelta) : fExpt(e), fTrueDelta(trueDelta) {}
685  virtual double ChiSq(osc::IOscCalculatorAdjustable* osc) const
686  {
687  const double dm = osc->GetDmsq32();
688  osc->SetDmsq32(+fabs(dm));
689  const double delta = osc->GetdCP();
690  osc->SetdCP(fTrueDelta);
691 
692  const double chi = fExpt->ChiSq(osc);
693 
694  osc->SetDmsq32(dm);
695  osc->SetdCP(delta);
696 
697  return chi;
698  }
700  {
701  return ChiSq(calc);
702  }
703 protected:
705  double fTrueDelta;
706 };
707 
708 void th23_delta_contour(osc::IOscCalculatorAdjustable* calc, IPrediction* predFHC, IPrediction* predRHC, NumuAnalysis* numuAna, double potFHC, double potRHC, TString suffix)
709 {
710  for(int kase = 0; kase < 6; ++kase){
711  TString kaseStr;
712  double ss2th23;
713  switch(kase){
714  case 0:
715  kaseStr = "nondegen_100";
716  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
717  calc->SetdCP(1.5*M_PI);
718  ss2th23 = 1;
719  break;
720  case 1:
721  kaseStr = "neardegen_100";
722  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
723  calc->SetdCP(0.5*M_PI);
724  ss2th23 = 1;
725  break;
726  case 2:
727  kaseStr = "nondegen_097";
728  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
729  calc->SetdCP(1.5*M_PI);
730  ss2th23 = .97;
731  break;
732  case 3:
733  kaseStr = "neardegen_097";
734  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
735  calc->SetdCP(0.5*M_PI);
736  ss2th23 = .97;
737  break;
738  case 4:
739  kaseStr = "nondegen_095";
740  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
741  calc->SetdCP(1.5*M_PI);
742  ss2th23 = .95;
743  break;
744  case 5:
745  kaseStr = "neardegen_095";
746  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
747  calc->SetdCP(0.5*M_PI);
748  ss2th23 = .95;
749  break;
750  default:
751  assert(0 && "Unknown case");
752  }
753 
754  double ssth23_lo = util::sqr(sin(asin(sqrt(ss2th23))/2));
755  double ssth23_hi = util::sqr(sin((M_PI-asin(sqrt(ss2th23)))/2));
756  if(ssth23_lo > ssth23_hi) std::swap(ssth23_lo, ssth23_hi);
757 
758  kFitSinSqTheta23.SetValue(calc, ssth23_hi);
759 
760  const Spectrum obsFHC = predFHC->Predict(calc).FakeData(potFHC);
761  const Spectrum obsRHC = predRHC->Predict(calc).FakeData(potRHC);
762 
763  SingleSampleExperiment exptFHC(predFHC, obsFHC);
764  SingleSampleExperiment exptRHC(predRHC, obsRHC);
765 
766  IChiSqExperiment* exptNumu = new HackExpt(numuAna->MakeExperiment(calc, potFHC), calc->GetdCP());
767 
768  // NOvANumuToy exptNumu;
769  // exptNumu.SetTrueValues(kFitSinSq2Theta23.GetValue(calc),
770  // fabs(calc->GetDmsq32()));
771 
772  MultiExperiment expts({&exptFHC, &exptRHC, exptNumu});
773 
774  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
775  FrequentistSurface surfNHFHC(&exptFHC, calc,
776  &kFitDeltaInPiUnits, 40, 0, 2,
777  &kFitSinSqTheta23, 80, 0, 1);
778 
779  FrequentistSurface surfNHRHC(&exptRHC, calc,
780  &kFitDeltaInPiUnits, 40, 0, 2,
781  &kFitSinSqTheta23, 80, 0, 1);
782 
783  FrequentistSurface surfNumuNH(exptNumu, calc,
784  &kFitDeltaInPiUnits, 40, 0, 2,
785  &kFitSinSqTheta23, 80, 0, 1);
786 
787  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
788  FrequentistSurface surfIHFHC(&exptFHC, calc,
789  &kFitDeltaInPiUnits, 40, 0, 2,
790  &kFitSinSqTheta23, 80, 0, 1);
791 
792  FrequentistSurface surfIHRHC(&exptRHC, calc,
793  &kFitDeltaInPiUnits, 40, 0, 2,
794  &kFitSinSqTheta23, 80, 0, 1);
795 
796  FrequentistSurface surfNumuIH(exptNumu, calc,
797  &kFitDeltaInPiUnits, 40, 0, 2,
798  &kFitSinSqTheta23, 80, 0, 1);
799 
800  new TCanvas;
801  surfNHFHC.DrawContour(Gaussian68Percent2D(surfNHFHC), kSolid, kBlue, 0);
802  surfNHRHC.DrawContour(Gaussian68Percent2D(surfNHFHC), kSolid, kCyan+2, 0);
803  surfIHFHC.DrawContour(Gaussian68Percent2D(surfIHFHC), kSolid, kRed, 0);
804  surfIHRHC.DrawContour(Gaussian68Percent2D(surfIHFHC), kSolid, kMagenta, 0);
805  surfNumuNH.DrawContour(Gaussian68Percent2D(surfNumuNH), kSolid, kBlack, 0);
806  surfNumuIH.DrawContour(Gaussian68Percent2D(surfNumuIH), kSolid, kGray+2, 0);
807 
808  gPad->Print("plots/th23_delta_sens_detail1sig_"+kaseStr+suffix+".eps");
809 
810  new TCanvas;
811  surfNHFHC.DrawContour(Gaussian2Sigma2D(surfNHFHC), kSolid, kBlue, 0);
812  surfNHRHC.DrawContour(Gaussian2Sigma2D(surfNHFHC), kSolid, kCyan+2, 0);
813  surfIHFHC.DrawContour(Gaussian2Sigma2D(surfIHFHC), kSolid, kRed, 0);
814  surfIHRHC.DrawContour(Gaussian2Sigma2D(surfIHFHC), kSolid, kMagenta, 0);
815  surfNumuNH.DrawContour(Gaussian2Sigma2D(surfNumuNH), kSolid, kBlack, 0);
816  surfNumuIH.DrawContour(Gaussian2Sigma2D(surfNumuIH), kSolid, kGray+2, 0);
817 
818  gPad->Print("plots/th23_delta_sens_detail2sig_"+kaseStr+suffix+".eps");
819 
820  new TCanvas;
821  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
822  FrequentistSurface surfNH(&expts, calc,
823  &kFitDeltaInPiUnits, 40, 0, 2,
824  &kFitSinSqTheta23, 40, 0.3, .7);
825  surfNH.DrawContour(Gaussian68Percent2D(surfNH), 7, kBlue);
826  surfNH.DrawContour(Gaussian2Sigma2D(surfNH), kSolid, kBlue);
827  surfNH.DrawBestFit(kBlue);
828 
829  surfNH.ToTH2(surfNH.BestLikelihood())->Write("th23_delta_sens_norm_"+kaseStr+suffix);
830 
831  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
832  FrequentistSurface surfIH(&expts, calc,
833  &kFitDeltaInPiUnits, 40, 0, 2,
834  &kFitSinSqTheta23, 40, 0.3, .7);
835  surfIH.DrawContour(Gaussian68Percent2D(surfIH), 7, kRed, surfNH.BestLikelihood());
836  surfIH.DrawContour(Gaussian2Sigma2D(surfIH), kSolid, kRed, surfNH.BestLikelihood());
837 
838  surfIH.ToTH2(surfNH.BestLikelihood())->Write("th23_delta_sens_inv_"+kaseStr+suffix);
839 
840  TGraph* g = new TGraph;
841  g->SetPoint(0, 0, .5);
842  g->SetPoint(1, 2, .5);
843  g->SetLineWidth(2);
844  g->SetLineStyle(7);
845  g->Draw("l");
846  gPad->SetGrid();
847  gPad->Print("plots/th23_delta_sens_"+kaseStr+suffix+".eps");
848  } // end for kase
849 }
850 
851 void octant(osc::IOscCalculatorAdjustable* calc, IPrediction* predFHC, IPrediction* predRHC, double potFHC, double potRHC, TString suffix)
852 {
853  for(int kase = 0; kase < 4; ++kase){
854  const double ss2th23 = kase%2 ? .95 : .97;
855  const char* thStr = kase%2 ? "095" : "097";
856  const bool upper = kase/2;
857  const char* octStr = upper ? "upper" : "lower";
858 
859  double ssth23_lo = util::sqr(sin(asin(sqrt(ss2th23))/2));
860  double ssth23_hi = util::sqr(sin((M_PI-asin(sqrt(ss2th23)))/2));
861  if(ssth23_lo > ssth23_hi) std::swap(ssth23_lo, ssth23_hi);
862 
863  const double ssth23_true = upper ? ssth23_hi : ssth23_lo;
864  const double ssth23_assum = upper ? ssth23_lo : ssth23_hi;
865 
866  TCanvas* c = new TCanvas;
867  TH2* axes = new TH2F("", ";#delta / #pi;significance of octant determination (#sigma)", 100, 0, 2, 100, 0, 4);
868  axes->Draw();
869  axes->GetXaxis()->CenterTitle();
870  axes->GetYaxis()->CenterTitle();
871  axes->SetTitle(TString::Format("sin^{2}2#theta_{23}=0.%d, %s octant", int(ss2th23*100+.5), octStr));
872 
873  TCanvas* cfrac = new TCanvas;
874  axes = new TH2F("", ";significance of octant determination (#sigma);Fraction of #delta", 100, 0, 4, 100, 0, 1);
875  axes->Draw();
876  axes->GetXaxis()->CenterTitle();
877  axes->GetYaxis()->CenterTitle();
878  axes->SetTitle(TString::Format("sin^{2}2#theta_{23}=0.%d, %s octant", int(ss2th23*100+.5), octStr));
879 
880  for(int trueHie = -1; trueHie <= +1; trueHie += 2){
881  TGraph* g = new TGraph;
882 
883  Progress prog("Scanning delta for octant plot");
884  for(double delta = 0; delta < 2.025; delta += .05){
885  prog.SetProgress(delta/2);
886  calc->SetdCP(delta*M_PI);
887 
888  kFitSinSqTheta23.SetValue(calc, ssth23_true);
889 
890  calc->SetDmsq32(trueHie*fabs(calc->GetDmsq32()));
891 
892  const Spectrum obsFHC = predFHC->Predict(calc).FakeData(potFHC);
893  const Spectrum obsRHC = predRHC->Predict(calc).FakeData(potRHC);
894 
895  SingleSampleExperiment exptFHC(predFHC, obsFHC);
896  SingleSampleExperiment exptRHC(predRHC, obsRHC);
897 
898  MultiExperiment expts({&exptFHC, &exptRHC});
899 
900 
901  kFitSinSqTheta23.SetValue(calc, ssth23_assum); // Assume wrong octant
902 
903  MinuitFitter fit(&expts, {&kFitDeltaInPiUnits});
904  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
905  double chi = fit.Fit(calc, IFitter::kQuiet);
906  // Try 180 degrees away. Seems to be enough to stabilize fit
907  calc->SetdCP(calc->GetdCP()+M_PI);
908  chi = std::min(chi, fit.Fit(calc, IFitter::kQuiet));
909 
910  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
911  chi = std::min(chi, fit.Fit(calc, IFitter::kQuiet));
912  calc->SetdCP(calc->GetdCP()+M_PI);
913  chi = std::min(chi, fit.Fit(calc, IFitter::kQuiet));
914  if(chi < 0) chi = 0;
915 
916  g->SetPoint(g->GetN(), delta, sqrt(chi));
917  } // end for delta
918  prog.Done();
919 
920  g->SetLineWidth(3);
921  if(trueHie == +1){
922  g->SetLineColor(kBlue);
923  }
924  else{
925  g->SetLineColor(kRed);
926  g->SetLineStyle(7);
927  }
928 
929  c->cd();
930  g->Draw("l same");
931  gOut->cd();
932  g->Write(TString::Format("octant_%s_%s_%s%s", thStr, octStr, (trueHie < 0) ? "inv" : "norm", suffix.Data()));
933 
934  // Fraction plot
935  std::vector<double> sigs;
936  for(int i = 0; i < g->GetN()-1; ++i){
937  double x0, y0, x1, y1;
938  g->GetPoint(i, x0, y0);
939  g->GetPoint(i+1, x1, y1);
940  // Interpolate to manufacture some more points
941  for(int n = 0; n < 10; ++n)
942  sigs.push_back(y0+n/10.*(y1-y0));
943  }
944  std::sort(sigs.begin(), sigs.end());
945 
946  TGraph* frac = new TGraph;
947  const unsigned int N = sigs.size();
948  frac->SetPoint(0, 0, 1);
949  for(unsigned int i = 0; i < N; ++i){
950  frac->SetPoint(frac->GetN(), sigs[i], 1-double(i+1)/N);
951  }
952  frac->SetPoint(frac->GetN(), sigs[N-1], 0);
953  frac->SetLineWidth(3);
954  if(trueHie > 0){
955  frac->SetLineColor(kBlue);
956  }
957  else{
958  frac->SetLineColor(kRed);
959  frac->SetLineStyle(7);
960  }
961  cfrac->cd();
962  frac->Draw("l same");
963  gOut->cd();
964  frac->Write(TString::Format("octantfrac_%s_%s_%s%s", thStr, octStr, (trueHie < 0) ? "inv" : "norm", suffix.Data()));
965  } // end for trueHie
966 
967  c->cd();
968  gPad->SetGrid();
969  gPad->Print(TString::Format("plots/octant_sens_%s_%s%s.eps", thStr, octStr, suffix.Data()));
970  cfrac->cd();
971  gPad->SetGrid();
972  gPad->Print(TString::Format("plots/octant_sens_frac_%s_%s%s.eps", thStr, octStr, suffix.Data()));
973  } // end for kase
974 }
975 
976 void delta_precision(osc::IOscCalculatorAdjustable* calc, IPrediction* predFHC, IPrediction* predRHC, double potFHC, double potRHC, TString suffix)
977 {
978  const int Nbins = 100;
979 
980  TH2* surfs[2];
981 
982  TGraph* one = new TGraph;
983  one->SetPoint(0, 0, 0);
984  one->SetPoint(1, 360, 0);
985  one->SetLineWidth(2);
986  one->SetLineStyle(7);
987 
988  TGraph* diag = new TGraph;
989  diag->SetPoint(0, 0, 0);
990  diag->SetPoint(1, 2, 2);
991  diag->SetLineWidth(2);
992  diag->SetLineStyle(7);
993 
994  TFile* f = new TFile("delta_err.root", "RECREATE");
995 
996  for(int trueHie = -1; trueHie <= +1; trueHie += 2){
997  const char* hieStr = (trueHie < 0) ? "inv" : "norm";
998  new TCanvas;
999 
1000  TH2* surfFit = new TH2F("", ";True value of #delta / #pi;Fit value of #delta / #pi", Nbins, 0, 2, Nbins, 0, 2);
1001  TH2* surfDiff = new TH2F("", ";True value of #delta (degrees);#delta_{reco}-#delta_{true} (degrees)", Nbins, 0, 360, Nbins-1, (-1+1./Nbins)*180, (+1-1./Nbins)*180);
1002 
1003  for(int ix = 1; ix <= Nbins; ++ix){
1004  const double x = surfFit->GetXaxis()->GetBinCenter(ix);
1005  kFitDeltaInPiUnits.SetValue(calc, x);
1006  calc->SetDmsq32(trueHie*fabs(calc->GetDmsq32()));
1007  const Spectrum obsFHC = predFHC->Predict(calc).FakeData(potFHC);
1008  const Spectrum obsRHC = predRHC->Predict(calc).FakeData(potRHC);
1009 
1010  SingleSampleExperiment exptFHC(predFHC, obsFHC);
1011  SingleSampleExperiment exptRHC(predRHC, obsRHC);
1012  MultiExperiment expts({&exptFHC, &exptRHC});
1013 
1014  calc->SetDmsq32(trueHie*fabs(calc->GetDmsq32()));
1015  TH1* slc = Slice(&expts, calc, &kFitDeltaInPiUnits, Nbins, 0, 2);
1016  slc->Draw("l");
1017  slc->GetYaxis()->SetRangeUser(0, 14);
1018  gPad->Update();
1019 
1020  for(int is = 1; is <= Nbins; ++is){
1021  const double y = slc->GetXaxis()->GetBinCenter(is);
1022  surfFit->Fill(x, y, slc->GetBinContent(is)+1e-10);
1023 
1024  double diff = y-x;
1025  if(diff < -1) diff += 2;
1026  if(diff > +1) diff -= 2;
1027 
1028  surfDiff->Fill(x*180, diff*180, slc->GetBinContent(is)+1e-10);
1029  }
1030 
1031  delete slc;
1032  }
1033 
1034  new TCanvas;
1035  surfFit->DrawClone("colz");
1036  gOut->cd();
1037  surfFit->Write(TString::Format("delta_err_fit_%s%s", hieStr, suffix.Data()));
1038 
1039  const double level = 1;
1040  surfFit->SetContour(1, &level);
1041  surfFit->SetLineWidth(2);
1042  surfFit->DrawClone("cont3 same");
1043  diag->Draw();
1044  gPad->Update();
1045  gPad->Print(TString::Format("plots/delta_err_fit_%s%s.eps", hieStr, suffix.Data()));
1046 
1047  new TCanvas;
1048  surfDiff->DrawClone("colz");
1049  gOut->cd();
1050  surfDiff->Write(TString::Format("delta_err_fit_%s%s", hieStr, suffix.Data()));
1051 
1052  surfDiff->SetContour(1, &level);
1053  surfDiff->SetLineWidth(2);
1054  surfDiff->DrawClone("cont3 same");
1055  one->Draw();
1056  gPad->Update();
1057  gPad->Print(TString::Format("plots/delta_err_diff_%s%s.eps", hieStr, suffix.Data()));
1058 
1059  surfs[(trueHie+1)/2] = surfDiff;
1060  }
1061 
1062  new TCanvas;
1063  TH2* axes = new TH2F("", ";True value of #delta (degrees);#delta_{reco}-#delta_{true} (degrees)", Nbins, 0, 360, Nbins, -180, +180);
1064  axes->Draw();
1065  axes->GetXaxis()->CenterTitle();
1066  axes->GetXaxis()->SetNdivisions(4, 6, 0, false);
1067  axes->GetYaxis()->CenterTitle();
1068  axes->GetYaxis()->SetNdivisions(4, 6, 0, false);
1069 
1070  surfs[0]->SetLineColor(kRed);
1071  surfs[0]->Draw("cont3 same");
1072  surfs[1]->SetLineColor(kBlue);
1073  surfs[1]->Draw("cont3 same");
1074  one->Draw("l");
1075 
1076  gOut->cd();
1077  surfs[0]->Write(TString::Format("delta_err_inv%s", suffix.Data()));
1078  surfs[1]->Write(TString::Format("delta_err_norm%s", suffix.Data()));
1079 
1080  gPad->Print(TString::Format("plots/delta_err_diff%s.eps", suffix.Data()));
1081 }
1082 
1083 void delta_slices(osc::IOscCalculatorAdjustable* calc, IPrediction* predFHC, IPrediction* predRHC, double potFHC, double potRHC, TString suffix)
1084 {
1085  const double ss2th23 = 0.95;
1086 
1087  for(int trueHie = -1; trueHie <= +1; trueHie += 2){
1088  for(int trueOct = -1; trueOct <= +1; trueOct += 2){
1089  for(int deltaIdx = 0; deltaIdx < 4; ++deltaIdx){
1090  calc->SetDmsq32(trueHie*fabs(calc->GetDmsq32()));
1091 
1092  kFitDeltaInPiUnits.SetValue(calc, .5*deltaIdx);
1093 
1094  double ssth23_lo = util::sqr(sin(asin(sqrt(ss2th23))/2));
1095  double ssth23_hi = util::sqr(sin((M_PI-asin(sqrt(ss2th23)))/2));
1096  if(ssth23_lo > ssth23_hi) std::swap(ssth23_lo, ssth23_hi);
1097 
1098  kFitSinSqTheta23.SetValue(calc, (trueOct < 0) ? ssth23_lo : ssth23_hi);
1099 
1100  const Spectrum obsFHC = predFHC->Predict(calc).FakeData(potFHC);
1101  const Spectrum obsRHC = predRHC->Predict(calc).FakeData(potRHC);
1102 
1103  SingleSampleExperiment exptFHC(predFHC, obsFHC);
1104  SingleSampleExperiment exptRHC(predRHC, obsRHC);
1105  MultiExperiment expts({&exptFHC, &exptRHC});
1106 
1107  for(int doSqrt = false; doSqrt <= true; ++doSqrt){
1108  const int Nbins = 100;
1109 
1110  new TCanvas;
1111  TH2* axes = new TH2F("", "3+3 years, sin^{2}2#theta_{13}=0.095, sin^{2}2#theta_{23}=0.95;#delta / #pi;#Delta#chi^{2}", 100, 0, 2, 100, 0, doSqrt ? 7 : 40);
1112  axes->GetXaxis()->CenterTitle();
1113  axes->GetYaxis()->CenterTitle();
1114 
1115  if(doSqrt) axes->GetYaxis()->SetTitle("Significance (#sigma)");
1116 
1117  axes->Draw();
1118 
1119  TLegend* leg = new TLegend(.3, .6, .5, .85);
1120  leg->SetFillStyle(0);
1121 
1122  for(int assumOct = -1; assumOct <= +1; assumOct += 2){
1123  kFitSinSqTheta23.SetValue(calc, (assumOct < 0) ? ssth23_lo : ssth23_hi);
1124  for(int assumHie = -1; assumHie <= +1; assumHie += 2){
1125  calc->SetDmsq32(assumHie*fabs(calc->GetDmsq32()));
1126  TH1* slc = (doSqrt ? SqrtSlice : Slice)(&expts, calc, &kFitDeltaInPiUnits, Nbins, 0, 2, 0, {});
1127 
1128  if(assumHie > 0) slc->SetLineColor(kBlue); else slc->SetLineColor(kRed);
1129  if(assumOct < 0) slc->SetLineStyle(7);
1130  slc->SetLineWidth(3);
1131 
1132  slc->Draw("l same");
1133  gOut->cd();
1134  slc->Write(TString::Format("delta_slice_%s_%s_%d%s_assum%s_assum%s%s", (trueHie < 0) ? "inv" : "norm", (trueOct < 0) ? "lo" : "up", deltaIdx, doSqrt ? "_sqrt" : "", (assumHie < 0) ? "inv" : "norm", (assumOct < 0) ? "lo" : "up", suffix.Data()));
1135 
1136  leg->AddEntry(slc, TString::Format("%s %s",
1137  (assumHie < 0) ? "IH, " : "NH,",
1138  (assumOct < 0) ? "#theta_{23}<#pi/4" : "#theta_{23}>#pi/4"), "l");
1139  }
1140  }
1141 
1142  leg->Draw();
1143  if(doSqrt) gPad->SetGrid();
1144  gPad->Print(TString::Format("plots/delta_%sslices_%s_%s_%gpi%s.eps", doSqrt ? "sqrt_" : "", (trueHie < 0) ? "inv" : "norm", (trueOct < 0) ? "lo" : "up", deltaIdx*.5, suffix.Data()));
1145  } // end for doSqrt
1146  } // end for deltaIdx
1147  } // end for trueOct
1148  } // end for trueHie
1149 }
1150 
1152 {
1154  resetCalc(&calc);
1155 
1156  biprob(&calc);
1157  resetCalc(&calc);
1158  // return;
1159 
1160  // toys(&calc);
1161  resetCalc(&calc);
1162 
1163 
1164  const std::string fnameMC = "$NOVA_DATA/mc/S13-06-18/hadd/fd_S13-06-18_fhc_nonswap.sim.caf.root";
1165 
1166  const std::string fnameSwap = "$NOVA_DATA/mc/S13-06-18/hadd/fd_S13-06-18_fhc_fluxswap.sim.caf.root";
1167 
1168  SpectrumLoader loader(fnameMC);
1169  SpectrumLoader loaderSwap(fnameSwap);
1170 
1171  PredictionNoExtrap pred(loader, loaderSwap,
1172  "LEM", Binning::Simple(20, .9, 1),
1173  kLEM, kPassesPresel);
1174 
1175  // PredictionNoExtrap predE(loader, loaderSwap,
1176  // "Calorimetric Energy", Binning::Simple(20, 0, 4),
1177  // kCaloE, kPassesPresel && kPassesLEM);
1178 
1179  NumuAnalysis numuAna(loader, loaderSwap);
1180 
1181  loader.Go();
1182  loaderSwap.Go();
1183 
1184  const std::string fnameRHC = "$NOVA_DATA/mc/S13-06-18/hadd/fd_S13-06-18_rhc_nonswap.sim.caf.root";
1185 
1186  const std::string fnameRHCSwap = "$NOVA_DATA/mc/S13-06-18/hadd/fd_S13-06-18_rhc_fluxswap.sim.caf.root";
1187 
1188  SpectrumLoader loaderRHC(fnameRHC);
1189  SpectrumLoader loaderRHCSwap(fnameRHCSwap);
1190 
1191  PredictionNoExtrap predRHC(loaderRHC, loaderRHCSwap,
1192  "LEM", Binning::Simple(20, .9, 1),
1193  kLEM, kPassesPresel);
1194 
1195  /// PredictionNoExtrap predERHC(loaderRHC, loaderRHCSwap,
1196  // "Calorimetric Energy", Binning::Simple(20, 0, 4),
1197  // kCaloE, kPassesPresel && kPassesLEM);
1198 
1199  loaderRHC.Go();
1200  loaderRHCSwap.Go();
1201 
1202  /*
1203  spects(&calc, &pred, &predRHC);
1204  resetCalc(&calc);
1205  spects(&calc, &predE, &predERHC);
1206  resetCalc(&calc);
1207  return;
1208  */
1209 
1210 
1211  gOut = new TFile("sensitivity.root", "RECREATE");
1212 
1213 
1214  // earlyAlternates(&calc, &pred, &predRHC);
1215  // return;
1216 
1217  // early(&calc, &pred, &predRHC);
1218  // resetCalc(&calc);
1219  // return;
1220 
1221  // early(&calc, &pred, &predRHC, 1.3, 1.3, "_fluxscale");
1222  // resetCalc(&calc);
1223 
1224 
1225  for(int exposureIdx = 0; exposureIdx < 6; ++exposureIdx){
1226  if(exposureIdx != 0 && exposureIdx != 4 && exposureIdx != 5) continue;
1227  double potFHC, potRHC;
1228  TString suffix;
1229  switch(exposureIdx){
1230  case 0: potFHC = potRHC = 18e20; suffix = ""; break;
1231  case 1: potFHC = potRHC = 18e20*1.3; suffix = "_fluxscale"; break;
1232  // case 1: potFHC = 18e20; potRHC = 18e20*1.35; suffix = "_rhcscale"; break;
1233  case 2: potFHC = potRHC = 18e20*5./3; suffix = "_5yr"; break;
1234  case 3: potFHC = potRHC = 18e20/3; suffix = "_1yr"; break;
1235  case 4: potFHC = potRHC = 18e20*2; suffix = "_2xexp"; break;
1236  case 5: potFHC = potRHC = 18e20*4; suffix = "_4xexp"; break;
1237  default: assert(0);
1238  }
1239 
1240  hierarchy(&calc, &pred, &predRHC, potFHC, potRHC, suffix);
1241  resetCalc(&calc);
1242 
1243  cpv(&calc, &pred, &predRHC, potFHC, potRHC, false, suffix);
1244  cpv(&calc, &pred, &predRHC, potFHC, potRHC, true, suffix+"_known");
1245  resetCalc(&calc);
1246 
1247  th23_delta_contour(&calc, &pred, &predRHC, &numuAna, potFHC, potRHC, suffix);
1248  resetCalc(&calc);
1249 
1250  octant(&calc, &pred, &predRHC, potFHC, potRHC, suffix);
1251  resetCalc(&calc);
1252 
1253  delta_precision(&calc, &pred, &predRHC, potFHC, potRHC, suffix);
1254  resetCalc(&calc);
1255 
1256  delta_slices(&calc, &pred, &predRHC, potFHC, potRHC, suffix);
1257  resetCalc(&calc);
1258  }
1259 }
T max(const caf::Proxy< T > &a, T b)
TGraph * Slice(const IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, MinuitFitter::FitOpts opts)
scan in one variable, holding all others constant
Definition: Fit.cxx:169
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:85
const Var kLEM
PID
Definition: Vars.cxx:24
virtual void SetL(double L)=0
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
cons_index_list< index_multi, nil_index_list > multi
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
HackExpt(IChiSqExperiment *e, double trueDelta)
Definition: sensitivity.C:684
Float_t y1[n_points_granero]
Definition: compare.C:5
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
double delta
Definition: runWimpSim.h:98
TFile * gOut
Definition: sensitivity.C:33
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
void hierarchy(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC, double potFHC, double potRHC, TString suffix)
Definition: sensitivity.C:446
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Float_t x1[n_points_granero]
Definition: compare.C:5
string fnameSwap
Definition: demo4.py:6
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:67
virtual void SetRho(double rho)=0
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
double ChiSq(osc::IOscCalculatorAdjustable *calc, const ana::SystShifts &) const
Definition: sensitivity.C:699
void delta_slices(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC, double potFHC, double potRHC, TString suffix)
Definition: sensitivity.C:1083
virtual void SetdCP(const T &dCP)=0
void th23_delta_contour(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC, NumuAnalysis *numuAna, double potFHC, double potRHC, TString suffix)
Definition: sensitivity.C:708
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
virtual void SetDmsq21(const T &dmsq21)=0
TGraphAsymmErrors * t2k
Definition: Xsec_final.C:30
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
sigs
#define M_PI
Definition: SbMath.h:34
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
void delta_precision(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC, double potFHC, double potRHC, TString suffix)
Definition: sensitivity.C:976
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
void biprob(osc::IOscCalculatorAdjustable *calc)
Definition: sensitivity.C:50
#define P(a, b, c, d, e, x)
Log-likelihood scan across two parameters.
Produce IChiSqExperiment objects encapsulating the analysis.
Definition: Numu.h:15
void resetCalc(osc::IOscCalculatorAdjustable *calc)
Definition: sensitivity.C:38
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:33
General interface to any calculator that lets you set the parameters.
void SetTruthParams(osc::IOscCalculatorAdjustable *osc)
Float_t E
Definition: plot.C:20
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:819
void SetPOT(double pot)
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:45
Base class defining interface for experiments.
T2K experiment, based on scaling published numbers.
void octant(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC, double potFHC, double potRHC, TString suffix)
Definition: sensitivity.C:851
const Cut kPassesPresel(PassesPreselFunc)
Does this event pass the nue preselection?
Definition: Cuts.h:14
Float_t d
Definition: plot.C:236
virtual void Go() override
Load all the registered spectra.
void cpv(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC, double potFHC, double potRHC, bool knownHie, TString suffix)
Definition: sensitivity.C:560
double frac(double x)
Fractional part.
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
loader
Definition: demo0.py:10
Very simple model allowing inclusion of reactor constraints.
Oscillation probability calculators.
Definition: Calcs.h:5
virtual T P(int flavBefore, int flavAfter, double E)=0
E in GeV; flavors as PDG codes (so, neg==>antinu)
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:170
TGraph * SqrtSlice(const IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, MinuitFitter::FitOpts opts)
Forward to Slice but sqrt the result for a crude significance.
Definition: Fit.cxx:180
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
virtual T GetDmsq32() const
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
void toys(osc::IOscCalculatorAdjustable *calc)
Definition: sensitivity.C:128
const FitSinSqTheta23 kFitSinSqTheta23
Definition: FitVars.cxx:15
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
double fTrueDelta
Definition: sensitivity.C:705
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void spects(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC)
Definition: sensitivity.C:100
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Prediction that just uses FD MC, with no extrapolation.
auto one()
Definition: PMNS.cxx:49
T min(const caf::Proxy< T > &a, T b)
void sensitivity()
Definition: sensitivity.C:1151
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
virtual double Fit(osc::IOscCalculatorAdjustable *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:56
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
void Done()
Call this when action is completed.
Definition: Progress.cxx:91
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
virtual double ChiSq(osc::IOscCalculatorAdjustable *osc) const
Definition: sensitivity.C:685
IChiSqExperiment * MakeExperiment(osc::IOscCalculator *calc, double pot)
Definition: Numu.cxx:25
double BestLikelihood() const
Definition: ISurface.h:29
static constexpr Double_t sr
Definition: Munits.h:164
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
gm Write()
IChiSqExperiment * fExpt
Definition: sensitivity.C:704
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:15