TEvtGen/EvtGenModels/EvtItgSimpsonIntegrator.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 //
   4 // Copyright Information: See EvtGen/COPYRIGHT
   5 //
   6 // Environment:
   7 //      This software is part of the EvtGen package developed jointly
   8 //      for the BaBar and CLEO collaborations.  If you use all or part
   9 //      of it, please give an appropriate acknowledgement.
  10 //
  11 // Module: EvtItgSimpsonIntegrator.hh
  12 //
  13 // Description:
  14 //      Abstraction of a generic function for use in integration methods elsewhere
  15 //      in this package. (Stolen and modified from
  16 //      the BaBar IntegrationUtils package - author: Phil Strother).
  17 //
  18 // Modification history:
  19 //
  20 //    Jane Tinslay                March 21, 2001       Module adapted for use in
  21 //                                                     EvtGen
  22 //
  23 //------------------------------------------------------------------------
  24 #include "EvtGenBase/EvtPatches.hh"
  25
  26 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
  27
  28 //-------------
  29 // C Headers --
  30 //-------------
  31 extern "C" {
  32 }
  33
  34 //---------------
  35 // C++ Headers --
  36 //---------------
  37
  38 #include <math.h>
  39
  40 //-------------------------------
  41 // Collaborating Class Headers --
  42 //-------------------------------
  43
  44 #include "EvtGenModels/EvtItgAbsFunction.hh"
  45 #include "EvtGenBase/EvtReport.hh"
  46 using std::endl;
  47
  48
  49 EvtItgSimpsonIntegrator::EvtItgSimpsonIntegrator(const EvtItgAbsFunction &theFunction, double precision, int maxLoop):
  50   EvtItgAbsIntegrator(theFunction),
  51   _precision(precision),
  52   _maxLoop(maxLoop)
  53 {}
  54
  55
  56 //--------------
  57 // Destructor --
  58 //--------------
  59
  60 EvtItgSimpsonIntegrator::~EvtItgSimpsonIntegrator()
  61 {}
  62
  63 double
  64 EvtItgSimpsonIntegrator::evaluateIt(double lower, double higher) const{
  65
  66   // report(INFO,"EvtGen")<<"in evaluate"<<endl;
  67   int j;
  68   double result(0.0);
  69   double s, st, ost(0.0);
  70   for (j=1;j<4;j++) {
  71     st = trapezoid(lower, higher, j, result);
  72     s = (4.0 * st - ost)/3.0;
  73     ost=st;
  74   }
  75
  76   double olds(s);
  77   st = trapezoid(lower, higher, j, result);
  78   s = (4.0 * st - ost)/3.0;
  79
  80   if (fabs(s - olds) < _precision*fabs(olds) || (s==0.0 && olds==0.0))     return s;
  81
  82   ost=st;
  83
  84   for (j=5;j<_maxLoop;j++){
  85
  86     st = trapezoid(lower, higher, j, result);
  87     s = (4.0 * st - ost)/3.0;
  88
  89     if (fabs(s - olds) < _precision*fabs(olds) || (s==0.0 && olds==0.0))    return s;
  90     olds=s;
  91     ost=st;
  92   }
  93
  94   report(ERROR,"EvtGen") << "Severe error in EvtItgSimpsonIntegrator.  Failed to converge after loop with 2**"
  95                  << _maxLoop << " calls to the integrand in." << endl;
  96
  97   return 0.0;
  98
  99 }