Spectrum.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/HistAxis.h"
6 #include "CAFAna/Core/Hist.h"
7 #include "CAFAna/Core/UtilsExt.h"
8 
9 #include <Eigen/Dense>
10 
11 #include "TAttLine.h"
12 
13 #include <set>
14 #include <string>
15 #include <vector>
16 
17 class TDirectory;
18 class TH1;
19 class TH2;
20 class TH3;
21 class TH1D;
22 
23 /// Oscillation analysis framework, runs over CAF files outside of ART
24 namespace ana
25 {
26  class Ratio;
27  class SpectrumLoaderBase;
28  template<class T> class _Var;
29  template<class T> class _Cut;
30  template<class T> class _MultiVar;
31 
32  template<class T> _Var<T> Unweighted();
33 
34  class SystShifts;
35  extern const SystShifts kNoShift;
36 
37  template<class T> class SpectrumSinkBase;
38 
39  /// Representation of a spectrum in any variable, with associated POT
40  class Spectrum
41  {
42  public:
43  friend class SpectrumLoaderBase;
44  friend class SpectrumSink;
45  friend class SpectrumSinkBase<Spectrum>;
46  friend class Ratio;
47 
49 
50  /// One constructor to rule them all
51  template<class T>
53  const _HistAxis<_Var<T>>& axis,
54  const _Cut<T>& cut,
55  const SystShifts& shift = kNoShift,
56  const _Var<T>& wei = Unweighted<T>(),
57  ESparse sparse = kDense);
58 
59  template<class T>
60  Spectrum(const std::string& label,
61  const Binning& bins,
62  SpectrumLoaderBase& loader,
63  const _Var<T>& var,
64  const _Cut<T>& cut,
65  const SystShifts& shift = kNoShift,
66  const _Var<T>& wei = Unweighted<T>(),
67  ESparse sparse = kDense);
68 
69  /// The only \ref MultiVar variant available
70  template<class T>
72  const _HistAxis<_MultiVar<T>>& axis,
73  const _Cut<T>& cut,
74  const SystShifts& shift = kNoShift,
75  const _Var<T>& wei = Unweighted<T>());
76 
77  /// Makes a spectrum from an eigen array
78  Spectrum(Eigen::ArrayXd&& h,
79  const LabelsAndBins& axis,
80  double pot, double livetime);
81 
82  /// Makes a spectrum from an eigen array of stan vars
84  const LabelsAndBins& axis,
85  double pot, double livetime);
86 
87  /// 2D Spectrum taking 2 HistAxis
88  template<class T>
90  const _HistAxis<_Var<T>>& xAxis,
91  const _HistAxis<_Var<T>>& yAxis,
92  const _Cut<T>& cut,
93  const SystShifts& shift = kNoShift,
94  const _Var<T>& wei = Unweighted<T>(),
95  ESparse sparse = kDense);
96 
97  /// 2D Spectrum of two Vars
98  template<class T>
99  Spectrum(const std::string& xLabel,
100  const std::string& yLabel,
101  SpectrumLoaderBase& loader,
102  const Binning& binsx, const _Var<T>& varx,
103  const Binning& binsy, const _Var<T>& vary,
104  const _Cut<T>& cut,
105  const SystShifts& shift = kNoShift,
106  const _Var<T>& wei = Unweighted<T>(),
107  ESparse sparse = kDense);
108 
109  /// 3D Spectrum taking 3 HistAxis
110  template<class T>
112  const _HistAxis<_Var<T>>& xAxis,
113  const _HistAxis<_Var<T>>& yAxis,
114  const _HistAxis<_Var<T>>& zAxis,
115  const _Cut<T>& cut,
116  const SystShifts& shift = kNoShift,
117  const _Var<T>& wei = Unweighted<T>(),
118  ESparse sparse = kDense);
119 
120  /// 3D Spectrum of three Vars
121  template<class T>
122  Spectrum(const std::string& xLabel,
123  const std::string& yLabel,
124  const std::string& zLabel,
125  SpectrumLoaderBase& loader,
126  const Binning& binsx, const _Var<T>& varx,
127  const Binning& binsy, const _Var<T>& vary,
128  const Binning& binsz, const _Var<T>& varz,
129  const _Cut<T>& cut,
130  const SystShifts& shift = kNoShift,
131  const _Var<T>& wei = Unweighted<T>(),
132  ESparse sparse = kDense);
133 
134  /// Expert constructor for ReweightableSpectrum et al
136  const LabelsAndBins& axis,
137  double pot,
138  double livetime)
139  : fHist(std::move(hist)), fPOT(pot), fLivetime(livetime), fAxis(axis)
140  {
141  }
142 
143  /// The only valid thing to do with such a spectrum is to assign something
144  /// else into it.
145  static Spectrum Uninitialized(){return Spectrum();}
146 
147  virtual ~Spectrum();
148 
149  Spectrum(const Spectrum& rhs);
150  Spectrum(Spectrum&& rhs);
151  Spectrum& operator=(const Spectrum& rhs);
152  Spectrum& operator=(Spectrum&& rhs);
153 
154  void Fill(double x, double w = 1);
155 
156  /// \brief Histogram made from this Spectrum, scaled to some exposure
157  ///
158  /// \param exposure POT or livetime (seconds)
159  /// \param col Histogram color (default black)
160  /// \param style Histogram line style (default solid)
161  /// \param expotype How to interpret exposure (kPOT (default) or kLivetime)
162  TH1D* ToTH1(double exposure,
163  Color_t col = kBlack,
164  Style_t style = kSolid,
165  EExposureType expotype = kPOT,
166  EBinType bintype = kBinContent) const;
167 
168  /// \brief Histogram made from this Spectrum, scaled to some exposure
169  ///
170  /// \param exposure POT or livetime (seconds)
171  /// \param expotype How to interpret exposure (kPOT (default) or kLivetime)
172  TH1D* ToTH1(double exposure,
173  EExposureType expotype,
174  EBinType bintype = kBinContent) const;
175 
176  /// Spectrum must be 2D to obtain TH2
177  TH2* ToTH2 (double exposure, EExposureType expotype = kPOT,
178  EBinType bintype = kBinContent) const;
179 
180  /// Spectrum must be 3D to obtain TH3
181  TH3* ToTH3 (double exposure, EExposureType expotype = kPOT,
182  EBinType bintype = kBinContent) const;
183 
184  TH1* ToTHX(double exposure, EExposureType expotype = kPOT,
185  EBinType bintype = kBinContent) const;
186 
187  bool HasStan() const {return fHist.HasStan();}
188  /// NB these don't have POT scaling. For expert high performance ops only!
189  const Eigen::ArrayXd& GetEigen() const {return fHist.GetEigen();}
191 
192  Eigen::ArrayXd GetEigen(double exposure, EExposureType expotype = kPOT) const;
193  Eigen::ArrayXstan GetEigenStan(double exposure, EExposureType expotype = kPOT) const;
194 
195  /// \brief Return total number of events scaled to \a pot
196  ///
197  /// \param exposure POT/livetime to scale to
198  /// \param err The statistical error on this total (optional)
199  /// \param expotype What the first parameter represents
200  double Integral(double exposure, double* err = 0,
201  EExposureType expotype = kPOT) const;
202 
203  /// \brief Return mean of 1D histogram
204  double Mean() const;
205 
206  /// \brief Mock data is \ref FakeData with Poisson fluctuations applied
207  ///
208  /// Use for low-budget MDCs, or just getting a sense of the expected scale
209  /// of statistical variation. NB seed = 0 is true random
210  Spectrum MockData(double pot, int seed = 0) const;
211 
212  /// \brief Asimov data is a MC spectrum scaled to the POT expected in the
213  /// data
214  ///
215  /// Use for sensitivity plots and testing fit convergence
216  Spectrum AsimovData(double pot) const;
217 
218  /// Use for sensitivity plots when fake cosmic data is needed.
219  /// Fake cosmic spectra can be added to FakeData by desired livetime.
220  Spectrum AsimovData(double pot, double livetime) const;
221 
222  /// Synonymous with AsimovData(). Retained for compatibility
223  Spectrum FakeData(double pot) const;
224  /// Synonymous with AsimovData(). Retained for compatibility
225  Spectrum FakeData(double pot, double livetime) const;
226 
227  double POT() const {return fPOT;}
228 
229  /// Seconds. For informational purposes only. No calculations use this.
230  double Livetime() const {return fLivetime;}
231 
232  /// DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
233  void OverridePOT(double newpot) {fPOT = newpot;}
234 
235  /// DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
236  void OverrideLivetime(double newlive) {fLivetime = newlive;}
237 
238  void Clear();
239 
240  /// Multiply this spectrum by a constant c
241  void Scale(double c);
242  void Scale(const stan::math::var& v);
243 
244  // Arithmetic operators are as if these are unlike samples, each a
245  // contribution to one total, not seperate sources of stats for the same
246  // sample.
247  Spectrum& operator+=(const Spectrum& rhs);
248  Spectrum operator+(const Spectrum& rhs) const;
249 
250  Spectrum& operator-=(const Spectrum& rhs);
251  Spectrum operator-(const Spectrum& rhs) const;
252 
253  Spectrum& operator*=(const Ratio& rhs);
254  Spectrum operator*(const Ratio& rhs) const;
255 
256  Spectrum& operator/=(const Ratio& rhs);
257  Spectrum operator/(const Ratio& rhs) const;
258 
259  void SaveTo(TDirectory* dir, const std::string& name) const;
260  static std::unique_ptr<Spectrum> LoadFrom(TDirectory* dir, const std::string& name);
261 
262  unsigned int NDimensions() const{return fAxis.NDimensions();}
263  const std::vector<std::string>& GetLabels() const {return fAxis.GetLabels();}
264  const std::vector<Binning>& GetBinnings() const {return fAxis.GetBinnings();}
265 
266  protected:
267  /// Constructor for Uninitialized()
269  : fHist(Hist::Uninitialized()),
270  fPOT(0), fLivetime(0),
271  fAxis(std::vector<std::string>(), std::vector<Binning>())
272  {
273  }
274 
275  /// Helper for constructors
276  Spectrum(const LabelsAndBins& axis, ESparse sparse = kDense);
277 
278  void RemoveLoader(Spectrum**);
279  void AddLoader(Spectrum**);
280 
281  /// Helper for operator+= and operator-=
282  Spectrum& PlusEqualsHelper(const Spectrum& rhs, int sign);
283 
285  double fPOT;
286  double fLivetime;
287 
288  /// Things that point at this Spectrum. Maintained by SpectrumLoader
289  std::set<Spectrum**> fReferences;
290 
292  };
293 
294  // Commutative
295  inline Spectrum operator*(const Ratio& lhs, const Spectrum& rhs){return rhs*lhs;}
296 }
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:264
const XML_Char * name
Definition: expat.h:151
const Eigen::ArrayXd & GetEigen() const
Definition: Hist.h:53
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Spectrum operator*(const Ratio &rhs) const
Definition: Spectrum.cxx:483
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
void AddLoader(Spectrum **)
Definition: Spectrum.cxx:373
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
LabelsAndBins fAxis
Definition: Spectrum.h:291
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:233
Spectrum & PlusEqualsHelper(const Spectrum &rhs, int sign)
Helper for operator+= and operator-=.
Definition: Spectrum.cxx:379
Spectrum & operator+=(const Spectrum &rhs)
Definition: Spectrum.cxx:448
void RemoveLoader(Spectrum **)
Definition: Spectrum.cxx:367
Definition: Hist.h:29
std::set< Spectrum ** > fReferences
Things that point at this Spectrum. Maintained by SpectrumLoader.
Definition: Spectrum.h:289
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:189
void Clear()
Definition: Spectrum.cxx:361
double fPOT
Definition: Spectrum.h:285
void Fill(double x, double w=1)
Definition: Spectrum.cxx:293
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
Spectrum & operator-=(const Spectrum &rhs)
Definition: Spectrum.cxx:462
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
TH1 * ToTHX(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Definition: Spectrum.cxx:211
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
Spectrum(Hist &&hist, const LabelsAndBins &axis, double pot, double livetime)
Expert constructor for ReweightableSpectrum et al.
Definition: Spectrum.h:135
EBinType
Definition: UtilsExt.h:28
static Spectrum Uninitialized()
Definition: Spectrum.h:145
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
unsigned int seed
Definition: runWimpSim.h:102
Spectrum & operator/=(const Ratio &rhs)
Definition: Spectrum.cxx:491
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
Spectrum operator/(const Ratio &rhs) const
Definition: Spectrum.cxx:498
Int_t col[ntarg]
Definition: Style.C:29
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: UtilsExt.h:35
Spectrum AsimovData(double pot) const
Asimov data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:318
#define pot
Spectrum & operator=(const Spectrum &rhs)
Definition: Spectrum.cxx:63
Spectrum()
Constructor for Uninitialized()
Definition: Spectrum.h:268
double fLivetime
Definition: Spectrum.h:286
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
void OverrideLivetime(double newlive)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:236
unsigned int NDimensions() const
Definition: Spectrum.h:262
loader
Definition: demo0.py:10
Regular histogram.
Definition: UtilsExt.h:30
double POT() const
Definition: Spectrum.h:227
Represent the ratio between two spectra.
Definition: Ratio.h:8
const SystShifts kNoShift
Definition: SystShifts.cxx:22
Base class for the various types of spectrum loader.
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Cut cut
Definition: exporter_fd.C:30
unsigned int NDimensions() const
Definition: LabelsAndBins.h:66
Spectrum operator+(const Spectrum &rhs) const
Definition: Spectrum.cxx:454
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:187
Collect information describing the x-axis of an analysis histogram.
Definition: HistAxis.h:18
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: StanUtils.h:7
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
_Var< T > Unweighted()
Definition: Var.h:93
virtual ~Spectrum()
Definition: Spectrum.cxx:78
Template for Cut and SpillCut.
Definition: Cut.h:15
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:263
Template for Var and SpillVar.
Float_t w
Definition: plot.C:20
bool HasStan() const
Definition: Spectrum.h:187
Spectrum operator-(const Spectrum &rhs) const
Definition: Spectrum.cxx:468
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
def sign(x)
Definition: canMan.py:197
Spectrum & operator*=(const Ratio &rhs)
Definition: Spectrum.cxx:476
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Spectrum.h:190
bool HasStan() const
Definition: Hist.h:52
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Hist.h:54
double Mean() const
Return mean of 1D histogram.
Definition: Spectrum.cxx:267
enum BeamMode string