Utilities.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <cassert>
4 #include <cmath>
5 #include <cxxabi.h>
6 #include <cfenv>
7 #include <map>
8 #include <memory>
9 #include <set>
10 #include <string>
11 #include <vector>
12 
13 #include "StandardRecord/SREnums.h"
14 #include "CAFAna/Core/FwdDeclare.h"
15 
16 template<class T> class TMatrixT;
18 template<class T> class TVectorT;
20 
21 class TArrayD;
22 class TDirectory;
23 class TH1;
24 class TH2;
25 class TH3;
26 class TF1;
27 class TH1D;
28 class TH1F;
29 class TH2F;
30 class TH2D;
31 class TH3D;
32 class TRandom;
33 class TVector3;
34 
35 namespace caf
36 {
37  class SRProxy;
38  class SRSpillProxy;
39 }
40 
41 namespace ana
42 {
43  class Binning;
44  class Ratio;
45 
46  enum EBinType
47  {
48  kBinContent, ///< Regular histogram
49  kBinDensity ///< Divide bin contents by bin widths
50  };
51 
52 
53  /// For use as an argument to \ref Spectrum::ToTH1
57  };
58 
59 
60  /// Return a different string each time, for creating histograms
62 
63  /// \brief Prevent histograms being added to the current directory
64  ///
65  /// Upon going out of scope, restores the previous setting
67  {
68  public:
71  protected:
72  bool fBackup;
73  };
74 
75  /// \brief ifdh calls between construction and destruction produce no output
76  ///
77  /// Upon going out of scope, restores the previous setting
78  class IFDHSilent
79  {
80  public:
81  IFDHSilent();
82  ~IFDHSilent();
83  protected:
84  bool fSet;
85  };
86 
87  /// \brief Alter floating-point exception flag
88  ///
89  /// Upon going out of scope, restores the previous setting
91  {
92  public:
93  FloatingExceptionOnNaN(bool enable = true);
95  protected:
96  fexcept_t fBackup;
97  };
98 
99  /** \brief Compute bin-to-bin covariance matrix from a collection of sets of bin contents.
100 
101  \param binSets Collection of sets of bins from which covariances should be calculated
102  \param firstBin The first bin that should be considered (inclusive)
103  \param lastBin The last bin that should be considered (inclusive). -1 means "last in set"
104 
105  \returns unique_ptr to TMatrixD containing computed covariance matrix unless binSets.size() < 2,
106  in which case the unique_ptr's target is nullptr.
107 
108  Note TH1D is a child class of TArrayD -- so you can pass a vector
109  of TH1D* to this method.
110  **/
111  std::unique_ptr<TMatrixD> CalcCovMx(const std::vector<TArrayD*> & binSets, int firstBin=0, int lastBin=-1);
112 
113 
114  /** \brief The log-likelihood formula from the PDG.
115 
116  \param exp The expected spectrum
117  \param obs The corresponding observed spectrum
118 
119  \returns The log-likelihood formula from the PDG
120  \f[ \chi^2=2\sum_i^{\rm bins}\left(e_i-o_i+o_i\ln\left({o_i\over e_i}\right)\right) \f]
121 
122  Includes underflow bin and an option for
123  overflow bin (off by default) and handles
124  zero observed or expected events correctly.
125  **/
126  double LogLikelihood(const TH1* exp, const TH1* obs, bool useOverflow = false);
127 
128  /** \brief The log-likelihood formula for a single bin
129 
130  \param exp Expected count
131  \param obs Observed count
132 
133  \returns The log-likelihood formula from the PDG
134  \f[ \chi^2=2\left(e-o+o\ln\left({o\over e}\right)\right) \f]
135 
136  Handles zero observed or expected events correctly.
137 
138  Templated so that it can handle usage with Stan vars and other numeric types.
139  (The horible third template parameter ensures this function can only be used
140  with types that accept conversion from double -- which means you can't pass
141  it a TH1* or a std::vector<stan::math::var>&, removing the ambiguity with
142  those other versions of LogLikelihood.)
143  **/
144  template <typename T, typename U,
146  auto LogLikelihood(T exp, U obs)
147  {
148  // With this value, negative expected events and one observed
149  // event gives a chisq from this one bin of 182.
150  const double minexp = 1e-40; // Don't let expectation go lower than this
151 
152  // this turns into a stan::var if any of its types are, otherwise it gives back the primitive type
153  decltype(exp - obs) chi = 0;
154 
155  assert(obs >= 0);
156  if(exp < minexp){
157  if(!obs) return chi;
158  exp = minexp;
159  }
160 
161  if(obs)
162  // This strange form is for numerical stability when e~o
163  chi += 2*obs*((exp-obs)/obs + log1p((obs-exp)/exp));
164  else
165  chi += 2*(exp-obs);
166 
167  return chi;
168  }
169 
170 
171  double LogLikelihoodDerivative(double e, double o, double dedx);
172 
173  double LogLikelihoodDerivative(const TH1D* eh, const TH1D* oh,
174  const std::vector<double>& dedx);
175 
176  /** \brief Chi-squared calculation using a covariance matrix.
177 
178  \param exp Expected bin counts
179  \param obs Observed bin counts
180  \param covmxinv Inverse of covariance matrix. Must have same dimensions as exp and obs
181 
182  \returns The chi^2 calculated according to the formula from the PDG:
183  \f[ \chi^2 = (\vec{o}-\vec{e})^{T} V^{-1} (\vec{o} - \vec{e}) \]
184 
185  Note that this implicitly assumes Gaussian statistics for the bin counts!
186  **/
187  double Chi2CovMx(const TVectorD* exp, const TVectorD* obs, const TMatrixD* covmxinv);
188 
189  /// Chi-squared calculation using covariance matrix (calls the TVectorD version internally).
190  double Chi2CovMx(const TH1* exp, const TH1* obs, const TMatrixD* covmxinv);
191 
192  /// \brief Internal helper for \ref Surface and \ref FCSurface
193  ///
194  /// Creates a histogram having bins \em centred at the min and max
195  /// coordinates
196  TH2F* ExpandedHistogram(const std::string& title,
197  int nbinsx, double xmin, double xmax,
198  int nbinsy, double ymin, double ymax);
199 
200  /// \brief Invert a symmetric matrix with possibly empty rows/columns.
201  ///
202  /// Invert a symmetric matrix that may have empty rows/columns,
203  /// which (strictly speaking) make it impossible to invert the matrix.
204  /// (This often arises when computing covariance matrices for predictions
205  /// which have empty bins; the covariance is 0 for the entire row/column
206  /// in that case.)
207  /// Since those rows/cols are not useful, we can sidestep the problem
208  /// by removing them (and the corresponding columns)
209  /// from the matrix, inverting that, then re-inserting
210  /// the null rows/columns.
211  std::unique_ptr<TMatrixD> SymmMxInverse(const TMatrixD& mx);
212 
213  /// Utility function to avoid need to switch on bins.IsSimple()
214  TH1D* MakeTH1D(const char* name, const char* title, const Binning& bins);
215  TH1F* MakeTH1F(const char* name, const char* title, const Binning& bins);
216  /// Utility function to avoid 4-way combinatorial explosion on the bin types
217  TH2D* MakeTH2D(const char* name, const char* title,
218  const Binning& binsx,
219  const Binning& binsy);
220  TH2F* MakeTH2F(const char* name, const char* title,
221  const Binning& binsx,
222  const Binning& binsy);
223  /// Utility function for 3D hist
224  // TODO: make bins more simple-like
225  TH3D* MakeTH3D(const char* name, const char* title,
226  const Binning& binsx,
227  const Binning& binsy,
228  const Binning& binsz);
229 
230  /// \brief For use with \ref Var2D
231  ///
232  /// Re-expand a histogram flattened by \ref Var2D into a 2D histogram for
233  /// plotting purposes. The binning scheme must match that used in the
234  /// original Var.
235  TH2* ToTH2(const Spectrum& s, double exposure, ana::EExposureType expotype,
236  const Binning& binsx, const Binning& binsy,
238 
239  /// Same as ToTH2, but with 3 dimensions
240  TH3* ToTH3(const Spectrum& s, double exposure, ana::EExposureType expotype,
241  const Binning& binsx, const Binning& binsy, const Binning& binsz,
243 
244  /// \brief For use with \ref Var2D
245  ///
246  /// Re-expand a flatenned histogram into a 2D histogram for
247  /// plotting purposes. The binning scheme must match that used in the
248  /// original Var.
249  TH2* ToTH2(const Ratio& r, const Binning& binsx, const Binning& binsy);
250 
251  /// Same as ToTH2, but with 3 dimensions
252  TH3* ToTH3(const Ratio& r, const Binning& binsx,
253  const Binning& binsy, const Binning& binsz);
254 
255  /// Helper for ana::ToTH2
256  TH2* ToTH2Helper(std::unique_ptr<TH1> h1,
257  const Binning& binsx,
258  const Binning& binsy,
260 
261  /// Helper for ana::ToTH3
262  TH3* ToTH3Helper(std::unique_ptr<TH1> h1,
263  const Binning& binsx,
264  const Binning& binsy,
265  const Binning& binsz,
267 
268  /// Find files matching a UNIX glob, plus expand environment variables
269  std::vector<std::string> Wildcard(const std::string& wildcardString);
270 
271  /// This is $SRT_PRIVATE_CONTEXT if you have used a test release,
272  /// otherwise $SRT_PUBLIC_CONTEXT
274 
275  /// This is $SRT_PRIVATE_CONTEXT if a CAFAna directory exists there,
276  /// otherwise $SRT_PUBLIC_CONTEXT
278 
279  /// Read list of input files from a text file, one per line
280  std::vector<std::string> LoadFileList(const std::string& listfile);
281 
282  /// \brief Extract map of metadata parameters from a CAF file
283  ///
284  /// \param dir The "meta" directory from the CAF file
285  /// \return A map from metadata field name to metadata value
286  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir);
287 
288  /// \brief \a base += \a add
289  ///
290  /// \param base The original source of strings, will be altered
291  /// \param add Strings to add to \a base if missing
292  /// \param mask Fields for which there was a conflict, will be altered
293  void CombineMetadata(std::map<std::string, std::string>& base,
294  const std::map<std::string, std::string>& add,
295  std::set<std::string>& mask);
296 
297  /// \brief Write map of metadata parameters into a CAF file
298  ///
299  /// \param dir The "meta" directory of the CAF file
300  /// \param meta Map from metadata field name to metadata value
301  void WriteCAFMetadata(TDirectory* dir,
302  const std::map<std::string, std::string>& meta);
303 
304 
305  /// \brief Average direction of NuMI neutrinos in a given detector
306  /// This function is a copy of
307  /// \ref geo::GeometryBase::NuMIBeamDirection ()
308  /// Any changes made to that function must be propagated to here.
309  TVector3 NuMIBeamDirection(caf::Det_t det);
310 
311  /// Is this a grid (condor) job?
312  bool RunningOnGrid();
313 
314  /// What's the process number for a grid job?
315  size_t JobNumber();
316 
317  bool SAMDefinitionExists(const std::string& def);
318 
319  bool AlmostEqual(double a, double b);
320 
321  /// utility method to figure out exactly what kind of object I am
322  template <typename T>
323  std::string DemangledTypeName(T* obj) { return abi::__cxa_demangle(typeid(*obj).name(), 0, 0, 0); }
324 
325  template <typename T>
326  constexpr char* DemangledTypeName(){ return abi::__cxa_demangle(typeid(T).name(), 0, 0, 0); }
327 
328  std::string pnfs2xrootd(std::string loc, bool unauth = false);
329 
330  // Calling this function will return a Fourier series, fit to the input
331  // histogram. Assumes x-axis covers one period
332  class FitToFourier
333  {
334  public:
335  FitToFourier(TH1* h, double xlo, double xhi, int NOsc);
336  ~FitToFourier();
337  TF1* Fit() const;
338  double operator()(double *x, double *par) const;
339  private:
340 
341  const TH1* fHist; // Histogram to fit
342  const double fxlo; // Lower bound
343  const double fxhi; // Upper bound - assumed to be 1 osc from the low end
344  const int fNOsc; // Highest harmonic to include
345 
346  };
347 
348  /// Generate a unique hash from a sequence of integers
349  std::size_t Hash(const std::vector<std::size_t> &vals);
350 
351 
352  /// Generate a unique hash for a SRSpill (run/subrun/event combination)
353  std::size_t Hash(const caf::SRSpillProxy * sr);
354 
355 
356  /// Generate a unique hash for a StandardRecord (run/subrun/cycle/event/slice combination)
357  std::size_t Hash(const caf::SRProxy * sr);
358 }
TH2D * MakeTH2D(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Definition: Utilities.cxx:345
Det_t
Which NOvA detector?
Definition: SREnums.h:7
const XML_Char * name
Definition: expat.h:151
size_t JobNumber()
Definition: Utilities.cxx:715
TH3D * MakeTH3D(const char *name, const char *title, const Binning &binsx, const Binning &binsy, const Binning &binsz)
Definition: Utilities.cxx:362
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
bool RunningOnGrid()
Definition: Utilities.cxx:708
TH3 * ToTH3Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Definition: Utilities.cxx:448
fVtxDx Fit("f")
TH2 * ToTH2Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:397
void CombineMetadata(std::map< std::string, std::string > &base, const std::map< std::string, std::string > &add, std::set< std::string > &mask)
Definition: Utilities.cxx:597
Divide bin contents by bin widths.
Definition: Utilities.h:49
std::vector< std::string > Wildcard(const std::string &wildcardString)
Definition: Utilities.cxx:502
Proxy for StandardRecord.
Definition: SRProxy.h:2237
std::unique_ptr< TMatrixD > SymmMxInverse(const TMatrixD &mx)
Definition: Utilities.cxx:219
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:741
ifdh calls between construction and destruction produce no output
Definition: Utilities.h:78
std::string FindCAFAnaDir()
Definition: Utilities.cxx:550
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
const Binning bins
Definition: sa_fd_arrays.h:120
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
::xsd::cxx::tree::type type
Definition: Database.h:110
Double_t ymax
Definition: plot.C:25
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Definition: Utilities.cxx:575
EBinType
Definition: Utilities.h:46
double LogLikelihoodDerivative(double e, double o, double dedx)
Definition: Utilities.cxx:154
const XML_Char * s
Definition: expat.h:262
const XML_Char int const XML_Char * value
Definition: expat.h:331
std::unique_ptr< TMatrixD > CalcCovMx(const std::vector< TArrayD * > &binSets, int firstBin, int lastBin)
Compute bin-to-bin covariance matrix from a collection of sets of bin contents.
Definition: Utilities.cxx:90
auto LogLikelihood(T exp, U obs)
The log-likelihood formula for a single bin.
Definition: Utilities.h:146
const double a
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: Utilities.h:54
const double kPOT
Definition: nue_pid_effs.C:26
Class for parsing *.fcl files for their metadata information.
Proxy for SRSpill.
Definition: SRProxy.h:1376
void WriteCAFMetadata(TDirectory *dir, const std::map< std::string, std::string > &meta)
Definition: Utilities.cxx:658
bool AlmostEqual(double a, double b)
Definition: Utilities.cxx:732
Regular histogram.
Definition: Utilities.h:48
std::string FindPackageDir(std::string Dir)
Definition: Utilities.cxx:530
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:374
Represent the ratio between two spectra.
Definition: Ratio.h:8
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Definition: Utilities.cxx:333
TH2F * MakeTH2F(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Definition: Utilities.cxx:353
TH1F * h1
UInt_t par
Definition: AnaPlotMaker.h:47
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::size_t Hash(const std::vector< std::size_t > &vals)
Definition: Utilities.cxx:772
TVectorT< double > TVectorD
Definition: Utilities.h:18
bool SAMDefinitionExists(const std::string &def)
Definition: Utilities.cxx:724
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
TVector3 NuMIBeamDirection(caf::Det_t det)
Definition: Utilities.cxx:676
Alter floating-point exception flag.
Definition: Utilities.h:90
Double_t ymin
Definition: plot.C:24
std::vector< std::string > LoadFileList(const std::string &listfile)
Definition: Utilities.cxx:556
double T
Definition: Xdiff_gwt.C:5
This module creates Common Analysis Files.
Definition: CSVMaker.h:8
TMatrixT< double > TMatrixD
Definition: Utilities.h:16
Prevent histograms being added to the current directory.
Definition: Utilities.h:66
Float_t e
Definition: plot.C:35
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Definition: Utilities.cxx:424
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TH1F * MakeTH1F(const char *name, const char *title, const Binning &bins)
Definition: Utilities.cxx:339
double Chi2CovMx(const TVectorD *e, const TVectorD *o, const TMatrixD *covmxinv)
Definition: Utilities.cxx:176
static constexpr Double_t sr
Definition: Munits.h:164
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax)
Definition: Utilities.cxx:198
h
Definition: demo3.py:41