NuSPlotFunctions.h

Go to the documentation of this file.

/*
  This file contains functions for plotting
  In this file:

  TLatex* DrawLatexLines(double x, double y, std::string sPOT)

  void DrawLineArrow(double cut, bool left, TH1* h)

  void SetHistOptions(TH1* h, double max, std::string title,
                      int ndiv, Color_t col, bool fill)

  void SetLegendOptions(TLegend* leg)

  void PlotSpectra(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, std::string det, int POT,
                   std::string option = "");

  void PlotSpectra(Spectrum spectraND[],
                   Spectrum spectraFD[],
                   Spectrum spectraSt[],
                   Spectrum spectraUn[],
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option = "");

  void PlotSpectra(Spectrum spectraND[],
                   Spectrum spectraFD[],
                   Spectrum spectraFS[],
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option = "");

  void PlotSpectra(IDecomp* decomp,
                   IPrediction* predFD, osc::IOscCalc* calc,
                   IPrediction* predSt, osc::IOscCalc* calcst,
                   Spectrum sCos,
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option = "");

  void PlotSpectra(IDecomp* decomp,
                   IPrediction* predFD,
                   Spectrum sCos,
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option = "");

  void PlotSpectra(IPrediction* predFD, Spectrum sCos, osc::IOscCalc* calc,
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option = "");

  void PlotPurEff(Spectrum pureff[], Spectrum effcyD[],
                  TDirectory* out, std::string name, std::string title);

  void PlotPurEff(Spectrum pureffND[], Spectrum effcyDND[],
                  Spectrum pureffFD[], Spectrum effcyDFD[],
                  Spectrum pureffFS[], Spectrum effcyDFS[],
                  TDirectory* out, std::string name, std::string title);

  Make2DFrom1D functions may be deprecated from 2D Spectrum functionality improvements
  I haven't audited those changes and how they affect these functions yet, however

  void Make2DFrom1D(Spectrum spectra[],
                    const Binning& binsX, const Binning& binsY,
                    TDirectory* rootOut,
                    std::string name, std::string det,
                    std::string labelX, std::string labelY,
                    std::string option = "");

  void Make2DFrom1D(Spectrum spectraND[], Spectrum spectraFD[], Spectrum spectraFS[],
                    const Binning& binsX, const Binning& binsY,
                    TDirectory* rootOut,
                    std::string name,
                    std::string labelX, std::string labelY,
                    std::string option = "");
*/
#pragma once

#include <cmath>
#include <iostream>
#include <stdio.h>
#include <string>

#include "TArrow.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2D.h"
#include "THStack.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMath.h"

#include "CAFAna/Analysis/Calcs.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Decomp/IDecomp.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Vars/Vars.h"

#include "OscLib/IOscCalc.h"
#include "OscLib/OscCalcPMNSOpt.h"

#include "Utilities/func/MathUtil.h"
#include "Utilities/rootlogon.C"

#define MAXSPEC 6
// Put a "NOvA Preliminary" tag in the corner

namespace ana
{
  std::vector<std::string> flavors{"All", "NC", "NuMu", "NuE", "Cos"};

  static double NCSCALE = 9e20;

  struct PlotOptions {
    int ndiv;
    double xL;
    double xR;
    double yB;
    double yT;
    double texX;
    double texY;

    double cutL;
    double cutR;
  };

  /// Draw common TLatex lines
  TLatex* DrawLatexLines(double x, double y, std::string sPOT)
  {
    TLatex *tex = new TLatex();
    tex->SetNDC();
    tex->SetTextFont(42);
    tex->SetTextSize(0.037);
    tex->SetTextAlign(11);

    std::string potStr = sPOT + " #times 10^{20} POT-equiv.";

    tex->DrawLatex(x, y, "#Deltam^{2}_{32} = 2.44x10^{-3} eV^{2}");
    tex->DrawLatex(x, y - 0.055, "#theta_{13} = 8.5#circ, ^{}#theta_{23} = 45#circ");
    tex->DrawLatex(x, y - 0.110, potStr.c_str());

    return tex;
  }

  //......................................................................

  /// Draw a dashed line and an arrow where a cut is made
  /// \param cut The point at which the cut is made
  /// \param left Bool to determine whether arrow should start at left and move right (true)
  ///             or start at right and move left (false)
  void DrawLineArrow(double cut, bool left, TH1* h)
  {
    double axLen = h->GetXaxis()->GetXmax() - h->GetXaxis()->GetXmin();
    double ayTop = h->GetMaximum();
    double delta = axLen/20.;
    if(!left) { delta *= -1.; }
    TLine* line = new TLine(cut, 0, cut, ayTop);
    TArrow* arrow = new TArrow(cut,         ayTop*0.95,
                               cut + delta, ayTop*0.95,
                               0.02, "|>");
    line->SetLineColor(kRed);
    line->SetLineStyle(2);
    line->SetLineWidth(3);
    line->Draw("same");

    arrow->SetFillColor(kRed);
    arrow->SetFillStyle(1001);
    arrow->SetLineColor(kRed);
    arrow->SetLineWidth(3);
    arrow->Draw();
  }

  //......................................................................

  /// Set common options for a TLegend
  void SetHistOptions(TH1* h, double max, std::string title,
                      int ndiv, Color_t col, bool fill)
  {
    CenterTitles(h);
    h->SetLineColor(col);
    if(fill) {
      FillWithDimColor(h);
    }
    h->SetMinimum(0);
    h->SetMaximum(max*1.05);
    if(ndiv) { h->GetXaxis()->SetNdivisions(ndiv); }
    h->SetTitle(title.c_str());
  }

  //......................................................................

  /// Set common options for a TLegend
  void SetLegendOptions(TLegend* leg)
  {
    leg->SetBorderSize(0);
    leg->SetFillColor(kWhite);
    leg->SetFillStyle(0);
    leg->SetFillStyle(0);
    leg->SetTextFont(42);
    leg->SetTextSize(0.037);
  }

  //......................................................................
  void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
                 std::string name, std::string det,
                 double POT, double potEquiv,
                 PlotOptions opt, Spectrum* dataspec)
  {
    char pot[4];
    sprintf(pot, "%4.2f", potEquiv);
    std::string spot = std::string(pot);
    std::string potStr = spot + " #times 10^{20} POT-equiv.";
    std::string NDpotStr = spot + " #times 10^{20} POT";

    std::cout << "POT: " << POT << std::endl;
    // Get each histogram from the Spectrum array
    TH1* hAll  = spectra[0].ToTH1(POT);
    TH1* hNC   = spectra[1].ToTH1(POT);
    TH1* hNumu = spectra[2].ToTH1(POT);
    TH1* hNue  = spectra[3].ToTH1(POT);
    TH1* hNutau = spectra[4].ToTH1(POT);
    TH1* hCos  = spectra[5].ToTH1(POT);
    // The default for these inputs are nullptr, only get a TH1 if they exist
    TH1* hData = (dataspec ? dataspec->ToTH1(dataspec->POT()) : nullptr);

    // Create the title and axis labels, which is an empty title,
    // the x label from an input Spectrum, and Events/POT
    std::string fullname = det + name;
    std::string xLabel(hNC->GetXaxis()->GetTitle());
    if(name.find("PTP") != std::string::npos) { xLabel = "p_{T}/p"; }
    if(name.find("CVN") != std::string::npos) { xLabel = "CVN NC Classifier"; }
    std::string yLabel = (det.compare("ND") == 0 ? "10^{3} Events": "Events");
    if(det.compare("FD") == 0){
      if(name.find("CalE") != std::string::npos ||
         name.find("ECal") != std::string::npos) { yLabel += " / 0.25 GeV"; }
    }

    if(det.compare("ND") == 0){
      if(name.find("CalE") != std::string::npos ||
         name.find("ECal") != std::string::npos) { yLabel += " / 0.25 GeV / "+NDpotStr; }
      else{
        yLabel += " / "+NDpotStr;
          }
    }
    std::string histTitle = ";" + xLabel + ";" + yLabel;

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

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

    int ndiv = opt.ndiv;

    // 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*)hNueStack->Clone();
    hNCStack->Add(hNC);

    // Manually set errors, and calculate the max value (including errors)
    double maxval = hAll->GetBinContent(hAll->GetMaximumBin());
    for(int i = 1, n = hNC->GetNbinsX(); i <= n; ++i) {
      double error = sqrt(hNC->GetBinContent(i));
      //hNC     ->SetBinError(i, error);
      //hNCStack->SetBinError(i, error);
      if(hData) {
        double data_value = hData->GetBinContent(i);
        double data_error = sqrt(data_value);
        //hData->SetBinError(i, data_error);
        maxval = std::max(maxval, data_value + data_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;
      if(hData) { hData->Scale(ndScale); }
    }

    // Stylize the unstacked plots
    const Color_t kCosmicBackgroundColor = kOrange+7;
    SetHistOptions(hNue,  maxval, histTitle, ndiv, kNueSignalColor,        false);
    SetHistOptions(hNumu, maxval, histTitle, ndiv, kNumuBackgroundColor,   false);
    SetHistOptions(hNC,   maxval, histTitle, ndiv, kNCBackgroundColor,     false);
    SetHistOptions(hCos,  maxval, histTitle, ndiv, kCosmicBackgroundColor, false);
    SetHistOptions(hAll,  maxval, histTitle, ndiv, kTotalMCColor,          false);
    if(hData) { SetHistOptions(hData, maxval, histTitle, ndiv, kBlack, false); }

    if(name.find("CosBDT") != std::string::npos) {
      hAll ->GetXaxis()->SetRangeUser(0.31, 0.7);
      hNC  ->GetXaxis()->SetRangeUser(0.31, 0.7);
      hNumu->GetXaxis()->SetRangeUser(0.31, 0.7);
      hNue ->GetXaxis()->SetRangeUser(0.31, 0.7);
      hCos ->GetXaxis()->SetRangeUser(0.31, 0.7);
      hNCStack  ->GetXaxis()->SetRangeUser(0.31, 0.7);
      hNueStack ->GetXaxis()->SetRangeUser(0.31, 0.7);
      hNumuStack->GetXaxis()->SetRangeUser(0.31, 0.7);
      hCosStack ->GetXaxis()->SetRangeUser(0.31, 0.7);
      if(hData) { hData->GetXaxis()->SetRangeUser(0.31, 0.7); }
    }

    if(hData) { hData->SetMarkerStyle(kFullCircle); }

    // Stylize the stacked plots
    SetHistOptions(hNueStack,  maxval, histTitle, ndiv, kNueSignalColor,        true);
    SetHistOptions(hNumuStack, maxval, histTitle, ndiv, kNumuBackgroundColor,   true);
    SetHistOptions(hNCStack,   maxval, histTitle, ndiv, kNCBackgroundColor,     true);
    SetHistOptions(hCosStack,  maxval, histTitle, ndiv, kCosmicBackgroundColor, true);

    std::string dataleg = ( (det.compare("ND") == 0) ? "ND Data" : "FD Data" );

    // Create a canvas for and draw the unstacked histograms
    TCanvas* c = new TCanvas(fullname.c_str(), fullname.c_str(), 800, 500);
    gPad->SetFillStyle(0);

    // Create a legend for the unstacked histograms
    double xL = opt.xL, xR = opt.xR, yB = opt.yB, yT = opt.yT;
    TLegend* leg = new TLegend(xL, yB, xR, yT);
    SetLegendOptions(leg);
    if(hData) { leg->AddEntry(hData, dataleg.c_str(), "lp"); }
    leg->AddEntry(hAll,  "Total Prediction", "l");
    leg->AddEntry(hNC,   "NC 3 Flavor Prediction", "l");
    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");
    }
    leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);
    if(std::abs(opt.texX - xL) < 0.2) { opt.texY = leg->GetY1() - 0.045; }

    hNue ->Draw("hist");
    hNumu->Draw("hist same");
    if(det.compare("FD") == 0) { hCos->Draw("hist same"); }
    hAll ->Draw("hist same");
    hNC  ->Draw("hist same");
    if(hData) { hData->Draw("same"); }
    if(opt.cutL != 1.) { DrawLineArrow(opt.cutL, true,  hNue); }
    if(opt.cutR != 0.) { DrawLineArrow(opt.cutR, false, hNue); }
    leg->Draw();
    if(det.compare("FD") == 0) { DrawLatexLines(opt.texX, opt.texY, spot); }

    gPad->RedrawAxis();
    if(hData) { Preliminary(); }
    else {      Simulation();  }

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

    // Create a canvas for and draw the stacked histograms
    fullname += "Stack";
    TCanvas* cStack = new TCanvas(fullname.c_str(), fullname.c_str(), 800, 500);
    gPad->SetFillStyle(0);

    TLegend* legStack = new TLegend(xL, yB, xR, yT);
    SetLegendOptions(legStack);
    if(hData) { legStack->AddEntry(hData, dataleg.c_str(), "lp"); }
    legStack->AddEntry(hNCStack,   "NC 3 Flavor Prediction", "f");
    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");
    }
    legStack->SetY1(legStack->GetY2() - legStack->GetNRows()*0.05);
    if(std::abs(opt.texX - xL) < 0.2) { opt.texY = leg->GetY1() - 0.045; }

    hNCStack  ->Draw("hist");
    hNueStack ->Draw("hist same");
    hNumuStack->Draw("hist same");
    if(det.compare("FD") == 0) { hCosStack->Draw("hist same"); }
    if(hData) { hData->Draw("same"); }
    if(opt.cutL != 1.) { DrawLineArrow(opt.cutL, true,  hNCStack); }
    if(opt.cutR != 0.) { DrawLineArrow(opt.cutR, false, hNCStack); }
    legStack->Draw();
    if(det.compare("FD") == 0) { DrawLatexLines(opt.texX, opt.texY, spot); }
    gPad->RedrawAxis();
    if(hData) { Preliminary(); }
    else {      Simulation();  }
    
    rootOut->WriteTObject(cStack);

    return;
  }

  //......................................................................

  /// Function for NuS group to plot spectra
  /// You probably won't call this version directly,
  /// but here are the parameters commmon to all versions:
  /// \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 option Additional options; see the function definition for its uses
  void PlotSpectra(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, std::string det, int POT,
                   std::string option)
  {
    // The input is in units of 1e20, turn this into the real scale
    double realPOT = (double)POT*1e20;
    double livetime = realPOT/2.5e18;

    char potBuff[2];
    sprintf(potBuff, "%d", POT);
    std::string spot = std::string(potBuff);
    std::string potStr = spot + " #times 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(livetime, kLivetime);

    // 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" : "Event");
    std::string histTitle = ";" + xLabel + ";" + yLabel;
    
    if(option.find("label") != std::string::npos) {
      fprintf(textOFS, "      %12.12s, %12.12s, %8.8s, %12.12s, %12.12s, %12.12s, %12.12s\n",
              "NC", "All", "Numu", "Nue", "Cosmic", "Purity", "FOM");
    }
    
    double nNC   = hNC ->Integral();
    double nAll  = hAll->Integral();
    double nNumu = hNumu->Integral();
    double nNue  = hNue ->Integral();
    double nCos  = hCos->Integral();
    double fom   = nNC/sqrt(nAll);
    fprintf(textOFS, "%2.2s -- %10E, %10E, %10E, %10E, %10E, %.6f, %10E\n",
            det.c_str(), nNC, nAll, nNumu, nNue, nCos, nNC/nAll, fom);

    // 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*)hNueStack->Clone();
    hNCStack->Add(hNC);

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

    if(det.compare("ND") == 0) {
      double ndScale = 1./1000.;
      hAll ->Scale(ndScale);
      hNC  ->Scale(ndScale);
      hNumu->Scale(ndScale);
      hNue ->Scale(ndScale);
      hCos ->Scale(ndScale);
      hNCStack  ->Scale(ndScale);
      hNueStack ->Scale(ndScale);
      hNumuStack->Scale(ndScale);
      maxval *= ndScale;
      maxvalstack *= ndScale;
    }

    if(option.find("logy") != std::string::npos) {
      maxval = maxval*2/1.05;
      maxvalstack = maxvalstack*2/1.05;
    }

    const Color_t kCosmicBackgroundColor = kOrange+7;
    SetHistOptions(hAll,  maxval, histTitle, -404, kTotalMCColor,          false);
    SetHistOptions(hNC,   maxval, histTitle, -404, kNCBackgroundColor,     false);
    SetHistOptions(hNue,  maxval, histTitle, -404, kNueSignalColor,        false);
    SetHistOptions(hNumu, maxval, histTitle, -404, kNumuBackgroundColor,   false);
    SetHistOptions(hCos,  maxval, histTitle, -404, kCosmicBackgroundColor, false);

    // Stylize the stacked plots
    SetHistOptions(hNCStack,   maxval, histTitle, -404, kNCBackgroundColor,     true);
    SetHistOptions(hNueStack,  maxval, histTitle, -404, kNueSignalColor,        true);
    SetHistOptions(hNumuStack, maxval, histTitle, -404, kNumuBackgroundColor,   true);
    SetHistOptions(hCosStack,  maxval, histTitle, -404, kCosmicBackgroundColor, true);

    if(option.find("logy") != std::string::npos) {
      hAll ->SetMinimum(0.001);
      hNC  ->SetMinimum(0.001);
      hNumu->SetMinimum(0.001);
      hNue ->SetMinimum(0.001);
      hCos ->SetMinimum(0.001);
      hNCStack  ->SetMinimum(0.001);
      hNumuStack->SetMinimum(0.001);
      hNueStack ->SetMinimum(0.001);
      hCosStack ->SetMinimum(0.001);
    }

    double xL = 0.5, xR = 0.7, yB = 0.64, yT = 0.84;
    
    if(option.find("area") != std::string::npos) {
      TLegend* leg = new TLegend(xL, yB, xR, yT);
      leg->AddEntry(hNC,   "NC Signal", "L");
      leg->AddEntry(hNue,  "#nu_{e} CC Background", "L");
      leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "L");
      if(det.compare("ND") != 0) { leg->AddEntry(hCos,  "Cosmic Background", "L"); }
      leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);

      hNC  ->Scale(1./nNC);
      hNumu->Scale(1./nNumu);
      hNue ->Scale(1./nNue);
      if(nCos > 0.) { hCos->Scale(1./nCos); }

      double maxval = 0.;
      maxval = std::max(maxval, hNC  ->GetBinContent(hNC  ->GetMaximumBin()));
      maxval = std::max(maxval, hNumu->GetBinContent(hNumu->GetMaximumBin()));
      maxval = std::max(maxval, hNue ->GetBinContent(hNue ->GetMaximumBin()));
      maxval = std::max(maxval, hCos ->GetBinContent(hCos ->GetMaximumBin()));

      maxval *= 1.1;
      hNC  ->SetMaximum(maxval);
      hNumu->SetMaximum(maxval);
      hNue ->SetMaximum(maxval);
      hCos ->SetMaximum(maxval);
      
      TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
      gPad->SetFillStyle(0);
      hNC  ->Draw("hist");
      hNumu->Draw("hist same");
      hNue ->Draw("hist same");
      if(det.compare("ND") != 0) { hCos ->Draw("hist same"); }
      leg->Draw();

      Simulation();
      
      rootOut->WriteTObject(c);
      
      return;
    }

    // Create a legend for the histograms
    TLegend* leg = new TLegend(xL, yB, xR, yT);
    SetLegendOptions(leg);
    leg->AddEntry(hAll,  "Total Prediction", "l");
    if(option.find("hist") != std::string::npos) {
      leg->AddEntry(hNC, "NC Prediction", "l");
    }
    else {
      leg->AddEntry(hNC, "NC Prediction", "le");
    }
    leg->AddEntry(hNue,  "#nu_{e} CC Background", "l");
    leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
    if(det.compare("ND") != 0) {
      leg->AddEntry(hCos,  "Cosmic Background", "l");
    }
    leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);

    TLegend* legStack = new TLegend(xL, yB, xR, yT);
    SetLegendOptions(legStack);
    if(option.find("hist") != std::string::npos) {
      legStack->AddEntry(hNCStack, "NC Prediction", "f");
    }
    else {
      legStack->AddEntry(hNCStack, "NC Prediction", "le");
    }
    legStack->AddEntry(hNueStack,  "#nu_{e} CC Background", "f");
    legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
    if(det.compare("ND") != 0) {
      legStack->AddEntry(hCosStack,  "Cosmic Background", "f");
    }
    legStack->SetY1(legStack->GetY2() - legStack->GetNRows()*0.05);

    // Create a canvas
    TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
    gPad->SetFillStyle(0);

    hNue ->Draw("hist");
    hNumu->Draw("hist same");
    if(det.compare("ND") != 0) { hCos ->Draw("hist same"); }

    if(option.find("hist") != std::string::npos) { hNC->Draw("hist same"); }
    else { hNC->Draw("same"); }

    hAll ->Draw("hist same");
    leg->Draw();
    if(det.compare("FD") == 0) { DrawLatexLines(xL + .025, leg->GetY1() - 0.045, spot); }
    if(option.find("logy") != std::string::npos) { c->SetLogy(); }

    Simulation();
    rootOut->WriteTObject(c);

    TCanvas* cStack = new TCanvas((name + "Stack").c_str(), title.c_str(), 800, 500);
    gPad->SetFillStyle(0);

    hNueStack ->Draw("hist");
    hNumuStack->Draw("hist same");
    if(det.compare("ND") != 0) { hCosStack ->Draw("hist same"); }

    if(option.find("hist") != std::string::npos) { hNCStack->Draw("hist same"); }
    else { hNCStack->Draw("same"); }

    legStack->Draw();
    if(det.compare("FD") == 0) { DrawLatexLines(xL + .025, legStack->GetY1() - 0.045, spot); }
    if(option.find("logy") != std::string::npos) { cStack->SetLogy(); }

    Simulation();

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

    return;
  }

  //......................................................................

  /// Make canvases for ND, FD with 3 flavor osc, FD with 4 flavor osc, FD with no osc
  void PlotSpectra(Spectrum spectraND[],
                   Spectrum spectraFD[],
                   Spectrum spectraSt[],
                   Spectrum spectraUn[],
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option)
  {
    fprintf(textOFS, "%s\n", name.c_str());

    std::string label = "ND " + title;
    std::string fullname = name + "ND";
    std::string fulloption = option + "label";
    PlotSpectra(spectraND, rootOut, textOFS, fullname, label, "ND", POT, fulloption);

    label = "FD " + title;
    fullname = name + "FD";
    PlotSpectra(spectraFD, rootOut, textOFS, fullname, label, "FD", POT, option);

    label = "FD (sterile) " + title;
    fullname = name + "St";
    PlotSpectra(spectraSt, rootOut, textOFS, fullname, label, "St", POT, option);

    label = "FD (Unoscillated) " + title;
    fullname = name + "Un";
    PlotSpectra(spectraUn, rootOut, textOFS, fullname, label, "Un", POT, option);

    fprintf(textOFS, "\n\n");

    return;
  }

  //......................................................................

  /// Make canvases for ND, FD with 3 flavor orsc, FD with 4 flavor osc
  void PlotSpectra(Spectrum spectraND[],
                   Spectrum spectraFD[],
                   Spectrum spectraFS[],
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option)
  {
    fprintf(textOFS, "%s\n", name.c_str());

    std::string label = "ND " + title;
    std::string fullname = name + "ND";
    std::string fulloption = option + "label";
    PlotSpectra(spectraND, rootOut, textOFS, fullname, label, "ND", POT, fulloption);

    label = "FD " + title;
    fullname = name + "FD";
    PlotSpectra(spectraFD, rootOut, textOFS, fullname, label, "FD", POT, option);

    label = "FD (swap) " + title;
    fullname = name + "FS";
    PlotSpectra(spectraFS, rootOut, textOFS, fullname, label, "FS", POT, option);

    fprintf(textOFS, "\n\n");

    return;
  }

  //......................................................................

  /// Make canvases for ND, FD with 3 flavor osc, FD with 4 flavor osc, FD with no osc
  /// Directly use the CAFAna object inputs
  void PlotSpectra(IDecomp* decomp,
                   IPrediction* predFD, osc::IOscCalc* calc,
                   IPrediction* predSt, osc::IOscCalc* calcst,
                   Spectrum sCos,
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option)
  {
    Spectrum empty = decomp->NCTotalComponent();
    empty.Clear();

    if(!sCos.POT()) {
      sCos.OverridePOT(sCos.Livetime()*2.5e18);
    }

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

    Spectrum spectraFD[MAXSPEC] = {
      predFD->Predict(calc) + sCos,
      predFD->PredictComponent(calc, Flavors::kAll,     Current::kNC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuTau,  Current::kCC, Sign::kBoth),
      sCos
    };

    Spectrum spectraSt[MAXSPEC] = {
      predSt->Predict(calcst) + sCos,
      predSt->PredictComponent(calcst, Flavors::kAll,     Current::kNC, Sign::kBoth),
      predSt->PredictComponent(calcst, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
      predSt->PredictComponent(calcst, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
      predSt->PredictComponent(calc, Flavors::kAllNuTau,  Current::kCC, Sign::kBoth),
      sCos
    };

    osc::NoOscillations noosc;

    Spectrum spectraUn[MAXSPEC] = {
      predFD->Predict(&noosc) + sCos,
      predFD->PredictComponent(&noosc, Flavors::kAll,     Current::kNC, Sign::kBoth),
      predFD->PredictComponent(&noosc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
      predFD->PredictComponent(&noosc, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
      predFD->PredictComponent(&noosc, Flavors::kAllNuTau,  Current::kCC, Sign::kBoth),
      sCos
    };

    PlotSpectra(spectraND, spectraFD, spectraSt, spectraUn,
                rootOut, textOFS, name, title, POT, option);

    return;
  }

  //......................................................................

  /// Make canvases for ND, FD with 3 flavor osc, FD with no osc
  /// Directly use the CAFAna object inputs
  void PlotSpectra(IDecomp* decomp,
                   IPrediction* predFD,
                   Spectrum sCos,
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT,
                   std::string option)
  {
    Spectrum empty = decomp->NCTotalComponent();
    empty.Clear();

    if(!sCos.POT()) {
      sCos.OverridePOT(sCos.Livetime()*2.5e18);
    }

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

    osc::NoOscillations noosc;
    osc::OscCalcPMNSOpt* calc = (osc::OscCalcPMNSOpt*)DefaultOscCalc();
    Spectrum spectraFD[MAXSPEC] = {
      predFD->Predict(calc) + sCos,
      predFD->PredictComponent(calc, Flavors::kAll,     Current::kNC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuTau,  Current::kCC, Sign::kBoth),
      sCos
    };

    Spectrum spectraUn[MAXSPEC] = {
      predFD->Predict(&noosc) + sCos,
      predFD->PredictComponent(&noosc, Flavors::kAll,     Current::kNC, Sign::kBoth),
      predFD->PredictComponent(&noosc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
      predFD->PredictComponent(&noosc, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuTau,  Current::kCC, Sign::kBoth),
      sCos
    };

    fprintf(textOFS, "%s\n", name.c_str());

    std::string label = "ND " + title;
    std::string fullname = name + "ND";
    std::string fulloption = option + "label";
    PlotSpectra(spectraND, rootOut, textOFS, fullname, label, "ND", POT, fulloption);

    label = "FD " + title;
    fullname = name + "FD";
    PlotSpectra(spectraFD, rootOut, textOFS, fullname, label, "FD", POT, option);

    label = "FD (Unoscillated) " + title;
    fullname = name + "Un";
    PlotSpectra(spectraUn, rootOut, textOFS, fullname, label, "Un", POT, option);

    fprintf(textOFS, "\n\n");

    return;
  }

  //......................................................................

  /// Make canvas for FD with 3 flavor osc directly from CAFAna object inputs
  void PlotSpectra(IPrediction* predFD, Spectrum sCos, osc::IOscCalc* calc,
                   TDirectory* rootOut, FILE* textOFS,
                   std::string name, std::string title, int POT, std::string option)
  {
    fprintf(textOFS, "%s\n", name.c_str());

    if(!sCos.POT()) {
      sCos.OverridePOT(sCos.Livetime()*2.5e18);
    }

    Spectrum spectraFD[MAXSPEC] = {
      predFD->Predict(calc) + sCos,
      predFD->PredictComponent(calc, Flavors::kAll,     Current::kNC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuE,  Current::kCC, Sign::kBoth),
      predFD->PredictComponent(calc, Flavors::kAllNuTau,  Current::kCC, Sign::kBoth),
      sCos
    };

    PlotSpectra(spectraFD, rootOut, textOFS, name, title, "FD", POT, option);

    return;
  }

  //......................................................................
  
  /// Function for NuS group to plot purity/efficiency
  /// Purity is pureff[1]/pureff[0], Efficiency is pureff[1]/effcyD[1]
  /// \param pureff[] array of spectra, should match the output of the SPECARR macro
  /// \param effcyD[] array of spectra, should match the output of the SPECARR macro
  /// \param out Already open file for root 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(Spectrum pureff[], Spectrum effcyD[],
                  TDirectory* out, std::string name, std::string title)
  {
    std::string fullname = name + "PurEff";
    std::string fulltitle = title + " Purity and Efficiency";

    TH1* hpur = pureff[1].ToTH1(NCSCALE);
    TH1* hpurD = pureff[0].ToTH1(NCSCALE);
    hpur->Divide(hpurD);
    TH1* heff = pureff[1].ToTH1(NCSCALE);
    TH1* heffD = effcyD[1].ToTH1(NCSCALE);
    heff->Divide(heffD);

    heff->SetLineColor(kRed);

    TH1* hnum = pureff[1].ToTH1(NCSCALE);
    double numerator = 0., efcncyD = 0., purityD = 0.;
    for(int i = 1, n = hnum->GetNbinsX(); i <= n; ++i) {
      numerator +=  hnum->GetBinContent(i);
      efcncyD   += heffD->GetBinContent(i);
      purityD   += hpurD->GetBinContent(i);
    }
    double efcncy = numerator/efcncyD;
    double purity = numerator/purityD;
    std::cout << title << " Efficiency, Purity: " << efcncy << ", " << purity << std::endl;

    TLegend* leg = new TLegend(0.5, 0.8, 0.8, 0.9);
    leg->AddEntry(hpur, "Purity",     "L");
    leg->AddEntry(heff, "Efficiency", "L");

    TCanvas* c = new TCanvas(fullname.c_str(), fulltitle.c_str(), 800, 500);
    gPad->SetFillStyle(0);
    if(hpur->GetBinContent(hpur->GetMaximumBin()) > heff->GetBinContent(heff->GetMaximumBin())) {
      hpur->SetTitle(fulltitle.c_str());
      hpur->Draw("hist");
      heff->Draw("hist same");
    }
    else {
      heff->SetTitle(fulltitle.c_str());
      heff->Draw("hist");
      hpur->Draw("hist same");
    }

    leg->Draw();
    out->WriteTObject(c);

    return;
  }

  //......................................................................

  /// Plot pur/eff for ND, FD assuming 3 flavor osc, FD assuming 4 flavor osc
  void PlotPurEff(Spectrum pureffND[], Spectrum effcyDND[],
                  Spectrum pureffFD[], Spectrum effcyDFD[],
                  Spectrum pureffFS[], Spectrum effcyDFS[],
                  TDirectory* out, std::string name, std::string title)
  {
    std::string label = "ND " + title;
    std::string fullname = name + "ND";
    PlotPurEff(pureffND, effcyDND, out, fullname, label);

    label = "FD " + title;
    fullname = name + "FD";
    PlotPurEff(pureffFD, effcyDFD, out, fullname, label);

    label = "FD (swap) " + title;
    fullname = name + "FS";
    PlotPurEff(pureffFS, effcyDFS, out, fullname, label);

    return;
  }

  //......................................................................

  /// Function for NuS group to create output plots for 'pseudo-2D spectra'
  /// spectra Array of spectra, should match the output of the SPECARR2D macro
  /// \param binsX Binning that describes the desired x axis
  /// \param binsY Binning that describes the desired y axis
  /// \param rootOut Already open file for root output
  /// \param name A name for the output, only part of the saved output object's name
  /// \param det The detector the spectra correspond to, should be 2 chars long
  /// \param labelX String for the x axis label
  /// \param labelY String for the y axis label
  /// \param option Additional options; see the function definition for its uses
  void Make2DFrom1D(Spectrum spectra[],
                    const Binning& binsX, const Binning& binsY,
                    TDirectory* rootOut,
                    std::string name, std::string det,
                    std::string labelX, std::string labelY,
                    std::string option)
  {
    spectra[0] += spectra[MAXSPEC - 1];

    int NX = binsX.NBins();
    int NY = binsY.NBins();

    std::string titlehelper = det + " " + labelY + " vs " + labelX + ", ";
    std::string labelhelper = ";" + labelX + ";" + labelY;
    std::string fulltitle = titlehelper;
    std::string namehelper = name + det;
    std::string fullname = namehelper;

    std::vector<TH2*> h2;

    for(unsigned int i = 0, n = MAXSPEC; i < n; ++i) {
      h2.push_back(spectra[i].ToTH2(NCSCALE));

      fulltitle = titlehelper + flavors[i] + labelhelper;
      fullname  = namehelper  + flavors[i];
      h2[i]->SetNameTitle(fullname.c_str(), fulltitle.c_str());

      if(option.find("norm") != std::string::npos) {
        for(int i_x = 1; i_x <= NX; ++i_x) {
          double column_norm = 0.;

          for(int i_y = 1; i_y <= NY; ++i_y) {
            column_norm += h2[i]->GetBinContent(i_x, i_y);
          }

          for(int i_y = 1; i_y <= NY; ++i_y) {
            double h2bin = ( (column_norm > 0.) ? h2[i]->GetBinContent(i_x, i_y)/column_norm : 0. );
            h2[i]->SetBinContent(i_x, i_y, h2bin);
          }
        } // end of for loop over x bins
      }
    } // end of loop over component plots

    if(option.find("avg") != std::string::npos) {
      std::vector<TH1*> hp;

      fulltitle = titlehelper + labelhelper;
      fullname  = namehelper;

      for(int i = 0, n = MAXSPEC; i < n; ++i) {
        hp.push_back(h2[i]->ProjectionX(fullname.c_str(), 1, NY));

        TAxis* ay = h2[i]->GetYaxis();
        for(int i_x = 1; i_x <= NX; ++i_x) {
          double weighted_val = 0.;

          for(int i_y = 1; i_y <= NY; ++i_y) {
            weighted_val += h2[i]->GetBinContent(i_x, i_y) * ay->GetBinCenter(i_y);
          }

          double val_norm = hp[i]->GetBinContent(i_x);
          hp[i]->SetBinContent(i_x, ((val_norm > 0.) ? weighted_val/val_norm : 0.));
        }

        hp[i]->SetTitle(fulltitle.c_str());
      }

      const Color_t kCosmicBackgroundColor = kOrange+1;
      hp[0]->SetLineColor(kTotalMCColor);
      hp[0]->SetLineStyle(2);
      hp[1]->SetLineColor(kNCBackgroundColor);
      hp[2]->SetLineColor(kNumuBackgroundColor);
      hp[3]->SetLineColor(kNueSignalColor);
      hp[4]->SetLineColor(kCosmicBackgroundColor);

      TLegend* leg = new TLegend(0.5, 0.7, 0.8, 0.9);
      double maxval = 0.;

      for(int i = 0, n = MAXSPEC; i < n; ++i) {
        leg->AddEntry(hp[i], flavors[i].c_str(), "L");

        if(hp[i]->GetBinContent(hp[i]->GetMaximumBin()) > maxval) {
          maxval = hp[i]->GetBinContent(hp[i]->GetMaximumBin());
        }
      }

      fulltitle = titlehelper;
      TCanvas* c = new TCanvas(fullname.c_str(), fulltitle.c_str(), 800, 500);
      gPad->SetFillStyle(0);

      maxval *= 1.1;
      hp[0]->SetMaximum(maxval);
      hp[0]->Draw("hist");

      for(int i = 1, n = MAXSPEC; i < n; ++i) {
        CenterTitles(hp[i]);
        hp[i]->SetMaximum(maxval);
        hp[i]->Draw("hist same");
      }

      leg->Draw();

      Simulation();
      rootOut->WriteTObject(c);
    }
    else {
      for(int i = 0, n = MAXSPEC; i < n; ++i) {
        fulltitle = titlehelper + flavors[i];
        fullname  = namehelper  + flavors[i];

        TCanvas* c = new TCanvas(fullname.c_str(), fulltitle.c_str(), 800, 500);
        gPad->SetFillStyle(0);

        CenterTitles(h2[i]);
        h2[i]->Draw("colz");
        Simulation();

        rootOut->WriteTObject(c);
      }
    }

    return;
  }

  //......................................................................

  /// Create output for 'pseudo-2D spectra' at ND, FD with 3 flavor osc, FD with 4 flavor osc
  void Make2DFrom1D(Spectrum spectraND[], Spectrum spectraFD[], Spectrum spectraFS[],
                    const Binning& binsX, const Binning& binsY,
                    TDirectory* rootOut,
                    std::string name,
                    std::string labelX, std::string labelY,
                    std::string option)
  {
    Make2DFrom1D(spectraND, binsX, binsY, rootOut, name, "ND", labelX, labelY, option);
    Make2DFrom1D(spectraFD, binsX, binsY, rootOut, name, "FD", labelX, labelY, option);
    Make2DFrom1D(spectraFS, binsX, binsY, rootOut, name, "FS", labelX, labelY, option);
  }
}
//#endif