00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include "EvtGenModels/EvtItgAbsIntegrator.hh"
00024
00025
00026
00027
00028 extern "C" {
00029 }
00030
00031 #include <math.h>
00032 #include <iostream>
00033
00034 #include "EvtGenBase/EvtReport.hh"
00035 #include "EvtGenModels/EvtItgAbsFunction.hh"
00036 using std::endl;
00037
00038
00039 EvtItgAbsIntegrator::EvtItgAbsIntegrator(const EvtItgAbsFunction &theFunction):
00040 _myFunction(theFunction)
00041 {}
00042
00043 EvtItgAbsIntegrator::~EvtItgAbsIntegrator()
00044 {}
00045
00046 double
00047 EvtItgAbsIntegrator::normalisation() const {
00048 return evaluateIt(_myFunction.lowerRange(), _myFunction.upperRange());
00049 }
00050
00051 double
00052 EvtItgAbsIntegrator::evaluate(double lower, double upper) const{
00053
00054 double newLower(lower), newUpper(upper);
00055
00056 boundsCheck(newLower, newUpper);
00057
00058 return evaluateIt(newLower, newUpper);
00059 }
00060
00061 double
00062 EvtItgAbsIntegrator::trapezoid(double lower, double higher, int n, double &result) const {
00063
00064 if (n==1) return 0.5*(higher-lower)*(_myFunction(lower) + _myFunction(higher));
00065
00066 int it, j;
00067
00068 for (it=1, j=1;j<n-1;j++) it <<=1;
00069
00070 double itDouble(it);
00071
00072 double sum(0.0);
00073
00074 double deltaX((higher - lower)/itDouble);
00075
00076 double x(lower + 0.5* deltaX);
00077
00078 for (j=1;j<=it;j++){
00079 sum+=_myFunction(x);
00080 x+=deltaX;
00081 }
00082
00083 result = 0.5*(result+(higher - lower)*sum/itDouble);
00084
00085 return result;
00086 }
00087
00088 void
00089 EvtItgAbsIntegrator::boundsCheck(double &lower, double &upper) const{
00090
00091 if (lower < _myFunction.lowerRange() ) {
00092 report(WARNING,"EvtGen") << "Warning in EvtItgAbsIntegrator::evaluate. Lower bound " << lower << " of integral "
00093 << " is less than lower bound " << _myFunction.lowerRange()
00094 << " of function. No contribution from this range will be counted." << endl;
00095 lower = _myFunction.lowerRange();
00096 }
00097
00098 if (upper > _myFunction.upperRange() ) {
00099 report(WARNING,"EvtGen") << "Warning in EvtItgAbsIntegrator::evaluate. Upper bound " << upper << " of integral "
00100 << " is greater than upper bound " << _myFunction.upperRange()
00101 << " of function. No contribution from this range will be counted." << endl;
00102 upper = _myFunction.upperRange();
00103 }
00104
00105 }