make_RHC_WrongSign_Numu_Pd7dData.C

Go to the documentation of this file.

// ----------------------------------------------------------------
// Make the Predictions at ND for numu, Data 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 "3FlavorAna/Vars/HistAxes.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_Pd7dData()
{
  // ND MC data set, Prod 4 RHC
  // nue
  std::string fname_numuRHC = "prod_caf_R17-11-14-prod4reco.k_nd_numi_rhc_epoch7d_v1_goodruns";
  SpectrumLoader ld_dataRHCNumu(fname_numuRHC);
  ld_dataRHCNumu.SetSpillCut(kStandardSpillCuts);

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

  // ------------------------------------------------------------------------------------------------
  // Numu DATA
  Spectrum sQ1NumuData_FSProtonRS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectRSNumuND0,kNoShift,kUnweighted);
  Spectrum sQ1NumuData_FSProtonWS(ld_dataRHCNumu,kNumuCCOptimisedAxis,kSelectWSNumuND0,kNoShift,kUnweighted);

  Spectrum sQ1NumuDataTot(ld_dataRHCNumu,kNumuCCOptimisedAxis,kNumuCutND2018,kNoShift,kUnweighted);
  // ------------------------------------------------------------------------------------------------

  // Press play
  ld_dataRHCNumu.Go();

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

  sQ1NumuData_FSProtonRS.SaveTo(outFile, "sQ1NumuData_FSProtonRS");
  sQ1NumuData_FSProtonWS.SaveTo(outFile, "sQ1NumuData_FSProtonWS");

  sQ1NumuDataTot.SaveTo(outFile, "sQ1NumuDataTot");

  outFile->Close();
}