NueAcceptSysts.cxx

Go to the documentation of this file.

#include "3FlavorAna/Systs/NueAcceptSysts.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "3FlavorAna/Vars/NueVars.h"

#include "CAFAna/Core/Utilities.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "TFile.h"
#include "TH1.h"
#include "TKey.h"
#include "TMath.h"

#include <iostream>

namespace ana
{
  /// FHC
  const NueAcceptSystBkg2018FHC kNueAcceptSystBkg2018FHC;
  const NueAcceptSystSignalKin2018FHC kNueAcceptSystSignalKin2018FHC;
  const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKin2020FHC;
  const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKinPtExtrap2020FHC(ana::kExtrapPt, "accept_signalkin_pTextrap_FHC_2020");


  /// RHC
  const NueAcceptSystBkg2018RHC kNueAcceptSystBkg2018RHC;
  const NueAcceptSystSignalKin2018RHC kNueAcceptSystSignalKin2018RHC;
  const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKin2020RHC;
  const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKinPtExtrap2020RHC(ana::kExtrapPt, "accept_signalkin_pTextrap_RHC_2020");


  /// FHC implementation
  void NueAcceptSystBkg2018FHC::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    // Nue acceptance systs of bkg from docdb-27935 ND subcomponents extrapolation
    // maximal core / peripheral FHC (RHC) difference: +- 1.5% / +-1.2% (+- 4.1% / +- 2.3%)

    /// FHC values
    const double core = .015;
    const double peri = .012;

    // Check if FHC, otherwise left unaltered
    if(kIsRHC(sr)) return;

    // Check if in FD, otherwise left unaltered
    if(sr->hdr.det != caf::kFARDET) return;

    // Check if neutrino, but not numu->nue signal, otherwise left unaltered
    if(sr->mc.nnu == 0) return;
    if(abs(sr->mc.nu[0].pdg) == 12 && abs(sr->mc.nu[0].pdgorig) == 14 && sr->mc.nu[0].iscc) return;
    // Left unaltered if tau neutrino
    if(abs(sr->mc.nu[0].pdg) == 16 && sr->mc.nu[0].iscc) return;

    if(sigma == 0) return;

    double kFDBkgCorr = core;

    // Check whether in peripheral or core sample
    if(kNue2018FDPeripheral(sr))
      kFDBkgCorr = peri;

    weight *= 1+kFDBkgCorr*sigma;
  }

  void NueAcceptSystSignalKin2018FHC::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    // Check if FHC, otherwise left unaltered
    if(kIsRHC(sr)) return;

    // Check if in FD, otherwise left unaltered
    if(sr->hdr.det != caf::kFARDET) return;

    // Check if signal or WS
    if(sr->mc.nnu == 0) return;
    if(
       abs(sr->mc.nu[0].pdg) != 12 ||
       abs(sr->mc.nu[0].pdgorig) != 14 ||
       !(sr->mc.nu[0].iscc)
      ) return;

    if(!fWeightHistFHC){

      // FD in file from docdb-27935 ExtrapSyst technote
      // Can be remade with 3FlavorAna/Ana2018/ExtrapSysts/ scripts
      // Using only correction of the weight with max impact, i.e. Cos weight for FHC and trueQ2 weight for RHC, see technote

      /// FHC
      const std::string kWeightsFnameFHC = FindCAFAnaDir()+"/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_FHC.root"; 
      const std::string kWeightNameFHC = "Cos_FD_KinematicsCorrection_FHC";

      TFile weightsFileFHC (kWeightsFnameFHC.c_str(),"read");

      if(weightsFileFHC.IsZombie()){
        std::cerr << "Warning: couldn't open " << kWeightsFnameFHC << std::endl;
        abort();
      }

      fWeightHistFHC = (TH1*) weightsFileFHC.Get(kWeightNameFHC.c_str())->Clone();

      // disassociate the histogram from files
      fWeightHistFHC -> SetDirectory(0);

      weightsFileFHC.Close();
    }

    if (kNue2018FDAllSamples(sr))
    {
      double kAnaBin = kNue2018AnaBin(sr);
      double kFDSignalWeight = 1;

      if(fWeightHistFHC->GetBinContent(kAnaBin) != 0)
        kFDSignalWeight = fWeightHistFHC->GetBinContent(kAnaBin);

      weight *= 1+(kFDSignalWeight-1)*sigma;
    }
  }

  void NueAcceptSystSignalKin2020FHC::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    // Check if FHC, otherwise left unaltered
    if(kIsRHC(sr)) return;

    // Check if in FD, otherwise left unaltered
    if(sr->hdr.det != caf::kFARDET) return;

    // Check if signal or WS
    if(sr->mc.nnu == 0) return;
    if(
       abs(sr->mc.nu[0].pdg) != 12 ||
       abs(sr->mc.nu[0].pdgorig) != 14 ||
       !(sr->mc.nu[0].iscc)
      ) return;

    if(!fWeightHistFHC){

      // FD in file from docdb-43917 ExtrapSyst technote
      // Can be remade with 3FlavorAna/Ana2020/AcceptSysts/ scripts
      // Using only correction of the weight with max impact, which is PtP weight for FHC with pT extrap and trueQ2 without pT extrap, see technote

      /// FHC
      std::string kWeightsFnameFHC;
      std::string kWeightNameFHC;

      if(extrapVar == ana::kExtrapEmpty){
        kWeightsFnameFHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_FHC.root";
        kWeightNameFHC = "trueQ2_FD_KinematicsCorrection_FHC";
      }
      else if(extrapVar == ana::kExtrapPt){
        kWeightsFnameFHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_pTExtrap_FHC.root";
        kWeightNameFHC = "PtP_FD_KinematicsCorrection_FHC";
      }
      else{
        std::cerr << "Error: Only pT extrap or no extrap options are available for acceptance syst!" << std::endl;
        abort();
      }

      TFile weightsFileFHC (kWeightsFnameFHC.c_str(),"read");

      if(weightsFileFHC.IsZombie()){
        std::cerr << "Warning: couldn't open " << kWeightsFnameFHC << std::endl;
        abort();
      }

      fWeightHistFHC = (TH1*) weightsFileFHC.Get(kWeightNameFHC.c_str())->Clone();

      // disassociate the histogram from files
      fWeightHistFHC -> SetDirectory(0);

      weightsFileFHC.Close();
    }

    if (kNue2020FDAllSamples(sr))
    {
      double kAnaBin = kNue2020AnaBin(sr) +1 ; //Need to add +1 to get the right bin# for GetBinContent unlike kNue2020Axis where we Fill the histogram
      double kFDSignalWeight = 1;

      if(fWeightHistFHC->GetBinContent(kAnaBin) != 0)
        kFDSignalWeight = fWeightHistFHC->GetBinContent(kAnaBin);

      weight *= 1+(kFDSignalWeight-1)*sigma;
    }
  }

  /// RHC implementation
  void NueAcceptSystBkg2018RHC::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    // Nue acceptance systs of bkg from docdb-27935 ND subcomponents extrapolation
    // maximal core / peripheral FHC (RHC) difference: +- 1.5% / +-1.2% (+- 4.1% / +- 2.3%)

    /// RHC values
    const double core = .041;
    const double peri = .023;

    // Check if RHC, otherwise left unaltered
    if(!kIsRHC(sr)) return;

    // Check if in FD, otherwise left unaltered
    if(sr->hdr.det != caf::kFARDET) return;

    // Check if neutrino, but not numu->nue signal, otherwise left unaltered
    if(sr->mc.nnu == 0) return;
    if(abs(sr->mc.nu[0].pdg) == 12 && abs(sr->mc.nu[0].pdgorig) == 14 && sr->mc.nu[0].iscc) return;
    // Left unaltered if tau neutrino
    if(abs(sr->mc.nu[0].pdg) == 16 && sr->mc.nu[0].iscc) return;

    if(sigma == 0) return;

    double kFDBkgCorr = core;

    // Check whether in peripheral or core sample
    if(kNue2018FDPeripheral(sr))
      kFDBkgCorr = peri;

    weight *= 1+kFDBkgCorr*sigma;
  }

  void NueAcceptSystSignalKin2018RHC::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    // Check if RHC, otherwise left unaltered
    if(!kIsRHC(sr)) return;
    
    // Check if in FD, otherwise left unaltered
    if(sr->hdr.det != caf::kFARDET) return;

    // Check if signal or WS
    if(sr->mc.nnu == 0) return;
    if(
       abs(sr->mc.nu[0].pdg) != 12 ||
       abs(sr->mc.nu[0].pdgorig) != 14 ||
       !(sr->mc.nu[0].iscc)
      ) return;

    if(!fWeightHistRHC){

      // FD in file from docdb-27935 ExtrapSyst technote
      // Can be remade with CAFAna/nue/Ana2018/ExtrapSysts/ scripts
      // Using only correction of the weight with max impact, i.e. Cos weight for FHC and trueQ2 weight for RHC, see technote

      /// RHC
      const std::string kWeightsFnameRHC = FindCAFAnaDir()+"/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_RHC.root"; 
      const std::string kWeightNameRHC = "trueQ2_FD_KinematicsCorrection_RHC";

      TFile weightsFileRHC (kWeightsFnameRHC.c_str(),"read");
      
      if(weightsFileRHC.IsZombie()){
        std::cerr << "Warning: couldn't open " << kWeightsFnameRHC << std::endl;
        abort();
      }
      
      fWeightHistRHC = (TH1*) weightsFileRHC.Get(kWeightNameRHC.c_str())->Clone();

      //disassociate the histogram from files
      fWeightHistRHC -> SetDirectory(0);
      
      weightsFileRHC.Close();
    }
    
    if (kNue2018FDAllSamples(sr))
      {
        double kAnaBin = kNue2018AnaBin(sr);
        double kFDSignalWeight = 1;

        if(fWeightHistRHC->GetBinContent(kAnaBin) != 0)
          kFDSignalWeight = fWeightHistRHC->GetBinContent(kAnaBin);
        
        weight *= 1+(kFDSignalWeight-1)*sigma;
      }
  }
  
  void NueAcceptSystSignalKin2020RHC::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    // Check if RHC, otherwise left unaltered
    if(!kIsRHC(sr)) return;

    // Check if in FD, otherwise left unaltered
    if(sr->hdr.det != caf::kFARDET) return;

    // Check if signal or WS
    if(sr->mc.nnu == 0) return;
    if(
       abs(sr->mc.nu[0].pdg) != 12 ||
       abs(sr->mc.nu[0].pdgorig) != 14 ||
       !(sr->mc.nu[0].iscc)
      ) return;

    if(!fWeightHistRHC){

      // FD in file from docdb-43917 ExtrapSyst technote
      // Can be remade with 3FlavorAna/Ana2020/AcceptSysts/ scripts
      // Using only correction of the weight with max impact, which is pT/p for RHC without pT extrap and cos for RHC with pT extrap, see technote

      /// RHC
      std::string kWeightsFnameRHC;
      std::string kWeightNameRHC;
      
      if(extrapVar == ana::kExtrapEmpty){
        kWeightsFnameRHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_RHC.root";
        kWeightNameRHC = "PtP_FD_KinematicsCorrection_RHC";
      }
      else if(extrapVar == ana::kExtrapPt){
        kWeightsFnameRHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_pTExtrap_RHC.root";
        kWeightNameRHC = "Cos_FD_KinematicsCorrection_RHC";
      }
      else{
        std::cerr << "Error: Only pT extrap and no extrap options are available for acceptance syst!" << std::endl;
        abort();
      }

      TFile weightsFileRHC (kWeightsFnameRHC.c_str(),"read");

      if(weightsFileRHC.IsZombie()){
        std::cerr << "Warning: couldn't open " << kWeightsFnameRHC << std::endl;
        abort();
      }

      fWeightHistRHC = (TH1*) weightsFileRHC.Get(kWeightNameRHC.c_str())->Clone();

      //disassociate the histogram from files
      fWeightHistRHC -> SetDirectory(0);

      weightsFileRHC.Close();
    }

    if (kNue2020FDAllSamples(sr))
    {
      double kAnaBin = kNue2020AnaBin(sr) + 1;  //Need to add +1 to get the right bin# for GetBinContent unlike kNue2020Axis where we Fill the histogram
      double kFDSignalWeight = 1;

      if(fWeightHistRHC->GetBinContent(kAnaBin) != 0)
        kFDSignalWeight = fWeightHistRHC->GetBinContent(kAnaBin);

      weight *= 1+(kFDSignalWeight-1)*sigma;
    }
  }
}