00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "EvtGenBase/EvtPatches.hh"
00025
00026 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
00027
00028
00029
00030
00031 extern "C" {
00032 }
00033
00034
00035
00036
00037
00038 #include <math.h>
00039
00040
00041
00042
00043
00044 #include "EvtGenModels/EvtItgAbsFunction.hh"
00045 #include "EvtGenBase/EvtReport.hh"
00046 using std::endl;
00047
00048
00049 EvtItgSimpsonIntegrator::EvtItgSimpsonIntegrator(const EvtItgAbsFunction &theFunction, double precision, int maxLoop):
00050 EvtItgAbsIntegrator(theFunction),
00051 _precision(precision),
00052 _maxLoop(maxLoop)
00053 {}
00054
00055
00056
00057
00058
00059
00060 EvtItgSimpsonIntegrator::~EvtItgSimpsonIntegrator()
00061 {}
00062
00063 double
00064 EvtItgSimpsonIntegrator::evaluateIt(double lower, double higher) const{
00065
00066
00067 int j;
00068 double result(0.0);
00069 double s, st, ost(0.0);
00070 for (j=1;j<4;j++) {
00071 st = trapezoid(lower, higher, j, result);
00072 s = (4.0 * st - ost)/3.0;
00073 ost=st;
00074 }
00075
00076 double olds(s);
00077 st = trapezoid(lower, higher, j, result);
00078 s = (4.0 * st - ost)/3.0;
00079
00080 if (fabs(s - olds) < _precision*fabs(olds) || (s==0.0 && olds==0.0)) return s;
00081
00082 ost=st;
00083
00084 for (j=5;j<_maxLoop;j++){
00085
00086 st = trapezoid(lower, higher, j, result);
00087 s = (4.0 * st - ost)/3.0;
00088
00089 if (fabs(s - olds) < _precision*fabs(olds) || (s==0.0 && olds==0.0)) return s;
00090 olds=s;
00091 ost=st;
00092 }
00093
00094 report(ERROR,"EvtGen") << "Severe error in EvtItgSimpsonIntegrator. Failed to converge after loop with 2**"
00095 << _maxLoop << " calls to the integrand in." << endl;
00096
00097 return 0.0;
00098
00099 }