genie_preds_make.C File Reference

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionNoOsc.h"
#include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "OscLib/IOscCalc.h"
#include "TFile.h"
#include "TString.h"
#include "3FlavorAna/Systs/GeniePCASyst.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Core/HistAxis.h"
#include "pca/genie_diag_utils.h"
#include <fstream>
#include <iostream>
#include <cmath>

Go to the source code of this file.

Functions
void	genie_preds_make (const int low, const int high)

Function Documentation

void genie_preds_make	(	const int	low,
		const int	high
	)

Definition at line 43 of file genie_preds_make.C.

References ana::assert(), ana::Energy, check_time_usage::float, ana::kNue2018FDPeripheral(), ana::kNueEnergy2018, caf::Proxy< caf::StandardRecord >::sel, and sr.

{
  
  const Var kNue2018CoreSelectionBin([](const caf::SRProxy *sr){ 
        bool isPeri = kNue2018FDPeripheral(sr);
        if (isPeri)
            return float(-5.0);
        if (sr->sel.cvnProd3Train.nueid > 0.75 && sr->sel.cvnProd3Train.nueid <= 0.97)
            return float(0.5);
        else if (sr->sel.cvnProd3Train.nueid > 0.97)
            return float(1.5);

        else return float(-5.0);
        });

  const Var kNue2018CoreAnaBin([&](const caf::SRProxy *sr){
        int selBin = kNue2018CoreSelectionBin(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 Binning kNue2018CoreBinning = Binning::Simple(18,0,18);
  const HistAxis kNue2018CoreAxis("NuE Energy / Analysis Bin",
                               kNue2018CoreBinning,kNue2018CoreAnaBin);
  const Cut kNueND = kNue2018NDCVNSsb;
  const Cut kNueFD = kNue2018FD;
  const HistAxis kNueAxis = kNue2018CoreAxis;
  const Cut kNumuND = kNumuCutND2018;
  const HistAxis kNumuAxis =  kNumuCCOptimisedAxis;
  const Cut kNumuFD = kNumuCutFD2018;
  
  const int N_PC = 10;
  std::vector<const ISyst*> systs;
  for(int i = 0; i < N_PC; i++)
    systs.push_back(GetGeniePrincipals2018Small(i));

  std::vector<Var> kPPFXFluxUnivWgt;
  
  Loaders loader_fd_fhc;
  loader_fd_fhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1",
                              caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
  loader_fd_fhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_full_v1",
                              caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
  loader_fd_fhc.SetLoaderPath("prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns",
                              caf::kNEARDET, Loaders::kData);
  loader_fd_fhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1",
                              caf::kNEARDET, Loaders::kMC);
  Loaders loader_fd_rhc;
  loader_fd_rhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1",
                              caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
  loader_fd_rhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.e_fd_genie_fluxswap_rhc_nova_v08_full_v1",
                              caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
  loader_fd_rhc.SetLoaderPath("prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns",
                              caf::kNEARDET, Loaders::kData);
  loader_fd_rhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1",
                              caf::kNEARDET, Loaders::kMC);

  const Var weight = kPPFXFluxCVWgt*kXSecCVWgt2018;
  
  std::map<std::string, PredictionInterp*> preds;
  auto calc = DefaultOscCalc();
  
  const NuePropExtrapGenerator kNuePropExtrap(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
  const NumuExtrapGenerator kNumuExtrap(kNumuAxis, kNumuFD, kNumuND, kNoShift, weight);
  PredictionInterp* pred_propnue_fhc = new PredictionInterp(systs, calc->Copy(), kNuePropExtrap, loader_fd_fhc);
  PredictionInterp* pred_propnue_rhc = new PredictionInterp(systs, calc->Copy(), kNuePropExtrap, loader_fd_rhc);
  PredictionInterp* pred_numu_fhc = new PredictionInterp(systs, calc->Copy(), kNumuExtrap, loader_fd_fhc);
  PredictionInterp* pred_numu_rhc = new PredictionInterp(systs, calc->Copy(), kNumuExtrap, loader_fd_rhc);
  
  preds.insert(std::make_pair("nue_prop_fhc", pred_propnue_fhc));
  preds.insert(std::make_pair("nue_prop_rhc", pred_propnue_rhc));
  preds.insert(std::make_pair("numu_fhc", pred_numu_fhc));
  preds.insert(std::make_pair("numu_rhc", pred_numu_rhc));

  auto allsysts = getAna2018SmallXsecSystsForPCA();
  const int totalsysts = allsysts.size();
  const int nuniverses = 1000;

  std::vector <std::pair <const ISyst*, SystMode> > systConfigs;
  for (int idx = 0; idx < totalsysts; ++idx)
    systConfigs.push_back({allsysts[idx], kSystGaussian});
  
  auto multiverse_shifts = GetSystShiftsMultiverse(systConfigs, nuniverses);
  assert(multiverse_shifts.size() == nuniverses);

  for(int UnivIdx = low; UnivIdx < high; UnivIdx++){
    const NuePropExtrapGenerator kNueUnivPropExtrap(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
    const NumuExtrapGenerator kNumuUnivExtrap(kNumuAxis, kNumuFD, kNumuND, kNoShift, weight);
    PredictionInterp* pred_propnue_univ_fhc = new PredictionInterp({}, calc->Copy(), kNueUnivPropExtrap, loader_fd_fhc,
                                                                   multiverse_shifts[UnivIdx]);
    PredictionInterp* pred_propnue_univ_rhc = new PredictionInterp({}, calc->Copy(), kNueUnivPropExtrap, loader_fd_rhc,
                                                                   multiverse_shifts[UnivIdx]);
    PredictionInterp* pred_numu_univ_fhc = new PredictionInterp({}, calc->Copy(), kNumuUnivExtrap, loader_fd_fhc,
                                                                multiverse_shifts[UnivIdx]);
    PredictionInterp* pred_numu_univ_rhc = new PredictionInterp({}, calc->Copy(), kNumuUnivExtrap, loader_fd_rhc,
                                                                multiverse_shifts[UnivIdx]);
    
    std::string predname_nue_prop_fhc = "nue_prop_univ_fhc_"+std::to_string(UnivIdx);
    preds.insert(std::make_pair(predname_nue_prop_fhc, pred_propnue_univ_fhc));
    std::string predname_nue_prop_rhc = "nue_prop_univ_rhc_"+std::to_string(UnivIdx);
    preds.insert(std::make_pair(predname_nue_prop_rhc, pred_propnue_univ_rhc));
    std::string predname_numu_fhc = "numu_univ_fhc_"+std::to_string(UnivIdx);
    preds.insert(std::make_pair(predname_numu_fhc, pred_numu_univ_fhc));
    std::string predname_numu_rhc = "numu_univ_rhc_"+std::to_string(UnivIdx);
    preds.insert(std::make_pair(predname_numu_rhc, pred_numu_univ_rhc));

  }

  loader_fd_fhc.SetSpillCut(kStandardSpillCuts);
  loader_fd_rhc.SetSpillCut(kStandardSpillCuts);

  loader_fd_fhc.Go();
  loader_fd_rhc.Go();

  const std::string filename = "genie_syst_pred.root";  
  TFile* file = new TFile(filename.c_str(), "recreate"); 
  std::for_each(preds.begin(), preds.end(),
                [&](std::pair<std::string, PredictionInterp*> pred){
                  pred.second->SaveTo(file, pred.first.c_str());
                });

}