EvtBtoXsgammaRootFinder.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
//
// Copyright Information: See EvtGen/COPYRIGHT
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Module: EvtBtoXsgammaRootFinder.cc
//
// Description:
//      Root finders for EvtBtoXsgammaKagan module.
//
// Modification history:
//
//      Jane Tinslay       March 21, 2001       Module created
//------------------------------------------------------------------------
#include "../EvtGenBase/EvtPatches.hh"
 
#include "../EvtGenBase/EvtReport.hh"
#include "EvtBtoXsgammaRootFinder.hh"
#include "EvtItgSimpsonIntegrator.hh"
#include "EvtItgTwoCoeffFcn.hh"
#include <math.h>
using std::endl;
 
//-------------
// C Headers --
//-------------
extern "C" {}
 
//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------
 
#define EVTITGROOTFINDER_MAXIT 100
#define EVTITGROOTFINDER_RELATIVEPRECISION 1.0e-16
 
EvtBtoXsgammaRootFinder::EvtBtoXsgammaRootFinder() {}
 
EvtBtoXsgammaRootFinder::~EvtBtoXsgammaRootFinder() {}
 
double EvtBtoXsgammaRootFinder::GetRootSingleFunc( const EvtItgAbsFunction* theFunc,
                                                   double functionValue, double lowerValue,
                                                   double upperValue, double precision ) {
 
  // Use the bisection to find the root.
  // Iterates until find root to the accuracy of precision
 
  double xLower = 0.0, xUpper = 0.0;
  double root = 0;
 
  double f1 = theFunc->value( lowerValue ) - functionValue;
  double f2 = theFunc->value( upperValue ) - functionValue;
 
  if ( f1 * f2 > 0.0 )
  {
    report( WARNING, "EvtGen" ) << "EvtBtoXsgammaRootFinder: No root in specified range !"
                                << endl;
    return 0;
  }
 
  // Already have root
  if ( fabs( f1 ) < precision )
  {
    root = lowerValue;
    return root;
  }
  if ( fabs( f2 ) < precision )
  {
    root = upperValue;
    return root;
  }
 
  // Orient search so that f(xLower) < 0
  if ( f1 < 0.0 )
  {
    xLower = lowerValue;
    xUpper = upperValue;
  }
  else
  {
    xLower = upperValue;
    xUpper = lowerValue;
  }
 
  double rootGuess = 0.5 * ( lowerValue + upperValue );
  double dxold     = fabs( upperValue - lowerValue );
  double dx        = dxold;
 
  double f = theFunc->value( rootGuess ) - functionValue;
 
  for ( int j = 0; j < EVTITGROOTFINDER_MAXIT; j++ )
  {
 
    dxold     = dx;
    dx        = 0.5 * ( xUpper - xLower );
    rootGuess = xLower + dx;
 
    // If change in root is negligible, take it as solution.
    if ( fabs( xLower - rootGuess ) < precision )
    {
      root = rootGuess;
      return root;
    }
 
    f = theFunc->value( rootGuess ) - functionValue;
 
    if ( f < 0.0 ) { xLower = rootGuess; }
    else { xUpper = rootGuess; }
  }
 
  report( WARNING, "EvtGen" ) << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
                              << "in EvtBtoXsgammaRootFinder::foundRoot exceeded!"
                              << " Returning false." << endl;
  return 0;
}
 
double EvtBtoXsgammaRootFinder::GetGaussIntegFcnRoot(
    EvtItgAbsFunction* theFunc1, EvtItgAbsFunction* theFunc2, double integ1Precision,
    double integ2Precision, int maxLoop1, int maxLoop2, double integLower, double integUpper,
    double lowerValue, double upperValue, double precision ) {
 
  // Use the bisection to find the root.
  // Iterates until find root to the accuracy of precision
 
  // Need to work with integrators
  EvtItgAbsIntegrator* func1Integ =
      new EvtItgSimpsonIntegrator( *theFunc1, integ1Precision, maxLoop1 );
  EvtItgAbsIntegrator* func2Integ =
      new EvtItgSimpsonIntegrator( *theFunc2, integ2Precision, maxLoop2 );
 
  // coefficient 1 of the integrators is the root to be found
  // need to set this to lower value to start off with
  theFunc1->setCoeff( 1, 0, lowerValue );
  theFunc2->setCoeff( 1, 0, lowerValue );
 
  double f1 = func1Integ->evaluate( integLower, integUpper ) -
              theFunc2->getCoeff( 1, 2 ) * func2Integ->evaluate( integLower, integUpper );
  theFunc1->setCoeff( 1, 0, upperValue );
  theFunc2->setCoeff( 1, 0, upperValue );
  double f2 = func1Integ->evaluate( integLower, integUpper ) -
              theFunc2->getCoeff( 1, 2 ) * func2Integ->evaluate( integLower, integUpper );
 
  double xLower = 0.0, xUpper = 0.0;
  double root = 0;
 
  if ( f1 * f2 > 0.0 )
  {
    report( WARNING, "EvtGen" ) << "EvtBtoXsgammaRootFinder: No root in specified range !"
                                << endl;
    return false;
  }
 
  // Already have root
  if ( fabs( f1 ) < precision )
  {
    root = lowerValue;
    return root;
  }
  if ( fabs( f2 ) < precision )
  {
    root = upperValue;
    return root;
  }
 
  // Orient search so that f(xLower) < 0
  if ( f1 < 0.0 )
  {
    xLower = lowerValue;
    xUpper = upperValue;
  }
  else
  {
    xLower = upperValue;
    xUpper = lowerValue;
  }
 
  double rootGuess = 0.5 * ( lowerValue + upperValue );
  double dxold     = fabs( upperValue - lowerValue );
  double dx        = dxold;
 
  theFunc1->setCoeff( 1, 0, rootGuess );
  theFunc2->setCoeff( 1, 0, rootGuess );
  double f = func1Integ->evaluate( integLower, integUpper ) -
             theFunc2->getCoeff( 1, 2 ) * func2Integ->evaluate( integLower, integUpper );
 
  for ( int j = 0; j < EVTITGROOTFINDER_MAXIT; j++ )
  {
 
    dxold     = dx;
    dx        = 0.5 * ( xUpper - xLower );
    rootGuess = xLower + dx;
 
    // If change in root is negligible, take it as solution.
    if ( fabs( xLower - rootGuess ) < precision )
    {
      root = rootGuess;
      return root;
    }
 
    theFunc1->setCoeff( 1, 0, rootGuess );
    theFunc2->setCoeff( 1, 0, rootGuess );
    f = func1Integ->evaluate( integLower, integUpper ) -
        theFunc2->getCoeff( 1, 2 ) * func2Integ->evaluate( integLower, integUpper );
 
    if ( f < 0.0 ) { xLower = rootGuess; }
    else { xUpper = rootGuess; }
  }
 
  report( WARNING, "EvtGen" ) << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
                              << "in EvtBtoXsgammaRootFinder::foundRoot exceeded!"
                              << " Returning false." << endl;
  return 0;
}