PredictionInterp.h
Go to the documentation of this file.
1 #pragma once
2 
5 
9 
10 #include <map>
11 #include <memory>
12 #include <unordered_map>
13 
14 #include "TMD5.h"
15 
16 class TH1;
17 
18 namespace ana
19 {
20  class Loaders;
21 
22  /// Implements systematic errors by interpolation between shifted templates
24  {
25  public:
26  enum EMode_t{
28  };
29 
30  /// \param systs What systematics we should be capable of interpolating
31  /// \param osc The oscillation point to expand around
32  /// \param predGen Construct an IPrediction from the following
33  /// information.
34  /// \param loaders The loaders to be passed on to the underlying prediction
35  /// \param shiftMC Underlying shift. Use with care. Mostly for
36  /// PredictionNumuFAHadE. Should not contain any of of
37  /// \a systs
38  /// \param mode In FHC the wrong-sign has bad stats and probably the
39  /// fits can't be split out reasonably. For RHC it's
40  /// important not to conflate them.
41  PredictionInterp(std::vector<const ISyst*> systs,
43  const IPredictionGenerator& predGen,
45  const SystShifts& shiftMC = kNoShift,
47 
48  virtual ~PredictionInterp();
49 
50 
51 
52  Spectrum Predict(osc::IOscCalc* calc) const override;
53  Spectrum Predict(osc::IOscCalcStan* calc) const override;
54 
55 
57  const SystShifts& shift) const override;
59  const SystShifts& shift) const override;
60 
62  Flavors::Flavors_t flav,
64  Sign::Sign_t sign) const override;
66  Flavors::Flavors_t flav,
68  Sign::Sign_t sign) const override;
69 
71  const SystShifts& shift,
72  Flavors::Flavors_t flav,
74  Sign::Sign_t sign) const override;
76  const SystShifts& shift,
77  Flavors::Flavors_t flav,
79  Sign::Sign_t sign) const override;
80 
81  virtual void SaveTo(TDirectory* dir, const std::string& name) const override;
82  static std::unique_ptr<PredictionInterp> LoadFrom(TDirectory* dir, const std::string& name);
83 
84  /// After calling this DebugPlots won't work fully and SaveTo won't work at
85  /// all.
86  void MinimizeMemory();
87 
88  void DebugPlot(const ISyst* syst,
92  Sign::Sign_t sign = Sign::kBoth) const;
93 
94  // If \a savePattern is not empty, print each pad. Must contain a "%s" to
95  // contain the name of the systematic.
97  const std::string& savePattern = "",
100  Sign::Sign_t sign = Sign::kBoth) const;
101 
102  void SetOscSeed(osc::IOscCalc* oscSeed);
103 
104  void DebugPlotColz(const ISyst* syst,
108  Sign::Sign_t sign = Sign::kBoth) const;
109 
111  const std::string& savePattern = "",
114  Sign::Sign_t sign = Sign::kBoth) const;
115 
116  bool SplitBySign() const {return fSplitBySign;}
119  kOther, ///< Taus, numu appearance
121  };
122 
123  PredictionInterp() : fOscOrigin(nullptr), fBinning(Spectrum::Uninitialized()), fSplitBySign(false) {}
124 
125  static void LoadFromBody(TDirectory* dir, PredictionInterp* ret,
126  std::vector<const ISyst*> veto = {});
127 
128  struct Coeffs{
129  Coeffs(double _a, double _b, double _c, double _d)
130  : a(_a), b(_b), c(_c), d(_d) {}
131  double a, b, c, d;
132  };
133 
134  /// Find coefficients describing this set of shifts
135  std::vector<std::vector<Coeffs>>
136  FitRatios(const std::vector<double>& shifts,
137  const std::vector<Eigen::ArrayXd>& ratios) const;
138 
139  /// Find coefficients describing the ratios from this component
140  std::vector<std::vector<Coeffs>>
141  FitComponent(const std::vector<double>& shifts,
142  const std::vector<IPrediction*>& preds,
143  Flavors::Flavors_t flav,
146  const std::string& systName) const;
147 
148  // What systs does this PredictionInterp know about?
149  std::vector<const ISyst*> GetAllSysts() const;
150 
153  bool nubar, // try to use fitsNubar if it exists
154  const SystShifts& shift) const;
155 
156  /// Helper for PredictComponentSyst
158  const TMD5* hash,
159  const SystShifts& shift,
160  Flavors::Flavors_t flav,
163  CoeffsType type) const;
164 
166  const TMD5* hash,
167  const SystShifts& shift,
168  Flavors::Flavors_t flav,
171  CoeffsType type) const;
172 
173  protected:
174  std::unique_ptr<IPrediction> fPredNom; ///< The nominal prediction
175 
177  {
178  double Stride() const {return shifts.size() > 1 ? shifts[1]-shifts[0] : 1;}
179 
180  std::string systName; ///< What systematic we're interpolating
181  std::vector<double> shifts; ///< Shift values sampled
182  std::vector<IPrediction*> preds;
183 
184  int nCoeffs; // Faster than calling size()
185 
186  /// Indices: [type][histogram bin][shift bin]
187  std::vector<std::vector<std::vector<Coeffs>>> fits;
188  /// Will be filled if signs are separated, otherwise not
189  std::vector<std::vector<std::vector<Coeffs>>> fitsNubar;
190 
191  void FillRemaps();
192 
193  // Same info as above but with more-easily-iterable index order
194  // [type][shift bin][histogram bin]. TODO this is ugly
195  std::vector<std::vector<std::vector<Coeffs>>> fitsRemap;
196  std::vector<std::vector<std::vector<Coeffs>>> fitsNubarRemap;
197  };
198 
199  mutable std::unordered_map<const ISyst*, ShiftedPreds> fPreds;
200  mutable std::unordered_map<const ISyst*, ShiftedPreds> fSumPreds;
201 
202  /// The oscillation values we assume when evaluating the coefficients
204 
205  mutable Spectrum fBinning; ///< Dummy spectrum to provide binning
206 
207  struct Key_t
208  {
212  bool operator<(const Key_t& rhs) const
213  {
214  return (std::make_tuple(flav, curr, sign) <
215  std::make_tuple(rhs.flav, rhs.curr, rhs.sign));
216  }
217  };
218  struct Val_t
219  {
220  TMD5 hash;
221  Spectrum nom; // todo: we can't cache stan::math::vars because they wind up getting invalidated when the Stan stack is cleared. but keeping only a <double> version around in this cache means that we're dumping the autodiff for the oscillation calculator part, which may mean Stan won't explore the space correctly. Not sure what to do here.
222  };
224 
226 
227  void InitFits() const;
228 
229  void InitFitsHelper(ShiftedPreds& sp,
230  std::vector<std::vector<std::vector<Coeffs>>>& fits,
231  Sign::Sign_t sign) const;
232 
233  /// Templated helper for \ref ShiftedComponent
234  template <typename T>
236  const TMD5* hash,
237  const SystShifts& shift,
238  Flavors::Flavors_t flav,
241  CoeffsType type) const;
242 
243  /// Templated helper for \ref PredictComponentSyst
244  template <typename T>
246  const SystShifts& shift,
247  Flavors::Flavors_t flav,
249  Sign::Sign_t sign) const;
250 
251  /// Helper for \ref ShiftSpectrum
252  template <typename T>
253  void ShiftBins(unsigned int N,
254  T* arr,
256  bool nubar,
257  const SystShifts& shift) const;
258  };
259 
260 }
void ShiftBins(unsigned int N, T *arr, CoeffsType type, bool nubar, const SystShifts &shift) const
Helper for ShiftSpectrum.
const XML_Char * name
Definition: expat.h:151
Implements systematic errors by interpolation between shifted templates.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::string systName
What systematic we&#39;re interpolating.
ThreadLocal< std::map< Key_t, Val_t > > fNomCache
std::vector< std::vector< std::vector< Coeffs > > > fitsNubarRemap
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
std::vector< std::vector< std::vector< Coeffs > > > fitsNubar
Will be filled if signs are separated, otherwise not.
void SetOscSeed(osc::IOscCalc *oscSeed)
void DebugPlotsColz(osc::IOscCalc *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::unique_ptr< IPrediction > fPredNom
The nominal prediction.
const XML_Char * s
Definition: expat.h:262
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Interactions of both types.
Definition: IPrediction.h:42
std::vector< double > shifts
Shift values sampled.
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::vector< std::vector< Coeffs > > FitComponent(const std::vector< double > &shifts, const std::vector< IPrediction * > &preds, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, const std::string &systName) const
Find coefficients describing the ratios from this component.
std::vector< IPrediction * > preds
Spectrum fBinning
Dummy spectrum to provide binning.
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::unordered_map< const ISyst *, ShiftedPreds > fSumPreds
Oscillation probability calculators.
Definition: Calcs.h:5
void DebugPlots(osc::IOscCalc *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
const SystShifts kNoShift
Definition: SystShifts.cxx:21
Taus, numu appearance.
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum _ShiftedComponent(osc::_IOscCalc< T > *calc, const TMD5 *hash, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, CoeffsType type) const
Templated helper for ShiftedComponent.
bool operator<(const Key_t &rhs) const
void DebugPlotColz(const ISyst *syst, osc::IOscCalc *calc, Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< std::vector< std::vector< Coeffs > > > fitsRemap
TDirectory * dir
Definition: macro.C:5
Coeffs(double _a, double _b, double _c, double _d)
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< std::string > fits
std::vector< const ISyst * > GetAllSysts() const
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void DebugPlot(const ISyst *syst, osc::IOscCalc *calc, Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
void InitFitsHelper(ShiftedPreds &sp, std::vector< std::vector< std::vector< Coeffs >>> &fits, Sign::Sign_t sign) const
double T
Definition: Xdiff_gwt.C:5
Given loaders and an MC shift, Generate() generates an IPrediction.
std::unordered_map< const ISyst *, ShiftedPreds > fPreds
osc::IOscCalc * fOscOrigin
The oscillation values we assume when evaluating the coefficients.
static void LoadFromBody(TDirectory *dir, PredictionInterp *ret, std::vector< const ISyst * > veto={})
All neutrinos, any flavor.
Definition: IPrediction.h:26
std::vector< std::vector< Coeffs > > FitRatios(const std::vector< double > &shifts, const std::vector< Eigen::ArrayXd > &ratios) const
Find coefficients describing this set of shifts.
Spectrum ShiftedComponent(osc::IOscCalc *calc, const TMD5 *hash, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign, CoeffsType type) const
Helper for PredictComponentSyst.
Spectrum _PredictComponentSyst(osc::_IOscCalc< T > *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Templated helper for PredictComponentSyst.
Spectrum ShiftSpectrum(const Spectrum &s, CoeffsType type, bool nubar, const SystShifts &shift) const
def sign(x)
Definition: canMan.py:197
std::vector< std::vector< std::vector< Coeffs > > > fits
Indices: [type][histogram bin][shift bin].