caf_numu_validation.C

Go to the documentation of this file.

// Fills spectra for systematic error band
#ifdef __CINT__
void caf_numu_validation(std::string kInputFileName, std::string kOutputFileName, bool isFD)
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NumuCuts.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 "CAFAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/NumuVars.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Vars/HistAxes.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 "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;

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

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;

};

const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
                         [](const caf::StandardRecord* sr)
                         {
                           if(sr->mc.nnu == 0) return false;
                           assert(sr->mc.nnu == 1);
                           return (sr->mc.nu[0].mode == caf::kMEC);
                         });

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

const TVector3 beamDirFD = NuMIBeamDirection(caf::kFARDET);
const TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);

const Var kNumuHadTrkE        = SIMPLEVAR(energy.numusimp.hadtrkE);
const Var kQePId              = SIMPLEVAR(sel.qepid.pid);
const Var kNonQeE             = SIMPLEVAR(energy.numusimp.trknonqeE);
const Var kQepidOffE          = SIMPLEVAR(sel.qepid.offE);
const Var kQepidEDiff         = SIMPLEVAR(sel.qepid.ediff);
const Var kQepidEDiffZ        = SIMPLEVAR(sel.qepid.ediffz);
const Var kQepidDedx          = SIMPLEVAR(sel.qepid.dedx);
const Var kQeE                = SIMPLEVAR(energy.numusimp.trkqeE);
const Var kQeAngleE           = SIMPLEVAR(energy.numusimp.angleE);

const Var kQeHadE({"energy.numusimp.trkqeE", "energy.mutrkE", "trk.nkalman",
                     "sel.remid.bestidx"},
            [](const caf::StandardRecord *sr)
            {
              if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
                return sr->energy.numusimp.trkqeE - sr->energy.mutrkE[sr->sel.remid.bestidx].E;
              return -5.f;
            });

const Var kNumuTrackE({"energy.mutrkE", "trk.nkalman", "sel.remid.bestidx"},
            [](const caf::StandardRecord *sr)
            {
              if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
                return sr->energy.mutrkE[sr->sel.remid.bestidx].E;
              return -5.f;
            });

const Var kNumuHadVisE({"energy.numusimp.hadcalE", "energy.numusimp.hadtrkE"}, 
            [](const caf::StandardRecord *sr)
            {
              return sr->energy.numusimp.hadcalE + sr->energy.numusimp.hadtrkE;
            });
const Var kSlcUnCalibNHit({"slc.nhit", "slc.ncalhit"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.nhit - sr->slc.ncalhit; 
            });

  const Var kSlcMeanTime({"slc.meantime"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.meantime / 1000.; // Nicer in \mu s
            });

  const Var kSlcStartTime({"slc.starttime"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.starttime / 1000.; // Nicer in \mu s
            });

  const Var kSlcEndTime({"slc.endtime"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.endtime / 1000.; // Nicer in \mu s
            });

  const Var kSlcMinX({"slc.boxmin.fX"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.boxmin.X() / 100.; // Nicer in m
            });
  
  const Var kSlcMinY({"slc.boxmin.fY"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.boxmin.Y() / 100.; // Nicer in m
            });

  const Var kSlcMinZ({"slc.boxmin.fZ"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.boxmin.Z() / 100.; // Nicer in m
            });

  const Var kSlcMaxX({"slc.boxmax.fX"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.boxmax.X() / 100.; // Nicer in m
            });

  const Var kSlcMaxY({"slc.boxmax.fY"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.boxmax.Y() / 100.; // Nicer in m
            });

  const Var kSlcMaxZ({"slc.boxmax.fZ"},
            [](const caf::StandardRecord *sr)
            {
              return sr->slc.boxmax.Z() / 100.; // Nicer in m
            });

  const Var kSlcExtentX({"slc.boxmax.fX", "slc.boxmin.fX"},
            [](const caf::StandardRecord *sr)
            {
              return (sr->slc.boxmax.X() - sr->slc.boxmin.X())/ 100.; // Nicer in m
            });

  const Var kSlcExtentY({"slc.boxmax.fY", "slc.boxmin.fY"},
            [](const caf::StandardRecord *sr)
            {
              return (sr->slc.boxmax.Y() - sr->slc.boxmin.Y())/ 100.; // Nicer in m
            });

  const Var kSlcExtentZ({"slc.boxmax.fZ", "slc.boxmin.fZ"},
            [](const caf::StandardRecord *sr)
            {
              return (sr->slc.boxmax.Z() - sr->slc.boxmin.Z())/ 100.; // Nicer in m
            });

  const Var kTrkNPlaneGap({"trk.kalman.nplanegap", "trk.nkalman", "sel.remid.bestidx"},
            [](const caf::StandardRecord *sr)
            {
              if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
                return sr->trk.kalman[sr->sel.remid.bestidx].nplanegap;
              return (unsigned short) 500;
            });
  const Var kCosNumi({"hdr.det", "sel.remid.bestidx", "trk.kalman.dir.fX", "trk.kalman.dir.fY",
                     "trk.kalman.dir.fZ"}, 
                    [](const caf::StandardRecord *sr)
                    {
                      if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999){
                        if (sr->hdr.det == 1){
                          return sr->trk.kalman[sr->sel.remid.bestidx].dir.Dot(beamDirND);
                        }
                        if (sr->hdr.det == 2){
                          return sr->trk.kalman[sr->sel.remid.bestidx].dir.Dot(beamDirFD);
                        }
                      }
                      return -5.; 
                    });
  const Var kDirX({"trk.nkalman", "trk.kalman.dir.fX", "sel.remid.bestidx"},
                  [](const caf::StandardRecord *sr)
                  {
                    if(sr->trk.nkalman < 1) return -5.;
                    return sr->trk.kalman[sr->sel.remid.bestidx].dir.X();
                  });

  const Var kDirY({"trk.nkalman", "trk.kalman.dir.fY", "sel.remid.bestidx"},
                  [](const caf::StandardRecord *sr)
                  {
                    if(sr->trk.nkalman < 1) return -5.;
                    return sr->trk.kalman[sr->sel.remid.bestidx].dir.Y();
                  });
  const Var kSlcCalEPerNHit({"slc.nhit", "slc.calE"}, 
            [](const caf::StandardRecord *sr)
            {
              if (sr->slc.nhit > 0) return sr->slc.calE / (1.78 * sr->slc.nhit);
              return -5.;
            });
  const Var kNonQeHadE({"energy.numusimp.trknonqeE", "energy.mutrkE", "trk.nkalman", 
                        "sel.remid.bestidx"},
            [](const caf::StandardRecord *sr)
            {
              if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999) 
        return sr->energy.numusimp.trknonqeE - sr->energy.mutrkE[sr->sel.remid.bestidx].E;
        return -5.f;
            });


class UsefulHist{

public:

  UsefulHist(Selection sel, TangibleAxis tanAxis, SpectrumLoader& loader, const Cut& bkg=kNoCut):
    fSel(sel), 
    fAxis(tanAxis),
    fHist(loader, fAxis.fObj, fSel.fObj),
    fName(tanAxis.fShortName + "-" + sel.fShortName){};

  Selection fSel;
  TangibleAxis fAxis; 
  Spectrum fHist;
  std::string fName;

};

void caf_numu_validation(std::string kInputFileName, std::string kOutputFileName, bool isFD)
{
  std::cout << "\nrun  : ---- Running CAF NuMU Validation.\n";
  std::cout << "run  :     input file name:    " << kInputFileName   << std::endl;
  std::cout << "run  :    output file name:    " << kOutputFileName  << std::endl;
  std::cout << "run  :     FD binning/cuts:    " << isFD             << std::endl;

  gStyle->SetMarkerStyle(kFullCircle);
  TGaxis::SetMaxDigits(3);
  
  SpectrumLoader loader(kInputFileName);

  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
  SpillVar spillRun({"det"}, [](const caf::SRSpill* spill) {return spill->det;});

  loader.SetSpillCut(kStandardSpillCuts);
  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);

  const Binning kXYBins = Binning::Simple(55,-2.20,2.20);
  const Binning kZBins = Binning::Simple(32,  0,16.00);
  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
  const Binning kFDZBins = Binning::Simple(55,0,60.00);
  const Binning kEnergyBinning = Binning::Simple(50,0,5);
  const Binning kHadronicEnergyBinning = Binning::Simple(50,0,5);
  
  std::vector<Selection> selections;

  if (isFD){
    selections.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej, "NumuContainFD",
                            "Selected events pass numu containment, cosmic rejection and slice quality cuts.  ");
    selections.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej && kNumuNCRej, "NumuFD", 
                            "Selected events pass numu containment, cosmic rejection, slice quality, and ReMId cuts.  ");
  }
  else{
    selections.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality, "NumuContainND",
                            "Selected events pass numu containment and slice quality cuts.  ");
    selections.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality && kNumuNCRej, "NumuND", 
                            "Selected events pass numu containment, slice quality, and ReMId cuts.  ");
  }

  std::vector<Selection> selectionsForQePId;

  if (isFD){
    selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej && kNumuNCRej 
                                    && SIMPLEVAR(trk.nkalman)==1, "NumuFD1trk",
                                    "Selected events pass numu containment, cosmic rejection, slice quality, ReMId cuts and have 1 track.  ");
    
    selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej && kNumuNCRej
                                    && SIMPLEVAR(trk.nkalman)==2, "NumuFD2trk",
                                    "Selected events pass numu containment, cosmic rejection, slice quality, ReMId cuts and have 2 tracks.  ");
  }
  else{
    selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality && kNumuNCRej 
                                    && SIMPLEVAR(trk.nkalman)==1, "NumuND1trk",
                                    "Selected events pass numu containment, slice quality, ReMId cuts and have 1 track.  ");
    
    selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality && kNumuNCRej
                                    && SIMPLEVAR(trk.nkalman)==2, "NumuND2trk",
                                    "Selected events pass numu containment, slice quality, ReMId cuts and have 2 tracks.  ");
  }

  std::vector<TangibleAxis> variables;

  variables.emplace_back(
    HistAxis("Reconstructed Neutrino Energy [GeV]", kEnergyBinning, kCCE), 
      "reco-numuE", "CC energy estimator.  ");

  variables.emplace_back(
    HistAxis("Slice N_{Hit}",  Binning::Simple(50, 0, 500), kNHit),
      "reco-slcNHit", "Number of hits in slice.  " );

  variables.emplace_back(
    HistAxis("Slice Calorimetric Energy [GeV]",  kEnergyBinning, kCaloE),
      "reco-calE", "Calorimetric energy of slice.  " );

  variables.emplace_back(
    HistAxis("Reconstructed Hadronic Energy [GeV]",  kHadronicEnergyBinning, kHadE), 
      "reco-hadE", "Hadronic enery, i.e. numu energy estimate minus muon track energy.  " );

  variables.emplace_back(
    HistAxis("Muon Track Energy [GeV]",  Binning::Simple(50, 0, 5), kNumuTrackE),
      "reco-numuTrackE", "Muon track Energy.  " );

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

  variables.emplace_back(
    HistAxis("Average Track Energy Per Hit [GeV]",  Binning::Simple(40, 0, 0.04), kTrkEPerNHit),
      "reco-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),
      "reco-hadNHit", "Number of hits in hadronic cluster.  ");

  variables.emplace_back(
    HistAxis("Off-track Calorimetric Hadronic Energy [GeV]", kHadronicEnergyBinning, kNumuHadCalE),
      "reco-hadCalE", "Sum of calibrated energy deposits in hadronic cluster.  " );

  variables.emplace_back(
    HistAxis("On-track Calorimetric Hadronic Energy [GeV]",  Binning::Simple(50, 0, 0.5), kNumuHadTrkE),
      "reco-hadTrkE", "Sum of calibrated hadronic energy deposits on the muon track.  " );

  variables.emplace_back(
    HistAxis("Calorimetric Hadronic Energy [GeV]",  kHadronicEnergyBinning, kNumuHadVisE),
      "reco-hadVisE", "Sum of calibrated hadronic energy, both on and off the muon track.  " );
  
  if(isFD){
    variables.emplace_back(
                           HistAxis("Slice Maximum Y [m]",  kFDXYBins, kSlcMaxY),
                           "reco-maxy", "Maximum Y position of slice.  " );
  }

  else{
    variables.emplace_back(
                           HistAxis("Slice Maximum Y [m]",  kXYBins, kSlcMaxY),
                           "reco-maxy", "Maximum Y position of slice.  " );
  }

  variables.emplace_back(
    HistAxis("Number of Tracks in Slice",  Binning::Simple(14, 1, 15), kNKalman),
      "reco-nkal", "Number of tracks in slice.  " );
  
  variables.emplace_back(
    HistAxis("Number of Hits in Slice",  Binning::Simple(50, 0, 500), kNHit),
      "reco-nhit", "Number of hits in slice.  " );
  
  if(isFD){
    variables.emplace_back(
                           HistAxis("Track Start X Position [m]",  kFDXYBins, kTrkStartX),
                           "reco-trkStartX", "Track start x position.  " );
    
    variables.emplace_back(
                           HistAxis("Track Start Y Position [m]",  kFDXYBins, kTrkStartY),
                           "reco-trkStartY", "Track start y position.  " );
    
    variables.emplace_back(
                           HistAxis("Track Start Z Position [m]",  kFDZBins, kTrkStartZ),
                           "reco-trkStartZ", "Track start z position.  " );
    
    variables.emplace_back(
                           HistAxis("Track End X Position [m]",  kFDXYBins, kTrkEndX),
                           "reco-trkEndX", "Track stop x position.  " );
    
    variables.emplace_back(
                           HistAxis("Track End Y Position [m]",  kFDXYBins, kTrkEndY),
                           "reco-trkEndY", "Track stop y position.  " );
    
    variables.emplace_back(
                           HistAxis("Track End Z Position [m]",  kFDZBins, kTrkEndZ),
                           "reco-trkEndZ", "Track stop z position.  " );
    
    variables.emplace_back(
                           HistAxis("Slice Duration [ns]",  Binning::Simple(50,0,3000), kSliceDuration),
                           "reco-sliceDuration", "Slice duration.  " );
  }

  else{
    variables.emplace_back(
                           HistAxis("Track Start X Position [m]",  kXYBins, kTrkStartX),
                           "reco-trkStartX", "Track start x position.  " );
    
    variables.emplace_back(
                           HistAxis("Track Start Y Position [m]",  kXYBins, kTrkStartY),
                           "reco-trkStartY", "Track start y position.  " );
    
    variables.emplace_back(
                           HistAxis("Track Start Z Position [m]",  kZBins, kTrkStartZ),
                           "reco-trkStartZ", "Track start z position.  " );
    
    variables.emplace_back(
                           HistAxis("Track End X Position [m]",  kXYBins, kTrkEndX),
                           "reco-trkEndX", "Track stop x position.  " );
    
    variables.emplace_back(
                           HistAxis("Track End Y Position [m]",  kXYBins, kTrkEndY),
                           "reco-trkEndY", "Track stop y position.  " );
    
    variables.emplace_back(
                           HistAxis("Track End Z Position [m]",  kZBins, kTrkEndZ),
                           "reco-trkEndZ", "Track stop z position.  " );
    
    variables.emplace_back(
                           HistAxis("Slice Duration [ns]",  Binning::Simple(50,0,500), kSliceDuration),
                           "reco-sliceDuration", "Slice duration.  " );
  }
  
  variables.emplace_back(
    HistAxis("Number of Hits on Primary Track",  Binning::Simple(50,0,500), kTrkNhits),
      "reco-trkNhits", "Number of hits on primary track.  ");
  
  variables.emplace_back(
    HistAxis("Primary Track Length [m]",  Binning::Simple(50,0,16), kTrkLength),
      "reco-trkLength", "Primary track length.  " );
  
  variables.emplace_back(
    HistAxis("ReMId Input: Scattering Log-likelihood",  Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
      "reco-scatLL", "ReMId scattering log log-likelihood for primary track.  " );
  
  variables.emplace_back(
    HistAxis("ReMId Input: dE/dx Log-likelihood",  Binning::Simple(50,-3,1), kReMIdDEDxLLH),
      "reco-dedxLL", "ReMId dE/dx log log-likelihood for primary track.  " );
  
  variables.emplace_back(
    HistAxis("ReMId Input: Non-hadronic Plane Fraction", Binning::Simple(50,0,1), kReMIdMeasFrac),
      "reco-nonHadPlaneFrac", "ReMId Non-hadronic plane fraction.  " );
  
  variables.emplace_back(
    HistAxis("ReMId",  kRemidBinning, kRemID),
      "reco-remid", "ReMId kNN score.  " );

  variables.emplace_back(
    HistAxis("Number of Uncalibrated Hits in Slice",  Binning::Simple(50, 0, 50), kSlcUnCalibNHit),
      "reco-slcUnCalibNHit", "Uncalibrated slice nhit.  " );

  variables.emplace_back(
    HistAxis("Slice Mean Time [#mus]",  Binning::Simple(60, -50, 550), kSlcMeanTime),
      "reco-slcMeanTime", "Slice mean time TNS in microseconds.  " );

  variables.emplace_back(
    HistAxis("Slice Start Time [#mus]",  Binning::Simple(60, -50, 550), kSlcStartTime),
      "reco-slcStartTime", "Slice start time TNS in microseconds.  " );

  variables.emplace_back(
    HistAxis("Slice End Time [#mus]",  Binning::Simple(60, -50, 550), kSlcEndTime),
      "reco-slcEndTime", "Slice end time TNS in microseconds.  " );

  if(isFD){
    variables.emplace_back(
                           HistAxis("Slice Minimum X [m]",  kFDXYBins, kSlcMinX),
                           "reco-slcMinX", "Slice Min X [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Minimum Y [m]",  kFDXYBins, kSlcMinY),
                           "reco-slcMinY", "Slice Min Y [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Minimum Z [m]",  kFDZBins, kSlcMinZ),
                           "reco-slcMinZ", "Slice Min Z [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Maximum X [m]",  kFDXYBins, kSlcMaxX),
                           "reco-slcMaxX", "Slice Max X [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Maximum Y [m]",  kFDXYBins, kSlcMaxY),
                           "reco-slcMaxY", "Slice Max Y [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Maximum Z [m]",  kFDZBins, kSlcMaxZ),
                           "reco-slcMaxZ", "Slice Max Z [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Extent X [m]",  Binning::Simple(55, 0, 16.0), kSlcExtentX),
                           "reco-slcExtentX", "Slice Extent X [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Extent Y [m]",  Binning::Simple(55, 0, 16.0), kSlcExtentY),
                           "reco-slcExtentY", "Slice Extent Y [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Extent Z [m]",  Binning::Simple(55, 0, 60.0), kSlcExtentZ),
                           "reco-slcExtentZ", "Slice Extent Z [m].  " );
  }
  
  else{
    variables.emplace_back(
                           HistAxis("Slice Minimum X [m]",  kXYBins, kSlcMinX),
                           "reco-slcMinX", "Slice Min X [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Minimum Y [m]",  kXYBins, kSlcMinY),
                           "reco-slcMinY", "Slice Min Y [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Minimum Z [m]",  kZBins, kSlcMinZ),
                           "reco-slcMinZ", "Slice Min Z [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Maximum X [m]",  kXYBins, kSlcMaxX),
                           "reco-slcMaxX", "Slice Max X [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Maximum Y [m]",  kXYBins, kSlcMaxY),
                           "reco-slcMaxY", "Slice Max Y [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Maximum Z [m]",  kZBins, kSlcMaxZ),
                           "reco-slcMaxZ", "Slice Max Z [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Extent X [m]",  Binning::Simple(55, 0, 4.4), kSlcExtentX),
                           "reco-slcExtentX", "Slice Extent X [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Extent Y [m]",  Binning::Simple(55, 0, 4.4), kSlcExtentY),
                           "reco-slcExtentY", "Slice Extent Y [m].  " );
    
    variables.emplace_back(
                           HistAxis("Slice Extent Z [m]",  Binning::Simple(18, 0, 16), kSlcExtentZ),
                           "reco-slcExtentZ", "Slice Extent Z [m].  " );
  }

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{X}",  Binning::Simple(100, -1, 1), kDirX),
      "reco-dirX", "X-direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{Y}",  Binning::Simple(100, -1, 1), kDirY),
      "reco-dirY", "Y-direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{Z}",  Binning::Simple(50, 0, 1), kDirZ),
      "reco-dirZ", "Z-direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Kalman Track Cos #theta_{NuMI}",  Binning::Simple(50, 0, 1), kCosNumi),
      "reco-cosNumi", "Beam direction of muon track.  " );

  variables.emplace_back(
    HistAxis("Number of Missing Planes in Primary Track",  Binning::Simple(50, 0, 50), 
                                                          kTrkNPlaneGap),
      "reco-trkNPlaneGap", "Track N Plane Gap.  " );

  variables.emplace_back(
    HistAxis("Visible Slice Energy Per Hit [GeV]",  Binning::Simple(40, 0, 0.04),
                                                          kSlcCalEPerNHit),
      "reco-slcCalEPerNHit", "Slice Energy Per Slice NHit.  " );

  /////////////////////////////////////////////////////
  std::vector<TangibleAxis> variablesForQePId;

  variablesForQePId.emplace_back(
    HistAxis("QePId",  Binning::Simple(51, -0.01, 1.01), kQePId),
      "reco-qepid", "QePId kNN score.  " );

  variablesForQePId.emplace_back(
    HistAxis("Non QE Neutrino Energy [GeV]",  Binning::Simple(50, 0, 5), kNonQeE),
      "reco-nonQeE", "Non QE Energy.  " );

  variablesForQePId.emplace_back(
    HistAxis("Hadronic Non-QE Energy [GeV]",  Binning::Simple(50, 0, 5), kNonQeHadE),
      "reco-nonQeHadE", "Hadronic Non QE Energy.  " );

  variablesForQePId.emplace_back(
    HistAxis("QePId Input: Off-track Energy Ratio",  Binning::Simple(100, 0, 1), kQepidOffE),
      "reco-QepidOffE", "Off-track Energy Ratio, QePId.  " );

  variablesForQePId.emplace_back(
    HistAxis("QePId Input: Fractional Energy Difference",  Binning::Simple(100, -1, 1), 
                                                           kQepidEDiff),
      "reco-QepidEDiff", "Fractional Energy Difference, QePId.  " );

  variablesForQePId.emplace_back(
    HistAxis("QePId Input: Fractional Energy Difference Z-test",  Binning::Simple(100, -5, 5),
                                                           kQepidEDiffZ),
      "reco-QepidEDiffZ", "Fractional Energy Difference Z-test, QePId.  " );

  variablesForQePId.emplace_back(
    HistAxis("QePId Input: dE/dx Ratio",  Binning::Simple(100, 0, 5), kQepidDedx),
      "reco-QepidDedx", "dE/dx Ratio, applied only to 2 trk events, QePId.  " );

  variablesForQePId.emplace_back(
    HistAxis("QE Neutrino Energy [GeV]",  Binning::Simple(50, 0, 5), kQeE),
      "reco-qeE", "QE energy.  " );

  variablesForQePId.emplace_back(
    HistAxis("Hadronic QE Energy [GeV]",  Binning::Simple(50, 0, 5), kQeHadE),
      "reco-qeHadE", "Hadronic energy for QE things.  " );

  variablesForQePId.emplace_back(
    HistAxis("QE Angle Neutrino Energy [GeV]",  Binning::Simple(50, 0, 5), kQeAngleE),
      "reco-qeAngleE", "QE angle formula energy.  " );

  variablesForQePId.emplace_back(
    HistAxis("Primary Track Energy [GeV]",  Binning::Simple(50, 0, 5), kNumuTrackE),
      "reco-numuTrackE", "Primary track Energy.  " );

  variablesForQePId.emplace_back(
    HistAxis("Visible Hadronic Energy [GeV]",  Binning::Simple(50, 0, 2), kNumuHadVisE),
      "reco-HadVisE", "Hadronic Energy .  " );


  std::vector<UsefulHist> hists;
  hists.reserve(selections.size() * variables.size() + selectionsForQePId.size() 
                                                          * variablesForQePId.size());
  
  for(const auto& sel:selections){
    for(const auto& variable:variables){
      hists.emplace_back(sel, variable, loader);
    }
  }

  for(const auto& sel:selectionsForQePId){
    for(const auto& variable:variablesForQePId){
      hists.emplace_back(sel, variable, loader);
    }
  }

  std::cout << "\nrun  : --- run loader.\n";
  loader.Go();
  std::cout << "\nrun  : --- done.\n";
  double pot = hists[0].fHist.POT();


  std::cout << "\nrun  : --- save output.\n";
  // TFile  inputfile(kInputFileName.c_str(),     "READ");
  TFile outputfile(kOutputFileName.c_str(), "RECREATE");

  for(const auto& hist:hists){
    TH1 *temp = hist.fHist.ToTH1(pot);
    std::string tempName = hist.fName;
    temp->Write(tempName.c_str());
  }

  std::cout << "\nrun  : --- save PoT.\n";

  TH1F *pothist = new TH1F("meta-TotalPOT",    "TotalPOT", 1, 0, 1);
  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1); 
  pothist->SetBinContent(1, pot);
  evthist->SetBinContent(1, spillHist->GetEntries());
  pothist->Write("meta-TotalPOT");
  evthist->Write("meta-TotalEvents");

  std::cout << "\nrun  : --- close output files.\n";
  // inputfile.Close();
  outputfile.Write();
  outputfile.Close();

  std::cout << "\nrun  : --- done.\n";
}

#endif