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  /// Value passed to --stride, or 1 if not specified
196  size_t Stride(bool allow_default = true);
197  /// Value passed to --offset, or 0 if not specified
198  size_t Offset(bool allow_default = true);
199  /// Value passed to --limit, or -1 if not specified
200  int Limit();
201 
202  /// What's the process number for a grid job?
203  size_t JobNumber();
204  size_t NumJobs();
205 
206  bool SAMDefinitionExists(const std::string& def);
207 
208  // Calling this function will return a Fourier series, fit to the input
209  // histogram. Assumes x-axis covers one period
211  {
212  public:
213  FitToFourier(TH1* h, double xlo, double xhi, int NOsc);
214  ~FitToFourier();
215  TF1* Fit() const;
216  double operator()(double *x, double *par) const;
217  private:
218 
219  const TH1* fHist; // Histogram to fit
220  const double fxlo; // Lower bound
221  const double fxhi; // Upper bound - assumed to be 1 osc from the low end
222  const int fNOsc; // Highest harmonic to include
223 
224  };
225 
226  /// Generate a unique hash from a sequence of integers
227  std::size_t Hash(const std::vector<std::size_t> &vals);
228 
229 
230  /// Generate a unique hash for a SRSpill (run/subrun/event combination)
231  std::size_t Hash(const caf::SRSpillProxy * sr);
232 
233 
234  /// Generate a unique hash for a StandardRecord (run/subrun/cycle/event/slice combination)
235  std::size_t Hash(const caf::SRProxy * sr);
236 }
const double fxhi
Definition: Utilities.h:221
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:437
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:358
fVtxDx Fit("f")
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
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:250
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:203
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:148
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:209
int Limit()
Value passed to –limit, or -1 if not specified.
Definition: Utilities.cxx:419
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:84
const double fxlo
Definition: Utilities.h:220
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:311
caf::StandardRecord * sr
size_t Stride(bool allow_default)
Value passed to –stride, or 1 if not specified.
Definition: Utilities.cxx:371
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
const TH1 * fHist
Definition: Utilities.h:219
std::string FindPackageDir(std::string Dir)
Definition: Utilities.cxx:183
Proxy for caf::SRSpill.
Definition: SRProxy.h:1346
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:459
size_t Offset(bool allow_default)
Value passed to –offset, or 0 if not specified.
Definition: Utilities.cxx:395
bool SAMDefinitionExists(const std::string &def)
Definition: UtilsExt.cxx:296
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
size_t NumJobs()
Definition: Utilities.cxx:448
std::map< std::string, std::string > GetCAFMetadata(TDirectory *dir)
Extract map of metadata parameters from a CAF file.
Definition: Utilities.cxx:228
assert(nhit_max >=nhit_nbins)
const int fNOsc
Definition: Utilities.h:222
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:328
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
enum BeamMode string