diff --git a/src/Statistical/StatUtils.h b/src/Statistical/StatUtils.h index 71672e2..fdf7bd9 100644 --- a/src/Statistical/StatUtils.h +++ b/src/Statistical/StatUtils.h @@ -1,294 +1,297 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef STATUTILS_H #define STATUTILS_H // C Includes #include "assert.h" #include #include #include #include #include #include #include #include // Root Includes #include "TDecompChol.h" #include "TDecompSVD.h" #include "TFile.h" #include "TGraphErrors.h" #include "TH1D.h" #include "TH2D.h" #include "TH2I.h" #include "TMath.h" #include "TMatrixDSym.h" #include "TRandom3.h" // Fit Includes #include "FitLogger.h" /*! * \addtogroup Utils * @{ */ //! Functions for handling statistics calculations namespace StatUtils { /* Chi2 Functions */ //! Get Chi2 using diagonal bin errors from the histogram. Masking applied //! before calculation if mask provided. Double_t GetChi2FromDiag(TH1D *data, TH1D *mc, TH1I *mask = NULL); //! Get Chi2 using diagonal bin errors from the histogram. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromDiag(TH2D *data, TH2D *mc, TH2I *map = NULL, TH2I *mask = NULL); //! Get Chi2 using an inverted covariance for the data Double_t GetChi2FromCov(TH1D *data, TH1D *mc, TMatrixDSym *invcov, TH1I *mask = NULL, double data_scale = 1, double covar_scale = 1E76, TH1D *outchi2perbin = NULL); //! Get Chi2 using an inverted covariance for the data //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromCov(TH2D *data, TH2D *mc, TMatrixDSym *invcov, TH2I *map = NULL, TH2I *mask = NULL, TH2D *outchi2perbin = NULL); //! Get Chi2 using an SVD method on the covariance before calculation. //! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work. Double_t GetChi2FromSVD(TH1D *data, TH1D *mc, TMatrixDSym *cov, TH1I *mask = NULL); //! Get Chi2 using an SVD method on the covariance before calculation. //! Method suggested by Rex at MiniBooNE. Shown that it doesn't actually work. //! Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromSVD(TH2D *data, TH2D *mc, TMatrixDSym *cov, TH2I *map = NULL, TH2I *mask = NULL); //! Get Chi2 using only the raw event rates given in each bin using a -2LL //! method. Double_t GetChi2FromEventRate(TH1D *data, TH1D *mc, TH1I *mask = NULL); //! Get Chi2 using only the raw event rates given in each bin using a -2LL //! method. Plots converted to 1D histograms before using 1D calculation. Double_t GetChi2FromEventRate(TH2D *data, TH2D *mc, TH2I *map = NULL, TH2I *mask = NULL); // Likelihood Functions //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromDiag(TH1D *data, TH1D *mc, TH1I *mask = NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromDiag(TH2D *data, TH2D *mc, TH2I *map = NULL, TH2I *mask = NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromCov(TH1D *data, TH1D *mc, TMatrixDSym *invcov, TH1I *mask = NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromCov(TH2D *data, TH2D *mc, TMatrixDSym *invcov, TH2I *map = NULL, TH2I *mask = NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromSVD(TH1D *data, TH1D *mc, TMatrixDSym *cov, TH1I *mask = NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromSVD(TH2D *data, TH2D *mc, TMatrixDSym *cov, TH2I *map = NULL, TH2I *mask = NULL); //! Placeholder for 1D binned likelihood method Double_t GetLikelihoodFromEventRate(TH1D *data, TH1D *mc, TH1I *mask = NULL); //! Placeholder for 2D binned likelihood method Double_t GetLikelihoodFromEventRate(TH2D *data, TH2D *mc, TH2I *map = NULL, TH2I *mask = NULL); /* NDOF Functions */ //! Return 1D Histogram NDOF considering masking and empty bins Int_t GetNDOF(TH1D *hist, TH1I *mask = NULL); //! Return 2D Histogram NDOF considering masking and empty bins Int_t GetNDOF(TH2D *hist, TH2I *map = NULL, TH2I *mask = NULL); /* Fake Data Functions */ //! Given a full covariance for a 1D data set throw the decomposition to //! generate fake data. throwdiag determines whether diagonal statistical errors //! are thrown. If no covariance is provided only statistical errors are thrown. TH1D *ThrowHistogram(TH1D *hist, TMatrixDSym *cov, bool throwdiag = true, TH1I *mask = NULL); //! Given a full covariance for a 2D data set throw the decomposition to //! generate fake data. Plots are converted to 1D histograms and the 1D //! ThrowHistogram is used, before being converted back to 2D histograms. TH2D *ThrowHistogram(TH2D *hist, TMatrixDSym *cov, TH2I *map = NULL, bool throwdiag = true, TH2I *mask = NULL); /* Masking Functions */ //! Given a mask histogram, mask out any bins in hist with non zero entries in //! mask. TH1D *ApplyHistogramMasking(TH1D *hist, TH1I *mask); //! Given a mask histogram, mask out any bins in hist with non zero entries in //! mask. TH2D *ApplyHistogramMasking(TH2D *hist, TH2I *mask); //! Given a mask histogram apply the masking procedure to each of the //! rows/columns in a covariance, before recalculating its inverse. TMatrixDSym *ApplyInvertedMatrixMasking(TMatrixDSym *mat, TH1I *mask); //! Given a mask histogram apply the masking procedure to each of the //! rows/columns in a covariance, before recalculating its inverse. Converts to //! 1D data before using the 1D ApplyInvertedMatrixMasking function and //! converting back to 2D. TMatrixDSym *ApplyInvertedMatrixMasking(TMatrixDSym *mat, TH2D *data, TH2I *mask, TH2I *map = NULL); //! Given a mask histogram apply the masking procedure to each of the //! rows/columns in a covariance TMatrixDSym *ApplyMatrixMasking(TMatrixDSym *mat, TH1I *mask); //! Given a mask histogram apply the masking procedure to each of the //! rows/columns in a covariance Converts to 1D data before using the 1D //! ApplyInvertedMatrixMasking function and converting back to 2D. TMatrixDSym *ApplyMatrixMasking(TMatrixDSym *mat, TH2D *data, TH2I *mask, TH2I *map = NULL); /* Covariance Handling Functions */ +//! Check if a matrix can be inverted with cholesky method +bool IsMatrixWellBehaved(TMatrixDSym* mat); + //! Return inverted matrix of TMatrixDSym TMatrixDSym *GetInvert(TMatrixDSym *mat, bool rescale = false); //! Return Cholesky Decomposed matrix of TMatrixDSym TMatrixDSym *GetDecomp(TMatrixDSym *mat); //! Return full covariances TMatrixDSym *GetCovarFromCorrel(TMatrixDSym *correl, TH1D *data); //! Given a normalisation factor for a dataset add in a new normalisation term //! to the covariance. void ForceNormIntoCovar(TMatrixDSym *&mat, TH1D *data, double norm); //! Given a normalisation factor for a dataset add in a new normalisation term //! to the covariance. Convertes 2D to 1D, before using 1D ForceNormIntoCovar void ForceNormIntoCovar(TMatrixDSym *mat, TH2D *data, double norm, TH2I *map = NULL); //! Given a dataset generate an uncorrelated covariance matrix using the bin //! errors. TMatrixDSym *MakeDiagonalCovarMatrix(TH1D *data, double scaleF = 1E38); //! Given a dataset generate an uncorrelated covariance matrix using the bin //! errors. TMatrixDSym *MakeDiagonalCovarMatrix(TH2D *data, TH2I *map = NULL, double scaleF = 1E38); //! Given a covariance set the errors in each bin on the data from the //! covariance diagonals. void SetDataErrorFromCov(TH1D *data, TMatrixDSym *cov, double scale = 1.0, bool ErrorCheck = true); //! Given a covariance set the errors in each bin on the data from the //! covariance diagonals. void SetDataErrorFromCov(TH2D *data, TMatrixDSym *cov, TH2I *map = NULL, double scale = 1.0, bool ErrorCheck = true); //! Given a covariance, extracts the shape-only matrix using the method from the //! MiniBooNE TN TMatrixDSym *ExtractShapeOnlyCovar(TMatrixDSym *full_covar, TH1D *data_hist, double data_scale = 1E38); TMatrixDSym *ExtractShapeOnlyCovar(TMatrixDSym *full_covar, TH2D *data_hist, TH2I *map = NULL, double data_scale = 1E38); /* Mapping Functions */ //! If no map is provided for the 2D histogram, generate one by counting up //! through the bins along x and y. TH2I *GenerateMap(TH2D *hist); //! Apply a map to a 2D histogram converting it into a 1D histogram. TH1D *MapToTH1D(TH2D *hist, TH2I *map); //! Apply a map to fill a TH2D from a TH1D; void MapFromTH1D(TH2 *fillhist, TH1 *fromhist, TH2I *map); //! Apply a map to a 2D mask convering it into a 1D mask. TH1I *MapToMask(TH2I *hist, TH2I *map); /// \brief Read TMatrixD from a text file /// /// - covfile = full path to text file /// - dimx = x dimensions of matrix /// - dimy = y dimensions of matrix /// /// Format of textfile should be: \n /// cov_11 cov_12 ... cov_1N \n /// cov_21 cov_22 ... cov_2N \n /// ... ... ... ... \n /// cov_N1 ... ... cov_NN \n /// /// If no dimensions are given, dimx and dimy are determined from rows/columns /// inside textfile. /// /// If only dimx is given a symmetric matrix is assumed. TMatrixD *GetMatrixFromTextFile(std::string covfile, int dimx = -1, int dimy = -1); /// \brief Read TMatrixD from a ROOT file /// /// - covfile = full path to root file (+';histogram') /// - histname = histogram name /// /// If no histogram name is given function assumes it has been appended /// covfile path as: \n /// 'covfile.root;histname' /// /// histname can point to a TMatrixD object, a TMatrixDSym object, or /// a TH2D object. TMatrixD *GetMatrixFromRootFile(std::string covfile, std::string histname = ""); /// \brief Calls GetMatrixFromTextFile and turns it into a TMatrixDSym TMatrixDSym *GetCovarFromTextFile(std::string covfile, int dim); /// \brief Calls GetMatrixFromRootFile and turns it into a TMatrixDSym TMatrixDSym *GetCovarFromRootFile(std::string covfile, std::string histname); }; // namespace StatUtils /*! @} */ #endif