root/Minuit2/MinimumError.h

0001 // @(#)root/minuit2:$Id$
0002 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005
0003
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
0007  *                                                                    *
0008  **********************************************************************/
0009
0010 #ifndef ROOT_Minuit2_MinimumError
0011 #define ROOT_Minuit2_MinimumError
0012
0013 #include "Minuit2/MnConfig.h"
0014 #include "Minuit2/MnMatrix.h"
0015 #include "Minuit2/MnPrint.h"
0016
0017 #include <memory>
0018
0019 namespace ROOT {
0020
0021 namespace Minuit2 {
0022
0023 /** MinimumError keeps the inv. 2nd derivative (inv. Hessian) used for
0024     calculating the Parameter step size (-V*g) and for the covariance Update
0025     (ErrorUpdator). The covariance matrix is equal to twice the inv. Hessian.
0026  */
0027 class MinimumError {
0028
0029 public:
0030    enum Status {
0031       MnUnset,
0032       MnPosDef,
0033       MnMadePosDef,
0034       MnNotPosDef,
0035       MnHesseFailed,
0036       MnInvertFailed,
0037       MnReachedCallLimit,
0038    };
0039
0040 public:
0041    MinimumError(unsigned int n) : fPtr{new Data{{n}, {0}, 1.0, MnUnset}} {}
0042
0043    MinimumError(const MnAlgebraicSymMatrix &mat, double dcov) : fPtr{new Data{mat, {0}, dcov, MnPosDef}} {}
0044
0045    MinimumError(const MnAlgebraicSymMatrix &mat, const MnAlgebraicSymMatrix &hess, double dcov) : fPtr{new Data{mat, hess, dcov, MnPosDef}} {}
0046
0047    MinimumError(const MnAlgebraicSymMatrix &mat, Status status) : fPtr{new Data{mat, {0}, 1.0, status}} {}
0048
0049    MnAlgebraicSymMatrix Matrix() const { return 2. * fPtr->fMatrix; } // why *2 ?
0050
0051    const MnAlgebraicSymMatrix &InvHessian() const { return fPtr->fMatrix; }
0052
0053    // calculate invert of matrix. Used to compute Hessian  by inverting matrix
0054    const MnAlgebraicSymMatrix & Hessian() const
0055    {
0056       if (fPtr->fHessian.size() == 0)
0057          fPtr->fHessian = InvertMatrix(fPtr->fMatrix);
0058       return fPtr->fHessian;
0059    }
0060
0061    static MnAlgebraicSymMatrix InvertMatrix(const MnAlgebraicSymMatrix & matrix, int & ifail) {
0062        // calculate inverse of given matrix
0063       MnAlgebraicSymMatrix tmp(matrix);
0064       ifail = ROOT::Minuit2::Invert(tmp);
0065       if (ifail != 0) {
0066          MnPrint print("MinimumError::Invert");
0067          print.Warn("Inversion fails; return diagonal matrix");
0068          for (unsigned int i = 0; i < matrix.Nrow(); ++i)
0069             for (unsigned int j = 0; j <= i; j++)
0070                tmp(i, j) = i == j ? 1. / matrix(i, i) : 0;
0071       }
0072       return tmp;
0073    }
0074    static MnAlgebraicSymMatrix InvertMatrix(const MnAlgebraicSymMatrix & matrix) {
0075       int ifail = 0;
0076       return InvertMatrix(matrix, ifail);
0077    }
0078
0079    double Dcovar() const { return fPtr->fDCovar; }
0080    Status GetStatus() const { return fPtr->fStatus; }
0081
0082    bool IsValid() const { return IsAvailable() && (IsPosDef() || IsMadePosDef() || IsNotPosDef()); }
0083    bool IsAccurate() const { return IsPosDef() && Dcovar() < 0.1; }
0084
0085    bool IsPosDef() const { return GetStatus() == MnPosDef; }
0086    bool IsMadePosDef() const { return GetStatus() == MnMadePosDef; }
0087    bool IsNotPosDef() const { return GetStatus() == MnNotPosDef; }
0088    bool HesseFailed() const { return GetStatus() == MnHesseFailed; }
0089    bool InvertFailed() const { return GetStatus() == MnInvertFailed; }
0090    bool HasReachedCallLimit() const { return GetStatus() == MnReachedCallLimit; }
0091    bool IsAvailable() const { return GetStatus() != MnUnset; }
0092
0093 private:
0094    struct Data {
0095       MnAlgebraicSymMatrix fMatrix;
0096       MnAlgebraicSymMatrix fHessian;  // optional stored also Hessian (used in Fumili)
0097       double fDCovar;
0098       Status fStatus;
0099    };
0100
0101    std::shared_ptr<Data> fPtr;
0102 };
0103
0104 } // namespace Minuit2
0105
0106 } // namespace ROOT
0107
0108 #endif // ROOT_Minuit2_MinimumError