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