Logo ROOT  
Reference Guide
MnMinos.cxx
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
10 #include "Minuit2/MnMinos.h"
12 #include "Minuit2/FCNBase.h"
14 #include "Minuit2/MnCross.h"
15 #include "Minuit2/MinosError.h"
16 #include "Minuit2/MnPrint.h"
17 
18 namespace ROOT {
19 
20 namespace Minuit2 {
21 
22 MnMinos::MnMinos(const FCNBase &fcn, const FunctionMinimum &min, unsigned int stra)
23  : fFCN(fcn), fMinimum(min), fStrategy(MnStrategy(stra))
24 {
25  MnPrint print("MnMinos");
26 
27  // construct from FCN + Minimum
28  // check if Error definition has been changed, in case re-update errors
29  if (fcn.Up() != min.Up()) {
30  print.Warn("MnMinos: UP value has changed, need to update FunctionMinimum class");
31  }
32 }
33 
34 MnMinos::MnMinos(const FCNBase &fcn, const FunctionMinimum &min, const MnStrategy &stra)
35  : fFCN(fcn), fMinimum(min), fStrategy(stra)
36 {
37  MnPrint print("MnMinos");
38 
39  // construct from FCN + Minimum
40  // check if Error definition has been changed, in case re-update errors
41  if (fcn.Up() != min.Up()) {
42  print.Warn("UP value has changed, need to update FunctionMinimum class");
43  }
44 }
45 
46 std::pair<double, double> MnMinos::operator()(unsigned int par, unsigned int maxcalls, double toler) const
47 {
48  // do Minos analysis given the parameter index returning a pair for (lower,upper) errors
49  MinosError mnerr = Minos(par, maxcalls, toler);
50  return mnerr();
51 }
52 
53 double MnMinos::Lower(unsigned int par, unsigned int maxcalls, double toler) const
54 {
55  // get lower error for parameter par
56 
57  MnCross aopt = Loval(par, maxcalls, toler);
58 
59  MinosError mnerr(par, fMinimum.UserState().Value(par), aopt, MnCross());
60 
61  return mnerr.Lower();
62 }
63 
64 double MnMinos::Upper(unsigned int par, unsigned int maxcalls, double toler) const
65 {
66  // upper error for parameter par
67 
68  MnCross aopt = Upval(par, maxcalls, toler);
69 
70  MinosError mnerr(par, fMinimum.UserState().Value(par), MnCross(), aopt);
71 
72  return mnerr.Upper();
73 }
74 
75 MinosError MnMinos::Minos(unsigned int par, unsigned int maxcalls, double toler) const
76 {
77  // do full minos error anlysis (lower + upper) for parameter par
78 
79  MnPrint print("MnMinos");
80 
81  MnCross up = Upval(par, maxcalls, toler);
82 
83  print.Debug("Function calls to find upper error", up.NFcn());
84 
85  MnCross lo = Loval(par, maxcalls, toler);
86 
87  print.Debug("Function calls to find lower error", lo.NFcn());
88 
89  print.Debug("return Minos error", lo.Value(), ",", up.Value());
90 
91  return MinosError(par, fMinimum.UserState().Value(par), lo, up);
92 }
93 
94 MnCross MnMinos::FindCrossValue(int direction, unsigned int par, unsigned int maxcalls, double toler) const
95 {
96  // get crossing value in the parameter direction :
97  // direction = + 1 upper value
98  // direction = -1 lower value
99  // pass now tolerance used for Migrad minimizations
100 
101  assert(direction == 1 || direction == -1);
102 
103  MnPrint print("MnMinos");
104 
105  print.Info("Determination of", direction == 1 ? "upper" : "lower", "Minos error for parameter", par);
106 
107  assert(fMinimum.IsValid());
108  assert(!fMinimum.UserState().Parameter(par).IsFixed());
109  assert(!fMinimum.UserState().Parameter(par).IsConst());
110 
111  if (maxcalls == 0) {
112  unsigned int nvar = fMinimum.UserState().VariableParameters();
113  maxcalls = 2 * (nvar + 1) * (200 + 100 * nvar + 5 * nvar * nvar);
114  }
115 
116  std::vector<unsigned int> para(1, par);
117 
119  double err = direction * upar.Error(par);
120  double val = upar.Value(par) + err;
121  // check if we do not cross limits
122  if (direction == 1 && upar.Parameter(par).HasUpperLimit()) {
123  val = std::min(val, upar.Parameter(par).UpperLimit());
124  }
125  if (direction == -1 && upar.Parameter(par).HasLowerLimit()) {
126  val = std::max(val, upar.Parameter(par).LowerLimit());
127  }
128  // recompute err in case it was truncated for the limit
129  err = val - upar.Value(par);
130  std::vector<double> xmid(1, val);
131  std::vector<double> xdir(1, err);
132 
133  double up = fFCN.Up();
134  unsigned int ind = upar.IntOfExt(par);
135  // get error matrix (methods return a copy)
137  // get internal parameters
138  const MnAlgebraicVector &xt = fMinimum.Parameters().Vec();
139  // LM: change to use err**2 (m(i,i) instead of err as in F77 version
140  double xunit = std::sqrt(up / m(ind, ind));
141  // LM (29/04/08) bug: change should be done in internal variables
142  // set the initial value for the other parmaeters that we are going to fit in MnCross
143  for (unsigned int i = 0; i < m.Nrow(); i++) {
144  if (i == ind)
145  continue;
146  double xdev = xunit * m(ind, i);
147  double xnew = xt(i) + direction * xdev;
148 
149  // transform to external values
150  unsigned int ext = upar.ExtOfInt(i);
151 
152  double unew = upar.Int2ext(i, xnew);
153 
154  // take into account limits
155  if (upar.Parameter(ext).HasUpperLimit()) {
156  unew = std::min(unew, upar.Parameter(ext).UpperLimit());
157  }
158  if (upar.Parameter(ext).HasLowerLimit()) {
159  unew = std::max(unew, upar.Parameter(ext).LowerLimit());
160  }
161 
162  print.Debug("Parameter", ext, "is set from", upar.Value(ext), "to", unew);
163 
164  upar.SetValue(ext, unew);
165  }
166 
167  upar.Fix(par);
168  upar.SetValue(par, val);
169 
170  print.Debug("Parameter", par, "is fixed and set from", fMinimum.UserState().Value(par), "to", val, "delta =", err);
171 
172  MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
173  MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
174 
175  print.Debug("aopt value found from MnFunctionCross =", aopt.Value());
176 
177  const char *par_name = upar.Name(par);
178  if (aopt.AtMaxFcn())
179  print.Warn("maximum number of function calls exceeded for Parameter", par_name);
180  if (aopt.NewMinimum())
181  print.Warn("new Minimum found while looking for Parameter", par_name);
182  if (direction == 1) {
183  if (aopt.AtLimit())
184  print.Warn("parameter", par_name, "is at Upper limit");
185  if (!aopt.IsValid())
186  print.Warn("could not find Upper Value for Parameter", par_name);
187  } else {
188  if (aopt.AtLimit())
189  print.Warn("parameter", par_name, "is at Lower limit");
190  if (!aopt.IsValid())
191  print.Warn("could not find Lower Value for Parameter", par_name);
192  }
193 
194  print.Info("end of Minos scan for", direction == 1 ? "up" : "low", "interval for parameter", upar.Name(par));
195 
196  return aopt;
197 }
198 
199 MnCross MnMinos::Upval(unsigned int par, unsigned int maxcalls, double toler) const
200 {
201  // return crossing in the lower parameter direction
202  return FindCrossValue(1, par, maxcalls, toler);
203 }
204 
205 MnCross MnMinos::Loval(unsigned int par, unsigned int maxcalls, double toler) const
206 {
207  // return crossing in the lower parameter direction
208  return FindCrossValue(-1, par, maxcalls, toler);
209 }
210 
211 } // namespace Minuit2
212 
213 } // namespace ROOT
ROOT::Minuit2::MnCross::Value
double Value() const
Definition: MnCross.h:89
ROOT::Minuit2::MnMinos::fFCN
const FCNBase & fFCN
Definition: MnMinos.h:65
m
auto * m
Definition: textangle.C:8
ROOT::Minuit2::MnMinos::fMinimum
const FunctionMinimum & fMinimum
Definition: MnMinos.h:66
ROOT::Minuit2::MnMinos::FindCrossValue
MnCross FindCrossValue(int dir, unsigned int, unsigned int maxcalls, double toler) const
internal method to get crossing value via MnFunctionCross
Definition: MnMinos.cxx:94
ROOT::Minuit2::MnFunctionCross
MnFunctionCross.
Definition: MnFunctionCross.h:29
ROOT::Minuit2::MnMinos::Lower
double Lower(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
calculate one side (negative or positive Error) of the Parameter give as input (optionally) maxcalls ...
Definition: MnMinos.cxx:53
ROOT::Minuit2::MnUserParameterState::Error
double Error(unsigned int) const
Definition: MnUserParameterState.cxx:409
ROOT::Minuit2::MinimumError::Matrix
MnAlgebraicSymMatrix Matrix() const
Definition: MinimumError.h:59
ROOT::Minuit2::FunctionMinimum::Up
double Up() const
Definition: FunctionMinimum.h:84
MnMinos.h
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::MnMinos::Loval
MnCross Loval(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:205
ROOT::Minuit2::MnPrint::Debug
void Debug(const Ts &... args)
Definition: MnPrint.h:138
ROOT::Minuit2::MinuitParameter::UpperLimit
double UpperLimit() const
Definition: MinuitParameter.h:158
ROOT::Minuit2::FunctionMinimum::IsValid
bool IsValid() const
Definition: FunctionMinimum.h:85
ROOT::Minuit2::MnCross::NewMinimum
bool NewMinimum() const
Definition: MnCross.h:94
FunctionMinimum.h
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
ROOT::Minuit2::MinimumParameters::Vec
const MnAlgebraicVector & Vec() const
Definition: MinimumParameters.h:38
ROOT::Minuit2::FCNBase
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:45
ROOT::Minuit2::MnCross::NFcn
unsigned int NFcn() const
Definition: MnCross.h:95
ROOT::Minuit2::MnMinos::fStrategy
MnStrategy fStrategy
Definition: MnMinos.h:67
MnFunctionCross.h
ROOT::Minuit2::MinuitParameter::HasUpperLimit
bool HasUpperLimit() const
Definition: MinuitParameter.h:156
MnCross.h
ROOT::Minuit2::MnUserParameterState::VariableParameters
unsigned int VariableParameters() const
Definition: MnUserParameterState.cxx:506
ROOT::Minuit2::FunctionMinimum::UserState
const MnUserParameterState & UserState() const
Definition: FunctionMinimum.h:71
ROOT::Minuit2::MinuitParameter::IsConst
bool IsConst() const
Definition: MinuitParameter.h:151
ROOT::Minuit2::FunctionMinimum::Error
const MinimumError & Error() const
Definition: FunctionMinimum.h:78
ROOT::Minuit2::FunctionMinimum::Fval
double Fval() const
Definition: FunctionMinimum.h:80
ROOT::Minuit2::MnPrint::Warn
void Warn(const Ts &... args)
Definition: MnPrint.h:126
ROOT::Minuit2::MnMinos::operator()
std::pair< double, double > operator()(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
returns the negative (pair.first) and the positive (pair.second) Minos Error of the Parameter
Definition: MnMinos.cxx:46
ROOT::Minuit2::MnCross
Definition: MnCross.h:19
ROOT::Minuit2::MnUserParameterState::Name
const char * Name(unsigned int) const
Definition: MnUserParameterState.cxx:473
ROOT::Minuit2::MinosError::Lower
double Lower() const
Definition: MinosError.h:54
ROOT::Minuit2::MinosError
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
ROOT::Minuit2::MnUserParameterState::SetValue
void SetValue(unsigned int, double)
Definition: MnUserParameterState.cxx:328
ROOT::Minuit2::FCNBase::Up
virtual double Up() const =0
Error definition of the function.
ROOT::Minuit2::MinuitParameter::IsFixed
bool IsFixed() const
Definition: MinuitParameter.h:152
ROOT::Minuit2::MnCross::IsValid
bool IsValid() const
Definition: MnCross.h:91
sqrt
double sqrt(double)
ROOT::Minuit2::MnUserParameterState::ExtOfInt
unsigned int ExtOfInt(unsigned int) const
Definition: MnUserParameterState.cxx:501
ROOT::Minuit2::MnMinos::Upper
double Upper(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:64
MinosError.h
ROOT::Minuit2::FunctionMinimum
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Definition: FunctionMinimum.h:33
ROOT::Minuit2::MnMinos::Upval
MnCross Upval(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:199
ROOT::Minuit2::MnUserParameterState
class which holds the external user and/or internal Minuit representation of the parameters and error...
Definition: MnUserParameterState.h:33
ROOT::Minuit2::MnUserParameterState::Parameter
const MinuitParameter & Parameter(unsigned int i) const
Definition: MnUserParameterState.cxx:230
ROOT::Minuit2::MnUserParameterState::IntOfExt
unsigned int IntOfExt(unsigned int) const
Definition: MnUserParameterState.cxx:496
ROOT::Minuit2::MnMinos::Minos
MinosError Minos(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
ask for MinosError (Lower + Upper) can be printed via std::cout
Definition: MnMinos.cxx:75
ROOT::Minuit2::MnStrategy
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
ROOT::Minuit2::MnMinos::MnMinos
MnMinos(const FCNBase &fcn, const FunctionMinimum &min, unsigned int stra=1)
construct from FCN + Minimum + strategy
Definition: MnMinos.cxx:22
ROOT::Minuit2::MinosError::Upper
double Upper() const
Definition: MinosError.h:63
ROOT::Minuit2::FunctionMinimum::Parameters
const MinimumParameters & Parameters() const
Definition: FunctionMinimum.h:77
FCNBase.h
ROOT::Minuit2::MinuitParameter::LowerLimit
double LowerLimit() const
Definition: MinuitParameter.h:157
ROOT::Minuit2::MnCross::AtLimit
bool AtLimit() const
Definition: MnCross.h:92
ROOT::Minuit2::MnUserParameterState::Value
double Value(unsigned int) const
Definition: MnUserParameterState.cxx:404
MnPrint.h
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
ROOT::Minuit2::MnUserParameterState::Int2ext
double Int2ext(unsigned int, double) const
Definition: MnUserParameterState.cxx:486
ROOT::Minuit2::MinuitParameter::HasLowerLimit
bool HasLowerLimit() const
Definition: MinuitParameter.h:155
ROOT::Minuit2::MnUserParameterState::Fix
void Fix(unsigned int)
Definition: MnUserParameterState.cxx:297
ROOT::Minuit2::MnPrint::Info
void Info(const Ts &... args)
Definition: MnPrint.h:132
ROOT::Minuit2::MnPrint
Definition: MnPrint.h:73
ROOT::Minuit2::MnCross::AtMaxFcn
bool AtMaxFcn() const
Definition: MnCross.h:93