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
18namespace ROOT {
19
20namespace Minuit2 {
21
22MnMinos::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
34MnMinos::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
46std::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
53double 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
64double 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
75MinosError 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
94MnCross 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
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
199MnCross 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
205MnCross 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
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:45
virtual double Up() const =0
Error definition of the function.
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
const MinimumParameters & Parameters() const
const MinimumError & Error() const
const MnUserParameterState & UserState() const
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
MnAlgebraicSymMatrix Matrix() const
Definition: MinimumError.h:48
const MnAlgebraicVector & Vec() const
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
double Upper() const
Definition: MinosError.h:63
double Lower() const
Definition: MinosError.h:54
bool AtLimit() const
Definition: MnCross.h:92
unsigned int NFcn() const
Definition: MnCross.h:95
bool AtMaxFcn() const
Definition: MnCross.h:93
bool IsValid() const
Definition: MnCross.h:91
bool NewMinimum() const
Definition: MnCross.h:94
double Value() const
Definition: MnCross.h:89
const FunctionMinimum & fMinimum
Definition: MnMinos.h:66
const FCNBase & fFCN
Definition: MnMinos.h:65
MnStrategy fStrategy
Definition: MnMinos.h:67
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
double Upper(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:64
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
MnCross Loval(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:205
MnCross FindCrossValue(int dir, unsigned int, unsigned int maxcalls, double toler) const
internal method to get crossing value via MnFunctionCross
Definition: MnMinos.cxx:94
MnMinos(const FCNBase &fcn, const FunctionMinimum &min, unsigned int stra=1)
construct from FCN + Minimum + strategy
Definition: MnMinos.cxx:22
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
MnCross Upval(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:199
void Debug(const Ts &... args)
Definition: MnPrint.h:147
void Info(const Ts &... args)
Definition: MnPrint.h:141
void Warn(const Ts &... args)
Definition: MnPrint.h:135
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
class which holds the external user and/or internal Minuit representation of the parameters and error...
double Int2ext(unsigned int, double) const
const MinuitParameter & Parameter(unsigned int i) const
unsigned int ExtOfInt(unsigned int) const
const char * Name(unsigned int) const
unsigned int IntOfExt(unsigned int) const
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
TMarker m
Definition: textangle.C:8