make_RHC_WrongSign_Numu_MC.C

Go to the documentation of this file.

// ----------------------------------------------------------------
// Make the Predictions at ND for numu, MC only,
//       to be combined with Matt's neutron study
// See docDB-29919 for tech note on WS study
// ----------------------------------------------------------------

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/WrongSignCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Vars/WrongSignBDTVars.h"

#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionNoOsc.h"
#include "OscLib/OscCalcPMNSOpt.h"

#include "TCanvas.h"
#include "TH2.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TFile.h"
#include "TProfile.h"
#include "TString.h"

using namespace ana;

void make_RHC_WrongSign_Numu_MC()
{
  // ND MC data set, Prod 4 RHC
  // nue
  std::string fname_numuRHC = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018";
  SpectrumLoader ld_numuRHC(fname_numuRHC);
  ld_numuRHC.SetSpillCut(kStandardSpillCuts);

  // Standard weights
  Var wei = kXSecCVWgt2018 * kPPFXFluxCVWgt;

  // Extra Cuts and Vars for analysis
  // -------------------------------------------------------------------------------------------------
  const Cut kCutHighCVNSSe = kCVNSSe >= 0.98; // Ana 2018 Hi CVN bin (docDB-27317)
  const Cut kCutLowCVNSSe  = kCVNSSe >= 0.89; // Ana 2018 Lo CVN bin has this as low edge

  const Cut kSelectRSNumuND0 = kNumuCutND2018 && kCVNFSProtonScore2018<0.49;
  const Cut kSelectRSNumuND1 = kNumuCutND2018 && !kCVNProngProtonCtNumu;
  const Cut kSelectRSNumuND2 = kNumuCutND2018 && kAntiNumuBDTCVN>0.11;

  const Cut kSelectWSNumuND0 = kNumuCutND2018 && kCVNFSProtonScore2018>=0.49;
  const Cut kSelectWSNumuND1 = kNumuCutND2018 && kCVNProngProtonCtNumu;
  const Cut kSelectWSNumuND2 = kNumuCutND2018 && kAntiNumuBDTCVN<=0.11;

  // Nue-style Binnings for BDT
  const Var kNueBDT3Bin([](const caf::SRProxy *sr)
                        {
                          if( kAntiNueBDTCVN(sr)>0. ) return float(0.5);
                          else return float(1.5);

                          return float(-1.5);
                        });

  const Var kNueBDT3AnaBinning([&](const caf::SRProxy *sr)
                               {
                                 int selBin = kNueBDT3Bin(sr);

                                 float nuE  = kNueEnergy2018(sr);
                                 int nuEBin = nuE/0.5;
                                 assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen");

                                 int anaBin = 9*selBin + nuEBin;
                                 return anaBin;
                               });

  const HistAxis kAxisFSProton("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueFSProtonAnaBinning);
  const HistAxis kAxisCVNProng("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueCVNProngAnaBinning);
  const HistAxis kAxisBDT3("NuE Energy / Analysis Bin",Binning::Simple(18,0,18),kNueBDT3AnaBinning);
  // -------------------------------------------------------------------------------------------------

  // Predictions/Spectra
  PredictionNoOsc pNumuFSProtonRS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectRSNumuND0,kNoShift,wei);
  PredictionNoOsc pNumuFSProtonWS(ld_numuRHC,kNumuCCOptimisedAxis,kSelectWSNumuND0,kNoShift,wei);
  PredictionNoOsc pNumuTot(ld_numuRHC,kNumuCCOptimisedAxis,kNumuCutND2018,kNoShift,wei);

  //--- Splitting into NC from nu and antinu using spectra for this analysis
  Spectrum sNumuFSProtonRSNC(ld_numuRHC,kNumuCCOptimisedAxis,
                               kSelectRSNumuND0 && kIsNC && !kIsAntiNu,kNoShift,wei);
  Spectrum sNumuFSProtonRSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
                                   kSelectRSNumuND0 && kIsNC && kIsAntiNu,kNoShift,wei);
  Spectrum sNumuFSProtonWSNC(ld_numuRHC,kNumuCCOptimisedAxis,
                               kSelectWSNumuND0 && kIsNC && !kIsAntiNu,kNoShift,wei);
  Spectrum sNumuFSProtonWSNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
                                   kSelectWSNumuND0 && kIsNC && kIsAntiNu,kNoShift,wei);

  Spectrum sNumuTotNC(ld_numuRHC,kNumuCCOptimisedAxis,
                        kNumuCutND2018 && kIsNC && !kIsAntiNu,kNoShift,wei);
  Spectrum sNumuTotNCAnti(ld_numuRHC,kNumuCCOptimisedAxis,
                            kNumuCutND2018 && kIsNC && kIsAntiNu,kNoShift,wei);

  // Press play
  ld_numuRHC.Go();

  // Save the predictions
  TFile *outFile = new TFile("./WrongSign-XCheck2019-NumuMC.root","RECREATE");
  outFile->cd();

  pNumuFSProtonRS.SaveTo(outFile, "pNumuFSProtonRS");
  pNumuFSProtonWS.SaveTo(outFile, "pNumuFSProtonWS");
  pNumuTot.SaveTo(outFile, "pNumuTot");

  sNumuFSProtonRSNC.SaveTo(outFile, "sNumuFSProtonRSNC");
  sNumuFSProtonRSNCAnti.SaveTo(outFile, "sNumuFSProtonRSNCAnti");
  sNumuFSProtonWSNC.SaveTo(outFile, "sNumuFSProtonWSNC");
  sNumuFSProtonWSNCAnti.SaveTo(outFile, "sNumuFSProtonWSNCAnti");
  sNumuTotNC.SaveTo(outFile, "sNumuTotNC");
  sNumuTotNCAnti.SaveTo(outFile, "sNumuTotNCAnti");

  outFile->Close();
}