PlotNus17Prediction.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;
};

void TexInfo();

TLegend* MakeLegend(double x1, double y1, double x2, double y2, double textsize = 0.03)
{
  TLegend* leg = new TLegend(x1, y1, x2, y2);
  leg->SetBorderSize(0);
  leg->SetFillColor(kWhite);
  leg->SetFillStyle(0);
  leg->SetFillStyle(0);
  leg->SetTextFont(42);
  leg->SetTextSize(textsize);
  return leg;
}

void SigmaLines(double xmin, double xmax);

void PlotText(const double xpos, const double start, const double step);

// Reset calculator values after each fit to preset defaults
void ResetAngles(osc::OscCalcSterile* calc);

// 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

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 = {});


void CompareSlices(TH1* stats, TH1* systs, std::string angle, TDirectory* rootOut, double xmin, double xmax);

// 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 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 = {});

// 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);



void PlotNus17Prediction(bool plotSensitivities = false)
{
  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 = "/nova/ana/steriles/Nus17/";
  std::string filenm = "predInterpSysts_nus17_v1.root";
  

  std::string loadLocation = folder + filenm + ".root";
  std::string saveLocation = filenm + "Ana.root";
  std::string textLocation = 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 + ":/sCos1").c_str());
  Spectrum sCosmic              = *Spectrum          ::LoadFrom(gDirectory);

  loadDir->cd((loadLocation + ":/pFDPENum1").c_str());
  PredictionExtrap pFDPurEffNum = *PredictionExtrap  ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/pFDPDen1").c_str());
  PredictionExtrap pFDPurityDen = *PredictionExtrap  ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/pFDEDen1").c_str());
  PredictionExtrap pFDEfcncyDen = *PredictionExtrap  ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sNDPENum1").c_str());
  Spectrum sNDPurEffNum         = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sNDPDen1").c_str());
  Spectrum sNDPurityDen         = *Spectrum          ::LoadFrom(gDirectory);
  loadDir->cd((loadLocation + ":/sNDEDen1").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);
  ResetSterileCalcToDefault(calc3f);

  Spectrum prediction9 = pred.Predict(calc3f);

  const double livetime9  = 440.;
  const double pot9       = 8.85e20;
  const double livetime24 = 440.*(24/9);
  const double pot24      = 24e20;

  sCosmic.Scale((440.)/sCosmic.Livetime());
  sCosmic.OverridePOT(9e20);
  
  Spectrum sCosmic9  = sCosmic.FakeData(9.49e20);
  Spectrum sCosmic24 = sCosmic.FakeData(24e20);

  Spectrum sFDdata9 = prediction9.FakeData(8.85e20) + sCosmic9;
  
  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  ResetAngles(calc4f);
  
  std::string labelEReco = "Energy Deposited in (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.0698132); // 4 degrees, 10 degrees=0.17453292
  //calc4f->SetAngle(3, 4, 0.349066); // 20 degrees, 35 degrees=0.61086524

  calc4f->SetAngle(2, 4, 0.3630285); // 20.8 degrees, 10 degrees=0.17453292
  calc4f->SetAngle(3, 4, 0.5445427); // 31.2 degrees, 35 degrees=0.61086524

  //calc4f->SetAngle(2, 4, 0.3630285); // 4 degrees, 10 degrees=0.17453292
  //calc4f->SetAngle(3, 4, 0.3979351); // 22.8 degrees, 35 degrees=0.61086524

  // Plot spectra for 9e20
  std::string titleHelper = " Spectra for 9e20 POT";
  PlotStack(&pred, sCosmic9, calc3f, calc4f, rootF, textF, "FDSpectra9", titleHelper, 9);
  PlotStack(&decompNC, rootF, textF, "NDSpectra9", titleHelper, 9);
  
  // Plot spectra for 24e20
  titleHelper = " Spectra for 24e20 POT";
  PlotStack(&pred, sCosmic24, calc3f, calc4f, rootF, textF, "FDSpectra24", titleHelper, 24);
  PlotStack(&decompNC, rootF, textF, "NDSpectra24", titleHelper, 24);

  // Plot purity and efficiency
  PlotPurEff(&pFDPurEffNum, &pFDPurityDen, &pFDEfcncyDen,
             sNDPurEffNum, sNDPurityDen, sNDEfcncyDen,
             calc3f, sCosmic9, rootF, textF, "PurEff", labelEReco);
  
  fclose(textF); // Text file is done now, close it
  
  ResetAngles(calc4f);

  if(plotSensitivities){
    // Create Experiment objects, one each for the 1 data yr and 3 data yr predictions
    GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.57, 0.23);
    
    CountingExperiment expt9c(&pred, sFDdata9, sCosmic9);
    MultiExperiment multi9c({&th23Constraint, &expt9c});
    
    SingleSampleExperiment expt9s(&pred, sFDdata9, sCosmic9);
    MultiExperiment multi9s({&th23Constraint, &expt9s});
    
    // 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);
    fitvars.push_back(&kFitDelta24InPiUnitsSterile);
    
    // Create fit objects
    MinuitFitter fit9c(&multi9c, fitvars);
    MinuitFitter fit9s(&multi9s, fitvars);
    
    // Run fits!
    fit9c.Fit(calc4f);
    ResetAngles(calc4f);
    fit9s.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.;
    
    std::string fullname, special;
    
    const Color_t kStatsFitCol = kBlack;
    
    const Color_t kFitColor = kRed;
    
    fullname = "34";
    special  = "statsC";
    // Plot 1D fits for 9e20 experiment: Counting
    TH1* h34C = Plot1DSlice(&multi9c, calc4f,
                             &kFitTheta34InDegreesSterile,
                             {&kFitTheta23InDegreesSterile,
                                 &kFitTheta24InDegreesSterile,
                                 &kFitDelta24InPiUnitsSterile},
                             avals.nbins34, avals.min34, avals.max34,
                             rootF, fullname, special, kStatsFitCol
                             );
    
    special = "statsS";
    // Plot 1D fits for 9e20 experiment: Sample
    TH1* h34S = Plot1DSlice(&multi9s, calc4f,
                             &kFitTheta34InDegreesSterile,
                             {&kFitTheta23InDegreesSterile,
                                 &kFitTheta24InDegreesSterile,
                                 &kFitDelta24InPiUnitsSterile},
                             avals.nbins34, avals.min34, avals.max34,
                             rootF, fullname, special, kStatsFitCol
                             );
    
    CompareSlices(h34C, h34S, "34", rootF,
                  avals.min34, avals.max34);
    
    // Lower the number of bins for surfaces
    avals.nbins23 = 50;
    avals.nbins34 = 50;
    avals.nbins24 = 45;
    
    
    special = "34vs24";
    // Plot the surfaces for 9e20 counting experiment
    Plot2DSlice(&multi9c, &multi9s, calc4f,
                &kFitTheta34InDegreesSterile,
                avals.nbins34, avals.min34, avals.max34,
                &kFitTheta24InDegreesSterile,
                avals.nbins24, avals.min24, avals.max24,
                rootF, special,
                {&kFitTheta23InDegreesSterile,
                    &kFitDelta24InPiUnitsSterile}
                 );
        
  } // if(plotSensitivities)
  rootF->Close(); // Close the output root file
}

//........................................................................
TH1* Plot1DSlice(IExperiment* expt, osc::OscCalcSterile* calc,
                 const IFitVar* var, std::vector<const IFitVar*> profVars,
                 int nbins, double min, double max,
                 TDirectory* rootOut, std::string indices,
                 std::string special, Color_t color,
                 std::vector<const ISyst*> systs)
{
  // 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...


  std::string fullname = indices + special;
  // Formatting
  h->SetLineColor(color);
  h->SetMaximum(5);
  h->SetMinimum(0);
  h->SetName(fullname.c_str()); // Object name for saving
  rootOut->WriteTObject(h);     // Save raw histogram


  std::cout << fullname << std::endl;
  std::string cname  = "c1D" + fullname;
  std::string ctitle = "Theta" + indices + "_ChiSq1DSlice_" + special;
  std::string htitle = "#theta_{"+indices+"} (deg.)";
  h->GetXaxis()->SetTitle(htitle.c_str());

  // Actually draw things together
  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
  h->GetXaxis()->SetRangeUser(min+0.25,max);
  h->Draw("l");
  TexInfo();
  SigmaLines(min, max);
  PlotText(0.135, 0.82, 0.055);
  TLegend* leg = 0;
  // Special case legend position for each angle slice
  if(indices.compare("34") == 0)
    leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
  if(indices.compare("24") == 0)
    leg = MakeLegend(0.600, 0.35, 0.900, 0.50, 0.037);
  if(indices.compare("23") == 0)
    leg = MakeLegend(0.350, 0.35, 0.650, 0.50, 0.037);
  if(special.compare("stats") == 0)
    leg->AddEntry(h,"stat. only","L");
  if(special.compare("systs") == 0)
    leg->AddEntry(h,"syst.","L");
  leg->Draw();

  Preliminary();
  CenterTitles(h);
  rootOut->WriteTObject(c);

  return h;
}


void CompareSlices(TH1* stats, TH1* systs, std::string indices, TDirectory* rootOut, double min, double max)
{

  std::string cname  = "c1DComp" + indices;
  std::string ctitle = "Theta" + indices + " ChiSq1Dcomp";
  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
  stats->GetXaxis()->SetRangeUser(min+0.25,max);
  systs->GetXaxis()->SetRangeUser(min+0.25,max);
  stats->Draw("l");
  systs->Draw("l same");
  TexInfo();
  SigmaLines(min, max);
  Preliminary();
  PlotText(0.135, 0.82, 0.055);
  TLegend* leg = 0;
  // Special case legend for angle slice
  if(indices.compare("34") == 0)
    leg = MakeLegend(0.135, 0.35, 0.435, 0.50, 0.037);
  if(indices.compare("24") == 0)
    leg = MakeLegend(0.600, 0.35, 0.900, 0.50, 0.037);
  if(indices.compare("23") == 0)
    leg = MakeLegend(0.350, 0.35, 0.650, 0.50, 0.037);
  leg->AddEntry(stats,"stat. [rate]","L");
  //leg->AddEntry(systs,"syst.","L");
  leg->AddEntry(systs,"stat. [shape]","L");
  leg->Draw();
  rootOut->WriteTObject(c);

  return;
}


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 indices,
                 std::vector<const IFitVar*> vars,
                 std::vector<const ISyst*> systs)
{

  ResetAngles(calc);
  std::string cname      = "c2D" + indices;
  std::string cnamestat  = "c2D" + indices + "statC";
  std::string cnamesyst  = "c2D" + indices + "statS";
  std::string hnamestat  = "h" + indices + "statC";
  std::string hnamesyst  = "h" + indices + "statS";
  std::string ctitle     = "Surface for Theta" + indices;
  std::string ctitlestat = "Surface for Theta" + indices + " statC";
  std::string ctitlesyst = "Surface for Theta" + indices + " systS";

  FrequentistSurface surfStat(expt1, calc,
                   varA, nbinsA, minA, maxA,
                   varB, nbinsB, minB, maxB,
                   vars);
  ResetAngles(calc);
  
  FrequentistSurface surfSyst(expt2, calc,
                   varA, nbinsA, minA, maxA,
                   varB, nbinsB, minB, maxB,
                   vars);
  ResetAngles(calc);

  TH1* hSurfStat = surfStat.ToTH2();     // Create a histogram from the surfaces
  hSurfStat->SetName(hnamestat.c_str()); // Set the name for saving
  rootOut->WriteTObject(hSurfStat);      // Save the raw histogram

  TH1* hSurfSyst = surfSyst.ToTH2();     // Create a histogram from the surfaces
  hSurfSyst->SetName(hnamesyst.c_str()); // Set the name for saving
  rootOut->WriteTObject(hSurfSyst);      // Save the raw histogram

  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str());

  surfStat.DrawContour(Gaussian68Percent1D(surfStat), kDashed, kAzure+2);
  surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid,  kAzure+2);
  surfSyst.DrawContour(Gaussian68Percent1D(surfSyst), kDashed, kBlack);
  surfSyst.DrawContour(Gaussian90Percent1D(surfSyst), kSolid,  kBlack);

  PlotText(0.135, 0.82, 0.055);
  TLegend* legcomp = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
  TH1* dummy   = new TH1F("","",1,0,1);

  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kDashed);
  dummy->SetLineWidth(2);
  legcomp->AddEntry(dummy->Clone(),"68% C.L. (stat.) [rate]","L");
  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kSolid);
  legcomp->AddEntry(dummy->Clone(),"90% C.L. (stat.) [rate]","L");
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kDashed);
  legcomp->AddEntry(dummy->Clone(),"68% C.L. (stat.) [shape]","L");
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kSolid);
  legcomp->AddEntry(dummy->Clone(),"90% C.L. (stat.) [shape]","L");
  legcomp->Draw();
  Preliminary();
  rootOut->WriteTObject(c);

  // Plot statisics only
  TCanvas* cStat = new TCanvas(cnamestat.c_str(), ctitlestat.c_str());

  surfStat.DrawContour(Gaussian68Percent1D(surfStat), kDashed, kAzure+2);
  surfStat.DrawContour(Gaussian90Percent1D(surfStat), kSolid,  kAzure+2);

  PlotText(0.135, 0.82, 0.055);
  TLegend* legStat = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);

  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kDashed);
  dummy->SetLineWidth(2);
  legStat->AddEntry(dummy->Clone(),"68% C.L. (stat.) [rate]","L");
  dummy->SetLineColor(kAzure+2);
  dummy->SetLineStyle(kSolid);
  legStat->AddEntry(dummy->Clone(),"90% C.L. (stat.) [rate]","L");
  legStat->Draw();
  Preliminary();
  rootOut->WriteTObject(cStat);

  // Plot stats. + systs. only
  TCanvas* cSyst = new TCanvas(cnamesyst.c_str(), ctitlesyst.c_str());

  surfSyst.DrawContour(Gaussian68Percent1D(surfSyst), kDashed, kBlack);
  surfSyst.DrawContour(Gaussian90Percent1D(surfSyst), kSolid,  kBlack);

  PlotText(0.135, 0.82, 0.055);
  TLegend* legSyst = MakeLegend(0.500, 0.62, 0.800, 0.82, 0.037);
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kDashed);
  dummy->SetLineWidth(2);
  legSyst->AddEntry(dummy->Clone(),"68% C.L. (stat.) [shape]","L");
  dummy->SetLineColor(kBlack);
  dummy->SetLineStyle(kSolid);
  legSyst->AddEntry(dummy->Clone(),"90% C.L. (stat.) [shape]","L");
  legSyst->Draw();
  Preliminary();
  rootOut->WriteTObject(cSyst);

  return;
}

void PlotText(const double xpos, const double start, const double step)
{
  TLatex* tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xpos, start, "NOvA 9.0#times10^{20} POT");
  tex->DrawLatex(xpos, start - 1*step, "#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
  tex->DrawLatex(xpos, start - 2*step, "#theta_{13} = 8.5#circ, ^{}#theta_{23} = 45#circ");
  tex->DrawLatex(xpos, start - 3*step, "#Deltam^{2}_{41} = 0.5 eV^{2}");
}

void TexInfo()
{
  const double xtex = 0.76;
  const double y68  = 0.28;
  const double y90  = 0.55;

  TLatex *tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(xtex, y68, "68% C.L."); // 1-sigma
  tex->DrawLatex(xtex, y90, "90% C.L."); // 2-sigma
}

void SigmaLines(double xmin, double xmax)
{
  TLine* line68 = new TLine();
  TLine* line90 = new TLine();
  // Fixes horizontal effects at max canvas value
  TLine* line05 = new TLine();
  // Dashed line for 1 and 2 sigma lines
  line68->SetLineStyle(kDashed);
  line90->SetLineStyle(kDashed);

  line68->SetLineWidth(2);
  line90->SetLineWidth(2);
  line05->SetLineWidth(2);
  line68->DrawLine(xmin, 1,    xmax, 1);
  line90->DrawLine(xmin, 2.71, xmax, 2.71);
  line05->DrawLine(xmin, 5,    xmax, 5);
}


//........................................................................
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* hNutau = spectra[4].ToTH1(realPOT);
  TH1* hCos  = spectra[5].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("Visible Energy (GeV)");
  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, %12.12s\n",
          "NC (3 Flav)", "All", "Numu", "Nue", "Nutau", "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 nNutau = hNutau->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, %10E\n\n",
          nNC3F, nAll, nNumu, nNue, nNutau, nCos, fom);

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

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

  const Color_t kCosmicBackgroundColor = kOrange+7;
  hAll ->SetLineColor(kTotalMCColor);
  hNC  ->SetLineColor(kNCBackgroundColor);
  hNumu->SetLineColor(kNumuBackgroundColor);
  hNutau->SetLineColor(kBeamNueBackgroundColor);
  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());
  hNutau ->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);

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

  // Stylize the stacked plots
  //CenterTitles(hNutauStack);
  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);
  //hNutauStack ->SetFillColor(kBeamNueBackgroundColor);
  //hNutauStack ->SetLineColor(kBeamNueBackgroundColor);
  hCosStack ->SetLineColor(kCosmicBackgroundColor);
  hCosStack ->SetFillColor(kCosmicBackgroundColor);

  hNCStack  ->SetMinimum(0);
  hNumuStack->SetMinimum(0);
  hNueStack ->SetMinimum(0);
  //hNutauStack ->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);
  if(det.compare("FD") == 0){
     hNutau->SetMaximum(maxval*1.1);
     hNutau->GetXaxis()->SetNdivisions(-405);
  }
  else{
    hNue->SetMaximum(maxval*1.1);
    hNue->GetXaxis()->SetNdivisions(-405);
  } 

  // 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(hNutau,"#nu_{#tau} CC Background", "l");
    leg->AddEntry(hCos,  "Cosmic Background", "l");
  }
  if(det.compare("FD") == 0)
    hNutau->Draw("hist");
  else
    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} = 4#circ, ^{}#theta_{34} = 10#circ");
    tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.44x10^{-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);
  //if(det.compare("FD") == 0){
  //hNutauStack->SetMaximum(maxval*1.1);
  //hNutauStack->GetXaxis()->SetNdivisions(-405);
  //}
  //else{
    hNueStack->SetMaximum(maxval*1.1);
    hNueStack->GetXaxis()->SetNdivisions(-405);
    //} 

  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(hNutauStack,  "#nu_{#tau} CC Background", "f");
    legStack->AddEntry(hCosStack,  "Cosmic Background", "f");
  }

  //if(det.compare("FD") == 0) {
  //hNutauStack->Draw("hist");
  //}
  //else{
    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} = 4#circ, ^{}#theta_{34} = 10#circ");
    tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.44x10^{-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,
    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),
    pred->PredictComponent(calc, Flavors::kAllNuTau,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(-405);
  
  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)
{
  ResetSterileCalcToDefault(calc);
  calc->SetDm(4, 0.5);          // #deltam41 = 0.5 eV^2
  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
  calc->SetAngle(3, 4, 0.);
  calc->SetAngle(2, 4, 0.);

  return;
}