genie_preds_make.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Loaders.h"
3 #include "CAFAna/Cuts/Cuts.h"
4 #include "CAFAna/Cuts/NueCutsSecondAna.h"
10 #include "CAFAna/Cuts/TruthCuts.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
17 #include "CAFAna/Systs/Systs.h"
18 #include "CAFAna/Systs/BeamSysts.h"
22 #include "OscLib/IOscCalc.h"
23 #include "TFile.h"
24 #include "TString.h"
25 
26 #include "CAFAna/Systs/Systs.h"
28 #include "CAFAna/Vars/Vars.h"
32 #include "CAFAna/Core/HistAxis.h"
33 
34 #include "pca/genie_diag_utils.h"
35 
36 #include <fstream>
37 #include <iostream>
38 #include <cmath>
39 
40 
41 using namespace ana;
42 //----------------------------------------------------------------------
43 void genie_preds_make(const int low, const int high)
44 {
45 
46  const Var kNue2018CoreSelectionBin([](const caf::SRProxy *sr){
47  bool isPeri = kNue2018FDPeripheral(sr);
48  if (isPeri)
49  return float(-5.0);
50  if (sr->sel.cvnProd3Train.nueid > 0.75 && sr->sel.cvnProd3Train.nueid <= 0.97)
51  return float(0.5);
52  else if (sr->sel.cvnProd3Train.nueid > 0.97)
53  return float(1.5);
54 
55  else return float(-5.0);
56  });
57 
58  const Var kNue2018CoreAnaBin([&](const caf::SRProxy *sr){
59  int selBin = kNue2018CoreSelectionBin(sr);
60 
61  float nuE = kNueEnergy2018(sr);
62  int nuEBin = nuE/0.5;
63 
64  assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen");
65 
66  int anaBin = 9*selBin + nuEBin;
67 
68  return anaBin;
69  });
70 
71  const Binning kNue2018CoreBinning = Binning::Simple(18,0,18);
72  const HistAxis kNue2018CoreAxis("NuE Energy / Analysis Bin",
73  kNue2018CoreBinning,kNue2018CoreAnaBin);
74  const Cut kNueND = kNue2018NDCVNSsb;
75  const Cut kNueFD = kNue2018FD;
76  const HistAxis kNueAxis = kNue2018CoreAxis;
77  const Cut kNumuND = kNumuCutND2018;
78  const HistAxis kNumuAxis = kNumuCCOptimisedAxis;
79  const Cut kNumuFD = kNumuCutFD2018;
80 
81  const int N_PC = 10;
82  std::vector<const ISyst*> systs;
83  for(int i = 0; i < N_PC; i++)
84  systs.push_back(GetGeniePrincipals2018Small(i));
85 
86  std::vector<Var> kPPFXFluxUnivWgt;
87 
88  Loaders loader_fd_fhc;
89  loader_fd_fhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1",
90  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
91  loader_fd_fhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_full_v1",
92  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
93  loader_fd_fhc.SetLoaderPath("prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns",
94  caf::kNEARDET, Loaders::kData);
95  loader_fd_fhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1",
96  caf::kNEARDET, Loaders::kMC);
97  Loaders loader_fd_rhc;
98  loader_fd_rhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1",
99  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
100  loader_fd_rhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.e_fd_genie_fluxswap_rhc_nova_v08_full_v1",
101  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
102  loader_fd_rhc.SetLoaderPath("prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns",
103  caf::kNEARDET, Loaders::kData);
104  loader_fd_rhc.SetLoaderPath("prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1",
105  caf::kNEARDET, Loaders::kMC);
106 
107  const Var weight = kPPFXFluxCVWgt*kXSecCVWgt2018;
108 
109  std::map<std::string, PredictionInterp*> preds;
110  auto calc = DefaultOscCalc();
111 
112  const NuePropExtrapGenerator kNuePropExtrap(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
113  const NumuExtrapGenerator kNumuExtrap(kNumuAxis, kNumuFD, kNumuND, kNoShift, weight);
114  PredictionInterp* pred_propnue_fhc = new PredictionInterp(systs, calc->Copy(), kNuePropExtrap, loader_fd_fhc);
115  PredictionInterp* pred_propnue_rhc = new PredictionInterp(systs, calc->Copy(), kNuePropExtrap, loader_fd_rhc);
116  PredictionInterp* pred_numu_fhc = new PredictionInterp(systs, calc->Copy(), kNumuExtrap, loader_fd_fhc);
117  PredictionInterp* pred_numu_rhc = new PredictionInterp(systs, calc->Copy(), kNumuExtrap, loader_fd_rhc);
118 
119  preds.insert(std::make_pair("nue_prop_fhc", pred_propnue_fhc));
120  preds.insert(std::make_pair("nue_prop_rhc", pred_propnue_rhc));
121  preds.insert(std::make_pair("numu_fhc", pred_numu_fhc));
122  preds.insert(std::make_pair("numu_rhc", pred_numu_rhc));
123 
124  auto allsysts = getAna2018SmallXsecSystsForPCA();
125  const int totalsysts = allsysts.size();
126  const int nuniverses = 1000;
127 
128  std::vector <std::pair <const ISyst*, SystMode> > systConfigs;
129  for (int idx = 0; idx < totalsysts; ++idx)
130  systConfigs.push_back({allsysts[idx], kSystGaussian});
131 
132  auto multiverse_shifts = GetSystShiftsMultiverse(systConfigs, nuniverses);
133  assert(multiverse_shifts.size() == nuniverses);
134 
135  for(int UnivIdx = low; UnivIdx < high; UnivIdx++){
136  const NuePropExtrapGenerator kNueUnivPropExtrap(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
137  const NumuExtrapGenerator kNumuUnivExtrap(kNumuAxis, kNumuFD, kNumuND, kNoShift, weight);
138  PredictionInterp* pred_propnue_univ_fhc = new PredictionInterp({}, calc->Copy(), kNueUnivPropExtrap, loader_fd_fhc,
139  multiverse_shifts[UnivIdx]);
140  PredictionInterp* pred_propnue_univ_rhc = new PredictionInterp({}, calc->Copy(), kNueUnivPropExtrap, loader_fd_rhc,
141  multiverse_shifts[UnivIdx]);
142  PredictionInterp* pred_numu_univ_fhc = new PredictionInterp({}, calc->Copy(), kNumuUnivExtrap, loader_fd_fhc,
143  multiverse_shifts[UnivIdx]);
144  PredictionInterp* pred_numu_univ_rhc = new PredictionInterp({}, calc->Copy(), kNumuUnivExtrap, loader_fd_rhc,
145  multiverse_shifts[UnivIdx]);
146 
147  std::string predname_nue_prop_fhc = "nue_prop_univ_fhc_"+std::to_string(UnivIdx);
148  preds.insert(std::make_pair(predname_nue_prop_fhc, pred_propnue_univ_fhc));
149  std::string predname_nue_prop_rhc = "nue_prop_univ_rhc_"+std::to_string(UnivIdx);
150  preds.insert(std::make_pair(predname_nue_prop_rhc, pred_propnue_univ_rhc));
151  std::string predname_numu_fhc = "numu_univ_fhc_"+std::to_string(UnivIdx);
152  preds.insert(std::make_pair(predname_numu_fhc, pred_numu_univ_fhc));
153  std::string predname_numu_rhc = "numu_univ_rhc_"+std::to_string(UnivIdx);
154  preds.insert(std::make_pair(predname_numu_rhc, pred_numu_univ_rhc));
155 
156  }
157 
158  loader_fd_fhc.SetSpillCut(kStandardSpillCuts);
159  loader_fd_rhc.SetSpillCut(kStandardSpillCuts);
160 
161  loader_fd_fhc.Go();
162  loader_fd_rhc.Go();
163 
164  const std::string filename = "genie_syst_pred.root";
165  TFile* file = new TFile(filename.c_str(), "recreate");
166  std::for_each(preds.begin(), preds.end(),
167  [&](std::pair<std::string, PredictionInterp*> pred){
168  pred.second->SaveTo(file, pred.first.c_str());
169  });
170 
171 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Var kNueEnergy2018([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC(sr);else return kNueEnergyFHC(sr);})
Definition: NueEnergy2018.h:25
const Cut Energy
caf::StandardRecord * sr
void genie_preds_make(const std::string ana="nue", const int low=0, const int high=1000)
const Cut kNue2018FDPeripheral(kNue2018FDPeripheralFunc)
assert(nhit_max >=nhit_nbins)
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141