Utilities.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/UtilsExt.h"
4 
5 #include <cassert>
6 #include <cmath>
7 #include <cxxabi.h>
8 #include <cfenv>
9 #include <map>
10 #include <memory>
11 #include <set>
12 #include <string>
13 #include <vector>
14 
15 #include "StandardRecord/SREnums.h"
16 #include "CAFAna/Core/FwdDeclare.h"
17 
18 template<class T> class TMatrixT;
20 template<class T> class TVectorT;
22 
23 class TArrayD;
24 class TDirectory;
25 class TH1;
26 class TH2;
27 class TH3;
28 class TF1;
29 class TH1D;
30 class TH1F;
31 class TH2F;
32 class TH2D;
33 class TH3D;
34 class TRandom;
35 class TVector3;
36 
38 
39 namespace ana
40 {
41  class Binning;
42  class Ratio;
43 
44  /// \brief ifdh calls between construction and destruction produce no output
45  ///
46  /// Upon going out of scope, restores the previous setting
47  class IFDHSilent
48  {
49  public:
50  IFDHSilent();
51  ~IFDHSilent();
52  protected:
53  bool fSet;
54  };
55 
56  /// \brief Alter floating-point exception flag
57  ///
58  /// Upon going out of scope, restores the previous setting
60  {
61  public:
62  FloatingExceptionOnNaN(bool enable = true);
64  protected:
65  fexcept_t fBackup;
66  };
67 
68  /** \brief Compute bin-to-bin covariance matrix from a collection of sets of bin contents.
69 
70  \param binSets Collection of sets of bins from which covariances should be calculated
71  \param firstBin The first bin that should be considered (inclusive)
72  \param lastBin The last bin that should be considered (inclusive). -1 means "last in set"
73 
74  \returns unique_ptr to TMatrixD containing computed covariance matrix unless binSets.size() < 2,
75  in which case the unique_ptr's target is nullptr.
76 
77  Note TH1D is a child class of TArrayD -- so you can pass a vector
78  of TH1D* to this method.
79  **/
80  std::unique_ptr<TMatrixD> CalcCovMx(const std::vector<TArrayD*> & binSets, int firstBin=0, int lastBin=-1);
81 
82  /** \brief The log-likelihood formula from the PDG.
83 
84  \param exp The expected spectrum
85  \param obs The corresponding observed spectrum
86 
87  \returns The log-likelihood formula from the PDG
88  \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]
89 
90  Includes underflow bin and an option for
91  overflow bin (off by default) and handles
92  zero observed or expected events correctly.
93  **/
94  double LogLikelihood(const TH1* exp, const TH1* obs, bool useOverflow = false);
95 
96  /** \brief The log-likelihood formula for a single bin
97 
98  \param exp Expected count
99  \param obs Observed count
100 
101  \returns The log-likelihood formula from the PDG
102  \f[ \chi^2=2\left(e-o+o\ln\left({o\over e}\right)\right) \f]
103 
104  Handles zero observed or expected events correctly.
105 
106  Templated so that it can handle usage with Stan vars and other numeric types.
107  (The horible third template parameter ensures this function can only be used
108  with types that accept conversion from double -- which means you can't pass
109  it a TH1* or a std::vector<stan::math::var>&, removing the ambiguity with
110  those other versions of LogLikelihood). The return type promotes to
111  stan::math::var if either T or U are.
112  **/
113  template <typename T, typename U,
114  typename std::enable_if_t<std::is_convertible_v<double, T> && std::is_convertible_v<double, U>, int> = 0>
115  decltype(T(0) - U(0)) LogLikelihood(T exp, U obs)
116  {
117  // With this value, negative expected events and one observed
118  // event gives a chisq from this one bin of 182.
119  const double minexp = 1e-40; // Don't let expectation go lower than this
120 
121  assert(obs >= 0);
122  if(exp < minexp){
123  if(!obs) return 0.;
124  exp = minexp;
125  }
126 
127  if(obs*1000 > exp){
128  // This strange form is for numerical stability when exp ~ obs
129  return 2*obs*((exp-obs)/obs + log1p((obs-exp)/exp));
130  }
131  else{
132  // But log1p doesn't like arguments near -1 (observation much smaller
133  // than expectation), and it's better to use the usual formula in that
134  // case.
135  if(obs){
136  return 2*(exp-obs + obs*log(obs/exp));
137  }
138  else{
139  return 2*exp;
140  }
141  }
142  }
143 
144  /// \brief Internal helper for \ref Surface and \ref FCSurface
145  ///
146  /// Creates a histogram having bins \em centred at the min and max
147  /// coordinates
148  TH2F* ExpandedHistogram(const std::string& title,
149  int nbinsx, double xmin, double xmax, bool xlog,
150  int nbinsy, double ymin, double ymax, bool ylog);
151 
152  /// This is $SRT_PRIVATE_CONTEXT if you have used a test release,
153  /// otherwise $SRT_PUBLIC_CONTEXT
155 
156  /// This is $SRT_PRIVATE_CONTEXT if a CAFAna directory exists there,
157  /// otherwise $SRT_PUBLIC_CONTEXT
159 
160  /// Read list of input files from a text file, one per line
161  std::vector<std::string> LoadFileList(const std::string& listfile);
162 
163  /// \brief Extract map of metadata parameters from a CAF file
164  ///
165  /// \param dir The "meta" directory from the CAF file
166  /// \return A map from metadata field name to metadata value
167  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir);
168 
169  /// \brief \a base += \a add
170  ///
171  /// \param base The original source of strings, will be altered
172  /// \param add Strings to add to \a base if missing
173  /// \param mask Fields for which there was a conflict, will be altered
174  void CombineMetadata(std::map<std::string, std::string>& base,
175  const std::map<std::string, std::string>& add,
176  std::set<std::string>& mask);
177 
178  /// \brief Write map of metadata parameters into a CAF file
179  ///
180  /// \param dir The "meta" directory of the CAF file
181  /// \param meta Map from metadata field name to metadata value
182  void WriteCAFMetadata(TDirectory* dir,
183  const std::map<std::string, std::string>& meta);
184 
185 
186  /// \brief Average direction of NuMI neutrinos in a given detector
187  /// This function is a copy of
188  /// \ref geo::GeometryBase::NuMIBeamDirection ()
189  /// Any changes made to that function must be propagated to here.
190  TVector3 NuMIBeamDirection(caf::Det_t det);
191 
192  /// Is this a grid (condor) job?
193  bool RunningOnGrid();
194 
195  /// What's the process number for a grid job?
196  size_t JobNumber();
197 
198  bool SAMDefinitionExists(const std::string& def);
199 
200  // Calling this function will return a Fourier series, fit to the input
201  // histogram. Assumes x-axis covers one period
203  {
204  public:
205  FitToFourier(TH1* h, double xlo, double xhi, int NOsc);
206  ~FitToFourier();
207  TF1* Fit() const;
208  double operator()(double *x, double *par) const;
209  private:
210 
211  const TH1* fHist; // Histogram to fit
212  const double fxlo; // Lower bound
213  const double fxhi; // Upper bound - assumed to be 1 osc from the low end
214  const int fNOsc; // Highest harmonic to include
215 
216  };
217 
218  /// Generate a unique hash from a sequence of integers
219  std::size_t Hash(const std::vector<std::size_t> &vals);
220 
221 
222  /// Generate a unique hash for a SRSpill (run/subrun/event combination)
223  std::size_t Hash(const caf::SRSpillProxy * sr);
224 
225 
226  /// Generate a unique hash for a StandardRecord (run/subrun/cycle/event/slice combination)
227  std::size_t Hash(const caf::SRProxy * sr);
228 }
const double fxhi
Definition: Utilities.h:213
Det_t
Which NOvA detector?
Definition: SREnums.h:7
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:370
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:363
fVtxDx Fit("f")
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
Int_t par
Definition: SimpleIterate.C:24
void CombineMetadata(std::map< std::string, std::string > &base, const std::map< std::string, std::string > &add, std::set< std::string > &mask)
base += add
Definition: Utilities.cxx:255
ifdh calls between construction and destruction produce no output
Definition: Utilities.h:47
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
std::string FindCAFAnaDir()
Definition: Utilities.cxx:208
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
TVectorT< double > TVectorD
Definition: Utilities.h:20
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, bool xlog, int nbinsy, double ymin, double ymax, bool ylog)
Internal helper for Surface and FCSurface.
Definition: Utilities.cxx:153
Double_t ymax
Definition: plot.C:25
std::vector< std::string > LoadFileList(const std::string &listfile)
Read list of input files from a text file, one per line.
Definition: Utilities.cxx:214
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:89
const double fxlo
Definition: Utilities.h:212
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)
Write map of metadata parameters into a CAF file.
Definition: Utilities.cxx:316
caf::StandardRecord * sr
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
const TH1 * fHist
Definition: Utilities.h:211
std::string FindPackageDir(std::string Dir)
Definition: Utilities.cxx:188
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:379
bool SAMDefinitionExists(const std::string &def)
Definition: UtilsExt.cxx:296
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Extract map of metadata parameters from a CAF file.
Definition: Utilities.cxx:233
assert(nhit_max >=nhit_nbins)
fvar< T > log1p(const fvar< T > &x)
Definition: log1p.hpp:11
const int fNOsc
Definition: Utilities.h:214
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:333
Alter floating-point exception flag.
Definition: Utilities.h:59
Double_t ymin
Definition: plot.C:24
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
add("abs", expr_type(int_type()), expr_type(int_type()))