genie_preds_make.C
Go to the documentation of this file.
1 #include "CAFAna/Analysis/Calcs.h"
2 #include "3FlavorAna/Systs/3FlavorAna2020Systs.h"
3 #include "CAFAna/Core/Loaders.h"
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Cuts/NueCutsSecondAna.h"
6 #include "3FlavorAna/Cuts/NueCuts2017.h"
7 #include "3FlavorAna/Cuts/NueCuts2018.h"
8 #include "3FlavorAna/Cuts/NueCuts2020.h"
9 #include "3FlavorAna/Cuts/NumuCuts.h"
10 #include "3FlavorAna/Cuts/NumuCuts2017.h"
11 #include "3FlavorAna/Cuts/NumuCuts2018.h"
12 #include "3FlavorAna/Cuts/NumuCuts2020.h"
13 #include "CAFAna/Cuts/TruthCuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/Prediction/IPrediction.h"
16 #include "CAFAna/Prediction/PredictionNoExtrap.h"
17 #include "CAFAna/Prediction/PredictionNoOsc.h"
18 #include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
19 #include "CAFAna/Prediction/PredictionInterp.h"
20 #include "CAFAna/Systs/Systs.h"
21 #include "CAFAna/Systs/BeamSysts.h"
22 #include "CAFAna/Vars/GenieWeights.h"
23 #include "CAFAna/Vars/PPFXWeights.h"
24 #include "3FlavorAna/Vars/HistAxes.h"
25 #include "OscLib/IOscCalc.h"
26 #include "TFile.h"
27 #include "TString.h"
28 
29 #include "CAFAna/Systs/Systs.h"
30 #include "3FlavorAna/Systs/GeniePCASyst.h"
31 #include "CAFAna/Vars/Vars.h"
32 #include "3FlavorAna/Vars/NueVarsExtra.h"
33 #include "3FlavorAna/Vars/NueVars.h"
34 #include "3FlavorAna/Vars/NumuVars.h"
35 #include "CAFAna/Core/HistAxis.h"
36 
37 #include "3FlavorAna/Ana2020/Systs/genie/genie_diag_utils.h"
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",
79  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
80  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1_numu2020",
81  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
82  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_numu2020",
83  caf::kNEARDET, Loaders::kData);
84  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_numu2020",
85  caf::kNEARDET, Loaders::kMC);
86 
87  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_numu2020",
88  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
89  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1_numu2020",
90  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
91  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_numu2020",
92  caf::kNEARDET, Loaders::kData);
93  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_numu2020",
94  caf::kNEARDET, Loaders::kMC);
95  } else {
96  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_nue2020",
97  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
98  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_fhc_nova_v08_full_v1_nue2020",
99  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
100  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_nue2020",
101  caf::kNEARDET, Loaders::kData);
102  loader_fd_fhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_nue2020",
103  caf::kNEARDET, Loaders::kMC);
104 
105  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_nue2020",
106  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kNonSwap);
107  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_fluxswap_rhc_nova_v08_full_v1_nue2020",
108  caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);
109  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_nue2020",
110  caf::kNEARDET, Loaders::kData);
111  loader_fd_rhc.SetLoaderPath("prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_nue2020",
112  caf::kNEARDET, Loaders::kMC);
113  }
114 
115  const Var weight = kPPFXFluxCVWgt*kXSecCVWgt2020;
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 }
ana::get3FlavorAna2020SmallXsecSysts
std::vector< const ISyst * > get3FlavorAna2020SmallXsecSysts(const EAnaType2020 ana)
Definition: 3FlavorAna2020Systs.cxx:281
caf::kNEARDET
Near Detector underground.
Definition: SREnums.h:10
plot_validation_datamc.pred
pred
Definition: plot_validation_datamc.py:175
ana::PredictionInterp
Implements systematic errors by interpolation between shifted templates.
Definition: PredictionInterp.h:23
caf::kFARDET
Far Detector at Ash River.
Definition: SREnums.h:11
ana::Binning
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
systs
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
ana::Loaders
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
ana::kNumu2020ND
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
shutoffs.filename
string filename
Definition: shutoffs.py:106
ana::kNumuFD
const Cut kNumuFD
Definition: NumuCuts.h:53
make_pair
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
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::DefaultOscCalc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
PandAna.cut.analysis_cuts.kNueFD
kNueFD
Definition: analysis_cuts.py:110
ana::kNue2020ND
const Cut kNue2020ND
Definition: NueCuts2020.h:178
compare_h5_caf.idx
int idx
Definition: compare_h5_caf.py:69
ana::NumuExtrapGenerator
Generates extrapolated Numu predictions.
Definition: PredictionGenerator3Flavor.h:20
ana::Loaders::Go
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
ana::Loaders::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
ana::NuePropExtrapGenerator
Generates extrapolated Nue predictions using ProportionalDecomp.
Definition: PredictionGenerator3Flavor.h:48
check_time_usage.float
float
Definition: check_time_usage.py:6
ana::kNumuND
const Cut kNumuND
Definition: NumuCuts.h:55
osc::OscCalcDumb::Copy
virtual IOscCalc * Copy() const
Definition: OscCalcDumb.h:19
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::GetSystShiftsMultiverse
std::vector< SystShifts > GetSystShiftsMultiverse(const std::vector< std::pair< const ISyst *, SystMode > > &systConfigs, const int nUniverses, const int seed)
Definition: SystShifts.cxx:282
genie_preds_make
void genie_preds_make(const std::string ana="nue", const int low=0, const int high=1000)
Definition: genie_preds_make.C:46
ana::kNumu2020FD
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
ana::kNue2020FD
const Cut kNue2020FD
Definition: NueCuts2020.h:65
ana::kNumuCCOptimisedAxis
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
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:21
genie_diag_utils.h
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ana::_HistAxis< Var >
ana::kNue2020AnaBin
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
file
TFile * file
Definition: cellShifts.C:17
ana::assert
assert(nhit_max >=nhit_nbins)
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
ana::_Cut< caf::SRProxy >
ana::_Var< caf::SRProxy >
nuniverses
const int nuniverses
Definition: genie_syst_pca.C:29
ana::Loaders::SetLoaderPath
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
ana::kNue2020FDPeripheral
const Cut kNue2020FDPeripheral(kNue2020FDPeripheralFunc)
ana::Binning::Simple
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
ana::kXSecCVWgt2020
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106
ana::GetGeniePrincipals2020Small
GeniePCASyst * GetGeniePrincipals2020Small(int PCIdx)
Definition: GeniePCASyst.cxx:41
Generated on Fri Oct 23 2020 05:11:28 for NOvA by   doxygen 1.8.11