make_predictions_systs.C
Go to the documentation of this file.
1 //analysis will pick up the list of variables and cuts, and extrapoltion types
2 //option will pick up dataset (e.g. concats, syst shifted), and list of systematics
3 
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/Systs/Systs.h"
17 #include "CAFAna/Systs/BeamSysts.h"
20 #include "CAFAna/Vars/HistAxes.h"
22 
23 #include <TFile.h>
24 #include <TString.h>
25 #include <vector>
26 #include <iostream>
27 
28 // TODO: change name of this header now
29 #include "CAFAna/shared/Ana2018/syst_variations.h"
30 
31 using namespace ana;
32 
34  const TString option, const TString analysis="allxp_nue",
35  const TString period="full", const TString outdir="./")
36 {
37 
38  // Get flux (FHC or RHC) from option
39  if (!(option.Contains("FHC") || option.Contains("RHC")))
40  {
41  std:: cerr << "Please mention FHC or RHC in option" << std::endl;
42  abort(); // unable to determine if FHC or RHC
43  }
44  auto flux = option.Contains("RHC") ? Loaders::kRHC : Loaders::kFHC;
45 
46  // Get Loaders
47  auto loaders = GetLoaders2018(option, period);
48  if (option.Contains("fakeNDData")) // swap ND data loaders to fake data
50  loaders->SetSpillCut(kStandardSpillCuts);
51 
52  // Get shifts
53  auto shifts = GetShifts2018(option);
54 
55  auto xpdefs = GetExtrapolationDefs(analysis, period, option);
56 
57  // Extrapolation and decomposition
58  // No extrap
59  bool run_nxp = analysis.Contains("nxp") || analysis.Contains("allxp");
60  // Proportional decomposition
61  bool run_xp_prop = analysis.Contains("xp_prop") ||
62  analysis.Contains("allxp_nue");
63  // Combo --> FHC only
64  bool run_xp_combo = (analysis.Contains("xp_combo") ||
65  analysis.Contains("allxp_nue")) &&
67  // Numu
68  bool run_xp_numu = analysis.Contains("xp_numu") ||
69  analysis.Contains ("allxp_numu");
70  uint NumExtrap = (uint)run_nxp + (uint)run_xp_prop +
71  (uint)run_xp_combo + (uint)run_xp_numu;
72 
73  std::cout << "NumExtrap: " << NumExtrap << std::endl;
74  if(analysis.Contains("allxp_nue") && analysis.Contains ("allxp_numu"))
75  {
76  std::cerr << "All nue and numu xp? Ambitious much?" << std::endl;
77  return;
78  }
79 
80  std::cout << analysis << std :: endl;
81 
83 
84  struct PredDef
85  {
86  const IPrediction * pred;
87  const TString xp_type;
88  const TString syst_name;
89  const TString sigma_name;
90  const TString edef_name;
91  };
92 
93  struct GenDef
94  {
95  const IPredictionGenerator *gen;
96  const TString xp_type;
97  const TString edef_name;
98  };
99 
100  std::vector <GenDef> gens;
101  std::vector < std::vector <PredDef> > preds;
102 
103  for(auto edef:xpdefs)
104  {
105  if(run_nxp)
106  gens.push_back(
107  {new NoExtrapGenerator(edef.axis, edef.fdcut, weight),
108  "pred_nxp", edef.name});
109  if(run_xp_prop)
110  {
111  if (flux == Loaders::kRHC) // use RHC extrap generator
112  {
113  gens.push_back({
115  edef.axis, axisNumuForNueSig, edef.fdcut,
116  edef.nuendcut, edef.numundcut, kNoShift, weight),
117  "pred_xp_prop", edef.name});
118  }
119  else
120  {
121  gens.push_back({
123  edef.axis, axisNumuForNueSig, edef.fdcut,
124  edef.nuendcut, edef.numundcut, kNoShift, weight),
125  "pred_xp_prop", edef.name});
126  }
127  }
128  if(run_xp_combo)
129  gens.push_back({
131  edef.axis, axisNumuForNueSig, edef.fdcut,
132  edef.nuendcut, edef.numundcut, kNoShift, weight),
133  "pred_xp_combo", edef.name});
134 
135  if(run_xp_numu)
136  gens.push_back({
138  edef.axis, edef.fdcut,
139  edef.numundcut, kNoShift, weight),
140  "pred_xp_numu", edef.name});
141  }//loop over xpdefs
142 
143  for (auto gen:gens)
144  {
145  std::vector <PredDef> temp_preds;
146  for (auto shift:shifts)
147  {
148  auto pred = gen.gen->Generate(*loaders, shift.shift);
149  temp_preds.push_back({pred.release(), gen.xp_type, shift.syst_name,
150  shift.sigma_name, gen.edef_name});
151  }//end shifts
152  preds.emplace_back(temp_preds);
153  }//end var/cut combos
154 
155  loaders->Go();
156 
157  //mkdir is so dumb
158  //hadd is the worst. had to split predictions a little
159  for(uint gIdx=0; gIdx < NumExtrap; ++gIdx)
160  {
161  // Configure what the tier is ==> Combine plotted variables and quantiles (if numu)
162  TString OutTier = gens[gIdx].xp_type;
163  std::cout << "\t\tOutTier is now " << OutTier << std::endl;
164  for (uint xIdx=0; xIdx < xpdefs.size(); ++xIdx)
165  {
166  int SInd = gIdx + ( xIdx * NumExtrap );
167  TString VarTier = gens[SInd].edef_name;
168  int VPos = gens[SInd].edef_name.Index('_');
169  if (VPos > 0)
170  VarTier = VarTier(0,VPos);
171  if (!OutTier.Contains(VarTier))
172  OutTier = OutTier + "_" + VarTier;
173  std::cout << "\t\tOutTier is now " << OutTier << std::endl;
174  }
175  if (analysis.Contains("numu")) OutTier = OutTier + "_AllQuantiles";
176  TString filename = (
177  outdir + OutTier + "_" + period + "_" + option + ".root");
178  std::cout << "\nMade a file called " << filename << std::endl;
179 
180  auto file = new TFile (filename,"recreate");
181  // Now loop through all of the things which want to go into this file...
182  for (uint xIdx=0; xIdx < xpdefs.size(); ++xIdx)
183  {
184  int SInd = gIdx + ( xIdx * NumExtrap );
185  for(uint sIdx = 0; sIdx < shifts.size(); ++sIdx)
186  {
187  if(preds[SInd][sIdx].syst_name != shifts[sIdx].syst_name ||
188  preds[SInd][sIdx].xp_type != gens[SInd].xp_type)
189  {
190  std::cerr << "ERROR: Loops - you're doing it wrong \n"
191  << "Expected gen" << SInd << " " << gens[SInd].xp_type
192  << " and syst " << sIdx << " " << shifts[sIdx].syst_name
193  << "; but got pred " << preds[SInd][sIdx].xp_type << " "
194  << preds[SInd][sIdx].syst_name << std::endl;
195  return;
196  }//end check
197  auto dir = file->GetDirectory(
198  preds[SInd][sIdx].xp_type+"_"+preds[SInd][sIdx].syst_name);
199  if(!dir)
200  dir = file->mkdir(
201  preds[SInd][sIdx].xp_type+"_"+preds[SInd][sIdx].syst_name);
202 
203  auto dir2 = dir->GetDirectory(preds[SInd][sIdx].sigma_name);
204  if(!dir2) dir2 = dir->mkdir(preds[SInd][sIdx].sigma_name);
205  TString prefix = (analysis.Contains("nue")?"nue_pred":"numu_pred");
206  preds[SInd][sIdx].pred->SaveTo(
207  dir2->mkdir(prefix + "_" + preds[SInd][sIdx].edef_name));
208  }//loop over shifts
209  }//loop over xpdefs
210  delete file;
211  std::cout << "Saved predictions to " << filename
212  << " (" << gIdx+1 << "/" << NumExtrap << ")" << std::endl;
213  }//loop over generators
214  std::cout << "Done" << std::endl;
215 }
const HistAxis axisNumuForNueSig
ana::Loaders * GetLoaders2018(const TString option, const TString period="full")
std::vector< ShiDef > GetShifts2018(const TString option)
std::map< TString, IPredictionGenerator * > gens
Definition: syst_header.h:387
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
Struct to hold prediction information.
OStream cerr
Definition: OStream.cxx:7
string filename
Definition: shutoffs.py:106
int gIdx
Definition: show_event.C:12
Generates FD-only predictions (no extrapolation)
Loaders::FluxType flux
Generates extrapolated Numu predictions.
Generates extrapolated Nue predictions using ProportionalDecomp.
std::vector< ExtrapDef > GetExtrapolationDefs(const TString analysis, const TString period)
const XML_Char * prefix
Definition: expat.h:380
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
TDirectory * dir
Definition: macro.C:5
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void make_predictions_systs(const TString option, const TString analysis="allxp_nue", const TString period="full", const TString outdir="./")
Generates extrapolated Nue predictions using Michel+BEN decomposition.
TFile * file
Definition: cellShifts.C:17
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Given loaders and an MC shift, Generate() generates an IPrediction.
void SwapNDDataLoader(Loaders *loaders, const TString option, const TString period="full")
const std::string outdir
unsigned int uint