Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
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
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
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
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 analysis (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/// Get crossing value in the parameter direction
95/// \param direction 1 upper value, -1 lower value
96/// \param par parameter index
97/// \param maxcalls maximum number of function calls, if 0, it will be replaced by
98/// `2 * (nvar + 1) * (200 + 100 * nvar + 5 * nvar * nvar` where `nvar` is number of variable parameters
99/// \param toler tolerance used for Migrad minimizations
100MnCross MnMinos::FindCrossValue(int direction, unsigned int par, unsigned int maxcalls, double toler) const
101{
102 assert(direction == 1 || direction == -1);
103
104 MnPrint print("MnMinos");
105
106 print.Info("Determination of", direction == 1 ? "upper" : "lower", "Minos error for parameter", par);
107
108 assert(fMinimum.IsValid());
109 assert(!fMinimum.UserState().Parameter(par).IsFixed());
110 assert(!fMinimum.UserState().Parameter(par).IsConst());
111
112 if (maxcalls == 0) {
113 unsigned int nvar = fMinimum.UserState().VariableParameters();
114 maxcalls = 2 * (nvar + 1) * (200 + 100 * nvar + 5 * nvar * nvar);
115 }
116
117 std::vector<unsigned int> para(1, par);
118
119 MnUserParameterState upar = fMinimum.UserState();
120 double err = direction * upar.Error(par);
121 double val = upar.Value(par) + err;
122 // check if we do not cross limits
123 if (direction == 1 && upar.Parameter(par).HasUpperLimit()) {
124 val = std::min(val, upar.Parameter(par).UpperLimit());
125 }
126 if (direction == -1 && upar.Parameter(par).HasLowerLimit()) {
127 val = std::max(val, upar.Parameter(par).LowerLimit());
128 }
129 // recompute err in case it was truncated for the limit
130 err = val - upar.Value(par);
131 std::vector<double> xmid(1, val);
132 std::vector<double> xdir(1, err);
133
134 double up = fFCN.Up();
135 unsigned int ind = upar.IntOfExt(par);
136 // get error matrix (methods return a copy)
137 MnAlgebraicSymMatrix m = fMinimum.Error().Matrix();
138 // get internal parameters
139 const MnAlgebraicVector &xt = fMinimum.Parameters().Vec();
140 // LM: change to use err**2 (m(i,i) instead of err as in F77 version
141 double xunit = std::sqrt(up / m(ind, ind));
142 // LM (29/04/08) bug: change should be done in internal variables
143 // set the initial value for the other parameters that we are going to fit in MnCross
144 for (unsigned int i = 0; i < m.Nrow(); i++) {
145 if (i == ind)
146 continue;
147 double xdev = xunit * m(ind, i);
148 double xnew = xt(i) + direction * xdev;
149
150 // transform to external values
151 unsigned int ext = upar.ExtOfInt(i);
152
153 double unew = upar.Int2ext(i, xnew);
154
155 // take into account limits
156 if (upar.Parameter(ext).HasUpperLimit()) {
157 unew = std::min(unew, upar.Parameter(ext).UpperLimit());
158 }
159 if (upar.Parameter(ext).HasLowerLimit()) {
160 unew = std::max(unew, upar.Parameter(ext).LowerLimit());
161 }
162
163 print.Debug("Parameter", ext, "is set from", upar.Value(ext), "to", unew);
164
165 upar.SetValue(ext, unew);
166 }
167
168 upar.Fix(par);
169 upar.SetValue(par, val);
170
171 print.Debug("Parameter", par, "is fixed and set from", fMinimum.UserState().Value(par), "to", val, "delta =", err);
172
173 MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
174 MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
175
176 print.Debug("aopt value found from MnFunctionCross =", aopt.Value());
177
178 const char *par_name = upar.Name(par);
179 if (aopt.AtMaxFcn())
180 print.Warn("maximum number of function calls exceeded for Parameter", par_name);
181 if (aopt.NewMinimum())
182 print.Warn("new Minimum found while looking for Parameter", par_name);
183 if (direction == 1) {
184 if (aopt.AtLimit())
185 print.Warn("parameter", par_name, "is at Upper limit");
186 if (!aopt.IsValid())
187 print.Warn("could not find Upper Value for Parameter", par_name);
188 } else {
189 if (aopt.AtLimit())
190 print.Warn("parameter", par_name, "is at Lower limit");
191 if (!aopt.IsValid())
192 print.Warn("could not find Lower Value for Parameter", par_name);
193 }
194
195 print.Info("end of Minos scan for", direction == 1 ? "up" : "low", "interval for parameter", upar.Name(par));
196
197 return aopt;
198}
199
200MnCross MnMinos::Upval(unsigned int par, unsigned int maxcalls, double toler) const
201{
202 // return crossing in the lower parameter direction
203 return FindCrossValue(1, par, maxcalls, toler);
204}
205
206MnCross MnMinos::Loval(unsigned int par, unsigned int maxcalls, double toler) const
207{
208 // return crossing in the lower parameter direction
209 return FindCrossValue(-1, par, maxcalls, toler);
210}
211
212} // namespace Minuit2
213
214} // namespace ROOT
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition FCNBase.h:49
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Class describing a symmetric matrix of size n.
Definition MnMatrix.h:438
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition MinosError.h:25
unsigned int NFcn() const
Definition MnCross.h:93
double Value() const
Definition MnCross.h:87
const FunctionMinimum & fMinimum
Definition MnMinos.h:64
const FCNBase & fFCN
Definition MnMinos.h:63
MnStrategy fStrategy
Definition MnMinos.h:65
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:206
MnCross FindCrossValue(int dir, unsigned int, unsigned int maxcalls, double toler) const
internal method to get crossing value via MnFunctionCross
Definition MnMinos.cxx:100
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:200
void Debug(const Ts &... args)
Definition MnPrint.h:135
void Info(const Ts &... args)
Definition MnPrint.h:129
void Warn(const Ts &... args)
Definition MnPrint.h:123
API class for defining four levels of strategies: low (0), medium (1), high (2), very high (>=3); act...
Definition MnStrategy.h:27
class which holds the external user and/or internal Minuit representation of the parameters and error...
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
TMarker m
Definition textangle.C:8