hyperon_nom_macro.C

Go to the documentation of this file.

#pragma once

#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/MultiVar.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SystShifts.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"

#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionGenerator.h"

#include "3FlavorAna/Systs/EnergySysts.h"

#include "Utilities/rootlogon.C"

#include "TCanvas.h"
#include "TH2.h"
#include "TProfile.h"
#include "TSystem.h"

using namespace ana;

void hyperon_nom_macro(bool IsFarDet = true)
{

  std::string sFD_ND   = ( IsFarDet == true ? "FD"  : "ND"   );
  std::string sOutFile = "nom_"+sFD_ND+".root";

  // --- Loaders ---
  std::string nom = "";
  if (IsFarDet) {
    nom = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_rhc_nova_v08_full_v1"; // 
  } 
  else {
    nom = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_rhc_nova_v08_full_v1";
  }

  SpectrumLoader ldr(nom);
  ldr.SetSpillCut(kStandardSpillCuts);

  // --- Weights ---
  const Var weight = kPPFXFluxCVWgt *  kXSecCVWgt2018;

  // --- Interaction types ---
  const std::vector<Cut>         GENIECuts = { kNoCut, kIsDytmanMEC, kIsDIS, kIsRes, kIsCoh, kIsQE };
  const std::vector<std::string> GENIEStr  = { "All", "MEC", "DIS", "RES", "COH", "QE"};
  const unsigned int nGENIE = GENIECuts.size();

  // --- Selection Cuts ---
  std::vector<std::string> cutNames; std::vector<Cut> kDetCuts;
  if (IsFarDet) {
    cutNames.emplace_back("kNoCut");        kDetCuts.emplace_back(kNoCut);
    cutNames.emplace_back("qual");          kDetCuts.emplace_back(kNumuQuality);
    cutNames.emplace_back("qual+cont");     kDetCuts.emplace_back(kNumuQuality && kNumuContainFD2017);
    cutNames.emplace_back("qual+cont+PID"); kDetCuts.emplace_back(kNumuQuality && kNumuContainFD2017 && kNumuPID2018);
    cutNames.emplace_back("full");          kDetCuts.emplace_back(kNumuCutFD2018);
  } 
  else {
    cutNames.emplace_back("kNoCut");        kDetCuts.emplace_back(kNoCut);
    cutNames.emplace_back("qual");          kDetCuts.emplace_back(kNumuQuality);
    cutNames.emplace_back("qual+cont");     kDetCuts.emplace_back(kNumuQuality && kNumuContainND2017);
    cutNames.emplace_back("full");          kDetCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2018);
  }
  unsigned int nDet = kDetCuts.size();

  // --- Quantiles ---
  auto cvmfs_dir = std::getenv("NUMUDATA_DIR");
  if (!cvmfs_dir) { std::cerr << "Couldn't find UPS dir for numu prediction!" << std::endl; return; }
  std::string fdspecfile = std::string(cvmfs_dir) + "/ana2018/Quantiles/quantiles__rhc_full__numu2018.root";

  TFile* inFile = TFile::Open(fdspecfile.c_str());
  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); 
  const int NHadEFracQuantiles = 4;
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
  HadEFracQuantCuts.push_back( kNoCut );
  std::vector<std::string> quantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };
  const unsigned int nQ = quantNames.size();

  // --- Additional Vars --- 
  // reco - true
  const Var kRmT([](const caf::SRProxy* sr){
    if (sr->mc.nnu == 0) return -5.;
    return kCCE(sr) - kTrueE(sr);
  });

  // (reco - true) / true 
  const Var kRmToT([](const caf::SRProxy* sr){
    if (sr->mc.nnu == 0) return -5.;
    if(kTrueE(sr) == 0) return -5.;
    return (kCCE(sr) - kTrueE(sr))/kTrueE(sr);
  });

  // number of prongs
  const Var kNprongs([](const caf::SRProxy* sr){
    double npng = -5;
    if( !sr->vtx.elastic.IsValid ) return npng;
    npng = sr->vtx.elastic.fuzzyk.npng;
    return npng;
  });

  // W
  const Var kW([](const caf::SRProxy* sr){
    double W = -5;
    if(sr->mc.nnu == 0) return W;
    if(sr->mc.nu[0].W2<0) return W;
    W = sqrt(sr->mc.nu[0].W2);
    return W;
  });

  // vertex X
  const Var kVx([](const caf::SRProxy* sr){
    double pos = -5;
    if(sr->mc.nnu == 0) return pos;
    pos = sr->mc.nu[0].vtx.x;
    return pos;
  });

  // vertex Y
  const Var kVy([](const caf::SRProxy* sr){
    double pos = -5;
    if(sr->mc.nnu == 0) return pos;
    pos = sr->mc.nu[0].vtx.y;
    return pos;
  });

  // vertex Z
  const Var kVz([](const caf::SRProxy* sr){
    double pos = -5;
    if(sr->mc.nnu == 0) return pos;
    pos = sr->mc.nu[0].vtx.z;
    return pos;
  });

  // --- Binning ---
  const Binning RmT_bins   = Binning::Simple(80, -2, 2);
  const Binning RmToT_bins = Binning::Simple(200, -1, 1);
  const Binning mu_bins    = Binning::Simple(100, 0, 5);
  const Binning had_bins   = Binning::Simple(40, 0, 2);
  const Binning png_bins   = Binning::Simple(10, 0, 9);
  const Binning W_bins     = Binning::Simple(60, 0, 3);
  const Binning id_bins    = Binning::Simple(120, -0.1, 1.1);
  const Binning xybins     = Binning::Simple(160, -800, 800);
  const Binning zbins      = Binning::Simple(610, 0, 6100);
  const Binning TLbins     = Binning::Simple(600,0.0,60.0);
  const Binning Scatbins   = Binning::Simple(200,-0.5,0.5); 
  const Binning dEdxbins   = Binning::Simple(200,-3.0,2.0); 
  const Binning MFbins     = Binning::Simple(120,-0.1,1.1); 

  // --- Spectra ---
  Spectrum *RmT_nom         [nDet][nQ][nGENIE]; 
  Spectrum *RmToT_nom       [nDet][nQ][nGENIE]; 
  Spectrum *trueHadE_nom    [nDet][nQ][nGENIE];
  Spectrum *recoHadE_nom    [nDet][nQ][nGENIE];
  Spectrum *png_nom         [nDet][nQ][nGENIE];
  Spectrum *trueNuE_nom     [nDet][nQ][nGENIE];
  Spectrum *recoNuE_nom     [nDet][nQ][nGENIE];
  Spectrum *trueMuonE_nom   [nDet][nQ][nGENIE];
  Spectrum *recoMuonE_nom   [nDet][nQ][nGENIE];
  Spectrum *remid_nom       [nDet][nQ][nGENIE];
  Spectrum *cvnprod3_nom    [nDet][nQ][nGENIE];
  Spectrum *cvn2017_nom     [nDet][nQ][nGENIE];
  Spectrum *Vx_nom          [nDet][nQ][nGENIE];
  Spectrum *Vy_nom          [nDet][nQ][nGENIE];
  Spectrum *Vz_nom          [nDet][nQ][nGENIE]; 
  Spectrum *trklen_nom      [nDet][nQ][nGENIE];
  Spectrum *scatLL_nom      [nDet][nQ][nGENIE];
  Spectrum *dedxLL_nom      [nDet][nQ][nGENIE];
  Spectrum *MF_nom          [nDet][nQ][nGENIE];

  for (unsigned int det=0; det<nDet; ++det) {
    for (unsigned int quant=0; quant<nQ; ++quant) {
      for (unsigned int gen=0; gen<nGENIE; ++gen){
        Cut nom_cut = GENIECuts[gen] && kDetCuts[det] && HadEFracQuantCuts[quant];
        std::string nomStr = " for Cut "+cutNames[det]+", "+quantNames[quant] +", "+GENIEStr[gen];

        RmT_nom         [det][quant][gen] = new Spectrum("Reco - True (GeV)"+nomStr, RmT_bins, ldr, kRmT, nom_cut, kNoShift, weight);
        RmToT_nom       [det][quant][gen] = new Spectrum("(Reco - True)/True"+nomStr, RmToT_bins, ldr, kRmToT, nom_cut, kNoShift, weight);
        trueHadE_nom    [det][quant][gen] = new Spectrum("true nu - reco muon"+nomStr, had_bins, ldr, kTrueE-kMuE, nom_cut);
        recoHadE_nom    [det][quant][gen] = new Spectrum("Hadronic Energy"+nomStr, had_bins, ldr, kHadE, nom_cut);
        png_nom         [det][quant][gen] = new Spectrum("Number of Prongs"+nomStr, png_bins, ldr, kNprongs, nom_cut);
        trueNuE_nom     [det][quant][gen] = new Spectrum("kTrueE"+nomStr, mu_bins, ldr, kTrueE, nom_cut);
        recoNuE_nom     [det][quant][gen] = new Spectrum("kCCE"+nomStr, mu_bins, ldr, kCCE, nom_cut);
        trueMuonE_nom   [det][quant][gen] = new Spectrum("kTrueMuonE"+nomStr, mu_bins, ldr, kTrueMuonE, nom_cut);
        recoMuonE_nom   [det][quant][gen] = new Spectrum("kMuE"+nomStr, mu_bins, ldr, kMuE, nom_cut);
        remid_nom       [det][quant][gen] = new Spectrum("remid"+nomStr, id_bins, ldr, SIMPLEVAR(sel.remid.pid), nom_cut);
        cvnprod3_nom    [det][quant][gen] = new Spectrum("cvnprod3"+nomStr, id_bins, ldr, SIMPLEVAR(sel.cvnProd3Train.numuid), nom_cut);
        cvn2017_nom     [det][quant][gen] = new Spectrum("cvn2017"+nomStr, id_bins, ldr, SIMPLEVAR(sel.cvn2017.numuid), nom_cut);
        Vx_nom          [det][quant][gen] = new Spectrum("true vertex x"+nomStr, xybins, ldr, kVx, nom_cut);
        Vy_nom          [det][quant][gen] = new Spectrum("true vertex y"+nomStr, xybins, ldr, kVy, nom_cut);
        Vz_nom          [det][quant][gen] = new Spectrum("true vertex z"+nomStr, zbins, ldr, kVz, nom_cut);
        trklen_nom      [det][quant][gen] = new Spectrum("trklen"+nomStr, TLbins, ldr, kTrkLength, nom_cut);
        scatLL_nom      [det][quant][gen] = new Spectrum("scatLL"+nomStr, Scatbins, ldr, kReMIdScatLLH, nom_cut);
        dedxLL_nom      [det][quant][gen] = new Spectrum("dEdxLL"+nomStr, dEdxbins, ldr, kReMIdDEDxLLH, nom_cut);
        MF_nom          [det][quant][gen] = new Spectrum("meas. frac."+nomStr, MFbins, ldr, kReMIdMeasFrac, nom_cut);
      }
    }
  }
  
  // --- Do all the things!
  ldr.Go();

  // --- Save all the things!
  TFile *OutFile = new TFile(sOutFile.c_str(),"RECREATE");
  for (unsigned int det=0; det<nDet; ++det) {
    for (unsigned int quant=0; quant<nQ; ++quant) {
      for (unsigned int gen=0; gen<nGENIE; ++gen) {
        std::string nomStr = cutNames[det]+"_"+quantNames[quant]+"_"+GENIEStr[gen]; 
        RmT_nom         [det][quant][gen] -> SaveTo(OutFile, TString("RmT_nom_")+TString(nomStr));
        RmToT_nom       [det][quant][gen] -> SaveTo(OutFile, TString("RmToT_nom_")+TString(nomStr));
        trueHadE_nom    [det][quant][gen] -> SaveTo(OutFile, TString("trueHadE_nom_")+TString(nomStr));
        recoHadE_nom    [det][quant][gen] -> SaveTo(OutFile, TString("recoHadE_nom_")+TString(nomStr));
        png_nom         [det][quant][gen] -> SaveTo(OutFile, TString("png_nom_")+TString(nomStr));
        trueNuE_nom     [det][quant][gen] -> SaveTo(OutFile, TString("trueNuE_nom_")+TString(nomStr));
        recoNuE_nom     [det][quant][gen] -> SaveTo(OutFile, TString("recoNuE_nom_")+TString(nomStr));
        trueMuonE_nom   [det][quant][gen] -> SaveTo(OutFile, TString("trueMuonE_nom_")+TString(nomStr));
        recoMuonE_nom   [det][quant][gen] -> SaveTo(OutFile, TString("recoMuonE_nom_")+TString(nomStr));
        remid_nom       [det][quant][gen] -> SaveTo(OutFile, TString("remid_nom_")+TString(nomStr));
        cvnprod3_nom    [det][quant][gen] -> SaveTo(OutFile, TString("cvnprod3_nom_")+TString(nomStr));
        cvn2017_nom     [det][quant][gen] -> SaveTo(OutFile, TString("cvn2017_nom_")+TString(nomStr));
        Vx_nom          [det][quant][gen] -> SaveTo(OutFile, TString("Vx_nom_")+TString(nomStr));
        Vy_nom          [det][quant][gen] -> SaveTo(OutFile, TString("Vy_nom_")+TString(nomStr));
        Vz_nom          [det][quant][gen] -> SaveTo(OutFile, TString("Vz_nom_")+TString(nomStr));
        trklen_nom      [det][quant][gen] -> SaveTo(OutFile, TString("trklen_nom_")+TString(nomStr));
        scatLL_nom      [det][quant][gen] -> SaveTo(OutFile, TString("scatLL_nom_")+TString(nomStr));
        dedxLL_nom      [det][quant][gen] -> SaveTo(OutFile, TString("dedxLL_nom_")+TString(nomStr));
        MF_nom          [det][quant][gen] -> SaveTo(OutFile, TString("MF_nom_")+TString(nomStr));            
      }
    }
  }
  OutFile->Close();
}