Utilities.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Core/UtilsExt.h"

#include <cassert>
#include <cmath>
#include <cxxabi.h>
#include <cfenv>
#include <map>
#include <memory>
#include <set>
#include <string>
#include <vector>

#include "StandardRecord/SREnums.h"
#include "CAFAna/Core/FwdDeclare.h"

template<class T> class TMatrixT;
typedef TMatrixT<double> TMatrixD;
template<class T> class TVectorT;
typedef TVectorT<double> TVectorD;

class TArrayD;
class TDirectory;
class TH1;
class TH2;
class TH3;
class TF1;
class TH1D;
class TH1F;
class TH2F;
class TH2D;
class TH3D;
class TRandom;
class TVector3;

#include "StandardRecord/Proxy/FwdDeclare.h"

namespace ana
{
  class Binning;
  class Ratio;

  /// \brief ifdh calls between construction and destruction produce no output
  ///
  /// Upon going out of scope, restores the previous setting
  class IFDHSilent
  {
  public:
    IFDHSilent();
    ~IFDHSilent();
  protected:
    bool fSet;
  };

  /// \brief Alter floating-point exception flag
  ///
  /// Upon going out of scope, restores the previous setting
  class FloatingExceptionOnNaN
  {
  public:
    FloatingExceptionOnNaN(bool enable = true);
    ~FloatingExceptionOnNaN();
  protected:
    fexcept_t fBackup;
  };

  /** \brief Compute bin-to-bin covariance matrix from a collection of sets of bin contents.

      \param binSets   Collection of sets of bins from which covariances should be calculated
      \param firstBin  The first bin that should be considered (inclusive)
      \param lastBin   The last bin that should be considered (inclusive).  -1 means "last in set"

      \returns  unique_ptr to TMatrixD containing computed covariance matrix unless binSets.size() < 2,
                in which case the unique_ptr's target is nullptr.

      Note TH1D is a child class of TArrayD -- so you can pass a vector
      of TH1D* to this method.
  **/
  std::unique_ptr<TMatrixD> CalcCovMx(const std::vector<TArrayD*> & binSets, int firstBin=0, int lastBin=-1);

  /** \brief The log-likelihood formula from the PDG.

      \param exp The expected spectrum
      \param obs The corresponding observed spectrum

      \returns The log-likelihood formula from the PDG
      \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]

      Includes underflow bin and an option for
      overflow bin (off by default) and handles
      zero observed or expected events correctly.
  **/
  double LogLikelihood(const TH1* exp, const TH1* obs, bool useOverflow = false);

  /** \brief The log-likelihood formula for a single bin

      \param exp Expected count
      \param obs Observed count

      \returns The log-likelihood formula from the PDG
      \f[ \chi^2=2\left(e-o+o\ln\left({o\over e}\right)\right) \f]

      Handles zero observed or expected events correctly.

      Templated so that it can handle usage with Stan vars and other numeric types.
      (The horible third template parameter ensures this function can only be used
       with types that accept conversion from double -- which means you can't pass
       it a TH1* or a std::vector<stan::math::var>&, removing the ambiguity with
       those other versions of LogLikelihood). The return type promotes to
       stan::math::var if either T or U are.
  **/
  template <typename T, typename U,
            typename std::enable_if_t<std::is_convertible_v<double, T> && std::is_convertible_v<double, U>, int> = 0>
  decltype(T(0) - U(0)) LogLikelihood(T exp, U obs)
  {
    // With this value, negative expected events and one observed
    // event gives a chisq from this one bin of 182.
    const double minexp = 1e-40; // Don't let expectation go lower than this

    assert(obs >= 0);
    if(exp < minexp){
      if(!obs) return 0.;
      exp = minexp;
    }

    if(obs*1000 > exp){
      // This strange form is for numerical stability when exp ~ obs
      return 2*obs*((exp-obs)/obs + log1p((obs-exp)/exp));
    }
    else{
      // But log1p doesn't like arguments near -1 (observation much smaller
      // than expectation), and it's better to use the usual formula in that
      // case.
      if(obs){
        return 2*(exp-obs + obs*log(obs/exp));
      }
      else{
        return 2*exp;
      }
    }
  }

  /// \brief Internal helper for \ref Surface and \ref FCSurface
  ///
  /// Creates a histogram having bins \em centred at the min and max
  /// coordinates
  TH2F* ExpandedHistogram(const std::string& title,
                          int nbinsx, double xmin, double xmax, bool xlog,
                          int nbinsy, double ymin, double ymax, bool ylog);

  /// This is $SRT_PRIVATE_CONTEXT if you have used a test release,
  /// otherwise $SRT_PUBLIC_CONTEXT
  std::string FindPackageDir( std::string Dir );

  /// This is $SRT_PRIVATE_CONTEXT if a CAFAna directory exists there,
  /// otherwise $SRT_PUBLIC_CONTEXT
  std::string FindCAFAnaDir();

  /// Read list of input files from a text file, one per line
  std::vector<std::string> LoadFileList(const std::string& listfile);

  /// \brief Extract map of metadata parameters from a CAF file
  ///
  /// \param dir The "meta" directory from the CAF file
  /// \return    A map from metadata field name to metadata value
  std::map<std::string, std::string> GetCAFMetadata(TDirectory* dir);

  /// \brief \a base += \a add
  ///
  /// \param base The original source of strings, will be altered
  /// \param add  Strings to add to \a base if missing
  /// \param mask Fields for which there was a conflict, will be altered
  void CombineMetadata(std::map<std::string, std::string>& base,
                       const std::map<std::string, std::string>& add,
                       std::set<std::string>& mask);

  /// \brief Write map of metadata parameters into a CAF file
  ///
  /// \param dir  The "meta" directory of the CAF file
  /// \param meta Map from metadata field name to metadata value
  void WriteCAFMetadata(TDirectory* dir,
                        const std::map<std::string, std::string>& meta);


  /// \brief Average direction of NuMI neutrinos in a given detector
  /// This function is a copy of
  /// \ref geo::GeometryBase::NuMIBeamDirection ()
  /// Any changes made to that function must be propagated to here.
  TVector3 NuMIBeamDirection(caf::Det_t det);

  /// Is this a grid (condor) job?
  bool RunningOnGrid();

  /// Value passed to --stride, or 1 if not specified
  size_t Stride(bool allow_default = true);
  /// Value passed to --offset, or 0 if not specified
  size_t Offset(bool allow_default = true);
  /// Value passed to --limit, or -1 if not specified
  int Limit();

  /// What's the process number for a grid job?
  size_t JobNumber();
  size_t NumJobs();

  bool SAMDefinitionExists(const std::string& def);

  // Calling this function will return a Fourier series, fit to the input
  // histogram.  Assumes x-axis covers one period
  class FitToFourier
  {
  public:
    FitToFourier(TH1* h, double xlo, double xhi, int NOsc);
    ~FitToFourier();
    TF1* Fit() const;
    double operator()(double *x, double *par) const;
  private:

    const TH1*   fHist; // Histogram to fit
    const double fxlo;  // Lower bound
    const double fxhi;  // Upper bound - assumed to be 1 osc from the low end
    const int    fNOsc; // Highest harmonic to include

  };

  /// Generate a unique hash from a sequence of integers
  std::size_t Hash(const std::vector<std::size_t> &vals);


  /// Generate a unique hash for a SRSpill (run/subrun/event combination)
  std::size_t Hash(const caf::SRSpillProxy * sr);


  /// Generate a unique hash for a StandardRecord (run/subrun/cycle/event/slice combination)
  std::size_t Hash(const caf::SRProxy * sr);
}