compare_pid_cut_sens.C
Go to the documentation of this file.
1 // This is the macro that creates the plots in docdb 13300. The input file
2 // "nue_first_ana_sens.root" is created by make_nue_firstana_extrap.C
3 
4 //scale real conditions pre+postshutdown to full exposure equivalent
5 //TODO update POT with final number
6 const double kPot = 3.52e20;
7 
8 const int kSurfBins = 50;
9 
10 
11 #include "CAFAna/Vars/FitVars.h"
20 #include "CAFAna/Fit/Fit.h"
22 #include "CAFAna/Core/Utilities.h"
24 #include "CAFAna/Systs/NueSysts.h"
25 #include "CAFAna/Analysis/Calcs.h"
26 using namespace ana;
27 
28 #include "Utilities/func/MathUtil.h"
29 
32 
33 #include "TCanvas.h"
34 #include "TH2.h"
35 #include "TMath.h"
36 
37 #include <iomanip>
38 #include <iostream>
39 
40 
42 {
43  osc::OscCalculatorDumb calcdumb;
44 
46 
47  Spectrum tot = pred->Predict(&calcdumb);
48 
49  std::cout << "sig: " << sig.Integral(kPot) << std::endl;
50  std::cout << "bkg: " << tot.Integral(kPot) - sig.Integral(kPot) << std::endl;
51  std::cout << "cos: " << cosmic->Integral(cosmic->POT()) << std::endl;
52 }
53 
56  const std::string& suffix,
57  nueSystPID systPID,
58  bool drawDayaBay,
59  std::vector<const IFitVar*> margVars = {})
60 {
61  new TCanvas;
62 
64  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
65 
66  FrequentistSurface surfNH(expt, calc,
69  margVars,
70  {getAllNueFirstAnaSysts(systPID)}, false);
71 
72  surfNH.SetTitle(suffix.c_str());
73 
74  surfNH.DrawContour(Gaussian68Percent1D(surfNH), 7, kRed);
75  surfNH.DrawContour(Gaussian90Percent1D(surfNH), kSolid, kRed);
76 
77  gPad->Update();
78 
80  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
81 
82  FrequentistSurface surfIH(expt, calc,
85  margVars,
86  {getAllNueFirstAnaSysts(systPID)}, false);
87 
88  surfIH.DrawContour(Gaussian68Percent1D(surfIH), 7, kBlue);
89  surfIH.DrawContour(Gaussian90Percent1D(surfIH), kSolid, kBlue);
90 
91  gPad->Update();
92 
93  //overlay dayabay best fit
95 
96  if (drawDayaBay){
97  const IChiSqExperiment* exptReact = DayaBayConstraint2014();
98  FrequentistSurface surfReact(exptReact, calc,
99  &kFitSinSq2Theta13, 100, 0, 0.2,
100  &kFitDeltaInPiUnits, 100, 0, 2);
101  surfReact.SetTitle("Daya Bay");
102  surfReact.DrawContour(Gaussian68Percent1D(surfReact), 7, kGreen);
103  surfReact.DrawContour(Gaussian90Percent1D(surfReact), 1, kGreen);
104  }
105 
106  gPad->Print(("conts_"+suffix+".eps").c_str());
107 }
108 
111  const std::string& suffix,
112  nueSystPID systPID,
113  std::vector<const IFitVar*> margVars = {})
114 {
115  new TCanvas;
116 
117  ResetOscCalcToDefault(calc);
118  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
119 
120  FrequentistSurface surfNH(expt, calc,
122  &kFitSinSqTheta23, kSurfBins, 0, 1,
123  margVars, {getAllNueFirstAnaSysts(systPID)}, false);
124 
125  surfNH.SetTitle(suffix.c_str());
126 
127  surfNH.DrawContour(Gaussian68Percent1D(surfNH), 7, kRed);
128  surfNH.DrawContour(Gaussian90Percent1D(surfNH), kSolid, kRed);
129 
130  gPad->Update();
131 
132  ResetOscCalcToDefault(calc);
133  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
134 
135  FrequentistSurface surfIH(expt, calc,
137  &kFitSinSqTheta23, kSurfBins, 0, 1,
138  margVars, {getAllNueFirstAnaSysts(systPID)}, false);
139 
140  surfIH.DrawContour(Gaussian68Percent1D(surfIH), 7, kBlue);
141  surfIH.DrawContour(Gaussian90Percent1D(surfIH), kSolid, kBlue);
142 
143  gPad->Print(("conts_th23_"+suffix+".eps").c_str());
144 }
145 
148  const std::string& suffix)
149 {
150  new TCanvas;
151 
152  ResetOscCalcToDefault(calc);
153  calc->SetDmsq32(+fabs(calc->GetDmsq32()));
154 
155  TH1* sliceNH = SqrtSlice(expt, calc, &kFitDeltaInPiUnits, 100, 0, 2, 0);
156 
157  ResetOscCalcToDefault(calc);
158  calc->SetDmsq32(-fabs(calc->GetDmsq32()));
159 
160  TH1* sliceIH = SqrtSlice(expt, calc, &kFitDeltaInPiUnits, 100, 0, 2, 0);
161 
162  sliceIH->SetTitle(suffix.c_str());
163  sliceIH->GetYaxis()->SetRangeUser(0, 2);
164 
165  sliceIH->SetLineColor(kBlue);
166  sliceIH->Draw("l");
167 
168  sliceNH->SetLineColor(kRed);
169  sliceNH->Draw("l same");
170 
171  gPad->Print(("slices_"+suffix+".eps").c_str());
172 }
173 
174 void MakeSignificanceTable(double n_h0, double n_ih, double n_nh,
175  double e_h0, double e_ih, double e_nh)
176 {
177  std::cout << std::fixed << std::setprecision(2) << std::endl;
178 
179  std::cout << "\\begin{tabular}{c|rrrr}" << std::endl;
180  std::cout << "N & $P=0$ & IH,${\\pi\\over2}$ & NH,${3\\pi\\over2}$ & Diff.\\\\" << std::endl;
181  std::cout << "\\hline" << std::endl;
182  for(int n = 0; n <= 14; ++n){
183  const double sig_h0 = CountingExperimentSigmaByLL(n, n_h0, e_h0);
184 
185  // "N or more"
186  const double p_ih = CountingExperimentPValueByLL(n, n_ih, e_ih);
187  const double p_nh = CountingExperimentPValueByLL(n, n_nh, e_nh);
188 
189  const double sig_ih = PValueToSigma(p_ih);
190  const double sig_nh = PValueToSigma(p_nh);
191 
192  const double diff = sqrt(fabs(sig_ih*sig_ih - sig_nh*sig_nh));
193 
194  std::cout << n << " & " << sig_h0 << " & " << sig_ih << " & " << sig_nh << " & " << diff << "\\\\" << std::endl;
195  }
196  std::cout << "\\end{tabular}" << std::endl;
197 }
198 
200 {
201  const std::string fname = "nue_first_ana_sens.root";
202 
204  nueSystPID systPID = kSystLID;
205  for(const std::string suffix: {"LID","LEM"}){
206  if (suffix == "LID") systPID = kSystLID;
207  if (suffix == "LEM") systPID = kSystLEM;
208  const PredictionInterp* pred = LoadFromFile<PredictionInterp>(fname, "prediction"+suffix).release();
209 
210  const Spectrum* cosmic = LoadFromFile<Spectrum>(fname, "fdCosmic"+suffix+"Scale").release();
211 
212 
213  std::cout << suffix << std::endl;
214  PrintEventCounts(pred, cosmic);
215 
216  ResetOscCalcToDefault(&calc);
217  calc.SetDmsq32(+fabs(calc.GetDmsq32()));
218 
219  Spectrum fake = pred->Predict(&calc).FakeData(kPot);
220 
221  // Adding cosmics in sucks. The event count is right but we have to do this
222  // to do this to force the POT to something that won't get incorrectly
223  // scaled.
224  Spectrum cosmic_add(cosmic->ToTH1(cosmic->POT()), fake.POT(), 0);
225 
226  fake += cosmic_add;
227 
228  CountingExperiment countExpt(pred, fake, *cosmic);
229  MultiExperiment expt13;
230  SolarConstraints solarExpt;
231  expt13.Add(&countExpt);
232  expt13.Add(T2KAtmConstraint());
233  expt13.Add(&solarExpt);
234  std::vector<const IFitVar*> margVars13{&kFitSinSq2Theta23,
235  &kFitDmSq32,
236  &kFitDmSq21,
238  Theta13DeltaContours(&expt13, &calc, suffix, systPID, true, margVars13);
239 
240  MultiExperiment expt23;
241  expt23.Add(&countExpt);
242  expt23.Add(&solarExpt);
243  std::vector<const IFitVar*> margVars23{&kFitDmSq21,
244  &kFitTanSqTheta12};
245  Theta23DeltaContours(&expt23, &calc, suffix, systPID, margVars23);
246 
247  // These are the parameters suggested by MINOS
248  ResetOscCalcToDefault(&calc);
249  calc.SetDmsq32(-fabs(calc.GetDmsq32()));
250  calc.SetdCP(TMath::Pi()/2);
251 
252  Spectrum fakeSlice1 = pred->Predict(&calc).FakeData(kPot);
253  fakeSlice1 += cosmic_add;
254 
255  CountingExperiment countSlice1(pred, fakeSlice1, *cosmic);
256 
257  DeltaSlices(&countSlice1, &calc, suffix+"_ih");
258 
259  //make delta slices with daya bay constraint
260  MultiExperiment exptSlice1;
261  exptSlice1.Add(&countSlice1);
262  exptSlice1.Add(T2KAtmConstraint());
263  exptSlice1.Add(&solarExpt);
264  exptSlice1.Add(DayaBayConstraint2014());
265  DeltaSlices(&exptSlice1, &calc, suffix+"_constrained_ih");
266 
267  // This is the other good case
268  ResetOscCalcToDefault(&calc);
269  calc.SetDmsq32(+fabs(calc.GetDmsq32()));
270  calc.SetdCP(3*TMath::Pi()/2);
271 
272  Spectrum fakeSlice2 = pred->Predict(&calc).FakeData(kPot);
273  fakeSlice2 += cosmic_add;
274 
275  CountingExperiment countSlice2(pred, fakeSlice2, *cosmic);
276 
277  DeltaSlices(&countSlice2, &calc, suffix+"_nh");
278 
279  //make delta slices with daya bay constraint
280  MultiExperiment exptSlice2;
281  exptSlice2.Add(&countSlice2);
282  exptSlice2.Add(T2KAtmConstraint());
283  exptSlice2.Add(&solarExpt);
284  exptSlice2.Add(DayaBayConstraint2014());
285  DeltaSlices(&exptSlice2, &calc, suffix+"_constrained_nh");
286 
287  // Doesn't matter much what osc parameters n_h0 is calculated at
288  osc::OscCalculatorDumb calcdumb;
289  const double n_h0 = (pred->Predict(&calcdumb) - pred->PredictComponent(&calcdumb, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth)).Integral(kPot) + cosmic->Integral(cosmic->POT());
290 
291  const double n_ih = fakeSlice1.Integral(fakeSlice1.POT());
292  const double n_nh = fakeSlice2.Integral(fakeSlice2.POT());
293  double e_h0 = CalcNueBkgFirstAnaSystError(systPID)*n_h0;
294  double e_ih = util::pythag(abs(n_ih-n_h0)*CalcNueSigFirstAnaSystError(systPID),e_h0);
295  double e_nh = util::pythag(abs(n_nh-n_h0)*CalcNueSigFirstAnaSystError(systPID),e_h0);
296 
297  cout<<"Significance table without errors\n";
298  MakeSignificanceTable(n_h0, n_ih, n_nh, 0, 0, 0);
299  cout<<"Significance table with systematic errors\n";
300  MakeSignificanceTable(n_h0, n_ih, n_nh, e_h0, e_ih, e_nh);
301  }
302 }
void Add(const IChiSqExperiment *expt)
TH2 * Gaussian90Percent1D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 1D in gaussian approximation.
Implements systematic errors by interpolation between shifted templates.
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void Theta23DeltaContours(IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *calc, const std::string &suffix, nueSystPID systPID, std::vector< const IFitVar * > margVars={})
Constraints on the parameters and from solar experiments.
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:16
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
const FitSinSq2Theta23 kFitSinSq2Theta23
double Integral(const Spectrum &s, const double pot, int cvnregion)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetTitle(const char *str)
Definition: ISurface.cxx:192
double CountingExperimentSigmaByLL(int N, double b0, double sigma)
void ResetOscCalcToDefault(osc::IOscCalculatorAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
const FitDmSq21 kFitDmSq21
Definition: FitVars.cxx:24
float abs(float number)
Definition: d0nt_math.hpp:39
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
void SetDmsq32(const T &dmsq32) override
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
void compare_pid_cut_sens()
Log-likelihood scan across two parameters.
double CalcNueSigFirstAnaSystError(NuESystPID type)
Definition: NueSysts.cxx:94
void SetdCP(const T &dCP) override
Charged-current interactions.
Definition: IPrediction.h:39
expt
Definition: demo5.py:34
General interface to any calculator that lets you set the parameters.
const AtmConstraint * T2KAtmConstraint()
Constraint from T2K&#39;s Oct 31, 2014 paper, arXiv 1410.1532v1.
Sum up livetimes from individual cosmic triggers.
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:819
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:45
Base class defining interface for experiments.
double CountingExperimentPValueByLL(int N, double b0, double sigma)
virtual void SetDmsq32(const T &dmsq32)=0
void DeltaSlices(IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *calc, const std::string &suffix)
void Theta13DeltaContours(IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *calc, const std::string &suffix, nueSystPID systPID, bool drawDayaBay, std::vector< const IFitVar * > margVars={})
double POT() const
Definition: Spectrum.h:263
OStream cout
Definition: OStream.cxx:6
const ReactorExperiment * DayaBayConstraint2014()
A ReactorExperiment initialized with the Nu2014 Daya Bay constraints.
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
virtual T GetDmsq32() const
std::vector< const ISyst * > getAllNueFirstAnaSysts(NuESystPID type)
Get a vector of all the nue group systs.
Definition: NueSysts.cxx:77
Spectrum Predict(osc::IOscCalculator *calc) const override
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
const FitSinSqTheta23 kFitSinSqTheta23
Definition: FitVars.cxx:15
void MakeSignificanceTable(double n_h0, double n_ih, double n_nh, double e_h0, double e_ih, double e_nh)
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
double CalcNueBkgFirstAnaSystError(NuESystPID type)
Definition: NueSysts.cxx:164
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
const double kPot
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
void PrintEventCounts(const IPrediction *pred, const Spectrum *cosmic)
Compare a data spectrum to MC expectation using only the event count.
TH2 * Gaussian68Percent1D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 1D in gaussian approximation.
const int kSurfBins
const FitTanSqTheta12 kFitTanSqTheta12
Definition: FitVars.cxx:22