genie_preds_make.C
Go to the documentation of this file.
3 #include "CAFAna/Core/Loaders.h"
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Cuts/NueCutsSecondAna.h"
13 #include "CAFAna/Cuts/TruthCuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Systs/Systs.h"
21 #include "CAFAna/Systs/BeamSysts.h"
25 #include "OscLib/IOscCalc.h"
26 #include "TFile.h"
27 #include "TString.h"
28 
29 #include "CAFAna/Systs/Systs.h"
31 #include "CAFAna/Vars/Vars.h"
35 #include "CAFAna/Core/HistAxis.h"
36 
38 
39 #include <fstream>
40 #include <iostream>
41 #include <cmath>
42 
43 
44 using namespace ana;
45 //----------------------------------------------------------------------
46 void genie_preds_make(const std::string ana = "nue", const int low = 0, const int high = 1000)
47 {
48 
49  const Var kNue2020CoreAnaBin([](const caf::SRProxy *sr){
50  bool isPeri = kNue2020FDPeripheral(sr);
51  if (isPeri)
52  return float(-5.0);
53  else return float(kNue2020AnaBin(sr));
54  });
55 
56  const Binning kNue2020CoreBinning = Binning::Simple(18,0,18);
57  const HistAxis kNue2020CoreAxis("NuE Energy / Analysis Bin",
58  kNue2020CoreBinning,kNue2020CoreAnaBin);
59  const Cut kNueND = kNue2020ND;
60  const Cut kNueFD = kNue2020FD;
61  const HistAxis kNueAxis = kNue2020CoreAxis;
62  const Cut kNumuND = kNumu2020ND;
63  const HistAxis kNumuAxis = kNumuCCOptimisedAxis;
64  const Cut kNumuFD = kNumu2020FD;
65 
66  const int N_PC = 20;
67  std::vector<const ISyst*> systs;
68  for(int i = 0; i < N_PC; i++)
69  systs.push_back(GetGeniePrincipals2020Small(i));
70 
71  std::vector<Var> kPPFXFluxUnivWgt;
72 
73  Loaders loader_fd_fhc;
74  Loaders loader_fd_rhc;
75 
76  if(ana == "numu")
77  {
78  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_numu2020",
80  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1_numu2020",
82  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_numu2020",
84  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_numu2020",
86 
87  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_numu2020",
89  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1_numu2020",
91  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_numu2020",
93  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_numu2020",
95  } else {
96  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_nue2020",
98  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1_nue2020",
100  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_nue2020",
102  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_nue2020",
104 
105  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_nue2020",
107  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1_nue2020",
109  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_nue2020",
111  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_nue2020",
113  }
114 
116 
117  std::map<std::string, PredictionInterp*> preds;
118  auto calc = DefaultOscCalc();
119 
120  if(low<0)
121  {
122  if(ana == "numu")
123  {
124  const NumuExtrapGenerator kNumuExtrap(kNumuAxis, kNumuFD, kNumuND, kNoShift, weight);
125 
126  PredictionInterp* pred_fhc = new PredictionInterp(systs, calc->Copy(), kNumuExtrap, loader_fd_fhc);
127  PredictionInterp* pred_rhc = new PredictionInterp(systs, calc->Copy(), kNumuExtrap, loader_fd_rhc);
128 
129  preds.insert(std::make_pair("numu_fhc", pred_fhc));
130  preds.insert(std::make_pair("numu_rhc", pred_rhc));
131 
132  } else {
133  const NuePropExtrapGenerator kNuePropExtrapFHC(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
134  const NuePropExtrapRHCGenerator kNuePropExtrapRHC(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
135 
136  PredictionInterp* pred_fhc = new PredictionInterp(systs, calc->Copy(), kNuePropExtrapFHC, loader_fd_fhc);
137  PredictionInterp* pred_rhc = new PredictionInterp(systs, calc->Copy(), kNuePropExtrapRHC, loader_fd_rhc);
138 
139  preds.insert(std::make_pair("nue_prop_fhc", pred_fhc));
140  preds.insert(std::make_pair("nue_prop_rhc", pred_rhc));
141  }
142  }
143 
144  auto allsysts = get3FlavorAna2020SmallXsecSysts();
145  const int totalsysts = allsysts.size();
146  const int nuniverses = 1000;
147 
148  std::vector <std::pair <const ISyst*, SystMode> > systConfigs;
149  for (int idx = 0; idx < totalsysts; ++idx)
150  systConfigs.push_back({allsysts[idx], kSystGaussian});
151 
152  auto multiverse_shifts = GetSystShiftsMultiverse(systConfigs, nuniverses);
153  assert(multiverse_shifts.size() == nuniverses);
154 
155  for(int UnivIdx = low; UnivIdx < high; UnivIdx++){
156  if(ana == "numu") {
157  const NumuExtrapGenerator kNumuUnivExtrap(kNumuAxis, kNumuFD, kNumuND, kNoShift, weight);
158 
159  PredictionInterp* pred_univ_fhc = new PredictionInterp({}, calc->Copy(), kNumuUnivExtrap, loader_fd_fhc, multiverse_shifts[UnivIdx]);
160  PredictionInterp* pred_univ_rhc = new PredictionInterp({}, calc->Copy(), kNumuUnivExtrap, loader_fd_rhc, multiverse_shifts[UnivIdx]);
161 
162  std::string predname_fhc = "numu_univ_fhc_"+std::to_string(UnivIdx);
163  preds.insert(std::make_pair(predname_fhc, pred_univ_fhc));
164  std::string predname_rhc = "numu_univ_rhc_"+std::to_string(UnivIdx);
165  preds.insert(std::make_pair(predname_rhc, pred_univ_rhc));
166  } else {
167  const NuePropExtrapGenerator kNueUnivPropExtrapFHC(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
168  const NuePropExtrapRHCGenerator kNueUnivPropExtrapRHC(kNueAxis, kNumuAxis, kNueFD, kNueND, kNumuND, kNoShift, weight);
169 
170  PredictionInterp* pred_univ_fhc = new PredictionInterp({}, calc->Copy(), kNueUnivPropExtrapFHC, loader_fd_fhc, multiverse_shifts[UnivIdx]);
171  PredictionInterp* pred_univ_rhc = new PredictionInterp({}, calc->Copy(), kNueUnivPropExtrapRHC, loader_fd_rhc, multiverse_shifts[UnivIdx]);
172 
173  std::string predname_fhc = "nue_prop_univ_fhc_"+std::to_string(UnivIdx);
174  preds.insert(std::make_pair(predname_fhc, pred_univ_fhc));
175  std::string predname_rhc = "nue_prop_univ_rhc_"+std::to_string(UnivIdx);
176  preds.insert(std::make_pair(predname_rhc, pred_univ_rhc));
177  }
178  }
179 
180  loader_fd_fhc.SetSpillCut(kStandardSpillCuts);
181  loader_fd_rhc.SetSpillCut(kStandardSpillCuts);
182 
183  loader_fd_fhc.Go();
184  loader_fd_rhc.Go();
185 
186  std::string suffix;
187  if(low<0)
188  suffix = "PCs";
189  else
190  suffix = std::to_string(low)+"_"+std::to_string(high);
191 
192  const std::string filename = "genie_syst_pred_"+ana+"_"+suffix+".root";
193  TFile* file = new TFile(filename.c_str(), "recreate");
194  std::for_each(preds.begin(), preds.end(),
195  [&](std::pair<std::string, PredictionInterp*> pred){
196  pred.second->SaveTo(file, pred.first.c_str());
197  });
198 }
std::vector< const ISyst * > get3FlavorAna2020SmallXsecSysts(const EAnaType2020 ana)
Near Detector underground.
Definition: SREnums.h:10
Implements systematic errors by interpolation between shifted templates.
Far Detector at Ash River.
Definition: SREnums.h:11
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
string filename
Definition: shutoffs.py:106
const Cut kNumuFD
Definition: NumuCuts.h:53
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const Cut kNue2020ND
Definition: NueCuts2020.h:178
Generates extrapolated Numu predictions.
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
Generates extrapolated Nue predictions using ProportionalDecomp.
const Cut kNumuND
Definition: NumuCuts.h:55
virtual IOscCalc * Copy() const
Definition: OscCalcDumb.h:19
caf::StandardRecord * sr
std::vector< SystShifts > GetSystShiftsMultiverse(const std::vector< std::pair< const ISyst *, SystMode > > &systConfigs, const int nUniverses, const int seed)
Definition: SystShifts.cxx:282
void genie_preds_make(const std::string ana="nue", const int low=0, const int high=1000)
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
const Cut kNue2020FD
Definition: NueCuts2020.h:65
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.cxx:21
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Var kNue2020AnaBin([](const caf::SRProxy *sr){int selBin=kNue2020SelectionBin(sr);float nuE=kNueEnergy2020(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;})
Use this Analysis Binning for Ana2020, official Binning.
Definition: NueCuts2020.h:191
TFile * file
Definition: cellShifts.C:17
assert(nhit_max >=nhit_nbins)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const int nuniverses
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
const Cut kNue2020FDPeripheral(kNue2020FDPeripheralFunc)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106
GeniePCASyst * GetGeniePrincipals2020Small(int PCIdx)