QuantileCuts2020.cxx

Go to the documentation of this file.

#include "3FlavorAna/Cuts/QuantileCuts2020.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Utilities.h"
#include "3FlavorAna/Vars/NueEnergy2020.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "3FlavorAna/Vars/HistAxes.h"

#include "TFile.h"
#include "TH2.h"

#include <iostream>

namespace ana
{

  //----------------------------------------------------------------------
  // 2020 Numu EhadFrac quantile cuts

  std::vector<Cut> GetNumuEhadFracQuantCuts2020( const bool isRHC, const unsigned int nquants )
  {
    if ( nquants <= 1 )
    {
      std::cout << "Wrong number of EhadFrac quantiles: " << nquants << std::endl;
      abort();
    }
    std::string beam = isRHC ? "rhc" : "fhc";
    std::string fname = "$NUMUDATA_LIB_PATH/ana2020/Quantiles/quantiles_" + beam + "_full_numu2020.root";
    std::cout << "\nGetting EhadFrac quantile hist from file:\n" << fname << std::endl;
    TFile* inFile = TFile::Open( fname.c_str() );
    TH2* quant_hist = (TH2*)inFile->FindObjectAny( "FDSpec2D" );
    quant_hist->SetDirectory(0);
    inFile->Close();
    delete inFile;
    return QuantileCutsFromTH2( quant_hist, kNumuCCOptimisedAxis2020, kHadEFracAxis2020, nquants );
  }

  //----------------------------------------------------------------------
  // Axes for constructing the Pt and Theta quantiles

  const HistAxis kAxisNumuPt( "Reco Muon #left|#vec{p}_{t}#right| (GeV)"   , Binning::Simple( 100, 0.,  5. ), kNumuMuonPt );
  const HistAxis kAxisNuePt ( "Reco Elecron #left|#vec{p}_{t}#right| (GeV)", Binning::Simple( 100, 0.,  5. ), kNueElectronPt2020 );

  const HistAxis kAxisNumuTheta( "Reco #theta_{#mu} (deg)", Binning::Simple( 180, 0., 180. ), kRecoThetaMu );
  const HistAxis kAxisNueTheta ( "Reco #theta_{e} (deg)"  , Binning::Simple( 180, 0., 180. ), kFirstProngTheta );

  const HistAxis kAxisNueBinsNumuEv( "Reco E_{#nu} (GeV)", Binning::Simple( 10, 0., 5. ), kNumuE2020 );
  const HistAxis kAxisNueBinsNueEv ( "Reco E_{#nu} (GeV)", Binning::Simple( 10, 0., 5. ), kNueEnergy2020 );

  //----------------------------------------------------------------------
  // 2020 Numu Pt quantile cuts

  std::vector<Cut> GetNumuPtQuantCuts2020( const bool isRHC, const unsigned int ehad_quant, const unsigned int nquants )
  {

    if ( ehad_quant < 1 || ehad_quant > 4 )
    {
      std::cout << "Wrong Ehad quant: " << ehad_quant << std::endl;
      abort();
    }
    if ( nquants <= 1 )
    {
      std::cout << "Wrong number of Numu Pt quantiles: " << nquants << std::endl;
      abort();
    }

    std::string beam = isRHC ? "rhc" : "fhc";
    std::string fname = "$NUMUDATA_LIB_PATH/ana2020/Quantiles/pt_quantile_hists_2020_nd_" + beam + ".root";
    std::cout << "\nRetrieving Numu Pt quant hists from file:\n" << fname << std::endl;
    TFile* inFile = TFile::Open( fname.c_str() );
    if ( inFile->IsZombie() )
    {
      std::cout << "Problem with file:\n" << fname << std::endl;
      abort();
    }

    TH2D* quant_hist = (TH2D*)inFile->Get( Form( "RecoPt_vs_RecoEv_NumuBins_Q%d", ehad_quant ) );
    quant_hist->SetDirectory(0);
    inFile->Close();
    delete inFile;

    std::cout << Form( "\nReturning Numu Pt quant cuts for %d quantiles in EhadFrac quantile %d\n", nquants, ehad_quant ) << std::endl;
    return QuantileCutsFromTH2( quant_hist, kNumuCCOptimisedAxis2020, kAxisNumuPt, nquants );
  }

  //----------------------------------------------------------------------
  // 2020 Nue Pt quantile cuts

  std::vector<Cut> GetNueQuantCuts2020( const bool isRHC, const caf::Det_t det, const unsigned int nquants, const ExtrapVar var )
  {

    if ( nquants <= 1 )
    {
      std::cout << "Wrong number of Nue extrap quantiles: " << nquants << std::endl;
      abort();
    }

    std::string varName;
    std::string histName;
    if ( var == kExtrapPt )
    {
      varName  = "Pt";
      histName = "RecoPt_vs_RecoEv_NueBins";
    }
    else if ( var == kExtrapTheta )
    {
      varName  = "Theta";
      histName = "RecoThMu_vs_RecoEv_NueBins";
    }
    else
    {
      std::cout << "Nue extrap variable must be kExtrapPt or kExtrapTheta" << std::endl;
      abort();
    }

    std::string beam = isRHC ? "rhc" : "fhc";
    std::string fname = "$NUMUDATA_LIB_PATH/ana2020/Quantiles/pt_quantile_hists_2020_nd_" + beam + ".root";
    std::cout << "\nRetrieving Nue quantile hists from file:\n" << fname << std::endl;
    TFile* inFile = TFile::Open( fname.c_str() );
    if ( inFile->IsZombie() )
    {
      std::cout << "Problem with file:\n" << fname << std::endl;
      abort();
    }

    std::cout << histName << std::endl;
    TH2D* quant_hist = (TH2D*)inFile->Get( histName.c_str() );
    quant_hist->SetDirectory(0);
    inFile->Close(); 
    delete inFile;

    std::vector<Cut> quant_cuts;

    if ( det == caf::kNEARDET )
    {
      std::cout << Form( "\nReturning cuts for %d ND Nue %s quantiles\n", nquants, varName.c_str() ) << std::endl;
      if      ( var == kExtrapPt    ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNumuEv, kAxisNumuPt   , nquants );
      else if ( var == kExtrapTheta ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNumuEv, kAxisNumuTheta, nquants );
    }
    else if ( det == caf::kFARDET )
    {
      std::cout << Form( "\nReturning cuts for %d FD Nue %s quantiles\n", nquants, varName.c_str() ) << std::endl;
      if      ( var == kExtrapPt    ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNueEv, kAxisNuePt   , nquants );
      else if ( var == kExtrapTheta ) quant_cuts = QuantileCutsFromTH2( quant_hist, kAxisNueBinsNueEv, kAxisNueTheta, nquants );
    }
    else
    {
      std::cout << "Detector must be caf::kNEARDET or caf::kFARDET" << std::endl;
      abort();
    }

    return quant_cuts;

  }

}