hadEFrac_nd_data_mc_systs.C

Go to the documentation of this file.

// Fills spectra for systematic error band.  Makes both pot norm with full syst band and area norm with shape only syst band
// have to manually make nd_syst_plots directory before running

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
//#include "Utilities/rootlogon.C"

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TGaxis.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TLatex.h"
#include "TStyle.h"
#include "THStack.h"
#include "TPaveText.h"
#include "TList.h"
#include "TGaxis.h"
#include "TAttLine.h"
#include "TAttMarker.h"

#include <cmath>
#include <iostream>
#include <vector>
#include <list>
#include <sstream>
#include <string>
#include <sstream>
#include <fstream>
#include <iomanip>



using namespace ana;


// weight bins by bin width for neutrino energy plots:
const Var kBinWidthWeight([](const caf::SRProxy *sr){
    const double Enu = kCCE(sr); // If want bin normalised
    //const double Enu = 6; // If don't want bin normalised.
    if(Enu < 0.75) return 0.1 / 0.75;
    if(Enu >= 0.75 && Enu < 1) return 0.1 / 0.25;
    if(Enu >= 1 && Enu < 2) return 0.1 / 0.1;
    if(Enu >= 2 && Enu < 3) return 0.1 / 0.25;
    if(Enu >= 3 && Enu < 4) return 0.1 / 0.5;
    if(Enu >= 4 && Enu < 5) return 0.1 / 1;
    else return 1.;
  });


// Put a "NOvA Preliminary" tag in the corner
void Preliminary()
{
  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
  prelim->SetTextColor(kBlue);
  prelim->SetNDC();
  prelim->SetTextSize(2/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}
/*
void CenterTitles(TH1* histo)
{
  histo->GetXaxis()->CenterTitle();
  histo->GetYaxis()->CenterTitle();
  histo->GetZaxis()->CenterTitle();
  histo->GetYaxis()->SetDecimals();
  histo->GetXaxis()->SetNoExponent();
}
*/
void CenterTitles()
{
  for(const auto& obj:*(gPad->GetListOfPrimitives()))
  {

    if(obj->InheritsFrom("TH1")) CenterTitles((TH1*)obj);
  }
}



void RedrawAxes()
{
  // --- Loop through the objects in the pad. 
  for(const auto& obj:*(gPad->GetListOfPrimitives())) {
    // --- Want to repaint all TH1's.
    if(obj->InheritsFrom("TH1")) {
      // --- Check if I have the Reco Energy canvas?      
      std::string XTit = ((TH1*)obj)->GetXaxis()->GetTitle();
      if ( XTit.find("Reconstructed Neutrino Energ") != std::string::npos ) {
        ((TH1*)obj)->GetYaxis()->SetTitle("Events / 0.1 GeV");
      }
      // Repaint the axes.
      ((TH1*)obj)->GetXaxis()->Paint();
      //((TH1*)obj)->GetYaxis()->SetTitleOffset(1.1);
      ((TH1*)obj)->GetYaxis()->Paint();
      ((TH1*)obj)->GetZaxis()->Paint();
    }
  }
}


TCanvas* LogClone(const TCanvas* can)
{

  TCanvas* clone = (TCanvas*)can->DrawClone();
  clone->cd();

  clone->SetCanvasSize(can->GetWw(), can->GetWh());

  float min =  std::numeric_limits<int>::max();
  float max =  0;

  for(const auto& obj:*(clone->GetListOfPrimitives()))
  {
    if(obj->InheritsFrom("TH1"))
    {
      TH1* h = (TH1*)obj;
      for(int i = 1; i <= h->GetNbinsX(); ++i)
      {
        int c = h->GetBinContent(i);
        if(c != 0 && c < min) min = c;
        if (c > max)          max = c;
      }
    }
  }
  min *= 0.8;
  float span = log10(max) - log10(min);
  span *= 1.5;
  max = pow(10, span + log10(min));
  for(const auto& obj:*(clone->GetListOfPrimitives()))
  {
    if(obj->InheritsFrom("TH1"))
    {
      ((TH1*)obj)->SetMinimum(min);
      ((TH1*)obj)->SetMaximum(max);
    }
  }



  TString name = TString(can->GetTitle()) + "_logy";
  clone->SetName(name);
  clone->SetTitle(name);
  clone->SetLogy();
  return clone;
}

template<class T> class Tangible
{

public:
  Tangible(const T& obj, const std::string&  shortName,
                                  const std::string&  blurb ):
  fObj(obj),
  fShortName(shortName),
  fBlurb(blurb)
  {};

  T           fObj;
  std::string fBlurb;
  std::string fShortName;

};


typedef Tangible<Cut> Selection;
typedef Tangible<HistAxis> TangibleAxis;



class DataMCPair
{
public:

  DataMCPair(Selection sel, TangibleAxis tanAxis,
             SpectrumLoader& loaderData, SpectrumLoader& loaderMC,
             std::vector<const ISyst*> systs, const Var BinWeight, const Cut& bkg=kNoCut ):
  fSel(sel),
  fAxis(tanAxis),
  //fData(loaderData, fAxis.fObj, fSel.fObj, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight),
  //fMC(loaderMC, fAxis.fObj, fSel.fObj, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight),
  //fMCBkg(loaderMC, fAxis.fObj, fSel.fObj && bkg, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight),
  fData (loaderData, fAxis.fObj, fSel.fObj       , kNoShift, BinWeight),
  fMC   (loaderMC  , fAxis.fObj, fSel.fObj       , kNoShift, BinWeight),
  fMCBkg(loaderMC  , fAxis.fObj, fSel.fObj && bkg, kNoShift, BinWeight),
  fShortName(fAxis.fShortName + "_" + fSel.fShortName)
  {
    fUps.reserve(systs.size());
    fDowns.reserve(systs.size());
    for(const auto& syst:systs){
      //fUps.emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, +1), kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight);
      //fDowns.emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, -1), kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight);
      fUps  .emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, +1), BinWeight);
      fDowns.emplace_back(loaderMC, fAxis.fObj, fSel.fObj, SystShifts(syst, -1), BinWeight);
    }
  };

  std::string ShortName() const
  {return fShortName;};
  const char* CName() const
  {return ShortName().c_str();};




  void DrawMCSyst(const int iSyst = -1) const
  {
    if(iSyst < 0)
      PlotWithSystErrorBand(fMC, fUps, fDowns, fData.POT(),
                kTotalMCColor, kTotalMCErrorBandColor, 1.5);
    else
    {
      assert(iSyst < fUps.size());
      PlotWithSystErrorBand(fMC, {fUps[iSyst]}, {fDowns[iSyst]}, fData.POT(),
                kTotalMCColor, kTotalMCErrorBandColor, 1.5);
    }

  }
  void OverlayDataMCSyst(const int iSyst = -1) const
  {
    DrawData();
    DrawMCSyst(iSyst);
    DrawData();
    DrawMCBkg();
    TLegend* leg = DrawLegendPOT();
    AestheticsPOT();

  };

  void DrawMCNormSyst(const int iSyst = -1) const
  {
        // New vectors for area normalized systs
    std::vector<Spectrum> normUps;
    std::vector<Spectrum> normDowns;

    Spectrum tempMC(fMC.ToTH1(fMC.POT()),fMC.POT(),fMC.Livetime());
    tempMC.OverridePOT(fMC.POT() * fMC.ToTH1(fData.POT())->GetSumOfWeights() / fData.ToTH1(fData.POT())->GetSumOfWeights());

    if (iSyst >= 0)
    {
      // We'll have one entry if there's just one syst to plot.
      normUps.reserve(1);
      normDowns.reserve(1);
    }
    else
    {
      // Otherwise we need all systs.
      normUps.reserve(fUps.size());
      normDowns.reserve(fUps.size());
    }
    for(size_t i = 0; i < fUps.size(); ++i)
    {
      // bail if we want a particular syst and it's not the one we want.
      if(iSyst >= 0 && (int)i != iSyst) continue;
      normUps.emplace_back(fUps[i]);
      Spectrum& norm = normUps.back();
      norm.OverridePOT(norm.POT() * norm.Integral(1) / tempMC.Integral(1));
    }
    for(size_t i = 0; i < fDowns.size(); ++i)
    {
      // bail if we want a particular syst and it's not the one we want.
      if(iSyst >= 0 && (int)i != iSyst) continue;
      normDowns.emplace_back(fDowns[i]);
      Spectrum& norm = normDowns.back();
      norm.OverridePOT(norm.POT() * norm.Integral(1) / tempMC.Integral(1));
    }
    PlotWithSystErrorBand(tempMC, normUps, normDowns, fData.POT(),
              kTotalMCColor, kTotalMCErrorBandColor, 1.5);

  }

  //// Weird thing, normalize the systs to area of the MC
  void OverlayDataMCSystNorm(const int iSyst = -1) const
  {

    const int tempcolor=kBlack;
    const int tempmarker=kFullCircle;
    TH1D* hDummy = fData.ToTH1(fData.POT());
    hDummy->GetYaxis()->SetTitle("Events / 0.1 GeV"); // only for cust bin neutrino plot
    //hDummy->GetYaxis()->SetTitleOffset(1.1);
    hDummy->SetMarkerColor(tempcolor);
    hDummy->SetMarkerStyle(tempmarker);
    hDummy->SetLineColor(tempcolor);
    hDummy->Draw();

    DrawMCNormSyst(iSyst);
    DrawData();
    DrawMCBkg();
    TLegend* leg = DrawLegendArea();
    AestheticsArea();
  };


  //// Weird thing, normalize the systs to area of the MC
  void OverlayDataMCSystExtraData(const DataMCPair& extraData,
                                      const std::string dataLegendTitle,
                                      const std::string extraLegendTitle,
                                      const int iSyst = -1,
                                      const int dataColor=kBlack,
                                      const int extraColor=kGray+1,
                                      const int dataMarker=kFullCircle,
                                      const int extraMarker=kFullCircle) const
  {

    DrawMCSyst(iSyst);
    DrawData(dataColor, dataMarker);
    extraData.DrawData(extraColor, extraMarker);
    DrawMCBkg();
    TLegend* leg = DrawLegendPOT(dataColor, dataLegendTitle, dataMarker);
    TH1F* colExtraData = new TH1F();
    colExtraData->SetMarkerColor(extraColor);
    colExtraData->SetMarkerStyle(extraMarker);
    colExtraData->SetLineColor(extraColor);
    leg->AddEntry((TObject*)colExtraData, extraLegendTitle.c_str(), "ple");
    AestheticsPOT();

  };



  //// Weird thing, normalize the systs to area of the MC
  void OverlayDataMCSystNormExtraData(const DataMCPair& extraData,
                                      const std::string dataLegendTitle,
                                      const std::string extraLegendTitle,
                                      const int iSyst = -1,
                                      const int dataColor=kBlack,
                                      const int extraColor=kGray+1,
                                      const int dataMarker=kFullCircle,
                                      const int extraMarker=kFullCircle) const
  {


    DrawMCNormSyst(iSyst);
    DrawData(dataColor, dataMarker);
    extraData.DrawData(extraColor);
    DrawMCBkg();
    TLegend* leg = DrawLegendArea(dataColor, dataLegendTitle, dataMarker);
    TH1F* colExtraData = new TH1F();
    colExtraData->SetMarkerColor(extraColor);
    colExtraData->SetMarkerStyle(extraMarker);
    colExtraData->SetLineColor(extraColor);
    leg->AddEntry((TObject*)colExtraData, extraLegendTitle.c_str(), "ple");
    AestheticsArea();

  };


  TLegend* DrawLegendArea(const int dataColor=kBlack,
                      std::string dataLegendTitle="Data",
                      const int dataMarker=kFullCircle) const
  {

    TLegend* leg = AutoPlaceLegend(0.36, 0.19, 0.68);

    TH1F* colMC = new TH1F();
    colMC->SetLineColor(kTotalMCColor);
    TH1F* colMCBkg = new TH1F();
    colMCBkg->SetLineColor(kBlue);
    TH1F* colData = new TH1F();
    colData->SetLineColor(dataColor);
    colData->SetMarkerColor(dataColor);
    colData->SetMarkerStyle(dataMarker);

    leg->AddEntry((TObject*)colMC, "Simulated Selected Events", "l");
    leg->AddEntry((TObject*)colMCBkg, "Simulated Background", "l");
    leg->AddEntry((TObject*)colData, dataLegendTitle.c_str(), "ple");
    TLegendEntry *entryns1=leg->AddEntry("error","Shape-only 1-#sigma syst. range","f");
    entryns1->SetFillStyle(1001);
    entryns1->SetFillColor(kRed-10);
    leg->SetTextSize(0.04);

    leg->SetFillColor(0);
    leg->SetFillStyle(0);

    for(const auto& obj:*leg->GetListOfPrimitives())
    {
      if(obj->InheritsFrom("TAttFill")){
        ((TAttFill*)obj)->SetFillStyle(0);
        //((TAttFill*)obj)->SetFillColor(0);
      }

    }
    leg->Draw();
    return leg;
  }

  TLegend* DrawLegendPOT(const int dataColor=kBlack,
                      std::string dataLegendTitle="Data",
                      const int dataMarker=kFullCircle) const
  {

    TLegend* leg = AutoPlaceLegend(0.36, 0.19, 0.68);

    TH1F* colMC = new TH1F();
    colMC->SetLineColor(kTotalMCColor);
    TH1F* colMCBkg = new TH1F();
    colMCBkg->SetLineColor(kBlue);
    TH1F* colData = new TH1F();
    colData->SetLineColor(dataColor);
    colData->SetMarkerColor(dataColor);
    colData->SetMarkerStyle(dataMarker);

    leg->AddEntry((TObject*)colMC, "Simulated Selected Events", "l");
    leg->AddEntry((TObject*)colMCBkg, "Simulated Background", "l");
    leg->AddEntry((TObject*)colData, dataLegendTitle.c_str(), "ple");
    TLegendEntry *entryns1=leg->AddEntry("error","Full 1-#sigma syst. range","f");
    entryns1->SetFillStyle(1001);
    entryns1->SetFillColor(kRed-10);
    leg->SetTextSize(0.04);

    leg->SetFillColor(0);
    leg->SetFillStyle(0);

    for(const auto& obj:*leg->GetListOfPrimitives())
    {
      if(obj->InheritsFrom("TAttFill")){
        ((TAttFill*)obj)->SetFillStyle(0);
        //((TAttFill*)obj)->SetFillColor(0);
      }

    }
    leg->Draw();
    return leg;
  }

  void AddExposureArea() const
  {
    std::stringstream expo;
    std::stringstream expoData;
    std::stringstream expoMC;
    float pot = fData.POT();
    float dataMean = (fData.ToTH1(pot)) -> GetMean();
    float dataMeanError = (fData.ToTH1(pot)) -> GetMeanError();
    float MCMean = (fMC.ToTH1(pot)) -> GetMean();
    float MCMeanError = (fMC.ToTH1(pot)) -> GetMeanError();
    expo.precision(2);
    expoData.precision(2);
    expoMC.precision(2);
    expo  << "ND area norm., "
          << TString::Format("%.02lf #times 10^{20}", pot/1e20) << " POT " <<endl;
    //expoData  << TString::Format("Data mean: %.02lf #pm %.02lf GeV  ", dataMean, dataMeanError) <<endl;
    //expoMC    << TString::Format("MC mean:   %.02lf #pm %.02lf GeV  ", MCMean, MCMeanError) <<endl;
    expoData  << TString::Format("Data mean: %.02lf GeV  ", dataMean) <<endl;
    expoMC    << TString::Format("MC mean:   %.02lf GeV  ", MCMean  ) <<endl;

    TLatex text;
    TLegend* place = AutoPlaceLegend(0.0001, 0.0001, 0.52);
    text.SetTextAlign(22);
    text.SetTextSize(0.04);

    // text.DrawLatexNDC(place->GetX1(), place->GetY1(), 
    //                (expo.str()).c_str() );

    text.DrawLatexNDC(place->GetX1(), place->GetY1(), 
                      ("#splitline{" + expo.str() + "}{ #splitline{" + 
                       expoData.str() + "}{" + expoMC.str() + "}}").c_str()
                      );

  }

  void AddExposurePOT() const
  {
    std::stringstream expo;
    float pot = fData.POT();
    expo.precision(2);
    expo  << "ND POT norm., "
          << TString::Format("%.02lf #times 10^{20}", pot/1e20) << " POT";
    TLatex text;
    TLegend* place = AutoPlaceLegend(0.0001, 0.0001, 0.56);
    text.SetTextAlign(22);
    text.SetTextSize(0.04);

    text.DrawLatexNDC(place->GetX1(), place->GetY1(), expo.str().c_str());

    if (ShortName().find("numuE") != std::string::npos) {
      std::cout << "\n" << ShortName() << "\t Data: " << fData.ToTH1(pot)->Integral()
                << "\t MC: " << fMC.ToTH1(pot)->Integral() << "\n";
    }
  }

  void AestheticsArea() const
  {
    CenterTitles();
    RedrawAxes();
    AddExposureArea();
  }

  void AestheticsPOT() const
  {
    CenterTitles();
    RedrawAxes();
    AddExposurePOT();
  }

  void DrawData(const int color=kBlack, const int marker=kFullCircle)const
  {
    TH1D* hData = fData.ToTH1(fData.POT());
    hData->GetYaxis()->SetTitle("Events / 0.1 GeV"); // only for cust bin neutrino plot
    hData->SetMarkerColor(color);
    hData->SetMarkerStyle(marker);
    hData->SetLineColor(color);

    hData->Draw("SAME");
  };

  void DrawMCBkg()const
  {
    TH1D* hMCBkg = fMCBkg.ToTH1(fData.POT());
    hMCBkg->GetYaxis()->SetTitle("Events / 0.1 GeV"); // only for cust bin neutrino plot
    hMCBkg->SetLineColor(kBlue);
    hMCBkg->Draw("HIST SAME");
  };

  float Purity()const
  {
    return 1 - fMCBkg.ToTH1(fData.POT())->GetEntries() /
                  fMC.ToTH1(fData.POT())->GetEntries();
        }

  Selection fSel;
  TangibleAxis fAxis;
  Spectrum fData;
  Spectrum fMC;
  Spectrum fMCBkg;
  std::vector<Spectrum> fUps;
  std::vector<Spectrum> fDowns;
  std::string fShortName;


};


//----------------------------------------------------------------------
  /// Get a vector of all the numu group systs
/*std::vector<const ISyst*> getAllNumuSystsLimited()
{
  std::vector<const ISyst*> systs =
    {&kNumuNormSyst,
     &kNumuRelNormSyst,
     &kNumuNCScaleSyst,
     &kNumuTauContaminationSyst,
     //&kNumuGEANTScaleSyst,
     //&kNumuGEANTNormSyst,
     &kFitHadEnergyScaleSyst,
     &kFDFitHadEnergyScaleSyst,
     &kMECScaleSyst,
     &kRPACCQESyst,
     &kFitMuEnergyScaleSyst,
     &kFDFitMuEnergyScaleSyst
    };
*/


  //if(includeBeam) systs.push_back(&kBeamAll);

  /*
    systs.push_back(&kNumuSummedSmallGENIESyst);
    std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
    largeGENIEParameters.insert(rwgt::fReweightMaNCEL);
    //    largeGENIEParameters.insert(rwgt::fReweightMaCCQE);
    largeGENIEParameters.insert(rwgt::fReweightMaCCQEshape);
    largeGENIEParameters.insert(rwgt::fReweightNormCCQE);
    largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
    largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
    largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
    largeGENIEParameters.insert(rwgt::fReweightMvNCRES);
    for (const auto & parameter:largeGENIEParameters)
    systs.push_back(new GenieRwgtTableSyst(parameter));
  */

//  return systs;
//}


// std::vector<const ISyst*> getSyst_AhtBY()
// {
//   std::vector<const ISyst*> systs =
//     {&kDIS_AhtBY_Syst};
//   return systs;
// }
// std::vector<const ISyst*> getSyst_BhtBY()
// {
//   std::vector<const ISyst*> systs =
//     {&kDIS_BhtBY_Syst};
//   return systs;
// }
// std::vector<const ISyst*> getSyst_CV1uBY()
// {
//   std::vector<const ISyst*> systs =
//     {&kDIS_CV1uBY_Syst};
//   return systs;
// }
// std::vector<const ISyst*> getSyst_CV2uBY()
// {
//   std::vector<const ISyst*> systs =
//     {&kDIS_CV2uBY_Syst};
//   return systs;
// }


/// Return an iterator if all keys are found in only one entry
//std::vector<DataMCPair>::iterator findPair(std::vector<DataMCPair>& list,
  DataMCPair* findPair(std::vector<DataMCPair>& list,
                                           std::vector<std::string> keys)
{
  int nFound = 0;

   std::vector<DataMCPair>::iterator it = list.end();

  for(std::vector<DataMCPair>::iterator entry =  list.begin();
                                       entry != list.end();
                                       ++entry)
  {
    /// Assume we're going to find it, innocent until guilty
    bool found = true;
    for(const auto& key:keys)
    {
      std::cout << entry->ShortName() << "  " << key << "  " << found << std::endl;
      std::cout << (entry->ShortName().find(key) == std::string::npos) << "  " << entry->ShortName().find(key) << "  " << std::string::npos << std::endl;
      // If any key is not present, this isn't the one.
      if(entry->ShortName().find(key) == std::string::npos)
      {

        found = false;
        std::cout << entry->ShortName() << "  " << key << "  " << found  << "  " <<entry->ShortName().find(key) << std::endl;
        break;

      }

    }
    if(found)
    {
     it = entry;
     nFound += 1;
    }

  }
  std::cout << "nFound: "<< nFound << std::endl;
  if(nFound == 0) throw "Failed to find match.";
  if(nFound > 1)  throw "Found too many matches.";

  return &(*it);
}


void WriteBlurb(const DataMCPair& pair, std::string dir, std::string name,
                std::string extra = "", std::string beginning = "")
{
  std::ofstream out(dir + name);
  out << pair.fAxis.fBlurb << pair.fSel.fBlurb << extra;
}

////////////////////////////////////////
// Making a text file for blessing:
void MakeTextFile( std::string OutName ) {
  // Determine output name, and caption.
  std::string FNa = "nd_syst_plots/"+OutName+".txt";
  // Replace underscores in the canvas name
  while (OutName.find("_") != std::string::npos) OutName.replace(OutName.find("_"),1," ");
  // Add a string about it being a cut level.
  if (OutName.find("Qual") != std::string::npos) OutName.insert(OutName.find(" Qual"), ", Cut level:");
  // Add a string about it being a cut level.
  if (OutName.find("pot") != std::string::npos) OutName.replace(OutName.find(" pot"),4,". With POT normalisation.");
  // Add a string about it being a cut level.
  if (OutName.find("area") != std::string::npos) OutName.replace(OutName.find(" area"),5,". With Area normalisation.");

  std::string Cap = "Plot showing the number of ND NuMi events in #nu_{#mu} Ana17 for variable: " + OutName;
  std::cout << "\nFile name: " << FNa << "\n\tCaption: " << Cap << std::endl;
  // Write to file.
  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
  TxtOut << Cap;
  TxtOut.close();
  // Done.
  return;
}

////////////////////////////////////////
// Main class function:
void hadEFrac_nd_data_mc_systs( bool RunBinNorm = false, bool IsTest=false )
{

  ////////////////////////////////////////
  // Get stuff for the hadEFrac Cut:
  ////////////////////////////////////////
  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used 

  std::string fdspecfile = "/nova/app/users/karlwarb/Workspace/NuMuSystematics/FDQuantileHists/final_full_FD_histo_for_quantile_cuts.root";
  TFile* inFile = TFile::Open( fdspecfile.c_str() );
  gDirectory->cd("dir_FDSpec2D");
  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);

  ////////////////////////////////////////
  // Define Axes and cuts.
  ////////////////////////////////////////
  
  std::vector<Cut> TieredCuts;
  std::vector<std::string> CutNames;
  //CutNames.emplace_back("Qual");               TieredCuts.emplace_back(kNumuQuality);
  CutNames.emplace_back("Qual_Cont");          TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017);
  CutNames.emplace_back("Qual_Cont_PID");      TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2017);
  //CutNames.emplace_back("Qual_Cont_PID_True"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2017 && kIsNumuCC);

  const HistAxis EnergyAxis  ( "Reconstructed Muon Energy (GeV)", Binning::Simple(50, 0, 5   ), kMuE        );
  const HistAxis LengthAxis  ( "Length of Primary Track (m)"    , Binning::Simple(50, 0, 16  ), kTrkLength  );
  const HistAxis HadEAxis    ( "Hadronic Energy (GeV)"          , Binning::Simple(50, 0, 3   ), kHadE       );
  const HistAxis HadEFracAxis( "Hadronic Energy Fraction"       , Binning::Simple(50, 0, 1.01), kHadEFrac   );
  const HistAxis CosNuMiAxis ( "Kalman Track Cos #theta_{NuMI}" , Binning::Simple(50, 0, 1.01), kCosNumi    );
  const HistAxis ReMIDAxis   ( "ReMID score"                    , Binning::Simple(50, 0, 1.01), kRemID      );
  const HistAxis CVNNuMuAxis ( "CVN NuMu score"                 , Binning::Simple(50, 0, 1.01), kCVNm       );

  ////////////////////////////////////////
  // Now make the plots....
  ////////////////////////////////////////
  gStyle->SetMarkerStyle(kFullCircle);
  TGaxis::SetMaxDigits(3);

  // --- Configure the definition to use.
  std::string fnameNDMC   = "prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
  std::string fnameNDData = "prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";
  if (!IsTest) {
    //fnameNDMC  = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017"; // Takes a VERY long time to run in one macro 
    //fnameNDMC  = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017_stride10";
    fnameNDMC  = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017_stride22";
    //fnameNDMC  = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2017_lim1";
    fnameNDData= "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_numu2017";
  }
  // --- Set the Loaders, along with spill cuts.
  SpectrumLoader loaderMC(fnameNDMC);
  SpectrumLoader loaderData(fnameNDData);
  loaderMC.SetSpillCut(kStandardSpillCuts);
  loaderData.SetSpillCut(kStandardSpillCuts);

  // --- Load the systematics
  std::vector<const ISyst*> systs = getAllDirectNumuSysts2017();
  //std::vector<const ISyst*> systs = getSyst_AhtBY();
  //std::vector<const ISyst*> systs = getSyst_BhtBY();
  //std::vector<const ISyst*> systs = getSyst_CV1uBY();
  //std::vector<const ISyst*> systs = getSyst_CV2uBY();
  //std::vector<const ISyst*> systs = getAllXsecSysts_2017();
  //std::vector<const ISyst*> systs = {&kMAQAGenieReducedSyst2017};

  // --- Make a vector of "Selections" which are a combination of Quantile && Tiered cuts.
  std::vector<Selection> selections;
  for (size_t cuts=0; cuts < TieredCuts.size(); ++cuts) {
    for (size_t quant=0; quant < HadEFracQuantCuts.size(); ++quant) {
      // --- Want no quantile cut selections too.
      if (quant == 0) {
        selections.emplace_back( TieredCuts[cuts], CutNames[cuts]+"_QuantAll",
                                 "Selected events pass "+CutNames[cuts]+" all quantiles");
      }
      selections.emplace_back( TieredCuts[cuts] && HadEFracQuantCuts[quant], CutNames[cuts]+"_Quant"+std::to_string(quant+1),
                               "Selected events pass "+CutNames[cuts]+" quantile "+std::to_string(quant+1)+".");
    }
  }

  // --- Make a vector of "TangibleAxis" which have the axis and variables for which I want to make spectra.
  std::vector<TangibleAxis> variables;

  if ( RunBinNorm ) {
    variables.emplace_back( kNumuCCOptimisedAxis, "NuMuE"     , "CC energy estimator.");
  } else {
    variables.emplace_back( EnergyAxis          , "MuonE"     , "Reconstructed energy of primary muon track.");
    variables.emplace_back( HadEAxis            , "HadE"      , "Hadronic enery, i.e. numu energy estimate minus muon track energy.  " );
    variables.emplace_back( HadEFracAxis        , "HadEFrac"  , "Hadronic enery fraction, i.e. ration of hadronic energy to numu energy.  " );
    variables.emplace_back( ReMIDAxis           , "ReMIDScore", "ReMId kNN score." );
    variables.emplace_back( CVNNuMuAxis         , "CVNNuMu"   , "CVN muon identification score." );
  }

  //variables.emplace_back( LengthAxis          , "TrkLength" , "Primary track length." );
  //variables.emplace_back( CosNuMiAxis         , "CosNuMi"   , "Beam direction of muon track.  " );


  // // // copied from sa script
  // // variables.emplace_back(HistAxis("Slice N_{Hit}",  Binning::Simple(50, 0, 500), kNHit),
  // //                          "slcNHit", "Number of hits in slice.  " );

  // // variables.emplace_back(HistAxis("Average Hadronic Energy Per Hit (GeV)", Binning::Simple(40, 0, 0.04), kHadEPerNHit),
  // //                          "hadEPerNHit", "Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit.  " );

  // // variables.emplace_back(HistAxis("Track Energy Per Hit (GeV)",  Binning::Simple(40, 0, 0.04), kTrkEPerNHit),
  // //                          "trkEPerNHit", "Average energy per hit on primary track, i.e. trk_E/trk_n_hit.  " );

  // // variables.emplace_back(HistAxis("Hadronic N_{Hit}",  Binning::Simple(50, 0, 100), kHadNHit),
  // //                          "hadNHit", "Number of hits in hadronic cluster.  ");


  // variables.emplace_back(HistAxis("Slice Maximum Y (m)",  kXYBins, kSlcMaxY),
  //                     "maxy", "Maximum Y position of slice.  " );
  // variables.emplace_back(HistAxis("Number of Tracks in Slice",  Binning::Simple(14, 1, 15), kNKalman),
  //                     "nkal", "Number of tracks in slice.  " );
  // // variables.emplace_back(HistAxis("Number of Hits in Slice",  Binning::Simple(50, 0, 500), kNHit),
  // //                          "nhit", "Number of hits in slice.  " );
  // variables.emplace_back(HistAxis("Muon Track cos(#theta_{Z})",  Binning::Simple(50, 0,1), kDirZ),
  //                     "dirZ", "Z-direction of muon track.  " );
  // variables.emplace_back(HistAxis("Track Start X Position (m)",  kXYBins, kTrkStartX),
  //                     "trkStartX", "Track start x position.  " );
  // variables.emplace_back(HistAxis("Track Start Y Position (m)",  kXYBins, kTrkStartY),
  //                     "trkStartY", "Track start y position.  " );
  // variables.emplace_back(
  //   HistAxis("Track Start Z Position (m)",  kZBins, kTrkStartZ),
  //     "trkStartZ", "Track start z position.  " );
  // variables.emplace_back(
  //   HistAxis("Track End X Position (m)",  kXYBins, kTrkEndX),
  //     "trkEndX", "Track stop x position.  " );
  // variables.emplace_back(
  //   HistAxis("Track End Y Position (m)",  kXYBins, kTrkEndY),
  //     "trkEndY", "Track stop y position.  " );
  // variables.emplace_back(
  //   HistAxis("Track End Z Position (m)",  kZBins, kTrkEndZ),
  //     "trkEndZ", "Track stop z position.  " );
  // // variables.emplace_back(
  // //   HistAxis("Slice Duration [ns]",  Binning::Simple(50,0,3000), kSliceDuration),
  // //     "sliceDuration", "Slice duration.  " );
  // // variables.emplace_back(
  // //   HistAxis("Number of Hits in Primary Track",  Binning::Simple(50,0,500), kTrkNhits),
  // //     "trkNhits", "Number of hits on primary track.  ");
  // // variables.emplace_back(
  // //   HistAxis("Length of Primary Track (m)",  Binning::Simple(50,0,16), kTrkLength),
  // //     "trkLength", "Primary track length.  " );

  // variables.emplace_back(
  //   HistAxis("Scattering Log-likelihood",  Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
  //     "scatLL", "ReMId scattering log log-likelihood for primary track.  " );
  // variables.emplace_back(
  //   HistAxis("dE/dx Log-likelihood",  Binning::Simple(50,-3,1), kReMIdDEDxLLH),
  //     "dedxLL", "ReMId dE/dx log log-likelihood for primary track.  " );
  // variables.emplace_back(
  //   HistAxis("Non-hadronic Plane Fraction",
  //                       Binning::Simple(50,0,1), kReMIdMeasFrac),
  //     "nonHadPlaneFrac", "ReMId Non-hadronic plane fraction.  " );

  // // finished copy from sa script

  // --- What is my weight going to be for the spectra?
  Var MyWeight             = kPPFXFluxCVWgt * kXSecCVWgt2017;
  if (RunBinNorm) MyWeight = kPPFXFluxCVWgt * kXSecCVWgt2017 * kBinWidthWeight;

  std::vector<DataMCPair> pairs;
    
  pairs.reserve(selections.size() * variables.size());
  for(const auto& sel:selections){
    for(const auto& variable:variables){
      pairs.emplace_back(sel, variable, loaderData, loaderMC, systs, MyWeight, !kIsNumuCC );
    }
  }
  
  loaderMC.Go();
  loaderData.Go();

  std::vector<TCanvas*> cans;

  for(const auto& pair:pairs) {
     
    std::cout << "\nNow looking pairs, " << pair.ShortName() << " - " << pair.CName() << ", it has purity " << pair.Purity() << std::endl;
      TString name(TString(pair.CName()) + "_allSysts_pot");

      //cans.push_back(new TCanvas(name,name,1200,1000));
      cans.push_back(new TCanvas(name,name));
      pair.OverlayDataMCSyst();
      WriteBlurb(pair, "nd_syst_plots/",
                 std::string(cans.back()->GetName())+".txt",
                 "Full systematic band shown.  ",
                 "Near detector data/MC comparison for numu first analysis.  " );
      /*
      cans.push_back(LogClone(cans.back()));
      WriteBlurb(pair, "nd_syst_plots/",
                 std::string(cans.back()->GetName())+".txt",
                 "Full systematic band shown.  Semi-log scale.",
                 "Near detector data/MC comparison for numu first analysis.  " );
      */

      // Other plot used by Vlad.
      //*
      TString nameNorm(TString(pair.CName()) + "_allSysts_area");
      //cans.push_back(new TCanvas(nameNorm,nameNorm,1200,1000));
      cans.push_back(new TCanvas(nameNorm,nameNorm));
      pair.OverlayDataMCSystNorm();
      WriteBlurb(pair, "nd_syst_plots/",
                 std::string(cans.back()->GetName())+".txt",
                 std::string("Each systematically shifted histogram (both up and down) ")
                 +std::string("is normalized to the area of the MC distribution, ")
                 +std::string("then summed in quadrature." ),
                 "Near detector data/MC comparison for numu first analysis.  " );

      //*/

      /*
      cans.push_back(LogClone(cans.back()));
      WriteBlurb(pair, "nd_syst_plots/",
                 std::string(cans.back()->GetName())+".txt",
                 std::string("Each systematically shifted histogram (both up and down) ")
                 +std::string("is normalized to the area of the MC distribution, ")
                 +std::string("then summed in quadrature.  Semi-log scale." ),
                 "Near detector data/MC comparison for numu first analysis.  ");
      */
    }
  
  // Open a file to save all of the canvases to.
  TFile *OutFile = new TFile("nd_syst_plots/NearDet_NuMiData_Systs_NuMuAna17.root", "RECREATE");
  OutFile->cd();

  for(const auto& can:cans) {
    can->cd();
    Preliminary();
    //can -> SetLeftMargin(0.13);
    can->SaveAs (TString("nd_syst_plots/") + can->GetName() + ".pdf");
    can->Write  ( can->GetName() );
    MakeTextFile( can->GetName() );
  }
}