PlotNus17Prediction.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TDirectory.h"
3 #include "TFile.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TLatex.h"
7 #include "TLegend.h"
8 #include "TLine.h"
9 #include "TPad.h"
10 
11 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Fit/Fit.h"
13 #include "CAFAna/Analysis/Style.h"
15 #include "CAFAna/Core/IFitVar.h"
17 #include "CAFAna/Cuts/TimingCuts.h"
28 
29 #include "OscLib/OscCalcSterile.h"
30 
32 
33 using namespace ana;
34 
35 // Definitions for organization
36 
37 // One struct to contain all values for angle fits
38 struct AngleValues{
39  int nbins23;
40  double min23;
41  double max23;
42 
43  int nbins34;
44  double min34;
45  double max34;
46 
47  int nbins24;
48  double min24;
49  double max24;
50 };
51 
52 void TexInfo();
53 
54 TLegend* MakeLegend(double x1, double y1, double x2, double y2, double textsize = 0.03)
55 {
56  TLegend* leg = new TLegend(x1, y1, x2, y2);
57  leg->SetBorderSize(0);
58  leg->SetFillColor(kWhite);
59  leg->SetFillStyle(0);
60  leg->SetFillStyle(0);
61  leg->SetTextFont(42);
62  leg->SetTextSize(textsize);
63  return leg;
64 }
65 
66 void SigmaLines(double xmin, double xmax);
67 
68 void PlotText(const double xpos, const double start, const double step);
69 
70 // Reset calculator values after each fit to preset defaults
72 
73 // Plot a single 1D slice canvas -- this function does the actual plotting
74 // The macro does not call this function directly in the main function
75 // \param expt The CAFAna experiment object, contains FD prediction and "data"
76 // \param calc Sterile calculator for oscillating the spectra
77 // \param var Parameter to fit for!
78 // \param profVars Vector of parameters to profile over
79 // \param rootOut Already open file for root output
80 // \param name Base of the string for creating the name for saving the root objects
81 
83  const IFitVar* var, std::vector<const IFitVar*> vars,
84  int nbins, double min, double max,
85  TDirectory* rootOut, std::string name, std::string special,
86  Color_t color, std::vector<const ISyst*> systs = {});
87 
88 
89 void CompareSlices(TH1* stats, TH1* systs, std::string angle, TDirectory* rootOut, double xmin, double xmax);
90 
91 // Plot 2D slice surfaces
92 // The macro does not call this function directly
93 // \param rootOut Already open file for root output
94 // \param prof String indicating whether surfaces are profiled or not (optional)
95 void Plot2DSlice(IExperiment* expt1,
96  IExperiment* expt2,
97  osc::OscCalcSterile* calc,
98  const IFitVar* varA,
99  int nbinsA, double minA, double maxA,
100  const IFitVar* varB,
101  int nbinsB, double minB, double maxB,
102  TDirectory* rootOut, std::string angles,
103  std::vector<const IFitVar*> vars,
104  std::vector<const ISyst*> systs = {});
105 
106 // Plot spectra and stack canvases -- the workhorse version
107 // The macro does not call this function directly
108 // \param spectra Array of spectra, should match the output of the SPECARR macro
109 // \param rootOut Already open file for root output
110 // \param textOFS Already open filestream for text output
111 // \param name A name for the output, only part of the saved output object's name
112 // \param title Title for the histograms. Typically this is the x axis label
113 // \param det The detector the spectra correspond to, should be 2 chars long for nice formatting
114 // \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
115 // \param sterile Pointer to a spectrum with 4 flavor osc prediction, includes all event types (optional)
116 // \param sterileNC Pointer to a spectrum with 4 flavor osc prediction, only includes NC events (optional)
117 void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
119  Spectrum* sterile = nullptr, Spectrum* sterileNC = nullptr);
120 
121 // Plot spectra and stack canvases for ND
122 // Sets up call to the workhorse version of PlotStack
123 // \param decomp The ND decomposition
124 // \param rootOut Already open file for root output
125 // \param textOFS Already open filestream for text output
126 // \param name A name for the output, only part of the saved output object's name
127 // \param title Title for the histograms. Typically this is the x axis label
128 // \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
129 void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
131 
132 // Plot spectra and stack canvases for FD
133 // Sets up call to the workhorse version of PlotStack
134 // \param pred The FD prediction object
135 // \param sCos The FD cosmic prediction spectrum
136 // \param calc Sterile calculator assuming 3 flavor oscillations
137 // \param calcSt Sterile calculator assuming 4 flavor oscillations
138 // \param rootOut Already open file for root output
139 // \param textOFS Already open filestream for text output
140 // \param name A name for the output, only part of the saved output object's name
141 // \param title Title for the histograms. Typically this is the x axis label
142 // \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
144  Spectrum& sCos,
145  osc::IOscCalc* calc,
146  osc::IOscCalc* calcSt,
147  TDirectory* rootOut, FILE* textOFS,
149 
150 // Function to plot purity and efficiency at both FD and ND
151 // \param pFDNum Prediction object for FD, used as nominator for FD purity and efficiency
152 // \param pFDPurDen Prediction object for FD, used as denominator for FD purity
153 // \param pFDPurDen Prediction object for FD, used as denominator for FD efficiency
154 // \param sNDNum Spectrum object for ND, used as nominator for ND purity and efficiency
155 // \param sNDPurDen Spectrum object for ND, used as denominator for ND purity
156 // \param sNDPurDen Spectrum object for ND, used as denominator for ND efficiency
157 // \param calc Sterile calculator assuming 3 flavor oscillations
158 // \param rootOut Already open file for root output
159 // \param textOFS Already open filestream for text output
160 // \param name A name for the output, only part of the saved output object's name
161 // \param title Title for the histograms. Typically this is the x axis label
162 void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
163  Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
164  osc::IOscCalc* calc, Spectrum& sCos,
165  TDirectory* rootOut, FILE* textOFS,
167 
168 
169 
170 void PlotNus17Prediction(bool plotSensitivities = false)
171 {
172  TH1::AddDirectory(0);
173 
174  // Set up path to file with filled CAFAna objects and for output of analysis
175  //std::string folder = "/nova/ana/steriles/Ana01/";
176  std::string folder = "/nova/ana/steriles/Nus17/";
177  std::string filenm = "predInterpSysts_nus17_v1.root";
178 
179 
180  std::string loadLocation = folder + filenm + ".root";
181  std::string saveLocation = filenm + "Ana.root";
182  std::string textLocation = filenm + "Ana.txt";
183 
184  // Load filled objects from file
185  TFile* rootL = new TFile(loadLocation.c_str(), "UPDATE");
186  TDirectory* tmpL = gDirectory;
187  TDirectory* loadDir = gDirectory;
188 
189  loadDir->cd((loadLocation + ":/decompNC").c_str());
190  ProportionalDecomp decompNC = *ProportionalDecomp::LoadFrom(gDirectory);
191  loadDir->cd((loadLocation + ":/prediction").c_str());
193  loadDir->cd((loadLocation + ":/sCos1").c_str());
194  Spectrum sCosmic = *Spectrum ::LoadFrom(gDirectory);
195 
196  loadDir->cd((loadLocation + ":/pFDPENum1").c_str());
197  PredictionExtrap pFDPurEffNum = *PredictionExtrap ::LoadFrom(gDirectory);
198  loadDir->cd((loadLocation + ":/pFDPDen1").c_str());
199  PredictionExtrap pFDPurityDen = *PredictionExtrap ::LoadFrom(gDirectory);
200  loadDir->cd((loadLocation + ":/pFDEDen1").c_str());
201  PredictionExtrap pFDEfcncyDen = *PredictionExtrap ::LoadFrom(gDirectory);
202  loadDir->cd((loadLocation + ":/sNDPENum1").c_str());
203  Spectrum sNDPurEffNum = *Spectrum ::LoadFrom(gDirectory);
204  loadDir->cd((loadLocation + ":/sNDPDen1").c_str());
205  Spectrum sNDPurityDen = *Spectrum ::LoadFrom(gDirectory);
206  loadDir->cd((loadLocation + ":/sNDEDen1").c_str());
207  Spectrum sNDEfcncyDen = *Spectrum ::LoadFrom(gDirectory);
208  tmpL->cd();
209  rootL->Close(); // Don't forget to close the file!
210 
211  // Set up oscillation calculator that uses default 3 flavor parameters
214 
215  Spectrum prediction9 = pred.Predict(calc3f);
216 
217  const double livetime9 = 440.;
218  const double pot9 = 8.85e20;
219  const double livetime24 = 440.*(24/9);
220  const double pot24 = 24e20;
221 
222  sCosmic.Scale((440.)/sCosmic.Livetime());
223  sCosmic.OverridePOT(9e20);
224 
225  Spectrum sCosmic9 = sCosmic.FakeData(9.49e20);
226  Spectrum sCosmic24 = sCosmic.FakeData(24e20);
227 
228  Spectrum sFDdata9 = prediction9.FakeData(8.85e20) + sCosmic9;
229 
230  // Set up oscillation calculator that uses default 3 flavor parameters
232  calc4f->SetNFlavors(4);
233  ResetAngles(calc4f);
234 
235  std::string labelEReco = "Energy Deposited in (GeV)";
236 
237  // Open files for saving analysis output
238  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
239  FILE* textF;
240  textF = fopen(textLocation.c_str(), "w");
241 
242  // Choose some angle values for making 4 flavor predictions
243  //calc4f->SetAngle(2, 4, 0.0698132); // 4 degrees, 10 degrees=0.17453292
244  //calc4f->SetAngle(3, 4, 0.349066); // 20 degrees, 35 degrees=0.61086524
245 
246  calc4f->SetAngle(2, 4, 0.3630285); // 20.8 degrees, 10 degrees=0.17453292
247  calc4f->SetAngle(3, 4, 0.5445427); // 31.2 degrees, 35 degrees=0.61086524
248 
249  //calc4f->SetAngle(2, 4, 0.3630285); // 4 degrees, 10 degrees=0.17453292
250  //calc4f->SetAngle(3, 4, 0.3979351); // 22.8 degrees, 35 degrees=0.61086524
251 
252  // Plot spectra for 9e20
253  std::string titleHelper = " Spectra for 9e20 POT";
254  PlotStack(&pred, sCosmic9, calc3f, calc4f, rootF, textF, "FDSpectra9", titleHelper, 9);
255  PlotStack(&decompNC, rootF, textF, "NDSpectra9", titleHelper, 9);
256 
257  // Plot spectra for 24e20
258  titleHelper = " Spectra for 24e20 POT";
259  PlotStack(&pred, sCosmic24, calc3f, calc4f, rootF, textF, "FDSpectra24", titleHelper, 24);
260  PlotStack(&decompNC, rootF, textF, "NDSpectra24", titleHelper, 24);
261 
262  // Plot purity and efficiency
263  PlotPurEff(&pFDPurEffNum, &pFDPurityDen, &pFDEfcncyDen,
264  sNDPurEffNum, sNDPurityDen, sNDEfcncyDen,
265  calc3f, sCosmic9, rootF, textF, "PurEff", labelEReco);
266 
267  fclose(textF); // Text file is done now, close it
268 
269  ResetAngles(calc4f);
270 
271  if(plotSensitivities){
272  // Create Experiment objects, one each for the 1 data yr and 3 data yr predictions
273  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.57, 0.23);
274 
275  CountingExperiment expt9c(&pred, sFDdata9, sCosmic9);
276  MultiExperiment multi9c({&th23Constraint, &expt9c});
277 
278  SingleSampleExperiment expt9s(&pred, sFDdata9, sCosmic9);
279  MultiExperiment multi9s({&th23Constraint, &expt9s});
280 
281  // Vector of parameters to fit: theta 23, 34, 24
282  std::vector<const IFitVar*> fitvars;
283  fitvars.push_back(&kFitTheta34InDegreesSterile);
284  fitvars.push_back(&kFitTheta24InDegreesSterile);
285  fitvars.push_back(&kFitTheta23InDegreesSterile);
286  fitvars.push_back(&kFitDelta24InPiUnitsSterile);
287 
288  // Create fit objects
289  MinuitFitter fit9c(&multi9c, fitvars);
290  MinuitFitter fit9s(&multi9s, fitvars);
291 
292  // Run fits!
293  fit9c.Fit(calc4f);
294  ResetAngles(calc4f);
295  fit9s.Fit(calc4f);
296  ResetAngles(calc4f);
297 
298  // Set up values for slice and surface plots
299  AngleValues avals;
300  avals.nbins23 = 200;
301  avals.min23 = 20.;
302  avals.max23 = 70.;
303  avals.nbins34 = 200;
304  avals.min34 = 0.;
305  avals.max34 = 50.;
306  avals.nbins24 = 180;
307  avals.min24 = 0.;
308  avals.max24 = 45.;
309 
310  std::string fullname, special;
311 
312  const Color_t kStatsFitCol = kBlack;
313 
314  const Color_t kFitColor = kRed;
315 
316  fullname = "34";
317  special = "statsC";
318  // Plot 1D fits for 9e20 experiment: Counting
319  TH1* h34C = Plot1DSlice(&multi9c, calc4f,
324  avals.nbins34, avals.min34, avals.max34,
325  rootF, fullname, special, kStatsFitCol
326  );
327 
328  special = "statsS";
329  // Plot 1D fits for 9e20 experiment: Sample
330  TH1* h34S = Plot1DSlice(&multi9s, calc4f,
334  &kFitDelta24InPiUnitsSterile},
335  avals.nbins34, avals.min34, avals.max34,
336  rootF, fullname, special, kStatsFitCol
337  );
338 
339  CompareSlices(h34C, h34S, "34", rootF,
340  avals.min34, avals.max34);
341 
342  // Lower the number of bins for surfaces
343  avals.nbins23 = 50;
344  avals.nbins34 = 50;
345  avals.nbins24 = 45;
346 
347 
348  special = "34vs24";
349  // Plot the surfaces for 9e20 counting experiment
350  Plot2DSlice(&multi9c, &multi9s, calc4f,
352  avals.nbins34, avals.min34, avals.max34,
354  avals.nbins24, avals.min24, avals.max24,
355  rootF, special,
356  {&kFitTheta23InDegreesSterile,
357  &kFitDelta24InPiUnitsSterile}
358  );
359 
360  } // if(plotSensitivities)
361  rootF->Close(); // Close the output root file
362 }
363 
364 //........................................................................
366  const IFitVar* var, std::vector<const IFitVar*> profVars,
367  int nbins, double min, double max,
368  TDirectory* rootOut, std::string indices,
369  std::string special, Color_t color,
370  std::vector<const ISyst*> systs)
371 {
372  // Create slice histogram with/without profiling based on input profVars vector size
373 
374  TH1* h = (profVars.size() > 0 ?
375  Profile(expt, calc, var, nbins, min, max, -1, profVars) :
376  Slice(expt, calc, var, nbins, min, max));
377 
378  ResetAngles(calc); // Don't forget to reset calculator...
379 
380 
381  std::string fullname = indices + special;
382  // Formatting
383  h->SetLineColor(color);
384  h->SetMaximum(5);
385  h->SetMinimum(0);
386  h->SetName(fullname.c_str()); // Object name for saving
387  rootOut->WriteTObject(h); // Save raw histogram
388 
389 
390  std::cout << fullname << std::endl;
391  std::string cname = "c1D" + fullname;
392  std::string ctitle = "Theta" + indices + "_ChiSq1DSlice_" + special;
393  std::string htitle = "#theta_{"+indices+"} (deg.)";
394  h->GetXaxis()->SetTitle(htitle.c_str());
395 
396  // Actually draw things together
397  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
398  h->GetXaxis()->SetRangeUser(min+0.25,max);
399  h->Draw("l");
400  TexInfo();
401  SigmaLines(min, max);
402  PlotText(0.135, 0.82, 0.055);
403  TLegend* leg = 0;
404  // Special case legend position for each angle slice
405  if(indices.compare("34") == 0)
406  leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
407  if(indices.compare("24") == 0)
408  leg = MakeLegend(0.600, 0.35, 0.900, 0.50, 0.037);
409  if(indices.compare("23") == 0)
410  leg = MakeLegend(0.350, 0.35, 0.650, 0.50, 0.037);
411  if(special.compare("stats") == 0)
412  leg->AddEntry(h,"stat. only","L");
413  if(special.compare("systs") == 0)
414  leg->AddEntry(h,"syst.","L");
415  leg->Draw();
416 
417  Preliminary();
418  CenterTitles(h);
419  rootOut->WriteTObject(c);
420 
421  return h;
422 }
423 
424 
425 void CompareSlices(TH1* stats, TH1* systs, std::string indices, TDirectory* rootOut, double min, double max)
426 {
427 
428  std::string cname = "c1DComp" + indices;
429  std::string ctitle = "Theta" + indices + " ChiSq1Dcomp";
430  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
431  stats->GetXaxis()->SetRangeUser(min+0.25,max);
432  systs->GetXaxis()->SetRangeUser(min+0.25,max);
433  stats->Draw("l");
434  systs->Draw("l same");
435  TexInfo();
436  SigmaLines(min, max);
437  Preliminary();
438  PlotText(0.135, 0.82, 0.055);
439  TLegend* leg = 0;
440  // Special case legend for angle slice
441  if(indices.compare("34") == 0)
442  leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
443  if(indices.compare("24") == 0)
444  leg = MakeLegend(0.600, 0.35, 0.900, 0.50, 0.037);
445  if(indices.compare("23") == 0)
446  leg = MakeLegend(0.350, 0.35, 0.650, 0.50, 0.037);
447  leg->AddEntry(stats,"stat. [rate]","L");
448  //leg->AddEntry(systs,"syst.","L");
449  leg->AddEntry(systs,"stat. [shape]","L");
450  leg->Draw();
451  rootOut->WriteTObject(c);
452 
453  return;
454 }
455 
456 
458  IExperiment* expt2,
459  osc::OscCalcSterile* calc,
460  const IFitVar* varA,
461  int nbinsA, double minA, double maxA,
462  const IFitVar* varB,
463  int nbinsB, double minB, double maxB,
464  TDirectory* rootOut, std::string indices,
465  std::vector<const IFitVar*> vars,
466  std::vector<const ISyst*> systs)
467 {
468 
469  ResetAngles(calc);
470  std::string cname = "c2D" + indices;
471  std::string cnamestat = "c2D" + indices + "statC";
472  std::string cnamesyst = "c2D" + indices + "statS";
473  std::string hnamestat = "h" + indices + "statC";
474  std::string hnamesyst = "h" + indices + "statS";
475  std::string ctitle = "Surface for Theta" + indices;
476  std::string ctitlestat = "Surface for Theta" + indices + " statC";
477  std::string ctitlesyst = "Surface for Theta" + indices + " systS";
478 
479  FrequentistSurface surfStat(expt1, calc,
480  varA, nbinsA, minA, maxA,
481  varB, nbinsB, minB, maxB,
482  vars);
483  ResetAngles(calc);
484 
485  FrequentistSurface surfSyst(expt2, calc,
486  varA, nbinsA, minA, maxA,
487  varB, nbinsB, minB, maxB,
488  vars);
489  ResetAngles(calc);
490 
491  TH1* hSurfStat = surfStat.ToTH2(); // Create a histogram from the surfaces
492  hSurfStat->SetName(hnamestat.c_str()); // Set the name for saving
493  rootOut->WriteTObject(hSurfStat); // Save the raw histogram
494 
495  TH1* hSurfSyst = surfSyst.ToTH2(); // Create a histogram from the surfaces
496  hSurfSyst->SetName(hnamesyst.c_str()); // Set the name for saving
497  rootOut->WriteTObject(hSurfSyst); // Save the raw histogram
498 
499  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str());
500 
501  surfStat.DrawContour(Gaussian68Percent1D(surfStat), kDashed, kAzure+2);
502  surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid, kAzure+2);
503  surfSyst.DrawContour(Gaussian68Percent1D(surfSyst), kDashed, kBlack);
504  surfSyst.DrawContour(Gaussian90Percent1D(surfSyst), kSolid, kBlack);
505 
506  PlotText(0.135, 0.82, 0.055);
507  TLegend* legcomp = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
508  TH1* dummy = new TH1F("","",1,0,1);
509 
510  dummy->SetLineColor(kAzure+2);
511  dummy->SetLineStyle(kDashed);
512  dummy->SetLineWidth(2);
513  legcomp->AddEntry(dummy->Clone(),"68% C.L. (stat.) [rate]","L");
514  dummy->SetLineColor(kAzure+2);
515  dummy->SetLineStyle(kSolid);
516  legcomp->AddEntry(dummy->Clone(),"90% C.L. (stat.) [rate]","L");
517  dummy->SetLineColor(kBlack);
518  dummy->SetLineStyle(kDashed);
519  legcomp->AddEntry(dummy->Clone(),"68% C.L. (stat.) [shape]","L");
520  dummy->SetLineColor(kBlack);
521  dummy->SetLineStyle(kSolid);
522  legcomp->AddEntry(dummy->Clone(),"90% C.L. (stat.) [shape]","L");
523  legcomp->Draw();
524  Preliminary();
525  rootOut->WriteTObject(c);
526 
527  // Plot statisics only
528  TCanvas* cStat = new TCanvas(cnamestat.c_str(), ctitlestat.c_str());
529 
530  surfStat.DrawContour(Gaussian68Percent1D(surfStat), kDashed, kAzure+2);
531  surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid, kAzure+2);
532 
533  PlotText(0.135, 0.82, 0.055);
534  TLegend* legStat = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
535 
536  dummy->SetLineColor(kAzure+2);
537  dummy->SetLineStyle(kDashed);
538  dummy->SetLineWidth(2);
539  legStat->AddEntry(dummy->Clone(),"68% C.L. (stat.) [rate]","L");
540  dummy->SetLineColor(kAzure+2);
541  dummy->SetLineStyle(kSolid);
542  legStat->AddEntry(dummy->Clone(),"90% C.L. (stat.) [rate]","L");
543  legStat->Draw();
544  Preliminary();
545  rootOut->WriteTObject(cStat);
546 
547  // Plot stats. + systs. only
548  TCanvas* cSyst = new TCanvas(cnamesyst.c_str(), ctitlesyst.c_str());
549 
550  surfSyst.DrawContour(Gaussian68Percent1D(surfSyst), kDashed, kBlack);
551  surfSyst.DrawContour(Gaussian90Percent1D(surfSyst), kSolid, kBlack);
552 
553  PlotText(0.135, 0.82, 0.055);
554  TLegend* legSyst = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
555  dummy->SetLineColor(kBlack);
556  dummy->SetLineStyle(kDashed);
557  dummy->SetLineWidth(2);
558  legSyst->AddEntry(dummy->Clone(),"68% C.L. (stat.) [shape]","L");
559  dummy->SetLineColor(kBlack);
560  dummy->SetLineStyle(kSolid);
561  legSyst->AddEntry(dummy->Clone(),"90% C.L. (stat.) [shape]","L");
562  legSyst->Draw();
563  Preliminary();
564  rootOut->WriteTObject(cSyst);
565 
566  return;
567 }
568 
569 void PlotText(const double xpos, const double start, const double step)
570 {
571  TLatex* tex = new TLatex();
572  tex->SetNDC();
573  tex->SetTextFont(42);
574  tex->SetTextSize(0.037);
575  tex->SetTextAlign(11);
576  tex->DrawLatex(xpos, start, "NOvA 9.0#times10^{20} POT");
577  tex->DrawLatex(xpos, start - 1*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
578  tex->DrawLatex(xpos, start - 2*step, "#theta_{13} = 8.5#circ, ^{}#theta_{23} = 45#circ");
579  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{41} = 0.5 eV^{2}");
580 }
581 
582 void TexInfo()
583 {
584  const double xtex = 0.76;
585  const double y68 = 0.28;
586  const double y90 = 0.55;
587 
588  TLatex *tex = new TLatex();
589  tex->SetNDC();
590  tex->SetTextFont(42);
591  tex->SetTextSize(0.037);
592  tex->SetTextAlign(11);
593  tex->DrawLatex(xtex, y68, "68% C.L."); // 1-sigma
594  tex->DrawLatex(xtex, y90, "90% C.L."); // 2-sigma
595 }
596 
597 void SigmaLines(double xmin, double xmax)
598 {
599  TLine* line68 = new TLine();
600  TLine* line90 = new TLine();
601  // Fixes horizontal effects at max canvas value
602  TLine* line05 = new TLine();
603  // Dashed line for 1 and 2 sigma lines
604  line68->SetLineStyle(kDashed);
605  line90->SetLineStyle(kDashed);
606 
607  line68->SetLineWidth(2);
608  line90->SetLineWidth(2);
609  line05->SetLineWidth(2);
610  line68->DrawLine(xmin, 1, xmax, 1);
611  line90->DrawLine(xmin, 2.71, xmax, 2.71);
612  line05->DrawLine(xmin, 5, xmax, 5);
613 }
614 
615 
616 //........................................................................
617 void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
619  Spectrum* sterile, Spectrum* sterileNC)
620 {
621  // The input is in units of 1e20, turn this into the real scale
622  double realPOT = (double)POT*1e20;
623 
624  char potBuff[2];
625  sprintf(potBuff, "%d", POT);
626  std::string potStr = std::string(potBuff) + " x 10^{20} POT";
627 
628  // Get each histogram from the Spectrum array
629  TH1* hAll = spectra[0].ToTH1(realPOT);
630  TH1* hNC = spectra[1].ToTH1(realPOT);
631  TH1* hNumu = spectra[2].ToTH1(realPOT);
632  TH1* hNue = spectra[3].ToTH1(realPOT);
633  TH1* hNutau = spectra[4].ToTH1(realPOT);
634  TH1* hCos = spectra[5].ToTH1(realPOT);
635  // The default for these inputs are nullptr, only get a TH1 if they exist
636  TH1* hStrl = (sterile ? sterile ->ToTH1(realPOT) : nullptr);
637  TH1* hStrlNC = (sterileNC ? sterileNC->ToTH1(realPOT) : nullptr);
638 
639  // Create the title and axis labels, which is an empty title,
640  // the x label from an input Spectrum, and Events/POT
641  std::string xLabel("Visible Energy (GeV)");
642  std::string yLabel = (det.compare("ND") == 0 ? "10^{3} Events" : "Events");
643  yLabel += "/" + potStr + " / 0.25 GeV";
644  std::string histTitle = ";" + xLabel + ";" + yLabel;
645 
646  // Print the column title and headers
647  fprintf(textOFS, "Event numbers at the %2.2s for %2de20 POT:\n", det.c_str(), POT);
648  if(hStrlNC) { fprintf(textOFS, "%12.12s, ", "NC (Sterile)"); }
649  else { fprintf(textOFS, "%14.14s", ""); }
650  fprintf(textOFS, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
651  "NC (3 Flav)", "All", "Numu", "Nue", "Nutau", "Cosmic", "FOM");
652 
653  // Get the integral number of events, then print them
654  double nNC3F = hNC ->Integral();
655  double nAll = hAll ->Integral();
656  double nNumu = hNumu->Integral();
657  double nNue = hNue ->Integral();
658  double nNutau = hNutau->Integral();
659  double nCos = hCos ->Integral();
660  double fom = nNC3F/sqrt(nAll);
661  if(hStrlNC) {
662  double nNCSt = hStrlNC->Integral();
663  fprintf(textOFS, "%10E, ", nNCSt);
664  }
665  else {
666  fprintf(textOFS, "%14.14s", "");
667  }
668  fprintf(textOFS, "%10E, %10E, %10E, %10E, %10E, %10E, %10E\n\n",
669  nNC3F, nAll, nNumu, nNue, nNutau, nCos, fom);
670 
671  // Stylize the unstacked plots
672  CenterTitles(hNue);
673  CenterTitles(hNumu);
674  CenterTitles(hNutau);
675  CenterTitles(hNC);
676  CenterTitles(hCos);
677  CenterTitles(hAll);
678 
679  hAll ->SetMinimum(0);
680  hNC ->SetMinimum(0);
681  hNumu->SetMinimum(0);
682  hNutau->SetMinimum(0);
683  hNue ->SetMinimum(0);
684  hCos ->SetMinimum(0);
685  if(hStrlNC) { hStrl->SetMinimum(0); }
686 
687  const Color_t kCosmicBackgroundColor = kOrange+7;
688  hAll ->SetLineColor(kTotalMCColor);
689  hNC ->SetLineColor(kNCBackgroundColor);
690  hNumu->SetLineColor(kNumuBackgroundColor);
691  hNutau->SetLineColor(kBeamNueBackgroundColor);
692  hNue ->SetLineColor(kNueSignalColor);
693  hCos ->SetLineColor(kCosmicBackgroundColor);
694  if(hStrlNC) {
695  hStrlNC->SetLineColor(kNCBackgroundColor);
696  hStrlNC->SetLineStyle(2);
697  }
698 
699  hAll ->SetTitle(histTitle.c_str());
700  hNC ->SetTitle(histTitle.c_str());
701  hNumu->SetTitle(histTitle.c_str());
702  hNue ->SetTitle(histTitle.c_str());
703  hNutau ->SetTitle(histTitle.c_str());
704  hCos ->SetTitle(histTitle.c_str());
705  if(hStrlNC) { hStrlNC->SetTitle(histTitle.c_str()); }
706 
707  // Create "stacked" copies of flavored histograms
708  // Cos Stack = Cos
709  TH1* hCosStack = (TH1*)hCos->Clone();
710 
711  // Numu Stack = Cos Stack + Numu = Cos + Numu
712  TH1* hNumuStack = (TH1*)hCosStack->Clone();
713  hNumuStack->Add(hNumu);
714 
715  // Nue Stack = Numu Stack + Nue = Cos + Numu + Nue
716  TH1* hNueStack = (TH1*)hNumuStack->Clone();
717  hNueStack->Add(hNue);
718 
719  //TH1* hNutauStack = (TH1*)hNueStack->Clone();
720  //hNutauStack->Add(hNutau);
721 
722  // NC Stack = Nue Stack + NC = Cos + Numu + Nue + NC
723  TH1* hNCStack = (TH1*)hNueStack->Clone();
724  hNCStack->Add(hNC);
725 
726  // Stylize the stacked plots
727  //CenterTitles(hNutauStack);
728  CenterTitles(hNueStack);
729  CenterTitles(hNumuStack);
730  CenterTitles(hNCStack);
731  CenterTitles(hCosStack);
732  if(hStrlNC) { CenterTitles(hStrlNC); }
733 
734  hNCStack ->SetLineColor(kNCBackgroundColor);
735  hNumuStack->SetLineColor(kNumuBackgroundColor);
736  hNumuStack->SetFillColor(kNumuBackgroundColor);
737  hNueStack ->SetLineColor(kNueSignalColor);
738  hNueStack ->SetFillColor(kNueSignalColor);
739  //hNutauStack ->SetFillColor(kBeamNueBackgroundColor);
740  //hNutauStack ->SetLineColor(kBeamNueBackgroundColor);
741  hCosStack ->SetLineColor(kCosmicBackgroundColor);
742  hCosStack ->SetFillColor(kCosmicBackgroundColor);
743 
744  hNCStack ->SetMinimum(0);
745  hNumuStack->SetMinimum(0);
746  hNueStack ->SetMinimum(0);
747  //hNutauStack ->SetMinimum(0);
748  hCosStack ->SetMinimum(0);
749 
750  // Manually set errors, and calculate the max value (including errors)
751  double maxval = 0.;
752  double maxvalstack = 0.;
753  for(int i = 1, n = hNC->GetNbinsX(); i <= n; ++i) {
754  double error = sqrt(hNC->GetBinContent(i));
755  // hNC ->SetBinError(i, error);
756  //hNCStack->SetBinError(i, error);
757  maxval = std::max(maxval, std::max(hAll->GetBinContent(i), hNC->GetBinContent(i) + error));
758  maxvalstack = std::max(maxvalstack, hNCStack->GetBinContent(i) + error);
759  }
760 
761  // Scale down the ND by 10e3 (don't forget the max values!)
762  if(det.compare("ND") == 0) {
763  double ndScale = 1./1000.;
764  hAll ->Scale(ndScale);
765  hNC ->Scale(ndScale);
766  hNue ->Scale(ndScale);
767  hNumu->Scale(ndScale);
768  hNCStack ->Scale(ndScale);
769  hNueStack ->Scale(ndScale);
770  hNumuStack->Scale(ndScale);
771  maxval *= ndScale;
772  maxvalstack *= ndScale;
773  }
774 
775  TLatex *tex = new TLatex();
776  tex->SetNDC();
777  tex->SetTextFont(42);
778  tex->SetTextSize(0.037);
779  tex->SetTextAlign(11);
780 
781  // Create a canvas for and draw the unstacked histograms
782  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
783  if(det.compare("FD") == 0){
784  hNutau->SetMaximum(maxval*1.1);
785  hNutau->GetXaxis()->SetNdivisions(-405);
786  }
787  else{
788  hNue->SetMaximum(maxval*1.1);
789  hNue->GetXaxis()->SetNdivisions(-405);
790  }
791 
792  // Create a legend for the unstacked histograms
793  double xL = 0.5, xR = 0.7, yB = 0.59, yT = 0.84;
794  if(det.compare("ND") == 0) { yB += 0.05; }
795  TLegend* leg = new TLegend(xL, yB, xR, yT);
796  leg->SetBorderSize(0);
797  leg->SetFillColor(kWhite);
798  leg->SetFillStyle(0);
799  leg->SetFillStyle(0);
800  leg->SetTextFont(42);
801  leg->SetTextSize(0.037);
802  leg->AddEntry(hAll, "Total Prediction", "l");
803  leg->AddEntry(hNC, "NC 3 Flavor Prediction", "le");
804  leg->AddEntry(hNue, "#nu_{e} CC Background", "l");
805  leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
806  if(det.compare("FD") == 0) {
807  leg->AddEntry(hNutau,"#nu_{#tau} CC Background", "l");
808  leg->AddEntry(hCos, "Cosmic Background", "l");
809  }
810  if(det.compare("FD") == 0)
811  hNutau->Draw("hist");
812  else
813  hNue ->Draw("hist");
814  hNumu->Draw("hist same");
815  if(det.compare("FD") == 0) { hCos->Draw("hist same"); }
816  hAll ->Draw("hist same");
817  hNC ->Draw("same");
818  leg->Draw();
819  if(det.compare("FD") == 0) {
820  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, ^{}#theta_{24} = 4#circ, ^{}#theta_{34} = 10#circ");
821  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.44x10^{-3} eV^{2}, ^{}#theta_{13} = 8.5#circ, ^{}#theta_{23} = 45#circ");
822  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
823  }
824 
825  gPad->RedrawAxis();
826  Simulation();
827 
828  // Save the canvas
829  rootOut->WriteTObject(c);
830 
831  // Create a canvas for and draw the stacked histograms
832  TCanvas* cStack = new TCanvas((name + "Stack").c_str(), title.c_str(), 800, 500);
833  //if(det.compare("FD") == 0){
834  //hNutauStack->SetMaximum(maxval*1.1);
835  //hNutauStack->GetXaxis()->SetNdivisions(-405);
836  //}
837  //else{
838  hNueStack->SetMaximum(maxval*1.1);
839  hNueStack->GetXaxis()->SetNdivisions(-405);
840  //}
841 
842  if(det.compare("ND") == 0) { yB += 0.05; }
843  TLegend* legStack = new TLegend(xL, yB, xR, yT);
844  legStack->SetBorderSize(0);
845  legStack->SetFillColor(kWhite);
846  legStack->SetFillStyle(0);
847  legStack->SetFillStyle(0);
848  legStack->SetTextFont(42);
849  legStack->SetTextSize(0.037);
850  legStack->AddEntry(hNCStack, "NC 3 Flavor Prediction", "le");
851  if(hStrlNC) { legStack->AddEntry(hStrlNC, "NC 4 Flavor Prediction", "l"); }
852  legStack->AddEntry(hNueStack, "#nu_{e} CC Background", "f");
853  legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
854  if(det.compare("FD") == 0) {
855  //legStack->AddEntry(hNutauStack, "#nu_{#tau} CC Background", "f");
856  legStack->AddEntry(hCosStack, "Cosmic Background", "f");
857  }
858 
859  //if(det.compare("FD") == 0) {
860  //hNutauStack->Draw("hist");
861  //}
862  //else{
863  hNueStack ->Draw("hist");
864  //}
865  hNumuStack->Draw("hist same");
866  if(det.compare("FD") == 0) { hCosStack->Draw("hist same"); }
867  hNCStack ->Draw("same");
868  if(hStrlNC) { hStrlNC->Draw("hist same"); }
869  legStack->Draw();
870  if(det.compare("FD") == 0) {
871  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, ^{}#theta_{24} = 4#circ, ^{}#theta_{34} = 10#circ");
872  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.44x10^{-3} eV^{2}, ^{}#theta_{13} = 8.5#circ, ^{}#theta_{23} = 45#circ");
873  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
874  }
875 
876  gPad->RedrawAxis();
877  Simulation();
878 
879  rootOut->WriteTObject(cStack);
880 
881  return;
882 }
883 
884 //........................................................................
885 void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
887 {
888  Spectrum empty = decomp->NCTotalComponent();
889  empty.Clear();
890 
891  Spectrum spectraND[MAXSPEC] = {
892  decomp->NCTotalComponent() +
893  decomp->NueComponent() + decomp->AntiNueComponent() +
894  decomp->NumuComponent() + decomp->AntiNumuComponent(),
895  decomp->NCTotalComponent(),
896  decomp->NumuComponent() + decomp->AntiNumuComponent(),
897  decomp->NueComponent() + decomp->AntiNueComponent(),
898  empty,
899  empty
900  };
901 
902  std::string det = "ND";
903  std::string fulltitle = det + title;
904  PlotStack(spectraND, rootOut, textOFS, name, fulltitle, det, POT);
905 
906  return;
907 }
908 
909 //........................................................................
910 void PlotStack(IPrediction* pred, Spectrum& sCos,
911  osc::IOscCalc* calc, osc::IOscCalc* calcSt,
912  TDirectory* rootOut, FILE* textOFS,
914 {
915  Spectrum spectraFD[MAXSPEC] = {
916  pred->Predict(calc) + sCos,
921  sCos
922  };
923 
924  Spectrum sterile = pred->Predict(calcSt) + sCos;
925  Spectrum sterileNC = pred->PredictComponent(calcSt, Flavors::kAll, Current::kNC, Sign::kBoth);
926 
927  std::string det = "FD";
928  std::string fulltitle = det + title;
929  PlotStack(spectraFD, rootOut, textOFS, name, fulltitle, det, POT, &sterile, &sterileNC);
930 
931  return;
932 }
933 
934 //........................................................................
935 void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
936  Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
937  osc::IOscCalc* calc, Spectrum& sCos,
938  TDirectory* rootOut, FILE* textOFS,
940 {
941  Spectrum sFDNum = pFDNum ->Predict(calc);
942  Spectrum sFDPurDen = pFDPurDen->Predict(calc) + sCos;
943  Spectrum sFDEffDen = pFDEffDen->Predict(calc);
944 
945  TH1* hFDNum = sFDNum .ToTH1(NCSCALE);
946  TH1* hFDPurDen = sFDPurDen.ToTH1(NCSCALE);
947  TH1* hFDEffDen = sFDEffDen.ToTH1(NCSCALE);
948 
949  TH1* hNDNum = sNDNum .ToTH1(NCSCALE);
950  TH1* hNDPurDen = sNDPurDen.ToTH1(NCSCALE);
951  TH1* hNDEffDen = sNDEffDen.ToTH1(NCSCALE);
952 
953  TH1* hFDEff = (TH1*)hFDNum->Clone();
954  hFDEff->Divide(hFDEffDen);
955 
956  TH1* hNDEff = (TH1*)hNDNum->Clone();
957  hNDEff->Divide(hNDEffDen);
958 
959  TH1* hFDPur = (TH1*)hFDNum->Clone();
960  hFDPur->Divide(hFDPurDen);
961 
962  TH1* hNDPur = (TH1*)hNDNum->Clone();
963  hNDPur->Divide(hNDPurDen);
964 
965  hFDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
966  hFDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
967  hNDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
968  hNDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
969 
970  hFDEff->SetLineColor(kBlue);
971  hNDEff->SetLineColor(kBlue);
972 
973  hFDPur->SetLineColor(kRed);
974  hNDPur->SetLineColor(kRed);
975 
976  hNDEff->SetLineStyle(2);
977  hNDPur->SetLineStyle(2);
978 
979  hFDEff->SetMaximum(1.);
980  hFDEff->SetMinimum(0.);
981  hFDEff->GetXaxis()->SetNdivisions(-405);
982 
983  hNDEff->SetMaximum(1.);
984  hNDEff->SetMinimum(0.);
985 
986  hFDPur->SetMaximum(1.);
987  hFDPur->SetMinimum(0.);
988 
989  hNDPur->SetMaximum(1.);
990  hNDPur->SetMinimum(0.);
991 
992  CenterTitles(hFDEff);
993  CenterTitles(hNDEff);
994  CenterTitles(hFDPur);
995  CenterTitles(hNDPur);
996 
997  TLegend* leg = new TLegend(0.6, 0.35, 0.85, 0.55);
998  leg->AddEntry(hFDEff, "FD MC Efficiency", "L");
999  leg->AddEntry(hNDEff, "ND MC Efficiency", "L");
1000  leg->AddEntry(hFDPur, "FD MC Purity", "L");
1001  leg->AddEntry(hNDPur, "ND MC Purity", "L");
1002 
1003  TCanvas* c = new TCanvas(name.c_str(), "Purity and Efficiency", 800, 500);
1004  hFDEff->Draw("hist");
1005  hNDEff->Draw("hist same");
1006  hFDPur->Draw("hist same");
1007  hNDPur->Draw("hist same");
1008  leg->Draw();
1009 
1010  gPad->RedrawAxis();
1011  Simulation();
1012 
1013  rootOut->WriteTObject(c);
1014 }
1015 
1016 
1017 //........................................................................
1019 {
1021  calc->SetDm(4, 0.5); // #deltam41 = 0.5 eV^2
1022  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
1023  calc->SetAngle(3, 4, 0.);
1024  calc->SetAngle(2, 4, 0.);
1025 
1026  return;
1027 }
void PlotText(const double xpos, const double start, const double step)
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string det, double POT, double potEquiv, PlotOptions opt, Spectrum *dataspec)
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Double_t angle
Definition: plot.C:86
TH2 * Gaussian90Percent1D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 1D in gaussian approximation.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
std::map< std::string, double > xmax
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual Spectrum AntiNueComponent() const =0
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Float_t y1[n_points_granero]
Definition: compare.C:5
void PlotPurEff(Spectrum pureff[], Spectrum effcyD[], TDirectory *out, std::string name, std::string title)
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
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:209
virtual Spectrum NumuComponent() const =0
#define MAXSPEC
Definition: NusLoadProd3.h:18
A simple Gaussian constraint on an arbitrary IFitVar.
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:237
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
static std::unique_ptr< PredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
void Plot2DSlice(IExperiment *expt1, IExperiment *expt2, osc::OscCalcSterile *calc, const IFitVar *varA, int nbinsA, double minA, double maxA, const IFitVar *varB, int nbinsB, double minB, double maxB, TDirectory *rootOut, std::string angles, std::vector< const IFitVar * > vars, std::vector< const ISyst * > systs={})
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
int stats(TString inputFilePath, Int_t firstRun, Int_t lastRun, Float_t thresh, TString myDet)
Definition: stats.C:13
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 Clear()
Definition: Spectrum.cxx:433
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
osc::OscCalcDumb calc
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
const Color_t kNumuBackgroundColor
Definition: Style.h:30
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:321
#define M_PI
Definition: SbMath.h:34
static std::unique_ptr< ProportionalDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const Color_t kNueSignalColor
Definition: Style.h:19
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
fclose(fg1)
Log-likelihood scan across two parameters.
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
void ResetAngles(osc::OscCalcSterile *calc)
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:64
const std::map< std::pair< std::string, std::string >, Variable > vars
void PlotNus17Prediction(bool plotSensitivities=false)
static std::unique_ptr< PredictionExtrap > LoadFrom(TDirectory *dir, const std::string &name)
TGraph * Slice(const IExperiment *expt, osc::IOscCalcAdjustable *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
virtual Spectrum AntiNumuComponent() const =0
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
void Preliminary()
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TH1 * Plot1DSlice(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var, std::vector< const IFitVar * > vars, int nbins, double min, double max, TDirectory *rootOut, std::string name, std::string special, Color_t color, std::vector< const ISyst * > systs={})
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Splits Data proportionally according to MC.
void SetDm(int i, double dm)
TLegend * MakeLegend(double x1, double y1, double x2, double y2, double textsize=0.03)
void TexInfo()
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
TLatex * tex
Definition: f2_nu.C:499
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
Base class defining interface for experiments.
Definition: IExperiment.h:14
Neutral-current interactions.
Definition: IPrediction.h:40
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Definition: IFitVar.h:16
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
A prediction object compatible with sterile oscillations.
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
static double NCSCALE
const Color_t kTotalMCColor
Definition: Style.h:16
Take the output of an extrapolation and oscillate it as required.
All neutrinos, any flavor.
Definition: IPrediction.h:26
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void SigmaLines(double xmin, double xmax)
const Color_t kNCBackgroundColor
Definition: Style.h:22
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
Compare a data spectrum to MC expectation using only the event count.
virtual Spectrum NueComponent() const =0
TH2 * Gaussian68Percent1D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 1D in gaussian approximation.
void CompareSlices(TH1 *stats, TH1 *systs, std::string angle, TDirectory *rootOut, double xmin, double xmax)
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17