Updated macro

[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / EquationSolver.icc

/*                                                                            \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