ROOT  6.06/09
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 
20 #include "Minuit2/MnPrint.h"
21 
22 
23 
24 namespace ROOT {
25 
26  namespace Minuit2 {
27 
28 
29 void PrintContourPoint(const std::pair<double,double> & point) {
30  std::cout << "\t x = " << point.first << " y = " << point.second << std::endl;
31 }
32 
33 std::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 
39 ContoursError 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 
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 
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;
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 
168 L300:
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
unsigned int NFcn() const
Definition: MnCross.h:65
void PrintContourPoint(const std::pair< double, double > &point)
Definition: MnContours.cxx:29
double par[1]
Definition: unuranDistr.cxx:38
#define MN_ERROR_MSG(str)
Definition: MnPrint.h:113
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
static int Level()
Definition: MnPrint.cxx:47
TVector3 cross(const TVector3 &v1, const TVector3 &v2)
Definition: CsgOps.cxx:816
#define assert(cond)
Definition: unittest.h:542
bool IsValid() const
Definition: MinosError.h:59
double Value() const
Definition: MnCross.h:59
const MnUserParameterState & UserState() const
const char * Name(unsigned int) const
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
API class for Minos Error analysis (asymmetric errors); minimization has to be done before and Minimu...
Definition: MnMinos.h:34
#define MN_ERROR_VAL2(loc, x)
Definition: MnPrint.h:132
bool IsValid() const
Definition: MnCross.h:61
bool minos
Definition: testMinim.cxx:33
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const FunctionMinimum & fMinimum
Definition: MnContours.h:65
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
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
class which holds the external user and/or internal Minuit representation of the parameters and error...
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:88
unsigned int Strategy() const
Definition: MnStrategy.h:39
unsigned int NFcn() const
Definition: MinosError.h:68
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
Double_t ey[n]
Definition: legend1.C:17
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
void SetValue(unsigned int, double)
double result[121]
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
Double_t ex[n]
Definition: legend1.C:17
const FCNBase & fFCN
Definition: MnContours.h:64