NumuCosmic.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"

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

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

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "TFile.h"

#include <iostream>
#include <cmath>
using namespace ana;

void NumuCosmic(bool isFHC )
{
  // What do my arguments mean?
  std::string beam = "rhc";
  if (isFHC)  beam = "fhc";
  
  // Set my axis.
  auto axis_numu = kNumuCCOptimisedAxis;
  
  // Set the default cuts to the combined quantile cuts for the Full FOM
  Cut cut_numufd = kNumu2020;

  // Set the string for the cosmic definitions.
  std::string dataset = "prod_caf_R19-11-18-prod5reco.i_fd_cosmic_"+beam+"_full_v1_goodruns";
  SpectrumLoader loaders(dataset);
  loaders.SetSpillCut(kStandardSpillCuts);

  //std::string infile_quant = (std::string)std::getenv("NUMUDATA_DIR")+"/lib/ana2020/Quantiles/quantiles_"+beam+"_full_numu2020.root";
  std::string infile_quant = "/cvmfs/nova.opensciencegrid.org/externals/numudata/v00.03/NULL/lib/ana2020/Quantiles/quantiles_"+beam+"_full_numu2020.root";
  TFile* infile = new TFile(infile_quant.c_str());

  TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D_LoosePTP");
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, axis_numu, kHadEFracAxis, 4);
  infile->Close();
  delete infile;
  
  // Make a vector of spectra.
  std::vector<Spectrum*> pred_cosmic_numu;

  // First, push back a combined quantile one.
  pred_cosmic_numu.push_back(new Spectrum(loaders, axis_numu, cut_numufd && kInCosmicTimingWindow,  kNoShift));
  // Then, loop through my Quantile cuts.
  for(int quant=0; quant<4; ++quant) {
    
    Cut ThisCut = cut_numufd && kInCosmicTimingWindow && HadEFracQuantCuts[quant];

    pred_cosmic_numu.push_back(new Spectrum(loaders, axis_numu, ThisCut,  kNoShift));
  }

  loaders.Go();

  TFile* outFile = new TFile(("fd_numu_cosmics_"+beam+".root").c_str(),"RECREATE");
  outFile->cd();
  for(int quant=0; quant<5; quant++){
      pred_cosmic_numu[quant]->SaveTo(outFile, Form("numu_cos_q%d", quant));
  }
  outFile->Close();
}