BlessedPlotsAna.C

Go to the documentation of this file.

#include "TCanvas.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TPad.h"

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Core/IFitVar.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "NuXAna/Vars/FitVarsSterile.h"

#include "OscLib/OscCalcSterile.h"

#include "NuXAna/macros/NuSPlotFunctions.h"

using namespace ana;

// Definitions for organization

// One struct to contain all values for angle fits
struct AngleValues{
  int nbins23;
  double min23;
  double max23;

  int nbins34;
  double min34;
  double max34;

  int nbins24;
  double min24;
  double max24;
};

// Plot a single 1D slice canvas -- this function does the actual plotting
// The macro does not call this function directly in the main function
// \param expt The CAFAna experiment object, contains FD prediction and "data"
// \param calc Sterile calculator for oscillating the spectra
// \param var Parameter to fit for!
// \param profVars Vector of parameters to profile over
// \param rootOut Already open file for root output
// \param name Base of the string for creating the name for saving the root objects
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);

// Plot 1D slice canvases
// This function sets up multiple calls to the single canvas version of Plot1DSlice
// \param expt The CAFAna experiment object, contains FD prediction and "data"
// \param calc Sterile calculator for oscillating the spectra
// \param rootOut Already open file for root output
// \param a Struct containing bin values for slice histograms
// \param years A string for the number of data years
void Plot1DSlices(IExperiment* expt, osc::OscCalcSterile* calc,
                  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
                  TDirectory* rootOut, AngleValues a,
                  std::string years, Color_t color);

// Plot 2D slice surfaces
// This function sets up multiple calls to the workhorse version of Plot2DSlices
// \param expt The CAFAna experiment object, contains FD prediction and "data"
// \param calc Sterile calculator for oscillating the spectra
// \param rootOut Already open file for root output
// \param a Struct containing bin values for slice histograms
// \param years A string for the number of data years
void Plot2DSlices(IExperiment* expt, osc::OscCalcSterile* calc,
                  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
                  TDirectory* rootOut, AngleValues a, std::string years);

// Plot 2D slice surfaces
// The macro does not call this function directly
// \param rootOut Already open file for root output
// \param prof String indicating whether surfaces are profiled or not (optional)
void Plot2DSlices(FrequentistSurface* s3424, FrequentistSurface* s2324, FrequentistSurface* s3423,
                  TDirectory* rootOut, std::string years, std::string prof = "");

// Plot spectra and stack canvases -- the workhorse version
// The macro does not call this function directly
// \param spectra Array of spectra, should match the output of the SPECARR macro
// \param rootOut Already open file for root output
// \param textOFS Already open filestream for text output
// \param name A name for the output, only part of the saved output object's name
// \param title Title for the histograms. Typically this is the x axis label
// \param det The detector the spectra correspond to, should be 2 chars long for nice formatting
// \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
// \param sterile Pointer to a spectrum with 4 flavor osc prediction, includes all event types (optional)
// \param sterileNC Pointer to a spectrum with 4 flavor osc prediction, only includes NC events (optional)
void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
               std::string name, std::string title, std::string det, int POT,
               Spectrum* sterile = nullptr, Spectrum* sterileNC = nullptr);

// Plot spectra and stack canvases for ND
// Sets up call to the workhorse version of PlotStack
// \param decomp The ND decomposition
// \param rootOut Already open file for root output
// \param textOFS Already open filestream for text output
// \param name A name for the output, only part of the saved output object's name
// \param title Title for the histograms. Typically this is the x axis label
// \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
               std::string name, std::string title, int POT);

// Plot spectra and stack canvases for FD
// Sets up call to the workhorse version of PlotStack
// \param pred The FD prediction object
// \param sCos The FD cosmic prediction spectrum
// \param calc Sterile calculator assuming 3 flavor oscillations
// \param calcSt Sterile calculator assuming 4 flavor oscillations
// \param rootOut Already open file for root output
// \param textOFS Already open filestream for text output
// \param name A name for the output, only part of the saved output object's name
// \param title Title for the histograms. Typically this is the x axis label
// \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
void PlotStack(IPrediction* pred,
               Spectrum& sCos,
               osc::IOscCalc* calc,
               osc::IOscCalc* calcSt,
               TDirectory* rootOut, FILE* textOFS,
               std::string name, std::string title, int POT);

// Function to plot purity and efficiency at both FD and ND
// \param pFDNum Prediction object for FD, used as nominator for FD purity and efficiency
// \param pFDPurDen Prediction object for FD, used as denominator for FD purity
// \param pFDPurDen Prediction object for FD, used as denominator for FD efficiency
// \param sNDNum Spectrum object for ND, used as nominator for ND purity and efficiency
// \param sNDPurDen Spectrum object for ND, used as denominator for ND purity
// \param sNDPurDen Spectrum object for ND, used as denominator for ND efficiency
// \param calc Sterile calculator assuming 3 flavor oscillations
// \param rootOut Already open file for root output
// \param textOFS Already open filestream for text output
// \param name A name for the output, only part of the saved output object's name
// \param title Title for the histograms. Typically this is the x axis label
void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
                Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
                osc::IOscCalc* calc, Spectrum& sCos,
                TDirectory* rootOut, FILE* textOFS,
                std::string name, std::string title);

// After running a fit, reset the values in a calculator to preset defaults
void ResetAngles(osc::OscCalcSterile* calc);

void BlessedPlotsAna()
{
  TH1::AddDirectory(0);

  // Set up path to file with filled CAFAna objects and for output of analysis
  //std::string folder = "/nova/ana/steriles/Ana01/";
  std::string folder = "";
  std::string filenm = "BlessedPlots14db";

  std::string loadLocation = folder + filenm + ".root";
  std::string saveLocation = folder + filenm + "Ana.root";
  std::string textLocation = folder + filenm + "Ana.txt";

  // Load filled objects from file
  TFile* rootL = new TFile(loadLocation.c_str(), "UPDATE");
  TDirectory* tmpL = gDirectory;
  TDirectory* loadDir = gDirectory;

  loadDir->cd((loadLocation + ":/decompNC").c_str());
  ProportionalDecomp decompNC   = *ProportionalDecomp::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/prediction").c_str());
  PredictionSterile pred        = *PredictionSterile ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sData1Yr").c_str());
  Spectrum sFDData1y            = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sData3Yr").c_str());
  Spectrum sFDData3y            = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sCos1Yr").c_str());
  Spectrum sCosmic1y            = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sCos3Yr").c_str());
  Spectrum sCosmic3y            = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/pFDPENum").c_str());
  PredictionExtrap pFDPurEffNum = *PredictionExtrap  ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/pFDPDen").c_str());
  PredictionExtrap pFDPurityDen = *PredictionExtrap  ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/pFDEDen").c_str());
  PredictionExtrap pFDEfcncyDen = *PredictionExtrap  ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sNDPENum").c_str());
  Spectrum sNDPurEffNum         = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sNDPDen").c_str());
  Spectrum sNDPurityDen         = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sNDEDen").c_str());
  Spectrum sNDEfcncyDen         = *Spectrum          ::LoadFrom(gDirectory);
  tmpL->cd();
  rootL->Close(); // Don't forget to close the file!

  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc3f = DefaultSterileCalc(3);

  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  calc4f->SetDm(4, 0.5); //Dm41 = 0.5 eV^2
  calc4f->SetAngle(2, 4, 0.);
  calc4f->SetAngle(3, 4, 0.);

  std::string labelEReco = "Calorimetric Energy (GeV)";

  // Open files for saving analysis output
  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
  FILE* textF;
  textF = fopen(textLocation.c_str(), "w");

  // Choose some angle values for making 4 flavor predictions
  calc4f->SetAngle(2, 4, 0.17453293); // 10 degrees
  calc4f->SetAngle(3, 4, 0.61086524); // 35 degrees

  // Plot spectra for 1 data yr
  std::string titleHelper = " Spectra for 6e20 POT";
  PlotStack(&pred, sCosmic1y, calc3f, calc4f, rootF, textF, "FDSpectra1yr", titleHelper, 6);
  PlotStack(&decompNC, rootF, textF, "NDSpectra1yr", titleHelper, 6);
  
  // Plot spectra for 3 data yr
  titleHelper = " Spectra for 18e20 POT";
  PlotStack(&pred, sCosmic3y, calc3f, calc4f, rootF, textF, "FDSpectra3yr", titleHelper, 18);
  PlotStack(&decompNC, rootF, textF, "NDSpectra3yr", titleHelper, 18);

  // Plot purity and efficiency
  PlotPurEff(&pFDPurEffNum, &pFDPurityDen, &pFDEfcncyDen,
             sNDPurEffNum, sNDPurityDen, sNDEfcncyDen,
             calc3f, sCosmic1y, rootF, textF, "PurEff", labelEReco);
  
  fclose(textF); // Text file is done now, close it
  ResetAngles(calc4f);
  
  // Create Experiment objects, one each for the 1 data yr and 3 data yr predictions
  GaussianConstraint th23Constraint1y(&kFitTheta23InDegreesSterile, 45., 8.73);
  CountingExperiment expt1y(&pred, sFDData1y, sCosmic1y);
  MultiExperiment multi1y({&th23Constraint1y, &expt1y});

  GaussianConstraint th23Constraint3y(&kFitTheta23InDegreesSterile, 45., 8.73);
  CountingExperiment expt3y(&pred, sFDData3y, sCosmic3y);
  MultiExperiment multi3y({&th23Constraint3y, &expt3y});

  // Vector of parameters to fit: theta 23, 34, 24
  std::vector<const IFitVar*> fitvars;
  fitvars.push_back(&kFitTheta34InDegreesSterile);
  fitvars.push_back(&kFitTheta24InDegreesSterile);
  fitvars.push_back(&kFitTheta23InDegreesSterile);

  // Create fit objects
  MinuitFitter fit1y(&multi1y, fitvars);
  MinuitFitter fit3y(&multi3y, fitvars);

  // Run fits!
  fit1y.Fit(calc4f);
  ResetAngles(calc4f);
  fit3y.Fit(calc4f);
  ResetAngles(calc4f);

  // Set up values for slice and surface plots
  AngleValues avals;
  avals.nbins23 = 200;
  avals.min23 = 20.;
  avals.max23 = 70.;
  avals.nbins34 = 200;
  avals.min34 = 0.;
  avals.max34 = 50.;
  avals.nbins24 = 180;
  avals.min24 = 0.;
  avals.max24 = 45.;

  const Color_t kFitColor = kRed;

  // Plot 1D fits for 1 data yr experiment
  Plot1DSlices(
    &multi1y, calc4f,
    &kFitTheta23InDegreesSterile,
    &kFitTheta34InDegreesSterile,
    &kFitTheta24InDegreesSterile,
    rootF, avals, "1", kFitColor
  );

  // Plot 1D fits for 3 data yr experiment
  Plot1DSlices(
    &multi3y, calc4f,
    &kFitTheta23InDegreesSterile,
    &kFitTheta34InDegreesSterile,
    &kFitTheta24InDegreesSterile,
    rootF, avals, "3", kFitColor
  );

  // Lower the number of bins for surfaces
  avals.nbins23 = 50;
  avals.nbins34 = 50;
  avals.nbins24 = 45;

  // Plot the surfaces for 1 data yr experiment
  Plot2DSlices(
    &multi1y, calc4f,
    &kFitTheta23InDegreesSterile,
    &kFitTheta34InDegreesSterile,
    &kFitTheta24InDegreesSterile,
    rootF, avals, "1"
  );

  // Plot the surfaces for 3 data yr experiment
  Plot2DSlices(
    &multi3y, calc4f,
    &kFitTheta23InDegreesSterile,
    &kFitTheta34InDegreesSterile,
    &kFitTheta24InDegreesSterile,
    rootF, avals, "3"
  );
  
  rootF->Close(); // Close the output root file
}

//........................................................................
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)
{
  // Create slice histogram with/without profiling based on input profVars vector size
  
  TH1* h = (profVars.size() > 0 ?
            Profile(expt, calc, var, nbins, min, max, -1, profVars) :
            Slice(expt, calc, var, nbins, min, max));

  ResetAngles(calc); // Don't forget to reset calculator...

  // Formatting:
  h->SetLineColor(color); // Set the line color
  h->SetMaximum(5); // Want the plot to go from 0 to 5
  h->SetMinimum(0);
  h->SetName(name.c_str()); // Set the object name for saving
  rootOut->WriteTObject(h); // Save the raw histogram

  // Create an object for writing LaTeX to canvas
  TLatex *tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);

  // Coordinates for LaTeX object
  double xTex = 0.92, y68 = 0.25, y90 = 0.52;

  // Dotted lines representing 1 and 2 sigma
  TLine* line68 = new TLine(min,  1,    max,  1);
  TLine* line90 = new TLine(min,  2.71, max,  2.71);

  // Sometimes histogram has horizontal artifact at max canvas value
  // This line hides it
  TLine* line5  = new TLine(min,  5,    max,  5);

  // Dashed line for 1 and 2 sigma lines
  line68->SetLineStyle(2);
  line90->SetLineStyle(2);

  // Formatting
  line68->SetLineWidth(2);
  line90->SetLineWidth(2);
  line5 ->SetLineWidth(2);

  // Input "name" should have "h" for histogram at the beginning
  // Make canvas name with input "name" but remove that "h"
  std::string cname = "c1DSlice" + name.substr(1, name.size()-1);

  // Theta subscripts should be at position 6, 7, get them
  std::string indices = name.substr(6, 2);

  // Create POT string based on whether "name" has 1 or 3 in it
  std::string cpotstr = ( name.find("1") != std::string::npos ?
                          "6e20 POT" : "18e20 POT" );

  // Create a title for canvas
  std::string ctitle = "Delta chi^2 for Theta" + indices;
  // Add info on profiling if applicable
  ctitle += ( name.find("Prof") != std::string::npos ?
              " (Profiling Over Other Angles) for " + cpotstr :
              " for " + cpotstr );

  // Create the canvas and draw!
  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
  h->Draw("]["); // Don't show vertical drop to 0 at edge of histogram
  line68->DrawClone("same");
  line90->DrawClone("same");
  line5 ->DrawClone("same");
  tex->DrawLatex(xTex, y68, "68%"); // Add label for 1 sigma line
  tex->DrawLatex(xTex, y90, "90%"); // Add label for 2 sigma line
  Simulation();
  rootOut->WriteTObject(c); // Save full canvas

  return;
}

//........................................................................
void Plot1DSlices(IExperiment* expt, osc::OscCalcSterile* calc,
                  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
                  TDirectory* rootOut, AngleValues a,
                  std::string years, Color_t color)
{
  // General form of the object name for saving
  std::string name = "h" + years + "yrTh";
  std::string fullname; // fullname adds 2 characters, the theta subscript indices

  // Plot unprofiled slices, for theta 23, then 34, then 24
  fullname = name + "23";
  Plot1DSlice(expt, calc, var23, {},
              rootOut, a.nbins23, a.min23, a.max23, fullname, color);
  fullname = name + "34";
  Plot1DSlice(expt, calc, var34, {},
              rootOut, a.nbins34, a.min34, a.max34, fullname, color);
  fullname = name + "24";
  Plot1DSlice(expt, calc, var24, {},
              rootOut, a.nbins24, a.min24, a.max24, fullname, color);

  
  // Plot profiled slices in same order as above
  // When slice is of one angle, the other two are profiled
  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
  fullname = name + "23Prof";
  Plot1DSlice(expt, calc, var23, {var34, var24},
              rootOut, a.nbins23, a.min23, a.max23, fullname, color);
  fullname = name + "34Prof";
  Plot1DSlice(expt, calc, var34, {var23, var24},
              rootOut, a.nbins34, a.min34, a.max34, fullname, color);
  fullname = name + "24Prof";
  Plot1DSlice(expt, calc, var24, {var23, var34},
              rootOut, a.nbins24, a.min24, a.max24, fullname, color);
  
  return;
}

//........................................................................
void Plot2DSlices(IExperiment* expt, osc::OscCalcSterile* calc,
                  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
                  TDirectory* rootOut, AngleValues a, std::string years)
{
  // Create unprofiled surfaces
  // 34 vs 24
  FrequentistSurface s3424(expt, calc,
                var34, a.nbins34, a.min34, a.max34,
                var24, a.nbins24, a.min24, a.max24);
  ResetAngles(calc);

  // 23 vs 24
  FrequentistSurface s2324(expt, calc,
                var23, a.nbins23, a.min23, a.max23,
                var24, a.nbins24, a.min24, a.max24);
  ResetAngles(calc);

  // 34 vs 23
  FrequentistSurface s3423(expt, calc,
                var34, a.nbins34, a.min34, a.max34,
                var23, a.nbins23, a.min23, a.max23);
  ResetAngles(calc);

  // Create profiled surfaces in same order as above
  // The third angle is the one that is profiled
  // I.e., the 34 vs 24 surface profiles 23
  FrequentistSurface s3424Prof(expt, calc,
                    var34, a.nbins34, a.min34, a.max34,
                    var24, a.nbins24, a.min24, a.max24,
                    {var23});
  ResetAngles(calc);

  FrequentistSurface s2324Prof(expt, calc,
                    var23, a.nbins23, a.min23, a.max23,
                    var24, a.nbins24, a.min24, a.max24,
                    {var34});
  ResetAngles(calc);

  FrequentistSurface s3423Prof(expt, calc,
                    var34, a.nbins34, a.min34, a.max34,
                    var23, a.nbins23, a.min23, a.max23,
                    {var24});
  ResetAngles(calc);

  // Plot all of the surfaces
  Plot2DSlices(&s3424,     &s2324,     &s3423,     rootOut, years);
  Plot2DSlices(&s3424Prof, &s2324Prof, &s3423Prof, rootOut, years, "Prof");

  return;
}

//........................................................................
void Plot2DSlices(FrequentistSurface* s3424, FrequentistSurface* s2324, FrequentistSurface* s3423,
                  TDirectory* rootOut, std::string years, std::string prof)
{
  // General format for the name for saving
  std::string name = "h" + years + "yrTh";
  std::string fullname;

  // Add indices for first theta, add second theta, and profile string if applicaable
  fullname = name + "34Th24" + prof;
  TH1* h3424 = s3424->ToTH2(); // Create a histogram from the surface
  h3424->SetName(fullname.c_str()); // Set the name for saving
  rootOut->WriteTObject(h3424); // Save the raw histogram

  fullname = name + "23Th24" + prof;
  TH1* h2324 = s2324->ToTH2();
  h2324->SetName(fullname.c_str());
  rootOut->WriteTObject(h2324);

  fullname = name + "34Th23" + prof;
  TH1* h3423 = s3423->ToTH2();
  h3423->SetName(fullname.c_str());
  rootOut->WriteTObject(h3423);

  // Create name for canvas
  std::string cname = "c2DSlice" + years + "yr" + prof;
  // Create string for POT based on whether name has 1 or 3 in it
  std::string potstr = ( name.find("1") != std::string::npos ?
                         "6e20 POT" : "18e20 POT" );
  // Create relevant title for the canvas
  std::string ctitle = (
    prof.size() > 0 ?
    "Surfaces (Profiling Over Other Angles) for " + potstr :
    "Surfaces for " + potstr
  );

  // Create triptych canvas and draw!
  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 800, 800);
  c->Divide(2, 2); // Split into 2 x 2 grid
  c->cd(1); // Go to top left
  s3424->Draw(); // Draw histogram
  s3424->DrawContour(Gaussian68Percent2D(*s3424), 7, kBlack); // Draw 1 sigma contour
  s3424->DrawContour(Gaussian90Percent2D(*s3424), 1, kBlack); // Draw 2 sigma contour
  s3424->DrawBestFit(kBlack); // Draw best fit point
  c->cd(2); // Go to top right
  s2324->Draw();
  s2324->DrawContour(Gaussian68Percent2D(*s2324), 7, kBlack);
  s2324->DrawContour(Gaussian90Percent2D(*s2324), 1, kBlack);
  s2324->DrawBestFit(kBlack);
  Simulation();
  c->cd(3); // Go to bottom left
  s3423->Draw();
  s3423->DrawContour(Gaussian68Percent2D(*s3423), 7, kBlack);
  s3423->DrawContour(Gaussian90Percent2D(*s3423), 1, kBlack);
  s3423->DrawBestFit(kBlack);
  rootOut->WriteTObject(c); // Save the root canvas

  return;
}

//........................................................................
void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
               std::string name, std::string title, std::string det, int POT,
               Spectrum* sterile, Spectrum* sterileNC)
{
  // The input is in units of 1e20, turn this into the real scale
  double realPOT = (double)POT*1e20;

  char potBuff[2];
  sprintf(potBuff, "%d", POT);
  std::string potStr = std::string(potBuff) + " x 10^{20} POT";

  // Get each histogram from the Spectrum array
  TH1* hAll  = spectra[0].ToTH1(realPOT);
  TH1* hNC   = spectra[1].ToTH1(realPOT);
  TH1* hNumu = spectra[2].ToTH1(realPOT);
  TH1* hNue  = spectra[3].ToTH1(realPOT);
  TH1* hCos  = spectra[4].ToTH1(realPOT);
  // The default for these inputs are nullptr, only get a TH1 if they exist
  TH1* hStrl   = (sterile   ? sterile  ->ToTH1(realPOT) : nullptr);
  TH1* hStrlNC = (sterileNC ? sterileNC->ToTH1(realPOT) : nullptr);

  // Create the title and axis labels, which is an empty title,
  // the x label from an input Spectrum, and Events/POT
  std::string xLabel(hAll->GetXaxis()->GetTitle());
  std::string yLabel = (det.compare("ND") == 0 ? "10^{3} Events" : "Events");
  yLabel += "/" + potStr + " / 0.25 GeV";
  std::string histTitle = ";" + xLabel + ";" + yLabel;

  // Print the column title and headers
  fprintf(textOFS, "Event numbers at the %2.2s for %2de20 POT:\n", det.c_str(), POT);
  if(hStrlNC) { fprintf(textOFS, "%12.12s, ", "NC (Sterile)"); }
  else      { fprintf(textOFS, "%14.14s", ""); }
  fprintf(textOFS, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
          "NC (3 Flav)", "All", "Numu", "Nue", "Cosmic", "FOM");

  // Get the integral number of events, then print them
  double nNC3F = hNC  ->Integral();
  double nAll  = hAll ->Integral();
  double nNumu = hNumu->Integral();
  double nNue  = hNue ->Integral();
  double nCos  = hCos ->Integral();
  double fom   = nNC3F/sqrt(nAll);
  if(hStrlNC) {
    double nNCSt = hStrlNC->Integral();
    fprintf(textOFS, "%10E, ", nNCSt);
  }
  else {
    fprintf(textOFS, "%14.14s", "");
  }
  fprintf(textOFS, "%10E, %10E, %10E, %10E, %10E, %10E\n\n",
          nNC3F, nAll, nNumu, nNue, nCos, fom);

  // Stylize the unstacked plots
  CenterTitles(hNue);
  CenterTitles(hNumu);
  CenterTitles(hNC);
  CenterTitles(hCos);
  CenterTitles(hAll);

  hAll ->SetMinimum(0);
  hNC  ->SetMinimum(0);
  hNumu->SetMinimum(0);
  hNue ->SetMinimum(0);
  hCos ->SetMinimum(0);
  if(hStrlNC) { hStrl->SetMinimum(0); }

  const Color_t kCosmicBackgroundColor = kOrange+1;
  hAll ->SetLineColor(kTotalMCColor);
  hNC  ->SetLineColor(kNCBackgroundColor);
  hNumu->SetLineColor(kNumuBackgroundColor);
  hNue ->SetLineColor(kNueSignalColor);
  hCos ->SetLineColor(kCosmicBackgroundColor);
  if(hStrlNC) {
    hStrlNC->SetLineColor(kNCBackgroundColor);
    hStrlNC->SetLineStyle(2);
  }

  hAll ->SetTitle(histTitle.c_str());
  hNC  ->SetTitle(histTitle.c_str());
  hNumu->SetTitle(histTitle.c_str());
  hNue ->SetTitle(histTitle.c_str());
  hCos ->SetTitle(histTitle.c_str());
  if(hStrlNC) { hStrlNC->SetTitle(histTitle.c_str()); }

  // Create "stacked" copies of flavored histograms
  // Cos Stack = Cos
  TH1* hCosStack  = (TH1*)hCos->Clone();

  // Numu Stack = Cos Stack + Numu = Cos + Numu
  TH1* hNumuStack = (TH1*)hCosStack->Clone();
  hNumuStack->Add(hNumu);

  // Nue Stack = Numu Stack + Nue = Cos + Numu + Nue
  TH1* hNueStack  = (TH1*)hNumuStack->Clone();
  hNueStack->Add(hNue);

  // NC Stack = Nue Stack + NC = Cos + Numu + Nue + NC
  TH1* hNCStack   = (TH1*)hNC->Clone();
  //hNCStack->Add(hNC);

  // Stylize the stacked plots
  CenterTitles(hNueStack);
  CenterTitles(hNumuStack);
  CenterTitles(hNCStack);
  CenterTitles(hCosStack);
  if(hStrlNC) { CenterTitles(hStrlNC); }

  hNCStack  ->SetLineColor(kNCBackgroundColor);
  hNumuStack->SetLineColor(kNumuBackgroundColor);
  hNumuStack->SetFillColor(kNumuBackgroundColor);
  hNueStack ->SetLineColor(kNueSignalColor);
  hNueStack ->SetFillColor(kNueSignalColor);
  hCosStack ->SetLineColor(kCosmicBackgroundColor);
  hCosStack ->SetFillColor(kCosmicBackgroundColor);

  hNCStack  ->SetMinimum(0);
  hNumuStack->SetMinimum(0);
  hNueStack ->SetMinimum(0);
  hCosStack ->SetMinimum(0);

  // Manually set errors, and calculate the max value (including errors)
  double maxval = 0.;
  double maxvalstack = 0.;
  for(int i = 1, n = hNC->GetNbinsX(); i <= n; ++i) {
    double error = sqrt(hNC->GetBinContent(i));
    hNC     ->SetBinError(i, error);
    hNCStack->SetBinError(i, error);
    maxval = std::max(maxval, std::max(hAll->GetBinContent(i), hNC->GetBinContent(i) + error));
    maxvalstack = std::max(maxvalstack, hNCStack->GetBinContent(i) + error);
  }

  // Scale down the ND by 10e3 (don't forget the max values!)
  if(det.compare("ND") == 0) {
    double ndScale = 1./1000.;
    hAll ->Scale(ndScale);
    hNC  ->Scale(ndScale);
    hNue ->Scale(ndScale);
    hNumu->Scale(ndScale);
    hNCStack  ->Scale(ndScale);
    hNueStack ->Scale(ndScale);
    hNumuStack->Scale(ndScale);
    maxval *= ndScale;
    maxvalstack *= ndScale;
  }

  TLatex *tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);

  // Create a canvas for and draw the unstacked histograms
  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
  hNue->SetMaximum(maxval*1.1);
  hNue->GetXaxis()->SetNdivisions(-406);

  // Create a legend for the unstacked histograms
  double xL = 0.5, xR = 0.7, yB = 0.59, yT = 0.84;
  if(det.compare("ND") == 0) { yB += 0.05; }
  TLegend* leg = new TLegend(xL, yB, xR, yT);
  leg->SetBorderSize(0);
  leg->SetFillColor(kWhite);
  leg->SetFillStyle(0);
  leg->SetFillStyle(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.037);
  leg->AddEntry(hAll,  "Total Prediction", "l");
  leg->AddEntry(hNC,   "NC 3 Flavor Prediction", "le");
  leg->AddEntry(hNue,  "#nu_{e} CC Background", "l");
  leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
  if(det.compare("FD") == 0) {
    leg->AddEntry(hCos,  "Cosmic Background", "l");
  }

  hNue ->Draw("hist");
  hNumu->Draw("hist same");
  if(det.compare("FD") == 0) { hCos->Draw("hist same"); }
  hAll ->Draw("hist same");
  hNC  ->Draw("same");
  leg->Draw();
  if(det.compare("FD") == 0) {
    tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
    tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
    tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
  }

  gPad->RedrawAxis();
  Simulation();

  // Save the canvas
  rootOut->WriteTObject(c);

  // Create a canvas for and draw the stacked histograms
  TCanvas* cStack = new TCanvas((name + "Stack").c_str(), title.c_str(), 800, 500);
  hNueStack->SetMaximum(maxvalstack*1.1);
  hNueStack->GetXaxis()->SetNdivisions(-406);

  if(det.compare("ND") == 0) { yB += 0.05; }
  TLegend* legStack = new TLegend(xL, yB, xR, yT);
  legStack->SetBorderSize(0);
  legStack->SetFillColor(kWhite);
  legStack->SetFillStyle(0);
  legStack->SetFillStyle(0);
  legStack->SetTextFont(42);
  legStack->SetTextSize(0.037);
  legStack->AddEntry(hNCStack,   "NC 3 Flavor Prediction", "le");
  if(hStrlNC) { legStack->AddEntry(hStrlNC, "NC 4 Flavor Prediction", "l"); }
  legStack->AddEntry(hNueStack,  "#nu_{e} CC Background", "f");
  legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
  if(det.compare("FD") == 0) {
    legStack->AddEntry(hCosStack,  "Cosmic Background", "f");
  }

  hNueStack ->Draw("hist");
  hNumuStack->Draw("hist same");
  if(det.compare("FD") == 0) { hCosStack->Draw("hist same"); }
  hNCStack  ->Draw("same");
  if(hStrlNC) { hStrlNC->Draw("hist same"); }
  legStack->Draw();
  if(det.compare("FD") == 0) {
    tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
    tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
    tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
  }

  gPad->RedrawAxis();
  Simulation();

  rootOut->WriteTObject(cStack);

  return;
}

//........................................................................
void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
               std::string name, std::string title, int POT)
{
  Spectrum empty = decomp->NCTotalComponent();
  empty.Clear();

  Spectrum spectraND[MAXSPEC] = {
    decomp->NCTotalComponent()   +
    decomp->NueComponent()  + decomp->AntiNueComponent() +
    decomp->NumuComponent() + decomp->AntiNumuComponent(),
    decomp->NCTotalComponent(),
    decomp->NumuComponent() + decomp->AntiNumuComponent(),
    decomp->NueComponent()  + decomp->AntiNueComponent(),
    empty
  };

  std::string det = "ND";
  std::string fulltitle = det + title;
  PlotStack(spectraND, rootOut, textOFS, name, fulltitle, det, POT);

  return;
}

//........................................................................
void PlotStack(IPrediction* pred, Spectrum& sCos,
               osc::IOscCalc* calc, osc::IOscCalc* calcSt,
               TDirectory* rootOut, FILE* textOFS,
               std::string name, std::string title, int POT)
{
  Spectrum spectraFD[MAXSPEC] = {
    pred->Predict(calc) + sCos,
    pred->PredictComponent(calc, Flavors::kAll,     Current::kNC, Sign::kBoth),
    pred->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
    pred->PredictComponent(calc, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
    sCos
  };

  Spectrum sterile   = pred->Predict(calcSt) + sCos;
  Spectrum sterileNC = pred->PredictComponent(calcSt, Flavors::kAll, Current::kNC, Sign::kBoth);

  std::string det = "FD";
  std::string fulltitle = det + title;
  PlotStack(spectraFD, rootOut, textOFS, name, fulltitle, det, POT, &sterile, &sterileNC);

  return;
}

//........................................................................
void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
                Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
                osc::IOscCalc* calc, Spectrum& sCos,
                TDirectory* rootOut, FILE* textOFS,
                std::string name, std::string title)
{
  Spectrum sFDNum    = pFDNum   ->Predict(calc);
  Spectrum sFDPurDen = pFDPurDen->Predict(calc) + sCos;
  Spectrum sFDEffDen = pFDEffDen->Predict(calc);

  TH1* hFDNum    = sFDNum   .ToTH1(NCSCALE);
  TH1* hFDPurDen = sFDPurDen.ToTH1(NCSCALE);
  TH1* hFDEffDen = sFDEffDen.ToTH1(NCSCALE);

  TH1* hNDNum    = sNDNum   .ToTH1(NCSCALE);
  TH1* hNDPurDen = sNDPurDen.ToTH1(NCSCALE);
  TH1* hNDEffDen = sNDEffDen.ToTH1(NCSCALE);

  TH1* hFDEff = (TH1*)hFDNum->Clone();
  hFDEff->Divide(hFDEffDen);

  TH1* hNDEff = (TH1*)hNDNum->Clone();
  hNDEff->Divide(hNDEffDen);

  TH1* hFDPur = (TH1*)hFDNum->Clone();
  hFDPur->Divide(hFDPurDen);

  TH1* hNDPur = (TH1*)hNDNum->Clone();
  hNDPur->Divide(hNDPurDen);

  hFDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
  hFDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
  hNDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
  hNDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());

  hFDEff->SetLineColor(kBlue);
  hNDEff->SetLineColor(kBlue);

  hFDPur->SetLineColor(kRed);
  hNDPur->SetLineColor(kRed);

  hNDEff->SetLineStyle(2);
  hNDPur->SetLineStyle(2);

  hFDEff->SetMaximum(1.);
  hFDEff->SetMinimum(0.);
  hFDEff->GetXaxis()->SetNdivisions(-406);
  
  hNDEff->SetMaximum(1.);
  hNDEff->SetMinimum(0.);

  hFDPur->SetMaximum(1.);
  hFDPur->SetMinimum(0.);

  hNDPur->SetMaximum(1.);
  hNDPur->SetMinimum(0.);

  CenterTitles(hFDEff);
  CenterTitles(hNDEff);
  CenterTitles(hFDPur);
  CenterTitles(hNDPur);

  TLegend* leg = new TLegend(0.6, 0.35, 0.85, 0.55);
  leg->AddEntry(hFDEff, "FD MC Efficiency", "L");
  leg->AddEntry(hNDEff, "ND MC Efficiency", "L");
  leg->AddEntry(hFDPur, "FD MC Purity",     "L");
  leg->AddEntry(hNDPur, "ND MC Purity",     "L");

  TCanvas* c = new TCanvas(name.c_str(), "Purity and Efficiency", 800, 500);
  hFDEff->Draw("hist");
  hNDEff->Draw("hist same");
  hFDPur->Draw("hist same");
  hNDPur->Draw("hist same");
  leg->Draw();

  gPad->RedrawAxis();
  Simulation();

  rootOut->WriteTObject(c);
}


//........................................................................
void ResetAngles(osc::OscCalcSterile* calc)
{
  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
  calc->SetAngle(3, 4, 0.);
  calc->SetAngle(2, 4, 0.);

  return;
}