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 #include "Minuit2/LaSum.h"
0017
0018 #include <memory>
0019
0020 namespace ROOT {
0021
0022 namespace Minuit2 {
0023
0024 /** MinimumError keeps the inv. 2nd derivative (inv. Hessian) used for
0025     calculating the Parameter step size (-V*g) and for the covariance Update
0026     (ErrorUpdator). The covariance matrix is equal to twice the inv. Hessian.
0027  */
0028 class MinimumError {
0029
0030 public:
0031    enum Status {
0032       MnUnset,
0033       MnPosDef,
0034       MnMadePosDef,
0035       MnNotPosDef,
0036       MnHesseFailed,
0037       MnInvertFailed,
0038       MnReachedCallLimit,
0039    };
0040
0041 public:
0042    MinimumError(unsigned int n) : fPtr{new Data{{n}, 1.0, MnUnset}} {}
0043
0044    MinimumError(const MnAlgebraicSymMatrix &mat, double dcov) : fPtr{new Data{mat, dcov, MnPosDef}} {}
0045
0046    MinimumError(const MnAlgebraicSymMatrix &mat, Status status) : fPtr{new Data{mat, 1.0, status}} {}
0047
0048    MnAlgebraicSymMatrix Matrix() const { return 2. * fPtr->fMatrix; } // why *2 ?
0049
0050    const MnAlgebraicSymMatrix &InvHessian() const { return fPtr->fMatrix; }
0051
0052    // calculate invert of matrix. Used to compute Hessian  by inverting matrix
0053    MnAlgebraicSymMatrix Hessian() const
0054    {
0055       return InvertMatrix(fPtr->fMatrix);
0056    }
0057
0058    static MnAlgebraicSymMatrix InvertMatrix(const MnAlgebraicSymMatrix & matrix, int & ifail) {
0059        // calculate inverse of given matrix
0060       MnAlgebraicSymMatrix tmp(matrix);
0061       ifail = ROOT::Minuit2::Invert(tmp);
0062       if (ifail != 0) {
0063          MnPrint print("MinimumError::Invert");
0064          print.Warn("Inversion fails; return diagonal matrix");
0065          for (unsigned int i = 0; i < matrix.Nrow(); ++i)
0066             for (unsigned int j = 0; j <= i; j++)
0067                tmp(i, j) = i == j ? 1. / matrix(i, i) : 0;
0068       }
0069       return tmp;
0070    }
0071    static MnAlgebraicSymMatrix InvertMatrix(const MnAlgebraicSymMatrix & matrix) {
0072       int ifail = 0;
0073       return InvertMatrix(matrix, ifail);
0074    }
0075
0076    double Dcovar() const { return fPtr->fDCovar; }
0077    Status GetStatus() const { return fPtr->fStatus; }
0078
0079    bool IsValid() const { return IsAvailable() && (IsPosDef() || IsMadePosDef() || IsNotPosDef()); }
0080    bool IsAccurate() const { return IsPosDef() && Dcovar() < 0.1; }
0081
0082    bool IsPosDef() const { return GetStatus() == MnPosDef; }
0083    bool IsMadePosDef() const { return GetStatus() == MnMadePosDef; }
0084    bool IsNotPosDef() const { return GetStatus() == MnNotPosDef; }
0085    bool HesseFailed() const { return GetStatus() == MnHesseFailed; }
0086    bool InvertFailed() const { return GetStatus() == MnInvertFailed; }
0087    bool HasReachedCallLimit() const { return GetStatus() == MnReachedCallLimit; }
0088    bool IsAvailable() const { return GetStatus() != MnUnset; }
0089
0090 private:
0091    struct Data {
0092       MnAlgebraicSymMatrix fMatrix;
0093       double fDCovar;
0094       Status fStatus;
0095    };
0096
0097    std::shared_ptr<Data> fPtr;
0098 };
0099
0100 } // namespace Minuit2
0101
0102 } // namespace ROOT
0103
0104 #endif // ROOT_Minuit2_MinimumError