Logo ROOT   6.10/09
Reference Guide
RooBrentRootFinder.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 /**
18 \file RooBrentRootFinder.cxx
19 \class RooBrentRootFinder
20 \ingroup Roofitcore
21 
22 Implement the abstract 1-dimensional root finding interface using
23 the Brent-Decker method. This implementation is based on the one
24 in the GNU scientific library (v0.99).
25 **/
26 
27 #include "RooFit.h"
28 
29 #include "RooBrentRootFinder.h"
30 #include "RooBrentRootFinder.h"
31 #include "RooAbsFunc.h"
32 #include <math.h>
33 #include "Riostream.h"
34 #include "RooMsgService.h"
35 
36 using namespace std;
37 
39 ;
40 
41 
42 ////////////////////////////////////////////////////////////////////////////////
43 /// Constructor taking function binding as input
44 
47  _tol(2.2204460492503131e-16)
48 {
49 }
50 
51 
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 /// Do the root finding using the Brent-Decker method. Returns a boolean status and
55 /// loads 'result' with our best guess at the root if true.
56 /// Prints a warning if the initial interval does not bracket a single
57 /// root or if the root is not found after a fixed number of iterations.
58 
60 {
61  _function->saveXVec() ;
62 
63  Double_t a(xlo),b(xhi);
64  Double_t fa= (*_function)(&a) - value;
65  Double_t fb= (*_function)(&b) - value;
66  if(fb*fa > 0) {
67  oocxcoutD((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): initial interval does not bracket a root: ("
68  << a << "," << b << "), value = " << value << " f[xlo] = " << fa << " f[xhi] = " << fb << endl;
69  return kFALSE;
70  }
71 
72  Bool_t ac_equal(kFALSE);
73  Double_t fc= fb;
74  Double_t c(0),d(0),e(0);
75  for(Int_t iter= 0; iter <= MaxIterations; iter++) {
76 
77  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
78  // Rename a,b,c and adjust bounding interval d
79  ac_equal = kTRUE;
80  c = a;
81  fc = fa;
82  d = b - a;
83  e = b - a;
84  }
85 
86  if (fabs (fc) < fabs (fb)) {
87  ac_equal = kTRUE;
88  a = b;
89  b = c;
90  c = a;
91  fa = fb;
92  fb = fc;
93  fc = fa;
94  }
95 
96  Double_t tol = 0.5 * _tol * fabs(b);
97  Double_t m = 0.5 * (c - b);
98 
99 
100  if (fb == 0 || fabs(m) <= tol) {
101  //cout << "RooBrentRootFinder: iter = " << iter << " m = " << m << " tol = " << tol << endl ;
102  result= b;
103  _function->restoreXVec() ;
104  return kTRUE;
105  }
106 
107  if (fabs (e) < tol || fabs (fa) <= fabs (fb)) {
108  // Bounds decreasing too slowly: use bisection
109  d = m;
110  e = m;
111  }
112  else {
113  // Attempt inverse cubic interpolation
114  Double_t p, q, r;
115  Double_t s = fb / fa;
116 
117  if (ac_equal) {
118  p = 2 * m * s;
119  q = 1 - s;
120  }
121  else {
122  q = fa / fc;
123  r = fb / fc;
124  p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
125  q = (q - 1) * (r - 1) * (s - 1);
126  }
127  // Check whether we are in bounds
128  if (p > 0) {
129  q = -q;
130  }
131  else {
132  p = -p;
133  }
134 
135  Double_t min1= 3 * m * q - fabs (tol * q);
136  Double_t min2= fabs (e * q);
137  if (2 * p < (min1 < min2 ? min1 : min2)) {
138  // Accept the interpolation
139  e = d;
140  d = p / q;
141  }
142  else {
143  // Interpolation failed: use bisection.
144  d = m;
145  e = m;
146  }
147  }
148  // Move last best guess to a
149  a = b;
150  fa = fb;
151  // Evaluate new trial root
152  if (fabs (d) > tol) {
153  b += d;
154  }
155  else {
156  b += (m > 0 ? +tol : -tol);
157  }
158  fb= (*_function)(&b) - value;
159 
160  }
161  // Return our best guess if we run out of iterations
162  oocoutE((TObject*)0,Eval) << "RooBrentRootFinder::findRoot(" << _function->getName() << "): maximum iterations exceeded." << endl;
163  result= b;
164 
166 
167  return kFALSE;
168 }
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
RooBrentRootFinder(const RooAbsFunc &function)
Constructor taking function binding as input.
STL namespace.
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:1956
#define oocoutE(o, a)
Definition: RooMsgService.h:47
#define oocxcoutD(o, a)
Definition: RooMsgService.h:81
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:146
virtual const char * getName() const
Definition: RooAbsFunc.h:59
virtual Bool_t findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value=0) const
Do the root finding using the Brent-Decker method.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const double tol
RooAbsRootFinder is the abstract interface for finding roots of real-valued 1-dimensional function th...
TRandom2 r(17)
TMarker * m
Definition: textangle.C:8
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method...
const RooAbsFunc * _function
const Bool_t kFALSE
Definition: RtypesCore.h:92
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
virtual void restoreXVec() const
Definition: RooAbsFunc.h:54
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
virtual void saveXVec() const
Definition: RooAbsFunc.h:51
Mother of all ROOT objects.
Definition: TObject.h:37
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double result[121]
float * q
Definition: THbookFile.cxx:87
const Bool_t kTRUE
Definition: RtypesCore.h:91
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23