Spectrum.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <memory>
4 #include <string>
5 #include <vector>
6 
7 #include "TAttLine.h"
8 
9 #include "CAFAna/Core/Var.h"
10 #include "CAFAna/Core/Cut.h"
11 #include "CAFAna/Core/HistAxis.h"
12 #include "CAFAna/Core/Utilities.h" // EExposureType et al
15 #include "CAFAna/Core/Var.h"
16 
17 class TDirectory;
18 class TH1;
19 class TH2;
20 class TH3;
21 class TH1D;
22 class TH1F;
23 
24 template<class T> class THnSparseT;
26 
27 /// Oscillation analysis framework, runs over CAF files outside of ART
28 namespace ana
29 {
30  class Binning;
31  class MultiVar;
32  class Ratio;
33 
34  class SpectrumStan;
35 
36  /// Representation of a spectrum in any variable, with associated POT
37  class Spectrum
38  {
39  public:
40  friend class SpectrumStan;
41  friend class SpectrumLoaderBase;
42  friend class SpectrumLoader;
43  friend class NullLoader;
44 
46 
47  Spectrum(const std::string& label, const Binning& bins,
49  const Var& var,
50  const Cut& cut,
51  const SystShifts& shift = kNoShift,
52  const Var& wei = kUnweighted);
53 
54  Spectrum(const std::string& label, const Binning& bins,
55  SpectrumLoaderBase& loader,
56  const MultiVar& var,
57  const Cut& cut,
58  const SystShifts& shift = kNoShift,
59  const Var& wei = kUnweighted);
60 
61  Spectrum(const std::string& label, const Binning& bins,
62  SpectrumLoaderBase& loader,
63  const NuTruthVar& var,
64  const NuTruthCut& cut,
65  const SystShifts& shifts = kNoShift,
66  const NuTruthVar& wei = kNuTruthUnweighted);
67 
68  Spectrum(const std::string& label, SpectrumLoaderBase& loader,
69  const Binning& binsx, const NuTruthVar& varx,
70  const Binning& binsy, const NuTruthVar& vary,
71  const NuTruthCut& cut,
72  const SystShifts& shifts = kNoShift,
73  const NuTruthVar& wei = kNuTruthUnweighted);
74 
76  const HistAxis& axis,
77  const Cut& cut,
78  const SystShifts& shift = kNoShift,
79  const Var& wei = kUnweighted);
80 
82  const MultiVarHistAxis& axis,
83  const Cut& cut,
84  const SystShifts& shift = kNoShift,
85  const Var& wei = kUnweighted);
86 
88  const NuTruthHistAxis& axis,
89  const NuTruthCut& cut,
90  const SystShifts& shifts = kNoShift,
92 
93  Spectrum(const std::string& label, const Binning& bins, ESparse sparse = kDense);
94  Spectrum(const std::string& label, double pot, double livetime, const Binning& bins);
95 
96  /// Copies \a h
97  Spectrum(TH1* h,
98  const std::vector<std::string>& labels,
99  const std::vector<Binning>& bins,
100  double pot, double livetime);
101 
102  // Special-case for 1D histogram
103  Spectrum(TH1* h, double pot, double livetime);
104 
105  /// Takes possession of \a h
106  Spectrum(std::unique_ptr<TH1D> h,
107  const std::vector<std::string>& labels,
108  const std::vector<Binning>& bins,
109  double pot, double livetime);
110 
111  Spectrum(std::unique_ptr<TH1F> h,
112  const std::vector<std::string>& labels,
113  const std::vector<Binning>& bins,
114  double pot, double livetime);
115 
116  /// 2D Spectrum of two Vars
117  Spectrum(const std::string& label, SpectrumLoaderBase& loader,
118  const Binning& binsx, const Var& varx,
119  const Binning& binsy, const Var& vary,
120  const Cut& cut,
121  const SystShifts& shift = kNoShift,
122  const Var& wei = kUnweighted);
123 
124  /// 2D Spectrum of two MultiVars
125  Spectrum(const std::string& label,
126  SpectrumLoaderBase& loader,
127  const Binning& binsx, const MultiVar& varx,
128  const Binning& binsy, const MultiVar& vary,
129  const Cut& cut,
130  const SystShifts& shift = kNoShift,
131  const Var& wei = kUnweighted);
132 
133  /// 2D Spectrum taking 2 HistAxis
135  const HistAxis& xAxis,
136  const HistAxis& yAxis,
137  const Cut& cut,
138  const SystShifts& shift = kNoShift,
139  const Var& wei = kUnweighted);
140 
141  Spectrum(const std::string& xLabel,
142  const std::string& yLabel,
143  SpectrumLoaderBase& loader,
144  const Binning& binsx, const Var& varx,
145  const Binning& binsy, const Var& vary,
146  const Cut& cut,
147  const SystShifts& shift = kNoShift,
148  const Var& wei = kUnweighted);
149 
150  Spectrum(const std::string& xLabel,
151  const std::string& yLabel,
152  SpectrumLoaderBase& loader,
153  const Binning& binsx, const MultiVar& varx,
154  const Binning& binsy, const MultiVar& vary,
155  const Cut& cut,
156  const SystShifts& shift = kNoShift,
157  const Var& wei = kUnweighted);
158 
159  /// 3D Spectrum of three Vars
160  Spectrum(const std::string& label, SpectrumLoaderBase& loader,
161  const Binning& binsx, const Var& varx,
162  const Binning& binsy, const Var& vary,
163  const Binning& binsz, const Var& varz,
164  const Cut& cut,
165  const SystShifts& shift = kNoShift,
166  const Var& wei = kUnweighted,
167  ESparse sparse = kDense);
168 
169  Spectrum(const std::string& xLabel,
170  const std::string& yLabel,
171  const std::string& zLabel,
172  SpectrumLoaderBase& loader,
173  const Binning& binsx, const Var& varx,
174  const Binning& binsy, const Var& vary,
175  const Binning& binsz, const Var& varz,
176  const Cut& cut,
177  const SystShifts& shift = kNoShift,
178  const Var& wei = kUnweighted,
179  ESparse sparse = kDense);
180 
181  /// 3D Spectrum taking 3 HistAxis
183  const HistAxis& xAxis,
184  const HistAxis& yAxis,
185  const HistAxis& zAxis,
186  const Cut& cut,
187  const SystShifts& shift = kNoShift,
188  const Var& wei = kUnweighted,
189  ESparse sparse = kDense);
190 
191  virtual ~Spectrum();
192 
193  Spectrum(const Spectrum& rhs);
194  Spectrum(Spectrum&& rhs);
195  Spectrum& operator=(const Spectrum& rhs);
196  Spectrum& operator=(Spectrum&& rhs);
197 
198  /// Copy conversion from SpectrumStan
199  Spectrum(const SpectrumStan& rhs);
200  Spectrum& operator=(const SpectrumStan& rhs);
201 
202 
203  void Fill(double x, double w = 1);
204 
205  /// \brief Histogram made from this Spectrum, scaled to some exposure
206  ///
207  /// \param exposure POT or livetime (seconds)
208  /// \param col Histogram color (default black)
209  /// \param style Histogram line style (default solid)
210  /// \param expotype How to interpret exposure (kPOT (default) or kLivetime)
211  TH1D* ToTH1(double exposure,
212  Color_t col = kBlack,
213  Style_t style = kSolid,
214  EExposureType expotype = kPOT,
215  EBinType bintype = kBinContent) const;
216 
217  /// \brief Histogram made from this Spectrum, scaled to some exposure
218  ///
219  /// \param exposure POT or livetime (seconds)
220  /// \param expotype How to interpret exposure (kPOT (default) or kLivetime)
221  TH1D* ToTH1(double exposure,
222  EExposureType expotype,
223  EBinType bintype = kBinContent) const
224  {
225  return ToTH1(exposure, kBlack, kSolid, expotype, bintype);
226  }
227 
228  /// Spectrum must be 2D to obtain TH2
229  TH2* ToTH2 (double exposure, EExposureType expotype = kPOT,
230  EBinType bintype = kBinContent) const;
231  /// Spectrum must be 2D to obtain TH2. Normalized to X axis.
232  TH2* ToTH2NormX(double exposure, EExposureType expotype = kPOT) const;
233 
234  /// Spectrum must be 3D to obtain TH3
235  TH3* ToTH3 (double exposure, EExposureType expotype = kPOT,
236  EBinType bintype = kBinContent) const;
237 
238  /// \brief Return total number of events scaled to \a pot
239  ///
240  /// \param exposure POT/livetime to scale to
241  /// \param err The statistical error on this total (optional)
242  /// \param expotype What the first parameter represents
243  double Integral(double exposure, double* err = 0,
244  EExposureType expotype = kPOT) const;
245 
246  /// \brief Return mean of 1D histogram
247  double Mean() const;
248 
249  /// \brief Mock data is \ref FakeData with Poisson fluctuations applied
250  ///
251  /// Use for low-budget MDCs, or just getting a sense of the expected scale
252  /// of statistical variation
253  Spectrum MockData(double pot, int idx = 0) const;
254  /// \brief Fake data is a MC spectrum scaled to the POT expected in the data
255  ///
256  /// Use for sensitivity plots and testing fit convergence
257  Spectrum FakeData(double pot) const;
258 
259  /// Use for sensitivity plots when fake cosmic data is needed.
260  /// Fake cosmic spectra can be added to FakeData by desired livetime.
261  Spectrum FakeData(double pot, double livetime) const;
262 
263  double POT() const {return fPOT;}
264 
265  /// Seconds. For informational purposes only. No calculations use this.
266  double Livetime() const {return fLivetime;}
267 
268  /// DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
269  void OverridePOT(double newpot) {fPOT = newpot;}
270 
271  /// DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
272  void OverrideLivetime(double newlive) {fLivetime = newlive;}
273 
274  void Clear();
275 
276  /// Multiply this spectrum by a constant c
277  void Scale(double c);
278 
279  // Arithmetic operators are as if these are unlike samples, each a
280  // contribution to one total, not seperate sources of stats for the same
281  // sample.
282  Spectrum& operator+=(const Spectrum& rhs);
283  Spectrum operator+(const Spectrum& rhs) const;
284 
285  Spectrum& operator-=(const Spectrum& rhs);
286  Spectrum operator-(const Spectrum& rhs) const;
287 
288  Spectrum& operator*=(const Ratio& rhs);
289  Spectrum operator*(const Ratio& rhs) const;
290 
291  Spectrum& operator/=(const Ratio& rhs);
292  Spectrum operator/(const Ratio& rhs) const;
293 
294  void SaveTo(TDirectory* dir) const;
295  static std::unique_ptr<Spectrum> LoadFrom(TDirectory* dir);
296 
297  unsigned int NDimensions() const{return fLabels.size();}
298  std::vector<std::string> GetLabels() const {return fLabels;}
299  std::vector<Binning> GetBinnings() const {return fBins;}
300 
301  protected:
302  Spectrum(const std::vector<std::string>& labels,
303  const std::vector<Binning>& bins,
304  ESparse sparse = kDense);
305 
306  Binning GetMultiDBinning() const;
307 
308  void ConstructHistogram(ESparse sparse = kDense);
309 
312 
313  /// Helper for operator+= and operator-=
314  Spectrum& PlusEqualsHelper(const Spectrum& rhs, int sign);
315 
316  TH1D* fHistD;
317  TH1F* fHistF;
319  double fPOT;
320  double fLivetime;
321 
322  /// This count is maintained by SpectrumLoader, as a sanity check
323  std::set<SpectrumLoaderBase*> fLoaderCount;
324 
325  std::vector<std::string> fLabels;
326  std::vector<Binning> fBins;
327  };
328 
329  // Commutative
330  inline Spectrum operator*(const Ratio& lhs, const Spectrum& rhs){return rhs*lhs;}
331 
332 }
TH1D * fHistD
Definition: Spectrum.h:316
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:652
Spectrum operator*(const Ratio &rhs) const
Definition: Spectrum.cxx:1016
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
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:564
const Var weight
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:269
Spectrum & PlusEqualsHelper(const Spectrum &rhs, int sign)
Helper for operator+= and operator-=.
Definition: Spectrum.cxx:905
Spectrum & operator+=(const Spectrum &rhs)
Definition: Spectrum.cxx:980
TH1F * fHistF
Definition: Spectrum.h:317
TH2 * ToTH2NormX(double exposure, EExposureType expotype=kPOT) const
Spectrum must be 2D to obtain TH2. Normalized to X axis.
Definition: Spectrum.cxx:682
void AddLoader(SpectrumLoaderBase *)
Definition: Spectrum.cxx:901
void ConstructHistogram(ESparse sparse=kDense)
Definition: Spectrum.cxx:524
void Clear()
Definition: Spectrum.cxx:889
double fPOT
Definition: Spectrum.h:319
void Fill(double x, double w=1)
Definition: Spectrum.cxx:785
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:753
Spectrum & operator-=(const Spectrum &rhs)
Definition: Spectrum.cxx:994
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:745
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const char * label
std::vector< std::string > fLabels
Definition: Spectrum.h:325
EBinType
Definition: Utilities.h:42
std::set< SpectrumLoaderBase * > fLoaderCount
This count is maintained by SpectrumLoader, as a sanity check.
Definition: Spectrum.h:323
THnSparseT< TArrayD > THnSparseD
Definition: Spectrum.h:24
Spectrum & operator/=(const Ratio &rhs)
Definition: Spectrum.cxx:1024
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:830
Spectrum operator/(const Ratio &rhs) const
Definition: Spectrum.cxx:1032
Int_t col[ntarg]
Definition: Style.C:29
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: Utilities.h:50
#define pot
Spectrum & operator=(const Spectrum &rhs)
Definition: Spectrum.cxx:426
double fLivetime
Definition: Spectrum.h:320
void OverrideLivetime(double newlive)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:272
unsigned int NDimensions() const
Definition: Spectrum.h:297
loader
Definition: demo0.py:10
void RemoveLoader(SpectrumLoaderBase *)
Definition: Spectrum.cxx:897
Regular histogram.
Definition: Utilities.h:44
double POT() const
Definition: Spectrum.h:263
THnSparseD * fHistSparse
Definition: Spectrum.h:318
Represent the ratio between two spectra.
Definition: Ratio.h:8
const SystShifts kNoShift
Definition: SystShifts.h:112
Base class for the various types of spectrum loader.
Binning GetMultiDBinning() const
Definition: Spectrum.cxx:510
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:798
Spectrum operator+(const Spectrum &rhs) const
Definition: Spectrum.cxx:986
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:713
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
std::vector< Binning > fBins
Definition: Spectrum.h:326
const NuTruthVar kNuTruthUnweighted
Definition: Var.h:104
Dummy loader that doesn&#39;t load any files.
TH1D * ToTH1(double exposure, EExposureType expotype, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.h:221
std::vector< std::string > GetLabels() const
Definition: Spectrum.h:298
virtual ~Spectrum()
Definition: Spectrum.cxx:324
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1066
Float_t w
Definition: plot.C:20
Spectrum operator-(const Spectrum &rhs) const
Definition: Spectrum.cxx:1000
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:100
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299
def sign(x)
Definition: canMan.py:197
Spectrum & operator*=(const Ratio &rhs)
Definition: Spectrum.cxx:1008
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1040
double Mean() const
Return mean of 1D histogram.
Definition: Spectrum.cxx:775
Spectrum(const std::string &label, const Binning &bins, SpectrumLoaderBase &loader, const Var &var, const Cut &cut, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
Definition: Spectrum.cxx:47