BlessedPlotsAnaByPeriod.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"
23 #include "CAFAna/Extrap/ExtrapSterile.h"
29 
30 #include "OscLib/OscCalcSterile.h"
31 
33 
34 using namespace ana;
35 
36 // Definitions for organization
37 
38 // One struct to contain all values for angle fits
39 struct AngleValues{
40  int nbins23;
41  double min23;
42  double max23;
43 
44  int nbins34;
45  double min34;
46  double max34;
47 
48  int nbins24;
49  double min24;
50  double max24;
51 };
52 
53 // Plot a single 1D slice canvas -- this function does the actual plotting
54 // The macro does not call this function directly in the main function
55 // \param expt The CAFAna experiment object, contains FD prediction and "data"
56 // \param calc Sterile calculator for oscillating the spectra
57 // \param var Parameter to fit for!
58 // \param profVars Vector of parameters to profile over
59 // \param rootOut Already open file for root output
60 // \param name Base of the string for creating the name for saving the root objects
62  const IFitVar* var, std::vector<const IFitVar*> profVars,
63  TDirectory* rootOut, int nbins, double min, double max,
64  std::string name, Color_t color);
65 
66 // Plot 1D slice canvases
67 // This function sets up multiple calls to the single canvas version of Plot1DSlice
68 // \param expt The CAFAna experiment object, contains FD prediction and "data"
69 // \param calc Sterile calculator for oscillating the spectra
70 // \param rootOut Already open file for root output
71 // \param a Struct containing bin values for slice histograms
73  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
74  TDirectory* rootOut, AngleValues a,
75  Color_t color);
76 
77 // Plot 2D slice surfaces
78 // This function sets up multiple calls to the workhorse version of Plot2DSlices
79 // \param expt The CAFAna experiment object, contains FD prediction and "data"
80 // \param calc Sterile calculator for oscillating the spectra
81 // \param rootOut Already open file for root output
82 // \param a Struct containing bin values for slice histograms
84  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
85  TDirectory* rootOut, AngleValues a);
86 
87 // Plot 2D slice surfaces
88 // The macro does not call this function directly
89 // \param rootOut Already open file for root output
90 // \param prof String indicating whether surfaces are profiled or not (optional)
92  TDirectory* rootOut, std::string prof = "");
93 
94 // Plot spectra and stack canvases -- the workhorse version
95 // The macro does not call this function directly
96 // \param spectra Array of spectra, should match the output of the SPECARR macro
97 // \param rootOut Already open file for root output
98 // \param textOFS Already open filestream for text output
99 // \param name A name for the output, only part of the saved output object's name
100 // \param title Title for the histograms. Typically this is the x axis label
101 // \param det The detector the spectra correspond to, should be 2 chars long for nice formatting
102 // \param sterile Pointer to a spectrum with 4 flavor osc prediction, includes all event types (optional)
103 // \param sterileNC Pointer to a spectrum with 4 flavor osc prediction, only includes NC events (optional)
104 void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
106  Spectrum* sterile = nullptr, Spectrum* sterileNC = nullptr);
107 
108 // Plot spectra and stack canvases for ND
109 // Sets up call to the workhorse version of PlotStack
110 // \param decomp The ND decomposition
111 // \param rootOut Already open file for root output
112 // \param textOFS Already open filestream for text output
113 // \param name A name for the output, only part of the saved output object's name
114 // \param title Title for the histograms. Typically this is the x axis label
115 void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
117 
118 // Plot spectra and stack canvases for FD
119 // Sets up call to the workhorse version of PlotStack
120 // \param pred The FD prediction object
121 // \param sCos The FD cosmic prediction spectrum
122 // \param calc Sterile calculator assuming 3 flavor oscillations
123 // \param calcSt Sterile calculator assuming 4 flavor oscillations
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
129  Spectrum& sCos,
131  osc::IOscCalc* calcSt,
132  TDirectory* rootOut, FILE* textOFS,
134 
135 // Function to plot purity and efficiency at both FD and ND
136 // \param pFDNum Prediction object for FD, used as nominator for FD purity and efficiency
137 // \param pFDPurDen Prediction object for FD, used as denominator for FD purity
138 // \param pFDPurDen Prediction object for FD, used as denominator for FD efficiency
139 // \param sNDNum Spectrum object for ND, used as nominator for ND purity and efficiency
140 // \param sNDPurDen Spectrum object for ND, used as denominator for ND purity
141 // \param sNDPurDen Spectrum object for ND, used as denominator for ND efficiency
142 // \param calc Sterile calculator assuming 3 flavor oscillations
143 // \param rootOut Already open file for root output
144 // \param textOFS Already open filestream for text output
145 // \param name A name for the output, only part of the saved output object's name
146 // \param title Title for the histograms. Typically this is the x axis label
147 void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
148  Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
149  osc::IOscCalc* calc, Spectrum& sCos,
150  TDirectory* rootOut, FILE* textOFS,
152 
153 // After running a fit, reset the values in a calculator to preset defaults
155 
157 {
158  TH1::AddDirectory(0);
159 
160  // Set up path to file with filled CAFAna objects and for output of analysis
161  //std::string folder = "/nova/ana/steriles/Ana01/";
162  std::string folder = "";
163  std::string filenm = "test_veto";
164 
165  std::string loadLocation = folder + filenm + ".root";
166  std::string saveLocation = folder + filenm + "Ana.root";
167  std::string textLocation = folder + filenm + "Ana.txt";
168 
169  // Load filled objects from file
170  TFile* rootL = new TFile(loadLocation.c_str(), "UPDATE");
171  TDirectory* tmpL = gDirectory;
172  TDirectory* loadDir = gDirectory;
173 
174  loadDir->cd((loadLocation + ":/decompNC").c_str());
175  ProportionalDecomp decompNC = *ProportionalDecomp::LoadFrom(gDirectory);
176 
177  loadDir->cd((loadLocation + ":/prediction").c_str());
179 
180  loadDir->cd((loadLocation + ":/sData").c_str());
181  Spectrum sFDData = *Spectrum::LoadFrom(gDirectory);
182 
183  loadDir->cd((loadLocation + ":/sCosSA").c_str());
184  Spectrum sCosmic = *Spectrum::LoadFrom(gDirectory);
185 
186  tmpL->cd();
187  rootL->Close(); // Don't forget to close the file!
188 
189  // Set up oscillation calculator that uses default 3 flavor parameters
191 
192  // Set up oscillation calculator that uses default 3 flavor parameters
194  calc4f->SetNFlavors(4);
195  calc4f->SetDm(4, 0.5); //Dm41 = 0.5 eV^2
196  calc4f->SetAngle(2, 4, 0.);
197  calc4f->SetAngle(3, 4, 0.);
198 
199  std::string labelEReco = "Calorimetric Energy (GeV)";
200 
201  // Open files for saving analysis output
202  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
203  FILE* textF;
204  textF = fopen(textLocation.c_str(), "w");
205 
206  // Choose some angle values for making 4 flavor predictions
207  calc4f->SetAngle(2, 4, 0.17453293); // 10 degrees
208  calc4f->SetAngle(3, 4, 0.61086524); // 35 degrees
209 
210  // Plot spectra for Ana01 prediction
211  std::string titleHelper = " Spectra for 6.05e20 POT-equiv.";
212  PlotStack(&pred, sCosmic, calc3f, calc4f, rootF, textF, "FDSpectra", titleHelper);
213  PlotStack(&decompNC, rootF, textF, "NDSpectra", titleHelper);
214 
215  fclose(textF); // Text file is done now, close it
216  ResetAngles(calc4f);
217 
218  // Create Experiment objects
219  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45., 8.73);
220  CountingExperiment expt(&pred, sFDData, sCosmic);
221  MultiExperiment multi({&th23Constraint, &expt});
222 
223  // Vector of parameters to fit: theta 23, 34, 24
224  std::vector<const IFitVar*> fitvars;
225  fitvars.push_back(&kFitTheta34InDegreesSterile);
226  fitvars.push_back(&kFitTheta24InDegreesSterile);
227  fitvars.push_back(&kFitTheta23InDegreesSterile);
228 
229  // Create fit objects
230  MinuitFitter fit(&multi, fitvars);
231 
232  // Run fits!
233  fit.Fit(calc4f);
234  ResetAngles(calc4f);
235 
236  // Set up values for slice and surface plots
237  AngleValues avals;
238  avals.nbins23 = 200;
239  avals.min23 = 20.;
240  avals.max23 = 70.;
241  avals.nbins34 = 200;
242  avals.min34 = 0.;
243  avals.max34 = 50.;
244  avals.nbins24 = 180;
245  avals.min24 = 0.;
246  avals.max24 = 45.;
247 
248  const Color_t kFitColor = kRed;
249 
250  // Plot 1D fits
251  Plot1DSlices(
252  &multi, calc4f,
256  rootF, avals, kFitColor
257  );
258 
259  // Lower the number of bins for surfaces
260  avals.nbins23 = 50;
261  avals.nbins34 = 50;
262  avals.nbins24 = 45;
263 
264  // Plot the surfaces
265  Plot2DSlices(
266  &multi, calc4f,
270  rootF, avals
271  );
272 
273  rootF->Close(); // Close the output root file
274 }
275 
276 //........................................................................
278  const IFitVar* var, std::vector<const IFitVar*> profVars,
279  TDirectory* rootOut, int nbins, double min, double max,
280  std::string name, Color_t color)
281 {
282  // Create slice histogram with/without profiling based on input profVars vector size
283 
284  TH1* h = (profVars.size() > 0 ?
285  Profile(expt, calc, var, nbins, min, max, -1, profVars) :
286  Slice(expt, calc, var, nbins, min, max));
287 
288  ResetAngles(calc); // Don't forget to reset calculator...
289 
290  // Formatting:
291  h->SetLineColor(color); // Set the line color
292  h->SetMaximum(5); // Want the plot to go from 0 to 5
293  h->SetMinimum(0);
294  h->SetName(name.c_str()); // Set the object name for saving
295  rootOut->WriteTObject(h); // Save the raw histogram
296 
297  // Create an object for writing LaTeX to canvas
298  TLatex *tex = new TLatex();
299  tex->SetNDC();
300  tex->SetTextFont(42);
301  tex->SetTextSize(0.037);
302  tex->SetTextAlign(11);
303 
304  // Coordinates for LaTeX object
305  double xTex = 0.92, y68 = 0.25, y90 = 0.52;
306 
307  // Dotted lines representing 1 and 2 sigma
308  TLine* line68 = new TLine(min, 1, max, 1);
309  TLine* line90 = new TLine(min, 2.71, max, 2.71);
310 
311  // Sometimes histogram has horizontal artifact at max canvas value
312  // This line hides it
313  TLine* line5 = new TLine(min, 5, max, 5);
314 
315  // Dashed line for 1 and 2 sigma lines
316  line68->SetLineStyle(2);
317  line90->SetLineStyle(2);
318 
319  // Formatting
320  line68->SetLineWidth(2);
321  line90->SetLineWidth(2);
322  line5 ->SetLineWidth(2);
323 
324  // Input "name" should have "h" for histogram at the beginning
325  // Make canvas name with input "name" but remove that "h"
326  std::cout << name << std::endl;
327  std::string cname = "c1DSlice" + name.substr(1, name.size()-1);
328 
329  // Theta subscripts should be at position 4, 5, get them
330  std::string indices = name.substr(4, 2);
331 
332  // Create POT string
333  std::string cpotstr = "6.05e20 POT-equiv.";
334 
335  // Create a title for canvas
336  std::string ctitle = "Delta chi^2 for Theta" + indices;
337 
338  // Add info on profiling if applicable
339  ctitle += " (Profiling Over Other Angles) for " + cpotstr;
340 
341  // Create the canvas and draw!
342  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
343  h->Draw("]["); // Don't show vertical drop to 0 at edge of histogram
344  line68->DrawClone("same");
345  line90->DrawClone("same");
346  line5 ->DrawClone("same");
347  tex->DrawLatex(xTex, y68, "68%"); // Add label for 1 sigma line
348  tex->DrawLatex(xTex, y90, "90%"); // Add label for 2 sigma line
349  Simulation();
350  rootOut->WriteTObject(c); // Save full canvas
351 
352  return;
353 }
354 
355 //........................................................................
357  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
358  TDirectory* rootOut, AngleValues a,
359  Color_t color)
360 {
361  // General form of the object name for saving
362  std::string name = "hTh";
363  std::string fullname; // fullname adds 2 characters, the theta subscript indices
364 
365  // Plot unprofiled slices, for theta 23, then 34, then 24
366  fullname = name + "23";
367  Plot1DSlice(expt, calc, var23, {},
368  rootOut, a.nbins23, a.min23, a.max23, fullname, color);
369  fullname = name + "34";
370  Plot1DSlice(expt, calc, var34, {},
371  rootOut, a.nbins34, a.min34, a.max34, fullname, color);
372  fullname = name + "24";
373  Plot1DSlice(expt, calc, var24, {},
374  rootOut, a.nbins24, a.min24, a.max24, fullname, color);
375 
376 
377  // Plot profiled slices in same order as above
378  // When slice is of one angle, the other two are profiled
379  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
380  fullname = name + "23Prof";
381  Plot1DSlice(expt, calc, var23, {var34, var24},
382  rootOut, a.nbins23, a.min23, a.max23, fullname, color);
383  fullname = name + "34Prof";
384  Plot1DSlice(expt, calc, var34, {var23, var24},
385  rootOut, a.nbins34, a.min34, a.max34, fullname, color);
386  fullname = name + "24Prof";
387  Plot1DSlice(expt, calc, var24, {var23, var34},
388  rootOut, a.nbins24, a.min24, a.max24, fullname, color);
389 
390  return;
391 }
392 
393 //........................................................................
395  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
396  TDirectory* rootOut, AngleValues a)
397 {
398  // Create unprofiled surfaces
399  // 34 vs 24
400  FrequentistSurface s3424(expt, calc,
401  var34, a.nbins34, a.min34, a.max34,
402  var24, a.nbins24, a.min24, a.max24);
403  ResetAngles(calc);
404 
405  // 23 vs 24
406  FrequentistSurface s2324(expt, calc,
407  var23, a.nbins23, a.min23, a.max23,
408  var24, a.nbins24, a.min24, a.max24);
409  ResetAngles(calc);
410 
411  // 34 vs 23
412  FrequentistSurface s3423(expt, calc,
413  var34, a.nbins34, a.min34, a.max34,
414  var23, a.nbins23, a.min23, a.max23);
415  ResetAngles(calc);
416 
417  // Create profiled surfaces in same order as above
418  // The third angle is the one that is profiled
419  // I.e., the 34 vs 24 surface profiles 23
420  FrequentistSurface s3424Prof(expt, calc,
421  var34, a.nbins34, a.min34, a.max34,
422  var24, a.nbins24, a.min24, a.max24,
423  {var23});
424  ResetAngles(calc);
425 
426  FrequentistSurface s2324Prof(expt, calc,
427  var23, a.nbins23, a.min23, a.max23,
428  var24, a.nbins24, a.min24, a.max24,
429  {var34});
430  ResetAngles(calc);
431 
432  FrequentistSurface s3423Prof(expt, calc,
433  var34, a.nbins34, a.min34, a.max34,
434  var23, a.nbins23, a.min23, a.max23,
435  {var24});
436  ResetAngles(calc);
437 
438  // Plot all of the surfaces
439  Plot2DSlices(&s3424, &s2324, &s3423, rootOut);
440  Plot2DSlices(&s3424Prof, &s2324Prof, &s3423Prof, rootOut, "Prof");
441 
442  return;
443 }
444 
445 //........................................................................
447  TDirectory* rootOut, std::string prof)
448 {
449  // General format for the name for saving
450  std::string name = "hTh";
451  std::string fullname;
452 
453  // Add indices for first theta, add second theta, and profile string if applicaable
454  fullname = name + "34Th24" + prof;
455  TH1* h3424 = s3424->ToTH2(); // Create a histogram from the surface
456  h3424->SetName(fullname.c_str()); // Set the name for saving
457  rootOut->WriteTObject(h3424); // Save the raw histogram
458 
459  fullname = name + "23Th24" + prof;
460  TH1* h2324 = s2324->ToTH2();
461  h2324->SetName(fullname.c_str());
462  rootOut->WriteTObject(h2324);
463 
464  fullname = name + "34Th23" + prof;
465  TH1* h3423 = s3423->ToTH2();
466  h3423->SetName(fullname.c_str());
467  rootOut->WriteTObject(h3423);
468 
469  // Create name for canvas
470  std::string cname = "c2DSlice" + prof;
471  // Create string for POT based on whether name has 1 or 3 in it
472  std::string potstr = "6.05e20 POT-equiv.";
473  // Create relevant title for the canvas
474  std::string ctitle = (
475  prof.size() > 0 ?
476  "Surfaces (Profiling Over Other Angles) for " + potstr :
477  "Surfaces for " + potstr
478  );
479 
480  // Create triptych canvas and draw!
481  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 800, 800);
482  c->Divide(2, 2); // Split into 2 x 2 grid
483  c->cd(1); // Go to top left
484  s3424->Draw(); // Draw histogram
485  s3424->DrawContour(Gaussian68Percent2D(*s3424), 7, kBlack); // Draw 1 sigma contour
486  s3424->DrawContour(Gaussian90Percent2D(*s3424), 1, kBlack); // Draw 2 sigma contour
487  s3424->DrawBestFit(kBlack); // Draw best fit point
488  c->cd(2); // Go to top right
489  s2324->Draw();
490  s2324->DrawContour(Gaussian68Percent2D(*s2324), 7, kBlack);
491  s2324->DrawContour(Gaussian90Percent2D(*s2324), 1, kBlack);
492  s2324->DrawBestFit(kBlack);
493  Simulation();
494  c->cd(3); // Go to bottom left
495  s3423->Draw();
496  s3423->DrawContour(Gaussian68Percent2D(*s3423), 7, kBlack);
497  s3423->DrawContour(Gaussian90Percent2D(*s3423), 1, kBlack);
498  s3423->DrawBestFit(kBlack);
499  rootOut->WriteTObject(c); // Save the root canvas
500 
501  return;
502 }
503 
504 //........................................................................
505 void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
507  Spectrum* sterile, Spectrum* sterileNC)
508 {
509 
510 
511  int POT = 6.05;
512  std::string potStr ="6.05e20 POT-equiv.";
513 
514  // Get each histogram from the Spectrum array
515  TH1* hAll = spectra[0].ToTH1(spectra[0].POT());
516  TH1* hNC = spectra[1].ToTH1(spectra[1].POT());
517  TH1* hNumu = spectra[2].ToTH1(spectra[2].POT());
518  TH1* hNue = spectra[3].ToTH1(spectra[3].POT());
519  TH1* hCos = spectra[4].ToTH1(spectra[4].POT());
520  // The default for these inputs are nullptr, only get a TH1 if they exist
521  TH1* hStrl = (sterile ? sterile ->ToTH1(sterile->POT()) : nullptr);
522  TH1* hStrlNC = (sterileNC ? sterileNC->ToTH1(sterileNC->POT()) : nullptr);
523 
524  // Create the title and axis labels, which is an empty title,
525  // the x label from an input Spectrum, and Events/POT
526  std::string xLabel(hAll->GetXaxis()->GetTitle());
527  std::string yLabel = (det.compare("ND") == 0 ? "10^{3} Events" : "Events");
528  yLabel += "/" + potStr + " / 0.25 GeV";
529  std::string histTitle = ";" + xLabel + ";" + yLabel;
530 
531  // Print the column title and headers
532  fprintf(textOFS, "Event numbers at the %2.2s for %2de20 POT-equiv.:\n", det.c_str(), POT);
533  if(hStrlNC) { fprintf(textOFS, "%12.12s, ", "NC (Sterile)"); }
534  else { fprintf(textOFS, "%14.14s", ""); }
535  fprintf(textOFS, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
536  "NC (3 Flav)", "All", "Numu", "Nue", "Cosmic", "FOM");
537 
538  // Get the integral number of events, then print them
539  double nNC3F = hNC ->Integral();
540  double nAll = hAll ->Integral();
541  double nNumu = hNumu->Integral();
542  double nNue = hNue ->Integral();
543  double nCos = hCos ->Integral();
544  double fom = nNC3F/sqrt(nAll);
545  if(hStrlNC) {
546  double nNCSt = hStrlNC->Integral();
547  fprintf(textOFS, "%10E, ", nNCSt);
548  }
549  else {
550  fprintf(textOFS, "%14.14s", "");
551  }
552  fprintf(textOFS, "%10E, %10E, %10E, %10E, %10E, %10E\n\n",
553  nNC3F, nAll, nNumu, nNue, nCos, fom);
554 
555  // Stylize the unstacked plots
556  CenterTitles(hNue);
557  CenterTitles(hNumu);
558  CenterTitles(hNC);
559  CenterTitles(hCos);
560  CenterTitles(hAll);
561 
562  hAll ->SetMinimum(0);
563  hNC ->SetMinimum(0);
564  hNumu->SetMinimum(0);
565  hNue ->SetMinimum(0);
566  hCos ->SetMinimum(0);
567  if(hStrlNC) { hStrl->SetMinimum(0); }
568 
569  const Color_t kCosmicBackgroundColor = kOrange+7;
570  hAll ->SetLineColor(kTotalMCColor);
571  hNC ->SetLineColor(kNCBackgroundColor);
572  hNumu->SetLineColor(kNumuBackgroundColor);
573  hNue ->SetLineColor(kNueSignalColor);
574  hCos ->SetLineColor(kCosmicBackgroundColor);
575  if(hStrlNC) {
576  hStrlNC->SetLineColor(kNCBackgroundColor);
577  hStrlNC->SetLineStyle(2);
578  }
579 
580  hAll ->SetTitle(histTitle.c_str());
581  hNC ->SetTitle(histTitle.c_str());
582  hNumu->SetTitle(histTitle.c_str());
583  hNue ->SetTitle(histTitle.c_str());
584  hCos ->SetTitle(histTitle.c_str());
585  if(hStrlNC) { hStrlNC->SetTitle(histTitle.c_str()); }
586 
587  // Create "stacked" copies of flavored histograms
588  // Cos Stack = Cos
589  TH1* hCosStack = (TH1*)hCos->Clone();
590 
591  // Numu Stack = Cos Stack + Numu = Cos + Numu
592  TH1* hNumuStack = (TH1*)hCosStack->Clone();
593  hNumuStack->Add(hNumu);
594 
595  // Nue Stack = Numu Stack + Nue = Cos + Numu + Nue
596  TH1* hNueStack = (TH1*)hNumuStack->Clone();
597  hNueStack->Add(hNue);
598 
599  // NC Stack = Nue Stack + NC = Cos + Numu + Nue + NC
600  TH1* hNCStack = (TH1*)hNC->Clone();
601  //hNCStack->Add(hNC);
602 
603  // Stylize the stacked plots
604  CenterTitles(hNueStack);
605  CenterTitles(hNumuStack);
606  CenterTitles(hNCStack);
607  CenterTitles(hCosStack);
608  if(hStrlNC) { CenterTitles(hStrlNC); }
609 
610  hNCStack ->SetLineColor(kNCBackgroundColor);
611  hNumuStack->SetLineColor(kNumuBackgroundColor);
612  hNumuStack->SetFillColor(kNumuBackgroundColor);
613  hNueStack ->SetLineColor(kNueSignalColor);
614  hNueStack ->SetFillColor(kNueSignalColor);
615  hCosStack ->SetLineColor(kCosmicBackgroundColor);
616  hCosStack ->SetFillColor(kCosmicBackgroundColor);
617 
618  hNCStack ->SetMinimum(0);
619  hNumuStack->SetMinimum(0);
620  hNueStack ->SetMinimum(0);
621  hCosStack ->SetMinimum(0);
622 
623  // Manually set errors, and calculate the max value (including errors)
624  double maxval = 0.;
625  double maxvalstack = 0.;
626  for(int i = 1, n = hNC->GetNbinsX(); i <= n; ++i) {
627  double error = sqrt(hNC->GetBinContent(i));
628  hNC ->SetBinError(i, error);
629  hNCStack->SetBinError(i, error);
630  maxval = std::max(maxval, std::max(hAll->GetBinContent(i), hNC->GetBinContent(i) + error));
631  maxvalstack = std::max(maxvalstack, hNCStack->GetBinContent(i) + error);
632  }
633 
634  // Scale down the ND by 10e3 (don't forget the max values!)
635  if(det.compare("ND") == 0) {
636  double ndScale = 1./1000.;
637  hAll ->Scale(ndScale);
638  hNC ->Scale(ndScale);
639  hNue ->Scale(ndScale);
640  hNumu->Scale(ndScale);
641  hNCStack ->Scale(ndScale);
642  hNueStack ->Scale(ndScale);
643  hNumuStack->Scale(ndScale);
644  maxval *= ndScale;
645  maxvalstack *= ndScale;
646  }
647 
648  TLatex *tex = new TLatex();
649  tex->SetNDC();
650  tex->SetTextFont(42);
651  tex->SetTextSize(0.037);
652  tex->SetTextAlign(11);
653 
654  // Create a canvas for and draw the unstacked histograms
655  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
656  hNue->SetMaximum(maxval*1.1);
657  hNue->GetXaxis()->SetNdivisions(-406);
658 
659  // Create a legend for the unstacked histograms
660  double xL = 0.5, xR = 0.7, yB = 0.59, yT = 0.84;
661  if(det.compare("ND") == 0) { yB += 0.05; }
662  TLegend* leg = new TLegend(xL, yB, xR, yT);
663  leg->SetBorderSize(0);
664  leg->SetFillColor(kWhite);
665  leg->SetFillStyle(0);
666  leg->SetFillStyle(0);
667  leg->SetTextFont(42);
668  leg->SetTextSize(0.037);
669  leg->AddEntry(hAll, "Total Prediction", "l");
670  leg->AddEntry(hNC, "NC 3 Flavor Prediction", "le");
671  leg->AddEntry(hNue, "#nu_{e} CC Background", "l");
672  leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
673  if(det.compare("FD") == 0) {
674  leg->AddEntry(hCos, "Cosmic Background", "l");
675  }
676 
677  hNue ->Draw("hist");
678  hNumu->Draw("hist same");
679  if(det.compare("FD") == 0) { hCos->Draw("hist same"); }
680  hAll ->Draw("hist same");
681  hNC ->Draw("same");
682  leg->Draw();
683  if(det.compare("FD") == 0) {
684  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
685  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
686  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
687  }
688 
689  gPad->RedrawAxis();
690  Simulation();
691 
692  // Save the canvas
693  rootOut->WriteTObject(c);
694 
695  // Create a canvas for and draw the stacked histograms
696  TCanvas* cStack = new TCanvas((name + "Stack").c_str(), title.c_str(), 800, 500);
697  hNueStack->SetMaximum(maxvalstack*1.1);
698  hNueStack->GetXaxis()->SetNdivisions(-406);
699 
700  if(det.compare("ND") == 0) { yB += 0.05; }
701  TLegend* legStack = new TLegend(xL, yB, xR, yT);
702  legStack->SetBorderSize(0);
703  legStack->SetFillColor(kWhite);
704  legStack->SetFillStyle(0);
705  legStack->SetFillStyle(0);
706  legStack->SetTextFont(42);
707  legStack->SetTextSize(0.037);
708  legStack->AddEntry(hNCStack, "NC 3 Flavor Prediction", "le");
709  if(hStrlNC) { legStack->AddEntry(hStrlNC, "NC 4 Flavor Prediction", "l"); }
710  legStack->AddEntry(hNueStack, "#nu_{e} CC Background", "f");
711  legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
712  if(det.compare("FD") == 0) {
713  legStack->AddEntry(hCosStack, "Cosmic Background", "f");
714  }
715 
716  hNueStack ->Draw("hist");
717  hNumuStack->Draw("hist same");
718  if(det.compare("FD") == 0) { hCosStack->Draw("hist same"); }
719  hNCStack ->Draw("same");
720  if(hStrlNC) { hStrlNC->Draw("hist same"); }
721  legStack->Draw();
722  if(det.compare("FD") == 0) {
723  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
724  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
725  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
726  }
727 
728  gPad->RedrawAxis();
729  Simulation();
730 
731  rootOut->WriteTObject(cStack);
732 
733  return;
734 }
735 
736 //........................................................................
737 void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
739 {
740  Spectrum empty = decomp->NCTotalComponent();
741  empty.Clear();
742 
743  Spectrum spectraND[MAXSPEC] = {
744  decomp->NCTotalComponent() +
745  decomp->NueComponent() + decomp->AntiNueComponent() +
746  decomp->NumuComponent() + decomp->AntiNumuComponent(),
747  decomp->NCTotalComponent(),
748  decomp->NumuComponent() + decomp->AntiNumuComponent(),
749  decomp->NueComponent() + decomp->AntiNueComponent(),
750  empty
751  };
752 
753  std::string det = "ND";
754  std::string fulltitle = det + title;
755  PlotStack(spectraND, rootOut, textOFS, name, fulltitle, det);
756 
757  return;
758 }
759 
760 //........................................................................
762  osc::IOscCalc* calc, osc::IOscCalc* calcSt,
763  TDirectory* rootOut, FILE* textOFS,
765 {
766  Spectrum spectraFD[MAXSPEC] = {
767  pred->Predict(calc) + sCos,
771  sCos
772  };
773 
774  Spectrum sterile = pred->Predict(calcSt) + sCos;
775  Spectrum sterileNC = pred->PredictComponent(calcSt, Flavors::kAll, Current::kNC, Sign::kBoth);
776 
777  std::string det = "FD";
778  std::string fulltitle = det + title;
779  PlotStack(spectraFD, rootOut, textOFS, name, fulltitle, det, &sterile, &sterileNC);
780 
781  return;
782 }
783 
784 //........................................................................
785 void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
786  Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
787  osc::IOscCalc* calc, Spectrum& sCos,
788  TDirectory* rootOut, FILE* textOFS,
790 {
791  Spectrum sFDNum = pFDNum ->Predict(calc);
792  Spectrum sFDPurDen = pFDPurDen->Predict(calc) + sCos;
793  Spectrum sFDEffDen = pFDEffDen->Predict(calc);
794 
795  TH1* hFDNum = sFDNum .ToTH1(NCSCALE);
796  TH1* hFDPurDen = sFDPurDen.ToTH1(NCSCALE);
797  TH1* hFDEffDen = sFDEffDen.ToTH1(NCSCALE);
798 
799  TH1* hNDNum = sNDNum .ToTH1(NCSCALE);
800  TH1* hNDPurDen = sNDPurDen.ToTH1(NCSCALE);
801  TH1* hNDEffDen = sNDEffDen.ToTH1(NCSCALE);
802 
803  TH1* hFDEff = (TH1*)hFDNum->Clone();
804  hFDEff->Divide(hFDEffDen);
805 
806  TH1* hNDEff = (TH1*)hNDNum->Clone();
807  hNDEff->Divide(hNDEffDen);
808 
809  TH1* hFDPur = (TH1*)hFDNum->Clone();
810  hFDPur->Divide(hFDPurDen);
811 
812  TH1* hNDPur = (TH1*)hNDNum->Clone();
813  hNDPur->Divide(hNDPurDen);
814 
815  hFDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
816  hFDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
817  hNDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
818  hNDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
819 
820  hFDEff->SetLineColor(kBlue);
821  hNDEff->SetLineColor(kBlue);
822 
823  hFDPur->SetLineColor(kRed);
824  hNDPur->SetLineColor(kRed);
825 
826  hNDEff->SetLineStyle(2);
827  hNDPur->SetLineStyle(2);
828 
829  hFDEff->SetMaximum(1.);
830  hFDEff->SetMinimum(0.);
831  hFDEff->GetXaxis()->SetNdivisions(-406);
832 
833  hNDEff->SetMaximum(1.);
834  hNDEff->SetMinimum(0.);
835 
836  hFDPur->SetMaximum(1.);
837  hFDPur->SetMinimum(0.);
838 
839  hNDPur->SetMaximum(1.);
840  hNDPur->SetMinimum(0.);
841 
842  CenterTitles(hFDEff);
843  CenterTitles(hNDEff);
844  CenterTitles(hFDPur);
845  CenterTitles(hNDPur);
846 
847  TLegend* leg = new TLegend(0.6, 0.35, 0.85, 0.55);
848  leg->AddEntry(hFDEff, "FD MC Efficiency", "L");
849  leg->AddEntry(hNDEff, "ND MC Efficiency", "L");
850  leg->AddEntry(hFDPur, "FD MC Purity", "L");
851  leg->AddEntry(hNDPur, "ND MC Purity", "L");
852 
853  TCanvas* c = new TCanvas(name.c_str(), "Purity and Efficiency", 800, 500);
854  hFDEff->Draw("hist");
855  hNDEff->Draw("hist same");
856  hFDPur->Draw("hist same");
857  hNDPur->Draw("hist same");
858  leg->Draw();
859 
860  gPad->RedrawAxis();
861  Simulation();
862 
863  rootOut->WriteTObject(c);
864 }
865 
866 
867 //........................................................................
869 {
870  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
871  calc->SetAngle(3, 4, 0.);
872  calc->SetAngle(2, 4, 0.);
873 
874  return;
875 }
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
enum BeamMode kOrange
void Plot1DSlice(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var, std::vector< const IFitVar * > profVars, TDirectory *rootOut, int nbins, double min, double max, std::string name, Color_t color)
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
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
cons_index_list< index_multi, nil_index_list > multi
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
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:148
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
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:361
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
osc::OscCalcDumb calc
const Color_t kNumuBackgroundColor
Definition: Style.h:30
#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:40
const Color_t kNueSignalColor
Definition: Style.h:19
fclose(fg1)
Log-likelihood scan across two parameters.
const int nbins
Definition: cellShifts.C:15
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
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)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
const double a
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
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:69
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 Plot1DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, Color_t color)
double POT() const
Definition: Spectrum.h:227
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
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
Splits Data proportionally according to MC.
void SetDm(int i, double dm)
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
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
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Sum MC predictions from different periods scaled according to data POT targets.
static double NCSCALE
void Plot2DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a)
const Color_t kTotalMCColor
Definition: Style.h:16
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
const Color_t kNCBackgroundColor
Definition: Style.h:22
enum BeamMode kBlue
void BlessedPlotsAnaByPeriod()
Compare a data spectrum to MC expectation using only the event count.
virtual Spectrum NueComponent() const =0
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string