-/* \r
- \r
- Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna\r
- amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru \r
- November. 2, 2006 \r
-\r
-*/ \r
-//This equation solver class is taken from GEANT4 and modified!!\r
-\r
-#include <iostream>\r
-#include <TMath.h>\r
-\r
-template <class Function>\r
-Bool_t EquationSolver<Function>::Brent(Function& theFunction) {\r
- const Double_t precision = 3.0e-8;\r
-\r
- // Check the interval before start\r
- if(fA > fB || TMath::Abs(fA-fB) <= fTolerance) {\r
- std::cerr << "EquationSolver::Brent: The interval must be properly set." << std::endl;\r
- return false;\r
- }\r
- Double_t fa = theFunction(fA);\r
- Double_t fb = theFunction(fB);\r
- if(fa*fb > 0.0) {\r
- std::cerr << "EquationSolver::Brent: The interval must include a root." << std::endl;\r
- return false;\r
- }\r
- \r
- Double_t c = fB;\r
- Double_t fc = fb;\r
- Double_t d = 0.0;\r
- //Double_t fd = 0.0;\r
- Double_t e = 0.0;\r
- //Double_t fe = 0.0;\r
- \r
- for(Int_t i=0; i < fMaxIter; i++) {\r
- // Rename a,b,c and adjust bounding interval d\r
- if(fb*fc > 0.0) {\r
- c = fA;\r
- fc = fa;\r
- d = fB - fA;\r
- e = d;\r
- }\r
- if(TMath::Abs(fc) < TMath::Abs(fb)) {\r
- fA = fB;\r
- fB = c;\r
- c = fA;\r
- fa = fb;\r
- fb = fc;\r
- fc = fa;\r
- }\r
- Double_t Tol1 = 2.0*precision*TMath::Abs(fB) + 0.5*fTolerance;\r
- Double_t xm = 0.5*(c-fB);\r
- if(TMath::Abs(xm) <= Tol1 || fb == 0.0) {\r
- fRoot = fB;\r
- return true;\r
- }\r
- // Inverse quadratic interpolation\r
- if(TMath::Abs(e) >= Tol1 && TMath::Abs(fa) > TMath::Abs(fb)) {\r
- Double_t s = fb/fa;\r
- Double_t p = 0.0;\r
- Double_t q = 0.0;\r
- if(fA == c) {\r
- p = 2.0*xm*s;\r
- q = 1.0 - s;\r
- } \r
- else {\r
- q = fa/fc;\r
- Double_t r = fb/fc;\r
- p = s*(2.0*xm*q*(q-r)-(fB-fA)*(r-1.0));\r
- q = (q-1.0)*(r-1.0)*(s-1.0);\r
- }\r
- // Check bounds\r
- if(p > 0.0) q = -q;\r
- p = TMath::Abs(p);\r
- Double_t min1 = 3.0*xm*q-TMath::Abs(Tol1*q);\r
- Double_t min2 = TMath::Abs(e*q);\r
- if (2.0*p < TMath::Min(min1,min2)) {\r
- // Interpolation\r
- e = d;\r
- d = p/q;\r
- } \r
- else {\r
- // Bisection\r
- d = xm;\r
- e = d;\r
- }\r
- } \r
- else {\r
- // Bounds decreasing too slowly, use bisection\r
- d = xm;\r
- e = d;\r
- }\r
- // Move last guess to a \r
- fA = fB;\r
- fa = fb;\r
- if (TMath::Abs(d) > Tol1) fB += d;\r
- else {\r
- if (xm >= 0.0) fB += TMath::Abs(Tol1);\r
- else fB -= TMath::Abs(Tol1);\r
- }\r
- fb = theFunction(fB);\r
- }\r
- std::cerr << "EquationSolver::Brent: Number of iterations exceeded." << std::endl;\r
- return false;\r
-}\r
-\r
-template <class Function> \r
-void EquationSolver<Function>::SetIntervalLimits(const Double_t Limit1, const Double_t Limit2) {\r
- if(TMath::Abs(Limit1-Limit2) <= fTolerance) {\r
- std::cerr << "EquationSolver::SetIntervalLimits: Interval must be wider than tolerance." << std::endl;\r
- return;\r
- }\r
- if(Limit1 < Limit2) {\r
- fA = Limit1;\r
- fB = Limit2;\r
- } \r
- else {\r
- fA = Limit2;\r
- fB = Limit1;\r
- }\r
- return;\r
-}\r
+/*
+
+ Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
+ amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
+ November. 2, 2006
+
+*/
+//This equation solver class is taken from GEANT4 and modified!!
+
+#include <iostream>
+#include <TMath.h>
+
+template <class Function>
+Bool_t EquationSolver<Function>::Brent(Function& theFunction) {
+ const Double_t precision = 3.0e-8;
+
+ // Check the interval before start
+ if(fA > fB || TMath::Abs(fA-fB) <= fTolerance) {
+ std::cerr << "EquationSolver::Brent: The interval must be properly set." << std::endl;
+ return false;
+ }
+ Double_t fa = theFunction(fA);
+ Double_t fb = theFunction(fB);
+ if(fa*fb > 0.0) {
+ std::cerr << "EquationSolver::Brent: The interval must include a root." << std::endl;
+ return false;
+ }
+
+ Double_t c = fB;
+ Double_t fc = fb;
+ Double_t d = 0.0;
+ //Double_t fd = 0.0;
+ Double_t e = 0.0;
+ //Double_t fe = 0.0;
+
+ for(Int_t i=0; i < fMaxIter; i++) {
+ // Rename a,b,c and adjust bounding interval d
+ if(fb*fc > 0.0) {
+ c = fA;
+ fc = fa;
+ d = fB - fA;
+ e = d;
+ }
+ if(TMath::Abs(fc) < TMath::Abs(fb)) {
+ fA = fB;
+ fB = c;
+ c = fA;
+ fa = fb;
+ fb = fc;
+ fc = fa;
+ }
+ Double_t Tol1 = 2.0*precision*TMath::Abs(fB) + 0.5*fTolerance;
+ Double_t xm = 0.5*(c-fB);
+ if(TMath::Abs(xm) <= Tol1 || fb == 0.0) {
+ fRoot = fB;
+ return true;
+ }
+ // Inverse quadratic interpolation
+ if(TMath::Abs(e) >= Tol1 && TMath::Abs(fa) > TMath::Abs(fb)) {
+ Double_t s = fb/fa;
+ Double_t p = 0.0;
+ Double_t q = 0.0;
+ if(fA == c) {
+ p = 2.0*xm*s;
+ q = 1.0 - s;
+ }
+ else {
+ q = fa/fc;
+ Double_t r = fb/fc;
+ p = s*(2.0*xm*q*(q-r)-(fB-fA)*(r-1.0));
+ q = (q-1.0)*(r-1.0)*(s-1.0);
+ }
+ // Check bounds
+ if(p > 0.0) q = -q;
+ p = TMath::Abs(p);
+ Double_t min1 = 3.0*xm*q-TMath::Abs(Tol1*q);
+ Double_t min2 = TMath::Abs(e*q);
+ if (2.0*p < TMath::Min(min1,min2)) {
+ // Interpolation
+ e = d;
+ d = p/q;
+ }
+ else {
+ // Bisection
+ d = xm;
+ e = d;
+ }
+ }
+ else {
+ // Bounds decreasing too slowly, use bisection
+ d = xm;
+ e = d;
+ }
+ // Move last guess to a
+ fA = fB;
+ fa = fb;
+ if (TMath::Abs(d) > Tol1) fB += d;
+ else {
+ if (xm >= 0.0) fB += TMath::Abs(Tol1);
+ else fB -= TMath::Abs(Tol1);
+ }
+ fb = theFunction(fB);
+ }
+ std::cerr << "EquationSolver::Brent: Number of iterations exceeded." << std::endl;
+ return false;
+}
+
+template <class Function>
+void EquationSolver<Function>::SetIntervalLimits(const Double_t Limit1, const Double_t Limit2) {
+ if(TMath::Abs(Limit1-Limit2) <= fTolerance) {
+ std::cerr << "EquationSolver::SetIntervalLimits: Interval must be wider than tolerance." << std::endl;
+ return;
+ }
+ if(Limit1 < Limit2) {
+ fA = Limit1;
+ fB = Limit2;
+ }
+ else {
+ fA = Limit2;
+ fB = Limit1;
+ }
+ return;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff