PredictionHandler.h
Go to the documentation of this file.
1 // PredictionHandler Class
2 // This is a refactoring of the Nus18PredictionHandler written by Mike Wallbank
3 //
4 // PredictionHandler acts as a base class with FDPredictionHandler and NDPredictionHandler
5 // acting as derived classes
6 //
7 // author : Adam Lister
8 // email : alister1@fnal.gov
9 
10 #include "CAFAna/Core/Loaders.h"
12 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Core/Binning.h"
14 #include "CAFAna/Core/ISyst.h"
15 #include "CAFAna/Core/SystShifts.h"
16 #include "CAFAna/Cuts/TimingCuts.h"
22 #include "CAFAna/Vars/HistAxes.h"
23 
25 
26 #include "NuXAna/Cuts/NusCuts18.h"
31 
32 #include "OscLib/OscCalcSterile.h"
33 
34 using namespace ana;
35 
36 // base class
38 {
39  public:
40  PredictionHandler(std::string predOpt, bool isData = false, Loaders* loaders = nullptr);
41  virtual ~PredictionHandler();
42 
43  virtual void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight, bool cutdata = true) = 0;
44 
45  virtual void AddVar(std::string name, std::string label, Var var, Binning binning, Cut cuts, SystShifts shift, Var weight, bool cutdata = true) = 0;
46 
47  virtual void AddVar(std::string name,
48  std::string labelx, std::string labely,
49  Var varx, Var vary,
50  Binning binningx, Binning binningy,
51  Cut cuts, SystShifts shift, Var weight, bool cutdata = true) = 0;
52 
53  virtual void AddVar(std::string name,
54  std::string labelx, std::string labely, std::string labelz,
55  Var varx, Var vary, Var varz,
56  Binning binningx, Binning binningy, Binning binningz,
57  Cut cuts, SystShifts shift, Var weight, bool cutdata = true) = 0;
58 
59 
60  virtual void Save(TFile* outFile) = 0;
61 
62  void AddLoaders(Loaders* loaders);
63 
64  void AddData(SpectrumLoader* loadersData);
65 
66  void AddSystematics(std::vector<const ISyst*> systs, osc::OscCalcSterile* calc);
67 
68  protected:
69  // Loaders
73 
74  // Systematics
75  std::vector<const ISyst*> fSysts;
77 
78  // Data product containers
79  std::map<std::string, IPrediction*> fPredictions;
80  std::map<std::string, Spectrum*> fSpectra;
81  std::map<std::string, Spectrum*> fCosSpec;
82  std::map<std::string, CheatDecomp*> fDecompositions;
83  std::map<std::string, std::map<std::string, std::pair<Spectrum*, Spectrum*> > > fShiftedSpectra;
84 
85  // options
87  bool fData;
88  bool fRunSysts;
89  bool fCosmics;
90 
91 };
92 
94  bool isData ,
95  Loaders* loaders) :
96  fPredOpt (predOpt),
97  fLoaders (loaders),
98  fData (isData) ,
99  fRunSysts (false) ,
100  fCosmics (false)
101 {
102 
103  if (fPredOpt != "Extrap" && fPredOpt != "FDExtrap" && fPredOpt != "NoExtrap" && fPredOpt != "Decomp")
104  {
105  std::cout << "PredictionHandler not set up for "
106  << fPredOpt
107  << "\nCurrent allowed options are:"
108  << "\n- FD: "
109  << "\n-- Extrap"
110  << "\n-- FDExtrap"
111  << "\n-- NoExtrap"
112  << "\n-- Decomp"
113  << std::endl;
114  std::abort();
115  }
116 }
117 
119 {
120  for (const auto& pred : fPredictions)
121  delete pred.second;
122  for (const auto& cos : fCosSpec)
123  delete cos.second;
124  for (const auto& spectrum : fSpectra)
125  delete spectrum.second;
126  for (const auto& decomp : fDecompositions)
127  delete decomp.second;
128  for (const auto& shift : fShiftedSpectra) {
129  for (const auto& syst : shift.second) {
130  delete syst.second.first;
131  delete syst.second.second;
132  }
133  }
134 }
135 
137  fLoaders = loaders;
138 }
139 
141  fData = true;
142  fLoadersData = loadersData;
143 }
144 
145 void PredictionHandler::AddSystematics(std::vector<const ISyst*> systs,
147 {
148  fRunSysts = true;
149  fSysts = systs;
150  fOscCalc = calc;
151 }
152 
153 // ND Prediction handler class
155 {
156  public:
157  NDPredictionHandler(std::string predOpt, bool isData = false, Loaders* loaders = nullptr)
158  : PredictionHandler(predOpt, isData, loaders)
159  {
160  }
161 
163  {
164  }
165 
166  void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
167  void AddVar(std::string name, std::string label, Var var, Binning binning, Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
168  void AddVar(std::string name,
169  std::string labelx, std::string labely,
170  Var varx, Var vary,
171  Binning binningx, Binning binningy,
172  Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
173 
174  void AddVar(std::string name,
175  std::string labelx, std::string labely, std::string labelz,
176  Var varx, Var vary, Var varz,
177  Binning binningx, Binning binningy, Binning binningz,
178  Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
179 
180  void Save(TFile* outFile);
181 };
182 
185  Var var,
186  Binning binning,
187  Cut cuts,
188  SystShifts shift,
189  Var weight,
190  bool cutdata)
191 {
192 
193  const HistAxis axis(label, binning, var);
194 
195  return AddVar(name, axis, cuts, shift, weight, cutdata);
196 
197 }
198 
200  std::string labelx,
201  std::string labely,
202  Var varx,
203  Var vary,
204  Binning binningx,
205  Binning binningy,
206  Cut cuts,
207  SystShifts shift,
208  Var weight,
209  bool cutdata)
210 {
211 
212  const HistAxis axis(labelx, binningx, varx, labely, binningy, vary);
213 
214  return AddVar(name, axis, cuts, shift, weight, cutdata);
215 
216 }
217 
219  std::string labelx, std::string labely, std::string labelz,
220  Var varx, Var vary, Var varz,
221  Binning binningx, Binning binningy, Binning binningz,
222  Cut cuts, SystShifts shift, Var weight, bool cutdata)
223 {
224  const HistAxis axis ({labelx, labely, labelz}, {binningx, binningy, binningz}, {varx, vary, varz});
225 
226  return AddVar(name, axis, cuts, shift, weight, cutdata);
227 }
228 
230  HistAxis axis,
231  Cut cuts,
232  SystShifts shift,
233  Var weight,
234  bool cutdata)
235 {
236 
237  CheatDecomp* decomp = nullptr;
238  IPrediction* pred = nullptr;
239  Spectrum* spec = nullptr;
240 
241  // Decomposition
242  if (fPredOpt == "Decomp")
244  axis, cuts, shift, weight);
245 
246  // Predictions
247  if (fPredOpt == "NDExtrap") {
248  const NDPredictionGenerator genND(axis, cuts, shift, weight);
249  pred = genND.Generate(*fLoaders).release();
250  }
251 
252  // Systematics
253  std::map<std::string, std::pair<Spectrum*, Spectrum*> > shifts;
254  if (fRunSysts) {
255  for (const auto& shift : fSysts) {
257  axis, cuts, SystShifts(shift, +1), weight);
259  axis, cuts, SystShifts(shift, -1), weight);
260  shifts[shift->ShortName()] = std::make_pair(shiftUp, shiftDown);
261  }
262  }
263 
264  // Make data spectrum
265  if (fData && cutdata)
266  spec = new Spectrum(*fLoadersData, axis, cuts);
267 
268  if (decomp)
270  if (pred)
272  if (shifts.size())
273  fShiftedSpectra[name] = shifts;
274  if (spec)
275  fSpectra[name] = spec;
276 
277  return;
278 
279 }
280 
282 {
283 
284  std::cout << "saving file..." << std::endl;
285 
286  // Save MC decompositions to file
287  for (std::map<std::string, CheatDecomp*>::const_iterator decompIt = fDecompositions.begin(); decompIt != fDecompositions.end(); ++decompIt)
288  decompIt->second->SaveTo(outFile->mkdir(Form("ND%sDecomp", decompIt->first.c_str())));
289 
290  // Save systematics to file
291  for (std::map<std::string, std::map<std::string, std::pair<Spectrum*, Spectrum*> > >::const_iterator shiftIt = fShiftedSpectra.begin();
292  shiftIt != fShiftedSpectra.end(); ++shiftIt) {
293  TString subdir = Form("ND%sDecompSysts", shiftIt->first.c_str());
294  outFile->cd();
295  outFile->mkdir(subdir);
296  for (std::map<std::string, std::pair<Spectrum*, Spectrum*> >::const_iterator systIt = shiftIt->second.begin();
297  systIt != shiftIt->second.end(); ++systIt) {
298  TString subsubdir = systIt->first.c_str();
299  outFile->cd(subdir);
300  gDirectory->mkdir(subsubdir);
301  gDirectory->cd(subsubdir);
302  systIt->second.first->SaveTo(gDirectory->mkdir("ShiftUp"));
303  systIt->second.second->SaveTo(gDirectory->mkdir("ShiftDown"));
304  }
305  }
306 
307  // Save data spectra to file
308  if (fData){
309  for (std::map<std::string, Spectrum*>::const_iterator specIt = fSpectra.begin(); specIt != fSpectra.end(); ++specIt){
310  specIt->second->SaveTo(outFile->mkdir(Form("ND%sSpec", specIt->first.c_str())));
311  }
312  }
313 
314 
315  return;
316 
317 }
318 
319 // FD Prediction handler class
321 {
322  public:
323  FDPredictionHandler(std::string predOpt, bool isData = false, Loaders* loaders = nullptr)
324  : PredictionHandler(predOpt, isData, loaders)
325  {
326  }
327 
329  {
330  }
331 
332  void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
333  void AddVar(std::string name, std::string label, Var var, Binning binning, Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
334  void AddVar(std::string name,
335  std::string labelx, std::string labely,
336  Var varx, Var vary,
337  Binning binningx, Binning binningy,
338  Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
339 
340  void AddVar(std::string name,
341  std::string labelx, std::string labely, std::string labelz,
342  Var varx, Var vary, Var varz,
343  Binning binningx, Binning binningy, Binning binningz,
344  Cut cuts, SystShifts shift, Var weight, bool cutdata = true);
345 
346  void Save(TFile* outFile);
347  void AddCosmics(SpectrumLoader* loadersCosmic);
348 };
349 
351 
352  const HistAxis axis(label, binning, var);
353 
354  return AddVar(name, axis, cuts, shift, weight, cutdata);
355 
356 }
357 
359  std::string labelx,
360  std::string labely,
361  Var varx,
362  Var vary,
363  Binning binningx,
364  Binning binningy,
365  Cut cuts,
366  SystShifts shift,
367  Var weight,
368  bool cutdata)
369 {
370 
371  const HistAxis axis(labelx, binningx, varx, labely, binningy, vary);
372 
373  return AddVar(name, axis, cuts, shift, weight, cutdata);
374 
375 }
376 
378  std::string labelx, std::string labely, std::string labelz,
379  Var varx, Var vary, Var varz,
380  Binning binningx, Binning binningy, Binning binningz,
381  Cut cuts, SystShifts shift, Var weight, bool cutdata)
382 {
383  const HistAxis axis ({labelx, labely, labelz}, {binningx, binningy, binningz}, {varx, vary, varz});
384 
385  return AddVar(name, axis, cuts, shift, weight, cutdata);
386 }
387 
389 
390  if (!fLoaders) {
391  std::cout << "Error: Loaders must be set before variables can be added." << std::endl;
392  abort();
393  }
394 
395  Spectrum* spec = nullptr;
396  Spectrum* cos = nullptr;
397  IPrediction* pred = nullptr;
398 
399  // Make the prediction generator
400  IPredictionGenerator* pred_gen = nullptr;
401  if (fPredOpt == "Extrap") {
403  cuts,
405  shift, weight);
406  pred_gen = genExtrap;
407  }
408  else if (fPredOpt == "FDExtrap") {
409  FDPredictionGenerator* genFD = new FDPredictionGenerator(axis, cuts, shift, weight);
410  pred_gen = genFD;
411  }
412  else if (fPredOpt == "NoExtrap") {
413  NoExtrapGenerator* genNoExtrap = new NoExtrapGenerator(axis, cuts, weight);
414  pred_gen = genNoExtrap;
415  }
416 
417  // Get the prediction
418  if (fRunSysts) {
419  PredictionInterp* pred_interp = new PredictionInterp(fSysts, fOscCalc, *pred_gen, *fLoaders);
420  pred = pred_interp;
421  }
422  else
423  pred = pred_gen->Generate(*fLoaders).release();
424 
425  // Make data spectrum
426  if (fData && cutdata)
427  spec = new Spectrum(*fLoadersData, axis, cuts && kInBeamSpill);
428 
429  // Make cosmic spectrum
430  if (fCosmics && cutdata)
431  cos = new Spectrum(*fLoadersCos, axis, cuts && kInCosmicTimingWindow);
432 
433  // Save the data products
434  if (pred)
436  if (spec)
437  fSpectra[name] = spec;
438  if (cos)
439  fCosSpec[name] = cos;
440 
441  return;
442 
443 }
444 
445 void FDPredictionHandler::Save(TFile* outFile) {
446 
447  // Save MC predictions to file
448  for (std::map<std::string, IPrediction*>::const_iterator predIt = fPredictions.begin();
449  predIt != fPredictions.end(); ++predIt)
450  predIt->second->SaveTo(outFile->mkdir(Form("FD%sPred", predIt->first.c_str())));
451 
452  // Save data spectra to file
453  if (fData)
454  for (std::map<std::string, Spectrum*>::const_iterator specIt = fSpectra.begin();
455  specIt != fSpectra.end(); ++specIt)
456  specIt->second->SaveTo(outFile->mkdir(Form("FD%sSpec", specIt->first.c_str())));
457 
458  // Save cosmics to file
459  if (fCosmics)
460  for (std::map<std::string, Spectrum*>::const_iterator cosIt = fCosSpec.begin();
461  cosIt != fCosSpec.end(); ++cosIt)
462  cosIt->second->SaveTo(outFile->mkdir(Form("FD%sCosSpec", cosIt->first.c_str())));
463 
464  return;
465 
466 }
467 
469  fCosmics = true;
470  fLoadersCos = loadersCosmic;
471 }
472 
473 
Near Detector underground.
Definition: SREnums.h:10
const XML_Char * name
Definition: expat.h:151
Implements systematic errors by interpolation between shifted templates.
void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight)
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 weight
subdir
Definition: cvnie.py:7
NDPredictionHandler(std::string predOpt, bool isData=false, Loaders *loaders=nullptr)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
void Save(TDirectory *outDir, PlotInfo info, osc::OscCalcSterile *calc, HistType ht, IPrediction *pred=nullptr, Spectrum *cosSpec=nullptr, Spectrum *data=nullptr)
const Color_t kMC
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
void Save(TFile *outFile)
std::map< std::string, Spectrum * > fSpectra
Adapt the PMNS_Sterile calculator to standard interface.
SpectrumLoader * fLoadersCos
void AddData(SpectrumLoader *loadersData)
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
PredictionHandler(std::string predOpt, bool isData=false, Loaders *loaders=nullptr)
Generates extrapolated NC predictions using ProportionalDecomp.
osc::OscCalcDumb calc
Generates FD-only predictions (no extrapolation)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
const Cut kInCosmicTimingWindow
Is the event far from the start and ends of the spill ? For FD cosmic selection.
Definition: TimingCuts.cxx:165
TFile * outFile
Definition: PlotXSec.C:135
virtual void Save(TFile *outFile)=0
osc::OscCalcSterile * fOscCalc
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
virtual void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight, bool cutdata=true)=0
void AddSystematics(std::vector< const ISyst * > systs, osc::OscCalcSterile *calc)
SpectrumLoader * fLoadersData
Class to handle ND predictions.
void AddCosmics(SpectrumLoader *loaders_cosmics)
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
FDPredictionHandler(std::string predOpt, bool isData=false, Loaders *loaders=nullptr)
std::vector< float > Spectrum
Definition: Constants.h:610
std::vector< const ISyst * > fSysts
void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight)
OStream cout
Definition: OStream.cxx:6
void Save(TFile *outFile)
std::map< std::string, CheatDecomp * > fDecompositions
std::map< std::string, Spectrum * > fCosSpec
void AddLoaders(Loaders *loaders)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
virtual ~PredictionHandler()
virtual std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const =0
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
T cos(T number)
Definition: d0nt_math.hpp:78
const HistAxis kNus18BinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNus18EnergyBinning, kCCE)
Definition: HistAxes.h:17
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
std::map< std::string, IPrediction * > fPredictions
Given loaders and an MC shift, Generate() generates an IPrediction.
Class to handle FD predictions.
Generates Near Detector predictions.
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
std::map< std::string, std::map< std::string, std::pair< Spectrum *, Spectrum * > > > fShiftedSpectra
enum BeamMode string