#include "RooFit.h"
#include "RooBrentRootFinder.h"
#include "RooBrentRootFinder.h"
#include "RooAbsFunc.h"
#include <math.h>
#include "Riostream.h"
ClassImp(RooBrentRootFinder)
;
RooBrentRootFinder::RooBrentRootFinder(const RooAbsFunc& function) :
RooAbsRootFinder(function)
{
}
Bool_t RooBrentRootFinder::findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value) const
{
Double_t a(xlo),b(xhi);
Double_t fa= (*_function)(&a) - value;
Double_t fb= (*_function)(&b) - value;
if(fb*fa > 0) {
cout << "RooBrentRootFinder::findRoot: initial interval does not bracket a root: ("
<< a << "," << b << "), value = " << value << endl;
return kFALSE;
}
Bool_t ac_equal(kFALSE);
Double_t fc= fb;
Double_t c(0),d(0),e(0);
for(Int_t iter= 0; iter <= MaxIterations; iter++) {
if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
ac_equal = kTRUE;
c = a;
fc = fa;
d = b - a;
e = b - a;
}
if (fabs (fc) < fabs (fb)) {
ac_equal = kTRUE;
a = b;
b = c;
c = a;
fa = fb;
fb = fc;
fc = fa;
}
Double_t tol = 0.5 * 2.2204460492503131e-16 * fabs(b);
Double_t m = 0.5 * (c - b);
if (fb == 0 || fabs(m) <= tol) {
result= b;
return kTRUE;
}
if (fabs (e) < tol || fabs (fa) <= fabs (fb)) {
d = m;
e = m;
}
else {
Double_t p, q, r;
Double_t s = fb / fa;
if (ac_equal) {
p = 2 * m * s;
q = 1 - s;
}
else {
q = fa / fc;
r = fb / fc;
p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
q = (q - 1) * (r - 1) * (s - 1);
}
if (p > 0) {
q = -q;
}
else {
p = -p;
}
Double_t min1= 3 * m * q - fabs (tol * q);
Double_t min2= fabs (e * q);
if (2 * p < (min1 < min2 ? min1 : min2)) {
e = d;
d = p / q;
}
else {
d = m;
e = m;
}
}
a = b;
fa = fb;
if (fabs (d) > tol) {
b += d;
}
else {
b += (m > 0 ? +tol : -tol);
}
fb= (*_function)(&b) - value;
}
cout << "RooBrentRootFinder::findRoot: maximum iterations exceeded." << endl;
result= b;
return kFALSE;
}
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.