Logo ROOT  
Reference Guide
MnContours.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/MnContours.h"
11#include "Minuit2/MnMinos.h"
12#include "Minuit2/MnMigrad.h"
15#include "Minuit2/FCNBase.h"
16#include "Minuit2/MnCross.h"
17#include "Minuit2/MinosError.h"
19
20#include "Minuit2/MnPrint.h"
21
22
23
24namespace ROOT {
25
26 namespace Minuit2 {
27
28
29void PrintContourPoint(const std::pair<double,double> & point) {
30 std::cout << "\t x = " << point.first << " y = " << point.second << std::endl;
31}
32
33std::vector<std::pair<double,double> > MnContours::operator()(unsigned int px, unsigned int py, unsigned int npoints) const {
34 // get contour as a pair of (x,y) points passing the parameter index (px, py) and the number of requested points (>=4)
35 ContoursError cont = Contour(px, py, npoints);
36 return cont();
37}
38
39ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int npoints) const {
40 // calculate the contour passing the parameter index (px, py) and the number of requested points (>=4)
41 // the fcn.UP() has to be set to the rquired value (see Minuit document on errors)
42 assert(npoints > 3);
43 unsigned int maxcalls = 100*(npoints+5)*(fMinimum.UserState().VariableParameters()+1);
44 unsigned int nfcn = 0;
45
46 std::vector<std::pair<double,double> > result; result.reserve(npoints);
47 std::vector<MnUserParameterState> states;
48 // double edmmax = 0.5*0.05*fFCN.Up()*1.e-3;
49
50 //double toler = 0.05;
51 double toler = 0.1; // use same defaut value as in Minos
52
53 //get first four points
54 // std::cout<<"MnContours: get first 4 params."<<std::endl;
56
57 double valx = fMinimum.UserState().Value(px);
58 double valy = fMinimum.UserState().Value(py);
59
60 MinosError mex = minos.Minos(px);
61 nfcn += mex.NFcn();
62 if(!mex.IsValid()) {
63 MN_ERROR_MSG("MnContours is unable to find first two points.");
64 return ContoursError(px, py, result, mex, mex, nfcn);
65 }
66 std::pair<double,double> ex = mex();
67
68 MinosError mey = minos.Minos(py);
69 nfcn += mey.NFcn();
70 if(!mey.IsValid()) {
71 MN_ERROR_MSG("MnContours is unable to find second two points.");
72 return ContoursError(px, py, result, mex, mey, nfcn);
73 }
74 std::pair<double,double> ey = mey();
75
76 MnMigrad migrad(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy()-1))));
77
78 migrad.Fix(px);
79 migrad.SetValue(px, valx + ex.second);
80 FunctionMinimum exy_up = migrad();
81 nfcn += exy_up.NFcn();
82 if(!exy_up.IsValid()) {
83 MN_ERROR_VAL2("MnContours: unable to find Upper y Value for x Parameter",px);
84 return ContoursError(px, py, result, mex, mey, nfcn);
85 }
86
87 migrad.SetValue(px, valx + ex.first);
88 FunctionMinimum exy_lo = migrad();
89 nfcn += exy_lo.NFcn();
90 if(!exy_lo.IsValid()) {
91 MN_ERROR_VAL2("MnContours: unable to find Lower y Value for x Parameter",px);
92 return ContoursError(px, py, result, mex, mey, nfcn);
93 }
94
95
96 MnMigrad migrad1(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy()-1))));
97 migrad1.Fix(py);
98 migrad1.SetValue(py, valy + ey.second);
99 FunctionMinimum eyx_up = migrad1();
100 nfcn += eyx_up.NFcn();
101 if(!eyx_up.IsValid()) {
102 MN_ERROR_VAL2("MnContours: unable to find Upper x Value for y Parameter",py);
103 return ContoursError(px, py, result, mex, mey, nfcn);
104 }
105
106 migrad1.SetValue(py, valy + ey.first);
107 FunctionMinimum eyx_lo = migrad1();
108 nfcn += eyx_lo.NFcn();
109 if(!eyx_lo.IsValid()) {
110 MN_ERROR_VAL2("MnContours: unable to find Lower x Value for y Parameter",py);
111 return ContoursError(px, py, result, mex, mey, nfcn);
112 }
113
114 double scalx = 1./(ex.second - ex.first);
115 double scaly = 1./(ey.second - ey.first);
116
117 result.push_back(std::pair<double,double>(valx + ex.first, exy_lo.UserState().Value(py)));
118 result.push_back(std::pair<double,double>(eyx_lo.UserState().Value(px), valy + ey.first));
119 result.push_back(std::pair<double,double>(valx + ex.second, exy_up.UserState().Value(py)));
120 result.push_back(std::pair<double,double>(eyx_up.UserState().Value(px), valy + ey.second));
121
122
124
125 // std::cout<<"MnContours: first 4 params finished."<<std::endl;
126 int printLevel = MnPrint::Level();
127
128 if (printLevel > 0 ) {
129 std::cout << "MnContour : List of found points " << std::endl;
130 std::cout << "\t Parameter x is " << upar.Name(px) << std::endl;
131 std::cout << "\t Parameter y is " << upar.Name(py) << std::endl;
132 }
133
134 if (printLevel > 0) {
135 for (unsigned int i = 0; i < 4; ++i)
136 PrintContourPoint(result[i] );
137 }
138
139 upar.Fix(px);
140 upar.Fix(py);
141
142 std::vector<unsigned int> par(2); par[0] = px; par[1] = py;
143 MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
144
145 for(unsigned int i = 4; i < npoints; i++) {
146
147 std::vector<std::pair<double,double> >::iterator idist1 = result.end()-1;
148 std::vector<std::pair<double,double> >::iterator idist2 = result.begin();
149 double dx = idist1->first - (idist2)->first;
150 double dy = idist1->second - (idist2)->second;
151 double bigdis = scalx*scalx*dx*dx + scaly*scaly*dy*dy;
152
153 for(std::vector<std::pair<double,double> >::iterator ipair = result.begin(); ipair != result.end()-1; ++ipair) {
154 double distx = ipair->first - (ipair+1)->first;
155 double disty = ipair->second - (ipair+1)->second;
156 double dist = scalx*scalx*distx*distx + scaly*scaly*disty*disty;
157 if(dist > bigdis) {
158 bigdis = dist;
159 idist1 = ipair;
160 idist2 = ipair+1;
161 }
162 }
163
164 double a1 = 0.5;
165 double a2 = 0.5;
166 double sca = 1.;
167
168L300:
169
170 if(nfcn > maxcalls) {
171 MN_ERROR_MSG("MnContours: maximum number of function calls exhausted.");
172 return ContoursError(px, py, result, mex, mey, nfcn);
173 }
174
175 double xmidcr = a1*idist1->first + a2*(idist2)->first;
176 double ymidcr = a1*idist1->second + a2*(idist2)->second;
177 double xdir = (idist2)->second - idist1->second;
178 double ydir = idist1->first - (idist2)->first;
179 double scalfac = sca*std::max(fabs(xdir*scalx), fabs(ydir*scaly));
180 double xdircr = xdir/scalfac;
181 double ydircr = ydir/scalfac;
182 std::vector<double> pmid(2); pmid[0] = xmidcr; pmid[1] = ymidcr;
183 std::vector<double> pdir(2); pdir[0] = xdircr; pdir[1] = ydircr;
184
185 MnCross opt = cross(par, pmid, pdir, toler, maxcalls);
186 nfcn += opt.NFcn();
187 if(!opt.IsValid()) {
188 // if(a1 > 0.5) {
189 if(sca < 0.) {
190 MN_ERROR_VAL2("MnContours : unable to find point on Contour",i+1);
191 MN_ERROR_VAL2("MnContours : found only i points",i);
192 return ContoursError(px, py, result, mex, mey, nfcn);
193 }
194 // a1 = 0.75;
195 // a2 = 0.25;
196 // std::cout<<"*****switch direction"<<std::endl;
197 sca = -1.;
198 goto L300;
199 }
200 double aopt = opt.Value();
201 if(idist2 == result.begin()) {
202 result.push_back(std::pair<double,double>(xmidcr+(aopt)*xdircr, ymidcr + (aopt)*ydircr));
203 if (printLevel > 0) PrintContourPoint( result.back() );
204 }
205 else {
206 result.insert(idist2, std::pair<double,double>(xmidcr+(aopt)*xdircr, ymidcr + (aopt)*ydircr));
207 if (printLevel > 0) PrintContourPoint( *idist2 );
208 }
209 }
210 if (printLevel >0)
211 std::cout << "MnContour: Number of contour points = " << result.size() << std::endl;
212
213 return ContoursError(px, py, result, mex, mey, nfcn);
214}
215
216
217 } // namespace Minuit2
218
219} // namespace ROOT
#define MN_ERROR_VAL2(loc, x)
Definition: MnPrint.h:132
#define MN_ERROR_MSG(str)
Definition: MnPrint.h:113
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
const MnUserParameterState & UserState() const
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
unsigned int NFcn() const
Definition: MinosError.h:78
void SetValue(unsigned int, double)
const FCNBase & fFCN
Definition: MnContours.h:64
std::vector< std::pair< double, double > > operator()(unsigned int, unsigned int, unsigned int npoints=20) const
ask for one Contour (points only) from number of points (>=4) and parameter indeces
Definition: MnContours.cxx:33
ContoursError Contour(unsigned int, unsigned int, unsigned int npoints=20) const
ask for one Contour ContoursError (MinosErrors + points) from number of points (>=4) and parameter in...
Definition: MnContours.cxx:39
const FunctionMinimum & fMinimum
Definition: MnContours.h:65
unsigned int NFcn() const
Definition: MnCross.h:67
bool IsValid() const
Definition: MnCross.h:63
double Value() const
Definition: MnCross.h:61
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
API class for Minos Error analysis (asymmetric errors); minimization has to be done before and Minimu...
Definition: MnMinos.h:34
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:84
static int Level()
Definition: MnPrint.cxx:47
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
unsigned int Strategy() const
Definition: MnStrategy.h:39
class which holds the external user and/or internal Minuit representation of the parameters and error...
const char * Name(unsigned int) const
Double_t ey[n]
Definition: legend1.C:17
Double_t ex[n]
Definition: legend1.C:17
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void PrintContourPoint(const std::pair< double, double > &point)
Definition: MnContours.cxx:29
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
static constexpr double second
Definition: first.py:1