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"
19 #include "CAFAna/Vars/HistAxes.h"
21 
22 #include <TFile.h>
23 #include <TString.h>
24 #include <vector>
25 #include <iostream>
26 
27 #include "CAFAna/shared/Ana2017/syst_variations.h"
28 
29 using namespace ana;
30 
31 
32 void make_predictions_systs(const TString option, const TString analysis = "allxp_nue", const TString period = "full", const TString outdir="./")
33 {
34 
35  auto loaders = GetLoaders2017(option, period);
36  if (option.Contains("fakeNDData")) SwapNDDataLoader(loaders, option, period);
37 
38  loaders->SetSpillCut(kStandardSpillCuts);
39 
40  auto shifts = GetShifts2017 (option);
41 
42  auto xpdefs = GetExtrapolationDefs(analysis, period);
43 
44  bool run_nxp = analysis.Contains("nxp_only") || analysis.Contains ("allxp");
45  bool run_xp_prop = analysis.Contains("xp_prop") || analysis.Contains ("allxp_nue");
46  bool run_xp_combo = analysis.Contains("xp_combo") || analysis.Contains ("allxp_nue");
47  bool run_xp_numu = analysis.Contains("xp_numu") || analysis.Contains ("allxp_numu");
48  uint NumExtrap = (uint)run_nxp + (uint)run_xp_prop + (uint)run_xp_combo + (uint)run_xp_numu;
49 
50  if(analysis.Contains("allxp_nue") && analysis.Contains ("allxp_numu")){
51  std::cerr << "All nue and numu xp? Ambitious much?" << std::endl;
52  return;
53  }
54 
55  std::cout << analysis << std :: endl;
56 
58 
59  if(option.Contains("RPA2018")) weight = kPPFXFluxCVWgt * kXSecCVWgt2018;
60 
61  struct PredDef{
62  const IPrediction * pred;
63  const TString xp_type;
64  const TString syst_name;
65  const TString sigma_name;
66  const TString edef_name;
67  };
68 
69  struct GenDef{
70  const IPredictionGenerator *gen;
71  const TString xp_type;
72  const TString edef_name;
73  };
74 
75  std::vector <GenDef> gens;
76  std::vector < std::vector <PredDef> > preds;
77 
78 
79  for(auto edef:xpdefs){
80 
81  if(run_nxp) gens.push_back ({new NoExtrapGenerator(edef.axis,
82  edef.fdcut,
83  weight),"pred_nxp", edef.name});
84  if(run_xp_prop)
85  gens.push_back({new NuePropExtrapGenerator(edef.axis,
87  edef.fdcut,
88  edef.nuendcut,
89  edef.numundcut,
90  kNoShift,
91  weight), "pred_xp_prop",edef.name});
92  if(run_xp_combo)
93  gens.push_back({new NueComboExtrapGenerator(edef.axis,
95  edef.fdcut,
96  edef.nuendcut,
97  edef.numundcut,
98  kNoShift,
99  weight) , "pred_xp_combo",edef.name});
100 
101  if(run_xp_numu)
102  gens.push_back({new NumuExtrapGenerator(edef.axis,
103  edef.fdcut,
104  edef.numundcut,
105  kNoShift,
106  weight),"pred_xp_numu",edef.name});
107 
108  }
109 
110  for (auto gen:gens){
111  std::vector <PredDef> temp_preds;
112  for (auto shift:shifts) {
113  auto pred = gen.gen->Generate(*loaders,
114  shift.shift);
115  temp_preds.push_back({pred.release(),gen.xp_type,
116  shift.syst_name,shift.sigma_name,gen.edef_name});
117  }//end shifts
118  preds.emplace_back(temp_preds);
119  }//end var/cut combos
120 
121  loaders->Go();
122 
123 
124  //mkdir is so dumb
125  //hadd is the worst. had to split predictions a little
126 
127  for(uint gIdx=0; gIdx < NumExtrap; ++gIdx){
128  // Configure what the tier is ==> Combine plotted variables and quantiles (if numu)
129  TString OutTier = gens[gIdx].xp_type;
130  for (uint xIdx=0; xIdx < xpdefs.size(); ++xIdx) {
131  int SInd = gIdx + ( xIdx * NumExtrap );
132  TString VarTier = gens[SInd].edef_name;
133  int VPos = gens[SInd].edef_name.Index('_');
134  if (VPos > 0)
135  VarTier = VarTier(0,VPos);
136  if (!OutTier.Contains(VarTier))
137  OutTier = OutTier + "_" + VarTier;
138  std::cout << "\t\tOutTier is now " << OutTier << std::endl;
139  }
140  if (analysis.Contains("numu")) OutTier = OutTier + "_AllQuantiles";
141  TString filename = (outdir + OutTier + "_" + period + "_" + option + ".root");
142  std::cout << "\nMade a file called " << filename << std::endl;
143 
144  auto file = new TFile (filename,"recreate");
145  // Now loop through all of the things which want to go into this file...
146  for (uint xIdx=0; xIdx < xpdefs.size(); ++xIdx) {
147  int SInd = gIdx + ( xIdx * NumExtrap );
148  for(uint sIdx = 0; sIdx < shifts.size(); ++sIdx){
149  if(preds[SInd][sIdx].syst_name != shifts[sIdx].syst_name ||
150  preds[SInd][sIdx].xp_type != gens[SInd].xp_type) {
151  std::cerr << "ERROR: Loops - you're doing it wrong \n"
152  << "Expected gen" << SInd << " " << gens[SInd].xp_type
153  << " and syst " << sIdx << " " << shifts[sIdx].syst_name
154  << "; but got pred " << preds[SInd][sIdx].xp_type << " " << preds[SInd][sIdx].syst_name
155  << std::endl;
156  return;
157  }//end check
158  auto dir = file->GetDirectory(preds[SInd][sIdx].xp_type+"_"+preds[SInd][sIdx].syst_name);
159  if(!dir) dir = file->mkdir(preds[SInd][sIdx].xp_type+"_"+preds[SInd][sIdx].syst_name);
160 
161  auto dir2 = dir->GetDirectory(preds[SInd][sIdx].sigma_name);
162  if(!dir2) dir2 = dir->mkdir(preds[SInd][sIdx].sigma_name);
163  TString prefix = (analysis.Contains("nue")?"nue_pred":"numu_pred");
164  preds[SInd][sIdx].pred->SaveTo(dir2, prefix + "_" + preds[SInd][sIdx].edef_name);
165  } //end shifts
166  }
167  delete file;
168  std::cout << "Saved predictions to " << filename
169  << " (" << gIdx+1 << "/" << NumExtrap << ")"
170  << std::endl;
171  } //end generators
172 
173  std::cout << "Done" << std::endl;
174 }
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
std::vector< ShiDef > GetShifts2017(const TString option)
int gIdx
Definition: show_event.C:12
Generates FD-only predictions (no extrapolation)
ana::Loaders * GetLoaders2017(const TString option, const TString period="full")
Generates extrapolated Numu predictions.
Generates extrapolated Nue predictions using ProportionalDecomp.
const HistAxis axis_numu_for_nuesig
std::vector< ExtrapDef > GetExtrapolationDefs(const TString analysis, const TString period)
const XML_Char * prefix
Definition: expat.h:380
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
TDirectory * dir
Definition: macro.C:5
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
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.
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
void SwapNDDataLoader(Loaders *loaders, const TString option, const TString period="full")
const std::string outdir
unsigned int uint