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"
18 #include "Minuit2/ContoursError.h"
19 #include "Minuit2/MnPrint.h"
20 
21 namespace ROOT {
22 
23 namespace Minuit2 {
24 
25 std::vector<std::pair<double, double>> MnContours::
26 operator()(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 
34 ContoursError 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;
54  MnMinos minos(fFCN, fMinimum, fStrategy);
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
ROOT::Minuit2::MnCross::Value
double Value() const
Definition: MnCross.h:89
ROOT::Minuit2::MnApplication::Fix
void Fix(unsigned int)
Definition: MnApplication.cxx:118
ex
Double_t ex[n]
Definition: legend1.C:17
ROOT::Minuit2::MnStrategy::Strategy
unsigned int Strategy() const
Definition: MnStrategy.h:38
ROOT::Minuit2::MnContours::operator()
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
first
Definition: first.py:1
ROOT::Minuit2::MnFunctionCross
MnFunctionCross.
Definition: MnFunctionCross.h:29
ROOT::Minuit2::MnMigrad
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:32
MnMinos.h
ROOT::Minuit2::MnContours::fFCN
const FCNBase & fFCN
Definition: MnContours.h:66
ROOT::Minuit2::FunctionMinimum::IsValid
bool IsValid() const
Definition: FunctionMinimum.h:85
FunctionMinimum.h
ROOT::Minuit2::MnCross::NFcn
unsigned int NFcn() const
Definition: MnCross.h:95
ROOT::Minuit2::ContoursError
Definition: ContoursError.h:23
MnFunctionCross.h
MnMigrad.h
MnCross.h
MnContours.h
ROOT::Minuit2::MnUserParameterState::VariableParameters
unsigned int VariableParameters() const
Definition: MnUserParameterState.cxx:506
ROOT::Minuit2::MnPrint::Error
void Error(const Ts &... args)
Definition: MnPrint.h:120
ROOT::Minuit2::FunctionMinimum::UserState
const MnUserParameterState & UserState() const
Definition: FunctionMinimum.h:71
ContoursError.h
ROOT::Minuit2::FunctionMinimum::Fval
double Fval() const
Definition: FunctionMinimum.h:80
ROOT::Math::gv_detail::dist
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Minuit2::MnCross
Definition: MnCross.h:19
ROOT::Minuit2::MnUserParameterState::Name
const char * Name(unsigned int) const
Definition: MnUserParameterState.cxx:473
ROOT::Minuit2::MinosError
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
ROOT::Minuit2::MnCross::IsValid
bool IsValid() const
Definition: MnCross.h:91
ey
Double_t ey[n]
Definition: legend1.C:17
ROOT::Minuit2::MnApplication::SetValue
void SetValue(unsigned int, double)
Definition: MnApplication.cxx:126
ROOT::Minuit2::MnContours::fMinimum
const FunctionMinimum & fMinimum
Definition: MnContours.h:67
ROOT::Minuit2::MnContours::Contour
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
MinosError.h
ROOT::Minuit2::FunctionMinimum
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Definition: FunctionMinimum.h:33
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::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::FunctionMinimum::NFcn
int NFcn() const
Definition: FunctionMinimum.h:82
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::MnContours::fStrategy
MnStrategy fStrategy
Definition: MnContours.h:68
ROOT::Minuit2::MnMinos
API class for Minos Error analysis (asymmetric errors); minimization has to be done before and Minimu...
Definition: MnMinos.h:33
FCNBase.h
ROOT::Minuit2::MinosError::NFcn
unsigned int NFcn() const
Definition: MinosError.h:84
ROOT::Minuit2::MnUserParameterState::Value
double Value(unsigned int) const
Definition: MnUserParameterState.cxx:404
MnPrint.h
TGeant4Unit::second
static constexpr double second
Definition: TGeant4SystemOfUnits.h:151
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::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::MinosError::IsValid
bool IsValid() const
Definition: MinosError.h:75