Nus18PredictionHandlers.h
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 // Nus18PredictionHandlers.h
3 //
4 // Mike Wallbank (wallbank@fnal.gov)
5 // May 2018
6 //
7 // General purpose classes to handle the generation, handling and
8 // writing of prediction objects.
9 // TODO: should probably be a base class from which FD and ND
10 // PredictionHandlers inherit.
11 //////////////////////////////////////////////////////////////////////////
12 
13 #include "CAFAna/Core/Loaders.h"
15 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/ISyst.h"
17 #include "CAFAna/Core/SystShifts.h"
18 #include "NuXAna/Cuts/NusCuts18.h"
20 #include "CAFAna/Cuts/TimingCuts.h"
30 #include "CAFAna/Vars/HistAxes.h"
31 #include "NuXAna/Vars/HistAxes.h"
32 #include "OscLib/OscCalcSterile.h"
33 
34 using namespace ana;
35 
36 /// Class to handle FD predictions
38 public:
39 
42 
43  void AddCosmics(SpectrumLoader* loaders_cosmics);
44  void AddData(SpectrumLoader* loaders_data);
45  void AddLoaders(Loaders* loaders);
46  void AddSystematics(std::vector<const ISyst*> systs, osc::OscCalcSterile* calc);
47  void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight);
48  void AddVar(std::string name, std::string label, Var var, Binning binning, Cut cuts, SystShifts shift, Var weight);
49  void Save(TFile* outFile);
50 
51 private:
52 
53  // Loaders
57 
58  // Systematics
59  std::vector<const ISyst*> fSysts;
61 
62  // Data product containers
63  std::map<std::string, IPrediction*> fPredictions;
64  std::map<std::string, Spectrum*> fSpectra;
65  std::map<std::string, Spectrum*> fCosSpec;
66 
67  // options
69  bool fData;
70  bool fRunSysts;
71  bool fCosmics;
72 
73 };
74 
76  fPredOpt(predOpt),
77  fLoaders(nullptr),
78  fData(false),
79  fRunSysts(false),
80  fCosmics(false) {
81 
82  // Make sure the Prediction Handler is set up to handle the prediction option
83  if (fPredOpt != "Extrap" and fPredOpt != "FDExtrap" and fPredOpt != "NoExtrap") {
84  std::cout << "The Prediction Handler is not set up to handle " << fPredOpt << " predictions." << std::endl
85  << "Currently options are Extrap (PredictionSterile), "
86  << "FDExtrap (FDPredictionSterile) or NoExtrap (PredictionNoExtrap)." << std::endl;
87  abort();
88  }
89 }
90 
92  for (const auto& pred : fPredictions)
93  delete pred.second;
94  for (const auto& spectrum : fSpectra)
95  delete spectrum.second;
96  for (const auto& cos : fCosSpec)
97  delete cos.second;
98 }
99 
101  fCosmics = true;
102  fLoadersCos = loaders_cosmics;
103 }
104 
106  fData = true;
107  fLoadersData = loaders_data;
108 }
109 
111  fLoaders = loaders;
112 }
113 
115  fRunSysts = true;
116  fSysts = systs;
117  fOscCalc = calc;
118 }
119 
121 
122  const HistAxis axis(label, binning, var);
123 
124  return AddVar(name, axis, cuts, shift, weight);
125 
126 }
127 
129 
130  if (!fLoaders) {
131  std::cout << "Error: Loaders must be set before variables can be added." << std::endl;
132  abort();
133  }
134 
135  Spectrum* spec = nullptr;
136  Spectrum* cos = nullptr;
137  IPrediction* pred = nullptr;
138 
139  // Make the prediction generator
140  IPredictionGenerator* pred_gen = nullptr;
141  if (fPredOpt == "Extrap") {
143  cuts, kNus18ND,
144  kNumuCutND2018, shift, weight);
145  pred_gen = genExtrap;
146  }
147  else if (fPredOpt == "FDExtrap") {
148  FDPredictionGenerator* genFD = new FDPredictionGenerator(axis, cuts, shift, weight);
149  pred_gen = genFD;
150  }
151  else if (fPredOpt == "NoExtrap") {
152  NoExtrapGenerator* genNoExtrap = new NoExtrapGenerator(axis, cuts, weight);
153  pred_gen = genNoExtrap;
154  }
155 
156  // Get the prediction
157  if (fRunSysts) {
158  PredictionInterp* pred_interp = new PredictionInterp(fSysts, fOscCalc, *pred_gen, *fLoaders);
159  pred = pred_interp;
160  }
161  else
162  pred = pred_gen->Generate(*fLoaders).release();
163 
164  // Make data spectrum
165  if (fData)
166  spec = new Spectrum(*fLoadersData, axis, cuts && kInBeamSpill);
167 
168  // Make cosmic spectrum
169  if (fCosmics)
170  cos = new Spectrum(*fLoadersCos, axis, cuts && kInCosmicTimingWindow);
171 
172  // Save the data products
173  if (pred)
175  if (spec)
176  fSpectra[name] = spec;
177  if (cos)
178  fCosSpec[name] = cos;
179 
180  return;
181 
182 }
183 
185 
186  // Save MC predictions to file
187  for (std::map<std::string, IPrediction*>::const_iterator predIt = fPredictions.begin();
188  predIt != fPredictions.end(); ++predIt)
189  predIt->second->SaveTo(outFile, Form("FD%sPred", predIt->first.c_str()));
190 
191  // Save data spectra to file
192  if (fData)
193  for (std::map<std::string, Spectrum*>::const_iterator specIt = fSpectra.begin();
194  specIt != fSpectra.end(); ++specIt)
195  specIt->second->SaveTo(outFile, Form("FD%sSpec", specIt->first.c_str()));
196 
197  // Save cosmics to file
198  if (fCosmics)
199  for (std::map<std::string, Spectrum*>::const_iterator cosIt = fCosSpec.begin();
200  cosIt != fCosSpec.end(); ++cosIt)
201  cosIt->second->SaveTo(outFile, Form("FD%sCosSpec", cosIt->first.c_str()));
202 
203  return;
204 
205 }
206 
207 /// Class to handle ND predictions
209 public:
210 
213 
214  void AddSystematics(std::vector<const ISyst*> systs, osc::OscCalcSterile* calc);
216  void AddVar(std::string name, std::string label, Var var, Binning binning, Cut cuts, SystShifts shift, Var weight);
217  void Save(TFile* outFile);
218 
219 private:
220 
221  // Loaders
223 
224  // Systematics
225  std::vector<const ISyst*> fSysts;
227 
228  // Data product containers
229  std::map<std::string, CheatDecomp*> fDecompositions;
230  std::map<std::string, IPrediction*> fPredictions;
231  std::map<std::string, Spectrum*> fSpectra;
232  std::map<std::string, std::map<std::string, std::pair<Spectrum*, Spectrum*> > > fShiftedSpectra;
233 
234  // options
236  bool fData;
237  bool fRunSysts;
238 
239 };
240 
242  fLoaders(loaders),
243  fPredOpt(predOpt),
244  fData(data),
245  fRunSysts(false) {
246 
247  // Make sure the Prediction Handler is set up to handle the prediction option
248  if (fPredOpt != "Decomp" and fPredOpt != "NDExtrap") {
249  std::cout << "The Prediction Handler is not set up to handle " << fPredOpt << " predictions." << std::endl;
250  std::cout << "Currently options are Decomp (CheatDecomp), or NDExtrap (NDPredictionSterile)." << std::endl;
251  abort();
252  }
253 }
254 
256  for (const auto& decomp : fDecompositions)
257  delete decomp.second;
258  for (const auto& shift : fShiftedSpectra) {
259  for (const auto& syst : shift.second) {
260  delete syst.second.first;
261  delete syst.second.second;
262  }
263  }
264  for (const auto& spectum : fSpectra)
265  delete spectum.second;
266 }
267 
269  fRunSysts = true;
270  fSysts = systs;
271  fOscCalc = calc;
272 }
273 
275 
276  const HistAxis axis(label, binning, var);
277 
278  return AddVar(name, axis, cuts, shift, weight);
279 
280 }
281 
283 
284  CheatDecomp* decomp = nullptr;
285  IPrediction* pred = nullptr;
286  Spectrum* spec = nullptr;
287 
288  // Decomposition
289  if (fPredOpt == "Decomp")
291  axis, cuts, shift, weight);
292 
293  // Predictions
294  if (fPredOpt == "NDExtrap") {
295  const NDPredictionGenerator genND(axis, cuts, shift, weight);
296  pred = genND.Generate(*fLoaders).release();
297  }
298 
299  // Systematics
300  std::map<std::string, std::pair<Spectrum*, Spectrum*> > shifts;
301  if (fRunSysts) {
302  for (const auto& shift : fSysts) {
304  axis, cuts, SystShifts(shift, +1), weight);
306  axis, cuts, SystShifts(shift, -1), weight);
307  shifts[shift->ShortName()] = std::make_pair(shiftUp, shiftDown);
308  }
309  }
310 
311  // Data
312  if (fData)
314  axis, cuts);
315 
316  if (decomp)
318  if (pred)
320  if (shifts.size())
321  fShiftedSpectra[name] = shifts;
322  if (spec)
323  fSpectra[name] = spec;
324 
325  return;
326 
327 }
328 
330 
331  // Save MC decompositions to file
332  for (std::map<std::string, CheatDecomp*>::const_iterator decompIt = fDecompositions.begin();
333  decompIt != fDecompositions.end(); ++decompIt)
334  decompIt->second->SaveTo(outFile, Form("ND%sDecomp", decompIt->first.c_str()));
335 
336  // Save systematics to file
337  for (std::map<std::string, std::map<std::string, std::pair<Spectrum*, Spectrum*> > >::const_iterator shiftIt = fShiftedSpectra.begin();
338  shiftIt != fShiftedSpectra.end(); ++shiftIt) {
339  TString subdir = Form("ND%sDecompSysts", shiftIt->first.c_str());
340  outFile->cd();
341  outFile->mkdir(subdir);
342  for (std::map<std::string, std::pair<Spectrum*, Spectrum*> >::const_iterator systIt = shiftIt->second.begin();
343  systIt != shiftIt->second.end(); ++systIt) {
344  TString subsubdir = systIt->first.c_str();
345  outFile->cd(subdir);
346  gDirectory->mkdir(subsubdir);
347  gDirectory->cd(subsubdir);
348  systIt->second.first->SaveTo(gDirectory, "ShiftUp");
349  systIt->second.second->SaveTo(gDirectory, "ShiftDown");
350  }
351  }
352 
353  // Save data spectra to file
354  if (fData)
355  for (std::map<std::string, Spectrum*>::const_iterator specIt = fSpectra.begin();
356  specIt != fSpectra.end(); ++specIt)
357  specIt->second->SaveTo(outFile, Form("ND%sSpec", specIt->first.c_str()));
358 
359  return;
360 
361 }
SpectrumLoader * fLoadersData
std::map< std::string, Spectrum * > fCosSpec
Near Detector underground.
Definition: SREnums.h:10
std::vector< const ISyst * > fSysts
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
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
std::map< std::string, IPrediction * > fPredictions
void Save(TFile *outFile)
Adapt the PMNS_Sterile calculator to standard interface.
std::map< std::string, Spectrum * > fSpectra
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
Generates extrapolated NC predictions using ProportionalDecomp.
osc::OscCalcDumb calc
std::map< std::string, CheatDecomp * > fDecompositions
Generates FD-only predictions (no extrapolation)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::map< std::string, Spectrum * > fSpectra
const char * label
osc::OscCalcSterile * fOscCalc
const XML_Char const XML_Char * data
Definition: expat.h:268
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
void AddLoaders(Loaders *loaders)
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
std::map< std::string, std::map< std::string, std::pair< Spectrum *, Spectrum * > > > fShiftedSpectra
Class to handle ND predictions.
void AddCosmics(SpectrumLoader *loaders_cosmics)
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
std::vector< float > Spectrum
Definition: Constants.h:610
void AddVar(std::string name, HistAxis axis, Cut cuts, SystShifts shift, Var weight)
OStream cout
Definition: OStream.cxx:6
FDPredictionHandler(std::string predOpt)
void Save(TFile *outFile)
void AddSystematics(std::vector< const ISyst * > systs, osc::OscCalcSterile *calc)
std::vector< const ISyst * > fSysts
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
virtual std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const =0
SpectrumLoader * fLoadersCos
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::map< std::string, IPrediction * > fPredictions
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
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
void AddSystematics(std::vector< const ISyst * > systs, osc::OscCalcSterile *calc)
NDPredictionHandler(Loaders *loaders, std::string predOpt, bool data)
osc::OscCalcSterile * fOscCalc
void AddData(SpectrumLoader *loaders_data)
enum BeamMode string