cafe_FD_predictions.C

Go to the documentation of this file.

// This macro creates unextrapolated predictions at the far detector for the
// nue selections (low, high, peripheral) and numu selections (Q1-Q4).  Both
// unoscillated and oscillated predictions are created and written into a ROOT file.
//
// The predictions from CAF are loaded up from SAM datasets to avoid having weights
// from a particular GENIE version baked into them.

#include <string>
#include <vector>
#include <set>

#include "TFile.h"
#include "TH2.h"

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/numu/Analysis2018/includes.h"
#include "DataQuality/TimePeak/Plots/Includes.h"
#include "DataQuality/TimePeak/Plots/VarsAndCuts.h"
#include "CAFAna/Vars/NueEnergy2018.h"
#include "NuXAna/Vars/NusVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"

using namespace ana;

// convenience enumeration for each selection type
// The CAFAna prediction files put all three nue selections
// into a single histogram
typedef enum SelType{
  lNue = 0,
  lNumuQ1,
  lNumuQ2,
  lNumuQ3,
  lNumuQ4
} SelType_t;

std::vector<std::string> gSelTypeNames({"NuE", "NuMuQ1", "NuMuQ2", "NuMuQ3", "NuMuQ4"});

//-----------------------------------------------------------------------------
// convenience struct to hold information about the selection, horn current,
// and prediction file name.
struct PredictionInfo{

  PredictionInfo(std::string    const& hornCur,
                 SelType_t      const& selType)
      : curName  (hornCur)
      , selection(selType)
  {
    fhc = (curName.find("FHC") != std::string::npos) ? true : false;
    std::string recolet("d");
    std::string lowerCur("fhc");
    std::string lowerSel("numu");

    if(fhc){
      pot = kAna2019FHCPOT;
    }
    else{
      recolet = "e";
      lowerCur = "rhc";
      pot = kAna2019RHCPOT;
    }

    // define cuts and vars we will need
    Cut     cuts   = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill};
    Var     energy = kCCE;
    Binning bins   = kNumuCCEOptimisedBinning;

    if(selection == lNue){
      bins     = kNue2018Binning;
      cuts     = {kNue2018FDAllSamples && kHasNeutrino && kInBeamSpill};
      energy   = kNue2018AnaBin;
      lowerSel = "nue";
    }
    else{
      std::string quantFile = "/cvmfs/nova.opensciencegrid.org/externals/numudata/v00.00/NULL/ana2018/Quantiles/";
      quantFile += (fhc) ? "quantiles__fhc_full__numu2018.root" : "quantiles__rhc_full__numu2018.root";
      TFile qtf(quantFile.c_str(), "READ");
      TH2F *FDSpec2D = dynamic_cast<TH2F*>(qtf.Get("FDSpec2D"));
      std::vector<Cut> numuQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);

      if(selection == lNumuQ1){
        cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[0]};
      }
      else if(selection == lNumuQ2){
        cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[1]};
      }
      else if(selection == lNumuQ3){
        cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[2]};
      }
      else if(selection == lNumuQ4){
        cuts = {kNumuCutFD2018 && kHasNeutrino && kInBeamSpill && numuQuantCuts[3]};
      }
    }

    // create Spectrum loaders from the concat datasets
    // TODO: no cosmic backgrounds or rock backgrounds can be loaded in this way
    std::string base = "prod_sumdecaf_R17-11-14-prod4reco." + recolet + "_fd_genie_";
    SpectrumLoader nonSwapLoader (base + "nonswap_"  + lowerCur + "_nova_v08_full_v1_" + lowerSel + "2018");
    SpectrumLoader fluxSwapLoader(base + "fluxswap_" + lowerCur + "_nova_v08_full_v1_" + lowerSel + "2018");
    SpectrumLoader tauSwapLoader (base + "tau_"      + lowerCur + "_nova_v08_full_v1_" + lowerSel + "2018");

    // now make the no extrapolation prediction thing
    cafPred = new PredictionNoExtrap(nonSwapLoader,
                                     fluxSwapLoader,
                                     tauSwapLoader,
                                     "E_{reco}  (GeV)",
                                     bins,
                                     energy,
                                     cuts,
                                     kNoShift,
                                     kXSecCVWgt2018 * kPPFXFluxCVWgt);

    // apparently we make the loaders do stuff after we give them to the prediction...
    nonSwapLoader .Go();
    fluxSwapLoader.Go();
    tauSwapLoader .Go();
  }

  //---------------------------------------------------------------
  std::string HistName(bool osc)
  {
    std::string totName = (gSelTypeNames[selection] + curName);
    totName += (osc) ? "Oscillated" : "Unoscillated";
    return totName;
  }

  //---------------------------------------------------------------
  TH1D* SpectrumFromPrediction(bool                     osc,
                               osc::IOscCalcAdjustable* calc)
  {
    TH1D* hist = nullptr;
    if(osc){
      hist = cafPred->Predict(calc).ToTH1(pot);
    }
    else
      hist = cafPred->PredictUnoscillated().ToTH1(pot);

    // set these histograms free from the input file by calling
    // SetDirectory(0)
    hist->SetDirectory(0);

    // the bool indicates if the histograms are oscillated or not
    hist->SetName((HistName(osc)).c_str());

    return hist;
  }

  //---------------------------------------------------------------

  bool                fhc;       ///< true for FHC, false for RHC
  double              pot;       ///< total POT for horn current
  SelType_t           selection; ///< numu + quantile or nue
  std::string         curName;   ///< string horn current value
  PredictionNoExtrap* cafPred;   ///< CAFAna prediction object
};

//-----------------------------------------------------------------------------
// Set the oscillation parameters.  Perhaps we want to set these on the command
// line at some point
void SetOscillationParameters(osc::IOscCalcAdjustable* calc)
{
  calc->SetdCP(0.17);
  calc->SetTh13(0.145426);
  calc->SetTh23(0.865744);
  calc->SetDmsq32(2.51e-3);
}

//-----------------------------------------------------------------------------
// main function to create the prediction histograms
// bool fhc: true if we want FHC events
// int  sel: the selection type according to SelType_t defined above
void cafe_FD_predictions(bool fhc,
                         int  sel)
{
  std::string curStr = (fhc) ? "FHC" : "RHC";

  // prediction information for the desired selection type
  // and horn current.
  PredictionInfo pi(curStr, (SelType_t)sel);

  // Set up the oscillation calculator
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  SetOscillationParameters(calc);

  // open the output file.
  std::string outFileName("caf_fd_predictions_");
  outFileName += (fhc) ? "fhc_" : "rhc_";
  outFileName += gSelTypeNames[sel];
  outFileName += ".root";
  TFile* outfile = new TFile(outFileName.c_str(), "RECREATE");

  // get and write prediction histograms for the oscillated and
  // unoscillated cases

  auto oscHist = pi.SpectrumFromPrediction(true, calc);
  oscHist->SetDirectory(outfile);
  oscHist->Write();

  auto noOscHist = pi.SpectrumFromPrediction(false, calc);
  noOscHist->SetDirectory(outfile);
  noOscHist->Write();

  // write and close the output file.
  outfile->Write();
  outfile->Close();
}