Page MenuHomeHEPForge
Files F19250584

EvtBtoXsgammaRootFinder.cpp
No OneTemporary
Actions

Size
5 KB
Referenced Files
None
Subscribers
None

EvtBtoXsgammaRootFinder.cpp
View Options

//--------------------------------------------------------------------------
//
//
// 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 "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
#include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
#include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
#include "EvtGenBase/EvtReport.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;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 5:42 AM (4 h, 8 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6494383
Default Alt Text
EvtBtoXsgammaRootFinder.cpp (5 KB)

Event Timeline

HEPForge.org · .