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