make_xsec_wgts_2018_hists.C

Go to the documentation of this file.

#ifdef __CINT__
void make_xsec_wgts_2018_hists
(
  const std::string beam,
  const std::string type,
  const std::string period = "full"
  //const std::string shift_non_mec = "",
)
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

//#include "StandardRecord/Proxy/SRProxy.h"

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Utilities.h"
#include "3FlavorAna/Core/HistAxes.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/TruthVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/HistAxes.h"
//#include "CAFAna/Vars/XsecTunes.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/XSecSysts.h"
//#include "MCReweight/func/HistWrapper.h"

#include "TTree.h"
#include "TH2.h"

#include <fstream>

using namespace ana;
using namespace caf;


TFile* file_mec_wgt_numu    = NULL;
TFile* file_mec_wgt_numubar = NULL;
TH2D* hist_mec_wgt_numu     = NULL;
TH2D* hist_mec_wgt_numubar  = NULL;

void LoadTunedMECWeights( const std::string shift_non_mec = "", const bool extra_stat_mc = false, const bool no_smoothing = false, const bool mec_refit = false )
{ 
  // "up" and "down" refer to shifting the total MC in the (QEL,RES) region up and down in (q0,q3)
  if ( !shift_non_mec.empty() && shift_non_mec != "Down" && shift_non_mec != "Up" ) throw std::runtime_error( Form( "Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );

  const bool finer_bins = true;
  std::string fit_tag = "";
  if ( !shift_non_mec.empty()  ) fit_tag += "_ShiftNonMEC" + shift_non_mec;
  if ( finer_bins              ) fit_tag += "_FinerBins";
  if ( extra_stat_mc           ) fit_tag += "_ExtraStatMC";

  std::string base_indir = mec_refit ? "/nova/app/users/mislivec/mec_weights" : "/nova/app/users/mislivec/MEC_rwgt/mec_fit_q0q3/test";
  std::string weights_file_name_numu    = Form( "%s/FHC_MECFit_2018%s/tuned_mec_weights.root", base_indir.c_str(), fit_tag.c_str() );
  std::string weights_file_name_numubar = Form( "%s/RHC_MECFit_2018%s/tuned_mec_weights.root", base_indir.c_str(), fit_tag.c_str() );

  if ( gSystem->AccessPathName( weights_file_name_numu.c_str() ) ) throw std::runtime_error( "Numu MEC weights file not found:\n" + weights_file_name_numu );
  if ( gSystem->AccessPathName( weights_file_name_numubar.c_str() ) ) throw std::runtime_error( "Numu MEC weights file not found:\n" + weights_file_name_numubar );

  file_mec_wgt_numu    = new TFile( weights_file_name_numu.c_str() );
  file_mec_wgt_numubar = new TFile( weights_file_name_numubar.c_str() );
  
  std::cout << "MEC weight files:" << std::endl;
  std::cout << file_mec_wgt_numu->GetName() << std::endl;
  std::cout << file_mec_wgt_numubar->GetName() << std::endl;

  std::string hist_name_numu    = no_smoothing ? "numu_mec_weights"    : "numu_mec_weights_smoothed";
  std::string hist_name_numubar = no_smoothing ? "numubar_mec_weights" : "numubar_mec_weights_smoothed";

  hist_mec_wgt_numu    = (TH2D*)file_mec_wgt_numu   ->Get( hist_name_numu.c_str() );
  hist_mec_wgt_numubar = (TH2D*)file_mec_wgt_numubar->Get( hist_name_numubar.c_str() ); 

  //hist_mec_wgt_numu    = (TH2D*)file_mec_wgt_numu   ->Get( "numu_mec_weights_smoothed" ); 
  //hist_mec_wgt_numubar = (TH2D*)file_mec_wgt_numubar->Get( "numubar_mec_weights_smoothed" );

}

const ana::Var GetTunedMECWeight( const bool interpolate = false )
{
  const Var kWeightTunedMEC( [ interpolate ]( const caf::SRProxy * sr )
  {
    if ( !ana::kIsDytmanMEC( sr ) ) return 1.0;

    double q0 = kTrueQ0( sr );
    double q3 = kTrueQ3( sr );

    double weight = 1.0;

    if ( sr->mc.nu[0].pdg > 0 )
    {
      if ( hist_mec_wgt_numu == NULL ) throw std::runtime_error( "Numu MEC weights histogram is NULL" );

      if ( q0 < hist_mec_wgt_numu->GetYaxis()->GetXmin() ||
           q0 > hist_mec_wgt_numu->GetYaxis()->GetXmax() ||
           q3 < hist_mec_wgt_numu->GetXaxis()->GetXmin() ||
           q3 > hist_mec_wgt_numu->GetXaxis()->GetXmax() ) return 1.0;

      weight = interpolate ?
               hist_mec_wgt_numu->Interpolate( q3, q0 ) :
               hist_mec_wgt_numu->GetBinContent( hist_mec_wgt_numu->FindFixBin( q3, q0 ) );
    }
    else
    {
      if ( hist_mec_wgt_numubar == NULL ) throw std::runtime_error( "Numubar MEC weights histogram is NULL" );
 
      if ( q0 < hist_mec_wgt_numubar->GetYaxis()->GetXmin() ||
           q0 > hist_mec_wgt_numubar->GetYaxis()->GetXmax() ||
           q3 < hist_mec_wgt_numubar->GetXaxis()->GetXmin() ||
           q3 > hist_mec_wgt_numubar->GetXaxis()->GetXmax() ) return 1.0;

      weight = interpolate ?
               hist_mec_wgt_numubar->Interpolate( q3, q0 ) :
               hist_mec_wgt_numubar->GetBinContent( hist_mec_wgt_numubar->FindFixBin( q3, q0 ) );
    } 
    if ( weight < 0 ) weight = 0.0;
    return weight;
  });
  return kWeightTunedMEC;
}

const Var kMinosRESSupp( []( const caf::SRProxy * sr )
{ 
  if ( !kIsRes( sr ) || !kIsNumuCC( sr ) ) return 1.0;
  if ( sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
  double Q2 = kTrueQ2( sr );
  return 1.010 / ( 1 + exp( 1 - sqrt( Q2 ) / 0.156 ) ); // From PhysRevD.91.012005
});

const Var kMINERvALowQ2RESSupp( []( const caf::SRProxy * sr )
{
  if ( !kIsRes( sr ) ) return 1.0;
  if ( sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen

  // arXiv:1903.01558
  // Joint FrAbs Tune
  double R1 = 0.32;
  double R2 = 0.50;

  double x1 = 0.0;
  double x2 = 0.35;
  double x3 = 0.7;

  double Q2 = kTrueQ2( sr );
  if ( Q2 > x3 ) return 1.0;

  double R = R2 * ( Q2 - x1 ) * ( Q2 - x3 ) / ( ( x2 - x1 ) * ( x2 - x3 ) ) + ( Q2 - x1 ) * ( Q2 - x2 ) / ( ( x3 - x1 ) * ( x3 - x2 ) );
  double weight = 1 - ( 1 - R1 ) * pow( 1 - R, 2 );
  if ( weight < 0 )
  {
    std::cout << "kMINERvALowQ2RESSupp = " << weight << std::endl;
    std::cout << "Q2 = " << Q2 << std::endl;
    weight = 0;
  }
  return weight;
});

/*
const Var kHitNucBin( []( const caf::SRProxy * sr )
{
  if ( !kIsDytmanMEC( sr ) ) return -1.0;
  else if ( kHitNuc( sr ) == 2000000200 ) return 0.0;
  else if ( kHitNuc( sr ) == 2000000201 ) return 1.0;
  else if ( kHitNuc( sr ) == 2000000202 ) return 2.0;
  else return 3.0;
});
*/

const Var kEmpiricalMECtoValenciaMECWgt_FixNP( []( const caf::SRProxy * sr )
{
  if ( !kIsDytmanMEC( sr ) ) return 1.0;

  double weight = kEmpiricalMECtoValenciaMECWgt( sr );

  const double np_frac_empirical = 0.89;
  const double np_frac_valencia  = 2.0 / 3.0;

  if ( ana::kHitNuc( sr ) == 2000000201 ) // np
    weight *= np_frac_valencia / np_frac_empirical;
  else
    weight *= ( 1 - np_frac_valencia ) / ( 1 - np_frac_empirical );
  
  return weight;
});

const Cut kCutDISAnomaly( []( const caf::SRProxy* sr )
{
  return sr->mc.nnu == 0;
});

const Cut kCutOnHydrogen( []( const caf::SRProxy* sr )
{ 
  return sr->mc.nu[0].tgtA == 1;
});

/*
const Var kTrueY( []( const caf::SRProxy* sr )
{
  return ( sr->mc.nnu == 0 ) ? -1.0 : float( sr->mc.nu[0].y );
});
*/

const Cut kIsOther = !kIsNumuCC || !( kIsQE || kIsRes || kIsDIS || kIsDytmanMEC );

const Cut kIsNumuCCMEC    = kIsNumuCC && kIsDytmanMEC && !kIsAntiNu;
const Cut kIsNumubarCCMEC = kIsNumuCC && kIsDytmanMEC &&  kIsAntiNu;

const Cut kWSelection( []( const caf::SRProxy* sr )
{
  double W = kRecoW( sr );
  return W > 1.2 && W < 1.5;
});

const Cut kIsMuonNeutrino( []( const caf::SRProxy* sr )
{
  return sr->hdr.ismc && sr->mc.nnu == 1 && sr->mc.nu[0].pdg == 14;
});

const Cut kIsMuonAntiNeutrino( []( const caf::SRProxy* sr )
{
  return sr->hdr.ismc && sr->mc.nnu == 1 && sr->mc.nu[0].pdg == -14;
});

/*
SystShifts kShiftsMECShapeUp;
kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixNu,     +1 );
kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixAntiNu, +1 );

SystShifts kShiftsMECShapeDown;
kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixNu,     -1 );
kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixAntiNu, -1 );
*/

//const Var kWeightTunedMECInterp = GetTunedMECWeight( true );

const NuTruthVar kTrueW_Eqn_NT( []( const caf::SRNeutrinoProxy * nu )
{
  const double M_p = 0.938272;
  double W2 = M_p * M_p + 2 * M_p * ( nu->y * nu->E ) - nu->q2;
  if ( W2 > 0 ) return sqrt( W2 );
  else          return -5.;
});
const Var kTrueW_Eqn = VarFromNuTruthVar( kTrueW_Eqn_NT, -5 );

const ana::Var GetWResolution( const bool fractional = false, const bool true_w_eqn = false )
{
  const Var kWResolution( [ fractional, true_w_eqn ]( const caf::SRProxy * sr )
  {
    const double trueW = true_w_eqn ? kTrueW_Eqn( sr ) : kTrueW( sr );
    const double recoW = kRecoW( sr );
    if ( trueW > 0 && recoW > 0 )
    {
      double wres = recoW - trueW;
      if ( fractional ) wres /= trueW;
      return wres;
    }
    return -10.;
  });
  return kWResolution;
}

class Distribution
{

  public :

    Distribution(
      std::string         name,
      bool                isMC,
      SpectrumLoaderBase& loader,
      const HistAxis&     axis,
      const Cut&          cut,
      const Var&          wgtMC = kUnweighted,
      const SystShifts&   shiftMC = kNoShift )
    {
      fName = name;

      if ( isMC )
      {
        fSpectra[ "TotMC"           ] = std::make_unique<Spectrum>( loader, axis, cut                             , shiftMC, wgtMC );
        fSpectra[ "QEL"             ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE       , shiftMC, wgtMC );
        fSpectra[ "MEC"             ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC );
        fSpectra[ "RES"             ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes      , shiftMC, wgtMC );
        fSpectra[ "DIS"             ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS      , shiftMC, wgtMC );
        fSpectra[ "Other"           ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther                 , shiftMC, wgtMC );

        //fSpectra[ "Numu"           ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonNeutrino           , shiftMC, wgtMC );
        //fSpectra[ "Numubar"        ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonAntiNeutrino       , shiftMC, wgtMC );

        /*
        fSpectra[ "QEL_WSelected"   ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE && kWSelection , shiftMC, wgtMC );
        fSpectra[ "MEC_WSelected"   ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC && kWSelection, shiftMC, wgtMC );
        fSpectra[ "RES_WSelected"   ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes && kWSelection, shiftMC, wgtMC );
        fSpectra[ "DIS_WSelected"   ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS && kWSelection, shiftMC, wgtMC );
        fSpectra[ "Other_WSelected" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther && kWSelection, shiftMC, wgtMC );
        */

        //fSpectra[ "TotMC_MECShapeUp"   ] = std::make_unique<Spectrum>( loader, axis, cut, kShiftsMECShapeUp  , wgtMC );
        //fSpectra[ "TotMC_MECShapeDown" ] = std::make_unique<Spectrum>( loader, axis, cut, kShiftsMECShapeDown, wgtMC );

        fSpectra[ "TotMC_Nominal"   ] = std::make_unique<Spectrum>( loader, axis, cut                             , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "QEL_Nominal"     ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE       , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "MEC_Nominal"     ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "RES_Nominal"     ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes      , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "DIS_Nominal"     ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS      , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "Other_Nominal"   ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther                 , shiftMC, kPPFXFluxCVWgt );

        /*
        fSpectra[ "Numu_Nominal"    ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonNeutrino          , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "Numubar_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonAntiNeutrino      , shiftMC, kPPFXFluxCVWgt );

        fSpectra[ "NonMEC_Nominal"         ] = std::make_unique<Spectrum>( loader, axis, cut && !kIsDytmanMEC                       , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "NonMEC_Nominal_Numu"    ] = std::make_unique<Spectrum>( loader, axis, cut && !kIsDytmanMEC && kIsMuonNeutrino    , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "NonMEC_Nominal_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && !kIsDytmanMEC && kIsMuonAntiNeutrino, shiftMC, kPPFXFluxCVWgt );
        */

        fSpectra[ "MEC_Numu"            ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC         , shiftMC, wgtMC );
        fSpectra[ "MEC_Numubar"         ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC      , shiftMC, wgtMC );
        fSpectra[ "MEC_Nominal_Numu"    ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC         , shiftMC, kPPFXFluxCVWgt );
        fSpectra[ "MEC_Nominal_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC      , shiftMC, kPPFXFluxCVWgt );

        fSpectra[ "MEC_Valencia"         ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt );
        fSpectra[ "MEC_Valencia_Numu"    ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC             , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt );
        fSpectra[ "MEC_Valencia_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC          , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt );

        fSpectra[ "MEC_Valencia_FixNP"         ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt_FixNP );
        fSpectra[ "MEC_Valencia_FixNP_Numu"    ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC             , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt_FixNP );
        fSpectra[ "MEC_Valencia_FixNP_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC          , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt_FixNP );

        fSpectra[ "MEC_MINERvA"          ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt * kMINERvA_Wgt_MEC );
        fSpectra[ "MEC_MINERvA_Numu"     ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC             , shiftMC, kPPFXFluxCVWgt * kMINERvA_Wgt_MEC );
        fSpectra[ "MEC_MINERvA_Numubar"  ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC          , shiftMC, kPPFXFluxCVWgt * kMINERvA_Wgt_MEC );

        fSpectra[ "RES_MINOS"          ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes, shiftMC, kPPFXFluxCVWgt * kMinosRESSupp );
        fSpectra[ "RES_MINERvA"        ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes, shiftMC, kPPFXFluxCVWgt * kMINERvALowQ2RESSupp );
      }
      else
      {
        fSpectra[ "Data"  ] = std::make_unique<Spectrum>( loader, axis, cut );
      }
    }

    void Write() const
    {
      for ( const auto & s : fSpectra )
      {
        double pot = s.second.get()->POT();
        std::string out_name = Form( "%s_%s", fName.c_str(), s.first.c_str() );
        int ndim = s.second.get()->NDimensions();
        if ( ndim == 1 )
        {
          s.second.get()->ToTH1( pot )->Write( out_name.c_str() );
        }
        else if ( ndim == 2 )
        {
          s.second.get()->ToTH2( pot )->Write( out_name.c_str() );
        }
        else throw std::runtime_error( Form( "Write operation for %d dimensions not defined", ndim ) );
        //s.second.get()->ToTH1( s.second.get()->POT() )->Write( Form( "%s_%s", fName.c_str(), s.first.c_str() ) );
      }
    }

    double POT() const { return fSpectra.begin()->second.get()->POT(); } 

  private :

    std::string fName;
    std::map<std::string,std::unique_ptr<Spectrum>> fSpectra;
};

void make_xsec_wgts_2018_hists
( 
  const std::string beam,
  const std::string type,
  const std::string period = "full"
  //const std::string shift_non_mec = "",
)
{
  if ( type != "data" && type != "mc" )
  {
    std::cout << "type must be \"data\" or \"mc\"" << std::endl;
    return;
  }
  const bool isMC = ( type == "mc" );

  /*
  const bool extra_stat_mc = false;
  const bool no_smoothing  = false;
  const bool mec_refit     = false;
  LoadTunedMECWeights( shift_non_mec, extra_stat_mc, no_smoothing, mec_refit );
  */

  std::string defn = "";
  if ( beam == "fhc" )
  {
    defn = isMC ?
           Form( "prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_%s_v1_numu2018", period.c_str() ) :
           Form( "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_%s_v1_addShortSimpleCVN_goodruns_numu2018", period.c_str() );
  }
  else if ( beam == "rhc" )
  {
    defn = isMC ?
           Form( "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_%s_v1_numu2018", period.c_str() ) :
           Form( "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_%s_v1_addShortSimpleCVN_goodruns_numu2018", period.c_str() );
  }
  if ( defn.empty() ) throw std::runtime_error( "SAM Definition is empty" );

  //defn = Form( "defname: %s with limit 1", defn.c_str() );
  std::cout << "defn = " << defn << std::endl;

  SpectrumLoader loader( defn );
  loader.SetSpillCut( kStandardSpillCuts );

  std::map< std::string, HistAxis > axes;
  axes.emplace( "RecoEhad"       , HistAxis( "Reco E_{had} (GeV)",       Binning::Simple( 50, 0., 2. ), kHadE ) );
  axes.emplace( "RecoEhadVis"    , HistAxis( "Visible E_{had} (GeV)",    Binning::Simple( 50, 0., 1. ), kNumuHadVisE ) );
  axes.emplace( "RecoEmu"        , HistAxis( "Reco E_{#mu} (GeV)",       Binning::Simple( 60, 0., 3. ), kMuE ) );
  axes.emplace( "RecoEnu"        , HistAxis( "Reco E_{#nu} (GeV)",       Binning::Simple( 40, 0., 4. ), kCCE ) );
  axes.emplace( "RecoQ2"         , HistAxis( "Reco Q^{2} (GeV^{2})",     Binning::Simple( 50, 0., 1. ), kRecoQ2 ) );
  axes.emplace( "RecoQ3"         , HistAxis( "Reco |#vec{q}| (GeV)",     Binning::Simple( 40, 0., 2. ), kRecoQmag ) );
  axes.emplace( "RecoW"          , HistAxis( "Reco W (GeV)",             Binning::Simple( 60, 0., 3. ), kRecoW ) );
  axes.emplace( "RecoEhadFrac"   , HistAxis( "Reco E_{had} Fraction",    Binning::Simple( 50, 0., 1. ), kHadEFrac ) );
  axes.emplace( "RecoCosThetaMu" , HistAxis( "Reco cos(#theta_{#mu})",   Binning::Simple( 50, 0., 1. ), kCosNumi ) );
  axes.emplace( "PrimTrkLen"     , HistAxis( "Primary Track Length (m)", Binning::Simple( 60, 0., 15.), kTrkLength ) );

  std::set< std::string > vars_with_quantiles;
  vars_with_quantiles.insert( "RecoEhad" );
  vars_with_quantiles.insert( "RecoEhadVis" );
  vars_with_quantiles.insert( "RecoEmu" );
  vars_with_quantiles.insert( "RecoEnu" );
  vars_with_quantiles.insert( "RecoQ2" );
  vars_with_quantiles.insert( "RecoQ3" );
  vars_with_quantiles.insert( "RecoW" );
  vars_with_quantiles.insert( "RecoCosThetaMu" );

  if ( isMC )
  {
    //axes.emplace( "NucleonPair", HistAxis( "Nucleon Pair"        , Binning::Simple(   5, -1, 4  ), kHitNucBin ) );
    axes.emplace( "TrueQ2"     , HistAxis( "True Q^{2} (GeV^{2})", Binning::Simple(  50, 0., 2. ), kTrueQ2 ) );
    //axes.emplace( "TrueQ0"     , HistAxis( "True q_{0} (GeV)"    , Binning::Simple( 200, 0., 2. ), kTrueQ0 ) );
    //axes.emplace( "TrueQ3"     , HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 200, 0., 2. ), kTrueQ3 ) );
    //axes.emplace( "TrueW"      , HistAxis( "True W (GeV)"        , Binning::Simple(  60, 0., 3. ), kTrueW  ) );
    //axes.emplace( "TrueY"      , HistAxis( "True y"              , Binning::Simple(  50, 0., 1. ), kTrueY  ) );
  }

  axes.emplace
  ( 
    "RecoEhadVis_vs_RecoQ3",
    HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag, "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE )
  );

  /*
  axes.emplace
  (
    "Fine_RecoEhadVis_vs_RecoQ3",
    HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 200, 0., 2. ), kRecoQmag, "Visible E_{had} (GeV)", Binning::Simple( 100, 0., 1. ), kNumuHadVisE )
  );
  */

  if ( isMC )
  {

    axes.emplace
    ( 
      "TrueQ0_vs_TrueQ3",
      HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 80, 0., 2. ), kTrueQ3, "True q_{0} (GeV)", Binning::Simple( 80,0,  2. ), kTrueQ0 )
    );

    /*
    axes.emplace
    ( 
      "Fine_TrueQ0_vs_TrueQ3",
      HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 160, 0., 2. ), kTrueQ3, "True q_{0} (GeV)", Binning::Simple( 160, 0., 2. ), kTrueQ0 )
    );
    */

    std::string title_wres_abs = "W_{reco} - W_{true} (GeV)";
    std::string title_wres_frac = "( W_{reco} - W_{true} ) / W_{true}";

    /*
    axes.emplace
    (
      "WRes_Abs",
      HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW    , title_wres_abs, Binning::Simple( 1000, -10., 10. ), GetWResolution( false, false ) )
    );

    axes.emplace
    (
      "WRes_Abs_Eqn",
      HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW_Eqn, title_wres_abs, Binning::Simple( 1000, -10., 10. ), GetWResolution( false, true ) )
    );

    axes.emplace
    (
      "WRes_Frac",
      HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW    , title_wres_frac, Binning::Simple( 1000, -10., 10. ), GetWResolution( true, false ) )
    );

    axes.emplace
    (
      "WRes_Frac_Eqn",
      HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW_Eqn, title_wres_frac, Binning::Simple( 1000, -10., 10. ), GetWResolution( true, true ) )
    );
    */

    /*
    axes.emplace
    (
      "RecoEhadVis_vs_TrueQ0",
      HistAxis( "True q_{0} (GeV)", Binning::Simple( 40,0,  1. ), kTrueQ0, "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE )
    );

    axes.emplace
    (
      "RecoQ3_vs_TrueQ3",
      HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kTrueQ3, "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag )
    );

    axes.emplace
    (
      "RecoEhad_vs_TrueQ0",
      HistAxis( "True q_{0} (GeV)", Binning::Simple( 40, 0, 2. ), kTrueQ0, "Reco E_{had} (GeV)", Binning::Simple( 40, 0., 2. ), kHadE )
    );
    */

  }

  // Hadronic energy fraction quantiles

  std::string InDirQuant = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
  std::string InFileQuant = InDirQuant + "/quantiles__" + beam + "_full__numu2018.root";
  std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
  TFile* inFile = TFile::Open( pnfs2xrootd(InFileQuant).c_str() );
  TH2* FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
  inFile->Close();
  delete inFile;

  const Cut kSelCuts  = kNumuCutND2018;
  const Var kWeightCV = kPPFXFluxCVWgt * kXSecCVWgt2018RPAFix_noDIS;

  SystShifts kShiftsMECShapeUp;
  kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixNu,     +1 );
  kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixAntiNu, +1 );

  SystShifts kShiftsMECShapeDown;
  kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixNu,     -1 );
  kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixAntiNu, -1 );

  SystShifts kShiftsRPACCQESuppUp;
  kShiftsRPACCQESuppUp.SetShift( &kRPACCQESuppSyst2019, +1 );

  SystShifts kShiftsRPACCQESuppDown;
  kShiftsRPACCQESuppDown.SetShift( &kRPACCQESuppSyst2019, -1 );

  SystShifts kShiftsRPARESUp;
  kShiftsRPARESUp.SetShift( &kRPARESSyst2019, +1 );

  SystShifts kShiftsRPARESDown;
  kShiftsRPARESDown.SetShift( &kRPARESSyst2019, -1 );

  std::vector<Distribution> dists;
  dists.reserve( 100 );

  for ( const auto & a : axes )
  {
    //if ( a.first != "RecoEhadVis" && a.first != "TrueQ0_vs_TrueQ3" ) continue;

    dists.emplace_back( a.first, isMC, loader, a.second, kSelCuts, kWeightCV, kNoShift );

    if ( isMC )
    {
      dists.emplace_back( a.first + "_MECShapeUp"  , isMC, loader, a.second, kSelCuts, kWeightCV, kShiftsMECShapeUp   );
      dists.emplace_back( a.first + "_MECShapeDown", isMC, loader, a.second, kSelCuts, kWeightCV, kShiftsMECShapeDown );
    }
    //if ( a.first == "RecoQ2" || a.first == "RecoEnu" )
    if ( vars_with_quantiles.find( a.first ) != vars_with_quantiles.end() )
    {
      for ( unsigned int icut = 0; icut < HadEFracQuantCuts.size(); ++icut )
      {
        std::string dist_quant_name = a.first + Form( "_EhadQuantile%d", icut + 1 );

        dists.emplace_back( dist_quant_name, isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kNoShift );
        /*
        dists.emplace_back( dist_quant_name + "_MECShapeUp"     , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsMECShapeUp   );
        dists.emplace_back( dist_quant_name + "_MECShapeDown"   , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsMECShapeDown );
        dists.emplace_back( dist_quant_name + "_RPACCQESuppUp"  , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPACCQESuppUp   );
        dists.emplace_back( dist_quant_name + "_RPACCQESuppDown", isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPACCQESuppDown );
        dists.emplace_back( dist_quant_name + "_RPARESUp"       , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPARESUp   );
        dists.emplace_back( dist_quant_name + "_RPARESDown"     , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPARESDown );
        */
      }
    }

    if ( a.first == "RecoQ2" )
    {
      dists.emplace_back( a.first + "_WSelected", isMC, loader, a.second, kSelCuts && kWSelection, kWeightCV, kNoShift );
    }
  }

  /*
  SystShifts kShiftsQELike;
  kShiftsQELike.SetShift( &kMAQEGenieReducedSyst2018, +1 );
  kShiftsQELike.SetShift( &kRPACCQESuppSyst2019, +1 );
  kShiftsQELike.SetShift( &kRPACCQEEnhSyst2019, +1 );
  kShiftsQELike.SetShift( GetGenieRwgtSyst( rwgt::fReweightCCQEPauliSupViaKF ), -1 );
  kShiftsQELike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMaCCRES ), -1 );
  kShiftsQELike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMvCCRES ), -1 );
  //kShiftsQELike.SetShift( &kDirectRelHadEScaleSyst2017, -1 );
  kShiftsQELike.SetShift( &kMECShapeSyst2018RPAFixNu, +1 );
  kShiftsQELike.SetShift( &kMECShapeSyst2018RPAFixAntiNu, +1 );

  SystShifts kShiftsRESLike;
  kShiftsRESLike.SetShift( &kMAQEGenieReducedSyst2018, -1 );
  kShiftsRESLike.SetShift( &kRPACCQESuppSyst2019, -1 );
  kShiftsRESLike.SetShift( &kRPACCQEEnhSyst2019, -1 );
  kShiftsRESLike.SetShift( GetGenieRwgtSyst( rwgt::fReweightCCQEPauliSupViaKF ), +1 );
  kShiftsRESLike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMaCCRES ), +1 );
  kShiftsRESLike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMvCCRES ), +1 );
  //kShiftsRESLike.SetShift( &kDirectRelHadEScaleSyst2017, +1 );
  kShiftsRESLike.SetShift( &kRPARESSyst2019, +1 ); // Turns off RES RPA
  kShiftsRESLike.SetShift( &kMECShapeSyst2018RPAFixNu, -1 );
  kShiftsRESLike.SetShift( &kMECShapeSyst2018RPAFixAntiNu, -1 );

  std::vector< std::string > shift_vars;
  shift_vars.push_back( "RecoEhadVis" );
  shift_vars.push_back( "RecoQ3" );
  shift_vars.push_back( "TrueQ0_vs_TrueQ3" );

  if ( isMC )
  {
    for ( const auto & v : shift_vars )
    {
      dists.emplace_back( v + "_ShiftsQELike" , isMC, loader, axes.at( v ), kSelCuts, kWeightCV, kShiftsQELike  );
      dists.emplace_back( v + "_ShiftsRESLike", isMC, loader, axes.at( v ), kSelCuts, kWeightCV, kShiftsRESLike );
    }
  }
  */

  loader.Go();

  double pot = dists.front().POT();
  
  TTree potTree( "pot", "pot" );
  potTree.Branch( Form( "pot_%s", type.c_str() ), &pot, Form( "pot_%s/D", type.c_str() ) );
  potTree.Fill();

  std::string out_dir = "hist_files";
  if( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );

  std::string out_file_name = Form( "hists_2018_%s_%s_%s", beam.c_str(), type.c_str(), period.c_str() );
  //if ( isMC && !shift_non_mec.empty() )
  //  out_file_name += "_ShiftNonMEC" + shift_non_mec;

  //TFile outputFile( Form( "%s/%s.root", out_dir.c_str(), out_file_name.c_str() ), "RECREATE" );
  TFile outputFile( Form( "%s.root", out_file_name.c_str() ), "RECREATE" );

  potTree.Write();

  for ( const auto & d : dists ) { d.Write(); }

  outputFile.Close();

}

#endif