make_xsec_tuning_hists_mp5.C

Go to the documentation of this file.

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/HistAxis.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 "3FlavorAna/Vars/HistAxes.h"
//#include "CAFAna/Vars/XsecTunes.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/XSecSysts.h"
//#include "MCReweight/func/HistWrapper.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "StandardRecord/Proxy/SRProxy.h"

#include "3FlavorAna/Cuts/NumuCuts2020.h"

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

#include <fstream>

using namespace ana;
using namespace caf;

const Var kMinervaMECGaussEnh( []( const caf::SRProxy* sr )
{
  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;

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

  double norm     = 10.58;
  double mean_q0  = 0.254;
  double mean_q3  = 0.508;
  double sigma_q0 = 0.057;
  double sigma_q3 = 0.129;
  double corr     = 0.875;

  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
             - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );

  return 1.0 + norm * exp( -0.5 * z / ( 1 - corr * corr ) );
});

const Var kDoubleGaussEnhMEC( []( const caf::SRProxy* sr )
{
  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;

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

  // FHC double gauss MEC fit from DocDB 42265-v2
  double norm_1     = 12.78;
  double mean_q0_1  = 0.33;
  double mean_q3_1  = 0.80;
  double sigma_q0_1 = 0.10;
  double sigma_q3_1 = 0.29;
  double corr_1     = 0.91;

  double norm_2     = 38.8;
  double mean_q0_2  = 0.036;
  double mean_q3_2  = 0.46;
  double sigma_q0_2 = 0.042;
  double sigma_q3_2 = 0.233;
  double corr_2     = 0.71;

  double z_1 = pow( ( q0 - mean_q0_1 ) / sigma_q0_1, 2 ) + pow( ( q3 - mean_q3_1 ) / sigma_q3_1, 2 )
             - 2 * corr_1 * ( q0 - mean_q0_1 ) * ( q3 - mean_q3_1 ) / ( sigma_q0_1 * sigma_q3_1 );

  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
             - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );

  return 1.0 + norm_1 * exp( -0.5 * z_1 / ( 1 - corr_1 * corr_1 ) ) + norm_2 * exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) );

});

const Var kMINOSLowQ2RESSupp( []( const caf::SRProxy * sr )
{
  if ( !kIsRes( 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 kVarFermiP( []( const caf::SRProxy* sr )
{
  if ( sr->mc.nnu != 1 ) return 0.;

  int A = sr->mc.nu[0].tgtA;
  int Z = sr->mc.nu[0].tgtZ;
  int N = A-Z;

  if ( A <= 1 ) return 0.;

  // From FermiMomentumForIsoscalarNucleonParametrization in $GENIE/src/Physics/NuclearState/NuclearUtils.cxx
  // Compute Fermi momentum for isoscalar nucleon (in GeV)
  double x = 1.0 / A;
  double P_Fermi = 0.27+x*(-1.12887857491+x*(9.72670908033-39.53095724456*x));

  // From $GENIE/src/Physics/Resonance/XSection/ReinSehgalRESPXSec.cxx
  bool is_p = sr->mc.nu[0].hitnuc == 2212;
  if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
  else     { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }

  return P_Fermi;
});

const Var kVarRESPauliBlock( []( const caf::SRProxy* sr )
{

  if ( !kIsRes( sr ) ) return 1.0;

  int hitnuc = sr->mc.nu[0].hitnuc;
  if ( hitnuc != 2212 && hitnuc != 2112 )
  {
    std::cout << "hitnuc = " << hitnuc << std::endl;
    return 1.0;
  }

  int A = sr->mc.nu[0].tgtA;
  if ( A < 6 ) return 1.0; // From $GENIE/src/Physics/Resonance/XSection/ReinSehgalRESPXSec.cxx

  double P_Fermi = kVarFermiP( sr );

  // From $GENIE/src/Framework/Conventions/Constants.h
  static const double kPionMass       =  1.3957018e-01;          // GeV
  static const double kPionMass2      =  TMath::Power(kPionMass,2);     // GeV^2
  static const double kProtonMass     =  9.38272081e-01;           // GeV
  static const double kNeutronMass    =  9.39565413e-01;           // GeV

  //bool is_p = hitnuc == 2212;
  //double Mnuc   = is_p ? kProtonMass : kNeutronMass;
  double Mnuc   = hitnuc == 2212 ? kProtonMass : kNeutronMass;
  double Mnuc2  = TMath::Power(Mnuc, 2);

  double W = kTrueW( sr );
  if ( !( W > 0 ) )
  {
    std::cout << "W = " << W << std::endl;
    return 1.0;
  }

  double W2 = TMath::Power(W,2);
  double Q2 = kTrueQ2( sr );
  double k  = 0.5 * (W2 - Mnuc2)/Mnuc;

  // From $GENIE/src/Physics/Resonance/XSection/ReinSehgalRESPXSec.cxx
  double FactorPauli_RES = 1.0;

  double k0 = 0., q = 0., q0 = 0.;

  if (P_Fermi > 0.)
  {
    k0 = (W2-Mnuc2-Q2)/(2*W);
    k = TMath::Sqrt(k0*k0+Q2);
    q0 = (W2-Mnuc2+kPionMass2)/(2*W);
    q = TMath::Sqrt(q0*q0-kPionMass2);
  }

  if (2*P_Fermi < k-q)
    FactorPauli_RES = 1.0;
  if (2*P_Fermi >= k+q)
    FactorPauli_RES = ((3*k*k+q*q)/(2*P_Fermi)-(5*TMath::Power(k,4)+TMath::Power(q,4)+10*k*k*q*q)/(40*TMath::Power(P_Fermi,3)))/(2*k);
  if (2*P_Fermi >= k-q && 2*P_Fermi <= k+q)
    FactorPauli_RES = ((q+k)*(q+k)-4*P_Fermi*P_Fermi/5-TMath::Power(k-q, 3)/(2*P_Fermi)+TMath::Power(k-q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*q*k);

  return FactorPauli_RES;

});



/*
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 Cut kIsOther = !kIsNumuCC || !( kIsQE || kIsRes || kIsDIS || kIsDytmanMEC );

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

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

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

class Decomp
{

  public :

    Decomp(
      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[ "RES_PauliBlock"  ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes      , shiftMC, wgtMC * kVarRESPauliBlock );

        fSpectra[ "MEC_MINERvA"             ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC * kMinervaMECGaussEnh );
        fSpectra[ "MEC_DoubleGauss"         ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC * kDoubleGaussEnhMEC );
        fSpectra[ "MEC_DoubleGauss_Numu"    ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC             , shiftMC, wgtMC * kDoubleGaussEnhMEC );
        fSpectra[ "MEC_DoubleGauss_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC          , shiftMC, wgtMC * kDoubleGaussEnhMEC );

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

      }
      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_tuning_hists_mp5
( 
  const std::string beam,
  const std::string type
)
{
  if ( type != "data" && type != "mc" )
  {
    std::cout << "type must be \"data\" or \"mc\"" << std::endl;
    return;
  }
  const bool isMC = ( type == "mc" );

  std::string defn = "";
  if ( beam == "fhc" )
  {
    defn = isMC ?
           "prod_caf_R19-02-23-miniprod5.m_nd_genie_G1810b0211a_nonswap_fhc_nova_v08_full_v1" :
           "prod_caf_R19-02-23-miniprod5.i_nd_numi_fhc_full_v1_goodruns_stride10";
  }
  else if ( beam == "rhc" )
  {
    defn = isMC ?
           "prod_caf_R19-02-23-miniprod5.m_nd_genie_G1810b0211a_nonswap_rhc_nova_v08_full_v1_nogeniereweight" :
           "prod_caf_R19-02-23-miniprod5.i_nd_numi_rhc_full_v1_goodruns_stride20";
  }
  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 ) );

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

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

  }

  // 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  = kNumuCutNDMiniprod5;
  const Var kWeightCV = kPPFXFluxCVWgt;

  std::vector<Decomp> decomps;
  decomps.reserve( 100 );

  for ( const auto & a : axes )
  {
    decomps.emplace_back( a.first, isMC, loader, a.second, kSelCuts, kWeightCV, kNoShift );
    decomps.emplace_back( a.first + "_WSelected", isMC, loader, a.second, kSelCuts && kWSelection, kWeightCV, kNoShift );

    if ( a.first != "RecoEhadFrac" )
    {
      for ( unsigned int icut = 0; icut < HadEFracQuantCuts.size(); ++icut )
      {
        std::string dist_quant_name = a.first + Form( "_EhadQuantile%d", icut + 1 );
        decomps.emplace_back( dist_quant_name, isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kNoShift );
      }
    }
  }

  loader.Go();

  double pot = decomps.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_file_name = Form( "hists_xec_tuning_mp5_%s_%s", beam.c_str(), type.c_str() );
  TFile outputFile( Form( "%s.root", out_file_name.c_str() ), "RECREATE" );

  potTree.Write();

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

  outputFile.Close();

}