ND_DataMC.C

Go to the documentation of this file.

// ND data/MC comparison plots for 3 flavor group.
// Makes POT norm with full syst band and area norm with shape only syst band

#include "CAFAna/Analysis/DataMCPair.h"
#include "3FlavorAna/Plotting/SharedPlotLists.h"
#include "3FlavorAna/Systs/3FlavorAna2020Systs.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "3FlavorAna/Cuts/QuantileCuts2020.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "3FlavorAna/Vars/DecompWeights.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "OscLib/IOscCalc.h"


#include "TFile.h"
#include "TSystem.h"

using namespace ana;

void ND_DataMC(TString sFHC_RHC, TString sNumu_Nue){
  bool isRHC = sFHC_RHC.Contains("rhc") ? true : false;

  // --- datasets --- 
  TString fnameMC = "prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_"+sNumu_Nue+"2020";
  TString fnameData = "prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_"+sNumu_Nue+"2020";

  if(isRHC){
    fnameMC = "prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_"+sNumu_Nue+"2020";
    fnameData = "prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_"+sNumu_Nue+"2020";
  }

  SpectrumLoader ldrMC((std::string)fnameMC);
  ldrMC.SetSpillCut(kStandardSpillCuts);

  SpectrumLoader ldrData((std::string)fnameData);
  ldrData.SetSpillCut(kStandardSpillCuts);

  // --- weights, shifts --- 
  Var weight = kPPFXFluxCVWgt*kXSecCVWgt2020;

  // --- output file --- 
  TString fileOut = sFHC_RHC+"_"+sNumu_Nue+".root";

  // --- numu-specific cuts, plots --- 
  if(sNumu_Nue.Contains("numu")){

    // selection cutflow
    std::vector<Cut> selCuts = {kNumuQuality&&kNumuContainND2020}; //, kNumu2020ND};
    std::vector<std::string> selNames = {"QualCont"}; //, "Full"};

//    std::vector<Cut> selCuts = {kNumu2020ND};
//    std::vector<std::string> selNames = {"Full"};
    unsigned int nSels = selNames.size();

    // quartiles 
    std::vector<Cut> hadEfracCuts = GetNumuEhadFracQuantCuts2020(isRHC, 4);
    std::vector<std::string> hadEfracNames = {"hadEfracQ1", "hadEfracQ2", "hadEfracQ3", "hadEfracQ4"};
    unsigned int nhadEfracQuants = hadEfracNames.size();
    
    std::vector<Cut> quantCuts;
    std::vector<std::string> quantNames;
    for(unsigned int q1 = 0; q1<nhadEfracQuants; ++q1){
      std::vector<Cut> ptCuts = GetNumuPtQuantCuts2020(isRHC, q1+1, 3);
      ptCuts.push_back(kNoCut);

      std::vector<std::string> ptNames = {"pTfracQ1", "pTfracQ2", "pTfracQ3", "pTfracAllQs"};
      unsigned int nptQuants = ptNames.size();

      for(unsigned int q2 = 0; q2<nptQuants; ++q2){
        quantCuts.push_back(hadEfracCuts[q1]&&ptCuts[q2]);
        quantNames.push_back(hadEfracNames[q1]+"_"+ptNames[q2]);
      }
    }

    quantCuts.push_back(kNoCut);
    quantNames.push_back("hadEfracAllQs_pTfracAllQs");
    unsigned int nQuants = quantNames.size();

    // variables 
    std::vector<HistDef> defs;
    AddHistDefNumuNDDataMC(defs);
    unsigned int nVars = defs.size(); 
    std::vector<std::string> varNames; std::vector<HistAxis> axis;
    for(unsigned int var = 0; var<nVars; ++var){
      varNames.emplace_back(defs[var].name);
      axis.emplace_back(defs[var].axis);
    }

    // systematics
    std::vector<const ISyst*> systs = get3FlavorAna2020AllSysts(EAnaType2020::kNumuAna, true, BeamType2020::kFHC);
    if(isRHC) systs = get3FlavorAna2020AllSysts(EAnaType2020::kNumuAna, true, BeamType2020::kRHC);
    DataMCPair* numuDataMC[nVars][nSels][nQuants];

    Cut WrongSign = kIsAntiNu;
    if(isRHC) WrongSign = kIsNu;

    for(unsigned int var = 0; var<nVars; ++var){
      for(unsigned int sel = 0; sel<nSels; ++sel){
        for(unsigned int quant = 0; quant<nQuants; ++quant){
          Cut thisCut = selCuts[sel] && quantCuts[quant];
          numuDataMC[var][sel][quant] = new DataMCPair(thisCut, axis[var], ldrData, ldrMC, weight, systs, {{"Total Background", !kIsNumuCC},{"Wrong Sign", WrongSign}});
        }
      }
    }

    ldrMC.Go();
    ldrData.Go();

    TFile* fout = new TFile(fileOut, "RECREATE");

    for(unsigned int var = 0; var<nVars; ++var){
      for(unsigned int sel = 0; sel<nSels; ++sel){
        for(unsigned int quant = 0; quant<nQuants; ++quant){
          TString name = varNames[var]+"_"+selNames[sel]+"_"+quantNames[quant];
          numuDataMC[var][sel][quant]->SaveTo(fout, name);
        }
      }
    }
  }

  // --- nue-specific cuts, plots --- 
  else{
    // selection cutflow
    double highedge;
    if(sFHC_RHC.Contains("fhc")) highedge = kNue2020CVNFHCHighEdge;
    else highedge = kNue2020CVNRHCHighEdge;
    Cut CVNeHigh = kCVNe_looseptp > highedge;
    std::vector<std::string> selNames = {"Presel", "CVNeLow", "CVNeHigh"};
    std::vector<Cut> selCuts = {kNue2020NDPresel, kNue2020ND, kNue2020ND&&CVNeHigh};
    unsigned int nSels = selNames.size();

    std::string decompFName = "nue2020_rhc_prop_decomp_ratios.root";
    std::string decompShortName = ""; //"prop";
    if(sFHC_RHC.Contains("fhc")) decompFName = "nue2020_fhc_combo_decomp_ratios.root";

    const Var kDecompWeight(DecompWeightFunc(FindCAFAnaDir()+
                                             "/nue/decomp/" + decompFName,
                                             decompShortName,
                                             kNue2020AnaBin));
    // variables
    std::vector<HistDef> defs;
    AddHistDefRecoND(defs);
    AddHistDefSlice(defs);
    AddHistDefNueND(defs);
    AddHistDefNueShower(defs);
    AddHistDefNuePID(defs);
    AddHistDefNueSelectionExtras(defs);
    AddHistDefProngCVN(defs);
    AddHistDefNueEnergy(defs);
    unsigned int nVars = defs.size(); 
    std::vector<std::string> varNames; std::vector<HistAxis> axis;
    for(unsigned int var = 0; var<nVars; ++var){
      varNames.emplace_back(defs[var].name);
      axis.emplace_back(defs[var].axis);
    }

    Spectrum* specData[nVars][nSels];
    IPrediction* predMC[nVars][nSels];
    IPrediction* predMCdecomp[nVars][nSels];

    for(unsigned int var = 0; var<nVars; ++var){
      for(unsigned int sel = 0; sel<nSels; ++sel){
        Cut thisCut = selCuts[sel];
        predMC[var][sel] = new PredictionNoExtrap(ldrMC, kNullLoader, kNullLoader, axis[var], thisCut, kNoShift, weight);
        predMCdecomp[var][sel] = new PredictionNoExtrap(ldrMC, kNullLoader, kNullLoader, axis[var], thisCut, kNoShift, weight*kDecompWeight);
        specData[var][sel] = new Spectrum(ldrData, axis[var], thisCut, kNoShift);
      }
    }

    ldrMC.Go();
    ldrData.Go();

    TString fileOut1(fileOut); fileOut1.ReplaceAll(".root", "_mc.root");
    TString fileOut2(fileOut); fileOut2.ReplaceAll(".root", "_data.root");
    TString fileOut3(fileOut); fileOut3.ReplaceAll(".root", "_mc_decomp.root");

    TFile* fout_mc = new TFile(fileOut1, "RECREATE");
    TFile* fout_data = new TFile(fileOut2, "RECREATE");
    TFile* fout_mcdecomp = new TFile(fileOut3, "RECREATE");

    std::vector<TString> compNames = {"nue", "numu", "antinue", "antinumu", "nc"};
    unsigned int nComp = compNames.size();

    for(unsigned int var = 0; var<nVars; ++var){
      for(unsigned int sel = 0; sel<nSels; ++sel){
        TString name = varNames[var]+"_"+selNames[sel];
        predMC[var][sel]->ComponentCC(12,12).Unoscillated().SaveTo(fout_mc, name+"_nue");
        predMC[var][sel]->ComponentCC(14,14).Unoscillated().SaveTo(fout_mc, name+"_numu");
        predMC[var][sel]->ComponentCC(-12,-12).Unoscillated().SaveTo(fout_mc, name+"_antinue");
        predMC[var][sel]->ComponentCC(-14,-14).Unoscillated().SaveTo(fout_mc, name+"_antinumu");
        predMC[var][sel]->ComponentNCTotal().SaveTo(fout_mc, name+"_nc");

        predMCdecomp[var][sel]->ComponentCC(12,12).Unoscillated().SaveTo(fout_mcdecomp, name+"_nue");
        predMCdecomp[var][sel]->ComponentCC(14,14).Unoscillated().SaveTo(fout_mcdecomp, name+"_numu");
        predMCdecomp[var][sel]->ComponentCC(-12,-12).Unoscillated().SaveTo(fout_mcdecomp, name+"_antinue");
        predMCdecomp[var][sel]->ComponentCC(-14,-14).Unoscillated().SaveTo(fout_mcdecomp, name+"_antinumu");
        predMCdecomp[var][sel]->ComponentNCTotal().SaveTo(fout_mcdecomp, name+"_nc");

        specData[var][sel]->SaveTo(fout_data, name);
      }
    }
  }
}