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#include "Minuit2/MnPrint.h"
20
21namespace ROOT {
22
23namespace Minuit2 {
24
25std::vector<std::pair<double, double>> MnContours::
26operator()(unsigned int px, unsigned int py, unsigned int npoints) const
27{
28 // get contour as a pair of (x,y) points passing the parameter index (px, py) and the number of requested points
29 // (>=4)
30 ContoursError cont = Contour(px, py, npoints);
31 return cont();
32}
33
34ContoursError MnContours::Contour(unsigned int px, unsigned int py, unsigned int npoints) const
35{
36 // calculate the contour passing the parameter index (px, py) and the number of requested points (>=4)
37 // the fcn.UP() has to be set to the rquired value (see Minuit document on errors)
38 assert(npoints > 3);
39 unsigned int maxcalls = 100 * (npoints + 5) * (fMinimum.UserState().VariableParameters() + 1);
40 unsigned int nfcn = 0;
41
42 MnPrint print("MnContours");
43
44 std::vector<std::pair<double, double>> result;
45 result.reserve(npoints);
46 std::vector<MnUserParameterState> states;
47 // double edmmax = 0.5*0.05*fFCN.Up()*1.e-3;
48
49 // double toler = 0.05;
50 double toler = 0.1; // use same defaut value as in Minos
51
52 // get first four points
53 // std::cout<<"MnContours: get first 4 params."<<std::endl;
55
56 double valx = fMinimum.UserState().Value(px);
57 double valy = fMinimum.UserState().Value(py);
58
59 MinosError mex = minos.Minos(px);
60 nfcn += mex.NFcn();
61 if (!mex.IsValid()) {
62 print.Error("unable to find first two points");
63 return ContoursError(px, py, result, mex, mex, nfcn);
64 }
65 std::pair<double, double> ex = mex();
66
67 MinosError mey = minos.Minos(py);
68 nfcn += mey.NFcn();
69 if (!mey.IsValid()) {
70 print.Error("unable to find second two points");
71 return ContoursError(px, py, result, mex, mey, nfcn);
72 }
73 std::pair<double, double> ey = mey();
74
75 MnMigrad migrad(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1))));
76
77 migrad.Fix(px);
78 migrad.SetValue(px, valx + ex.second);
79 FunctionMinimum exy_up = migrad();
80 nfcn += exy_up.NFcn();
81 if (!exy_up.IsValid()) {
82 print.Error("unable to find Upper y Value for x Parameter", px);
83 return ContoursError(px, py, result, mex, mey, nfcn);
84 }
85
86 migrad.SetValue(px, valx + ex.first);
87 FunctionMinimum exy_lo = migrad();
88 nfcn += exy_lo.NFcn();
89 if (!exy_lo.IsValid()) {
90 print.Error("unable to find Lower y Value for x Parameter", px);
91 return ContoursError(px, py, result, mex, mey, nfcn);
92 }
93
94 MnMigrad migrad1(fFCN, fMinimum.UserState(), MnStrategy(std::max(0, int(fStrategy.Strategy() - 1))));
95 migrad1.Fix(py);
96 migrad1.SetValue(py, valy + ey.second);
97 FunctionMinimum eyx_up = migrad1();
98 nfcn += eyx_up.NFcn();
99 if (!eyx_up.IsValid()) {
100 print.Error("unable to find Upper x Value for y Parameter", py);
101 return ContoursError(px, py, result, mex, mey, nfcn);
102 }
103
104 migrad1.SetValue(py, valy + ey.first);
105 FunctionMinimum eyx_lo = migrad1();
106 nfcn += eyx_lo.NFcn();
107 if (!eyx_lo.IsValid()) {
108 print.Error("unable to find Lower x Value for y Parameter", py);
109 return ContoursError(px, py, result, mex, mey, nfcn);
110 }
111
112 double scalx = 1. / (ex.second - ex.first);
113 double scaly = 1. / (ey.second - ey.first);
114
115 result.emplace_back(valx + ex.first, exy_lo.UserState().Value(py));
116 result.emplace_back(eyx_lo.UserState().Value(px), valy + ey.first);
117 result.emplace_back(valx + ex.second, exy_up.UserState().Value(py));
118 result.emplace_back(eyx_up.UserState().Value(px), valy + ey.second);
119
121
122 print.Info("List of found points", '\n', " Parameter x is", upar.Name(px), '\n', " Parameter y is", upar.Name(py),
123 '\n', result[0], '\n', result[1], '\n', result[2], '\n', result[3]);
124
125 upar.Fix(px);
126 upar.Fix(py);
127
128 std::vector<unsigned int> par(2);
129 par[0] = px;
130 par[1] = py;
131 MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
132
133 for (unsigned int i = 4; i < npoints; i++) {
134
135 auto idist1 = result.end() - 1;
136 auto idist2 = result.begin();
137 double dx = idist1->first - (idist2)->first;
138 double dy = idist1->second - (idist2)->second;
139 double bigdis = scalx * scalx * dx * dx + scaly * scaly * dy * dy;
140
141 for (auto ipair = result.begin(); ipair != result.end() - 1; ++ipair) {
142 double distx = ipair->first - (ipair + 1)->first;
143 double disty = ipair->second - (ipair + 1)->second;
144 double dist = scalx * scalx * distx * distx + scaly * scaly * disty * disty;
145 if (dist > bigdis) {
146 bigdis = dist;
147 idist1 = ipair;
148 idist2 = ipair + 1;
149 }
150 }
151
152 double a1 = 0.5;
153 double a2 = 0.5;
154 double sca = 1.;
155
156 L300:
157
158 if (nfcn > maxcalls) {
159 print.Error("maximum number of function calls exhausted");
160 return ContoursError(px, py, result, mex, mey, nfcn);
161 }
162
163 double xmidcr = a1 * idist1->first + a2 * (idist2)->first;
164 double ymidcr = a1 * idist1->second + a2 * (idist2)->second;
165 double xdir = (idist2)->second - idist1->second;
166 double ydir = idist1->first - (idist2)->first;
167 double scalfac = sca * std::max(std::fabs(xdir * scalx), std::fabs(ydir * scaly));
168 double xdircr = xdir / scalfac;
169 double ydircr = ydir / scalfac;
170 std::vector<double> pmid(2);
171 pmid[0] = xmidcr;
172 pmid[1] = ymidcr;
173 std::vector<double> pdir(2);
174 pdir[0] = xdircr;
175 pdir[1] = ydircr;
176
177 MnCross opt = cross(par, pmid, pdir, toler, maxcalls);
178 nfcn += opt.NFcn();
179 if (!opt.IsValid()) {
180 // if(a1 > 0.5) {
181 if (sca < 0.) {
182 print.Error("unable to find point on Contour", i + 1, '\n', "found only", i, "points");
183 return ContoursError(px, py, result, mex, mey, nfcn);
184 }
185 // a1 = 0.75;
186 // a2 = 0.25;
187 // std::cout<<"*****switch direction"<<std::endl;
188 sca = -1.;
189 goto L300;
190 }
191 double aopt = opt.Value();
192 if (idist2 == result.begin()) {
193 result.emplace_back(xmidcr + (aopt)*xdircr, ymidcr + (aopt)*ydircr);
194 print.Info(result.back());
195 } else {
196 result.insert(idist2, {xmidcr + (aopt)*xdircr, ymidcr + (aopt)*ydircr});
197 print.Info(*idist2);
198 }
199 }
200
201 print.Info("Number of contour points =", result.size());
202
203 return ContoursError(px, py, result, mex, mey, nfcn);
204}
205
206} // namespace Minuit2
207
208} // namespace ROOT
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
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:84
void SetValue(unsigned int, double)
const FCNBase & fFCN
Definition: MnContours.h:66
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:26
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:34
const FunctionMinimum & fMinimum
Definition: MnContours.h:67
unsigned int NFcn() const
Definition: MnCross.h:95
bool IsValid() const
Definition: MnCross.h:91
double Value() const
Definition: MnCross.h:89
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:32
API class for Minos Error analysis (asymmetric errors); minimization has to be done before and Minimu...
Definition: MnMinos.h:33
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
void Error(const Ts &... args)
Definition: MnPrint.h:129
void Info(const Ts &... args)
Definition: MnPrint.h:141
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:38
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)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double second
Definition: first.py:1