TEvtGen/EvtGenModels/EvtItgAbsIntegrator.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Copyright Information: See EvtGen/COPYRIGHT
   4 //
   5 // Environment:
   6 //      This software is part of the EvtGen package developed jointly
   7 //      for the BaBar and CLEO collaborations.  If you use all or part
   8 //      of it, please give an appropriate acknowledgement.
   9 //
  10 // Module: EvtItgIntegrator.cc
  11 //
  12 // Description:
  13 //      Simpson integrator (Stolen and modified from
  14 //      the BaBar IntegrationUtils package - author: Phil Strother).
  15 //
  16 // Modification history:
  17 //
  18 //    Jane Tinslay                March 21, 2001       Module adapted for use in
  19 //                                                     EvtGen
  20 //
  21 //------------------------------------------------------------------------
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include "EvtGenModels/EvtItgAbsIntegrator.hh"
  24
  25 //-------------
  26 // C Headers --
  27 //-------------
  28 extern "C" {
  29 }
  30
  31 #include <math.h>
  32 #include <iostream>
  33
  34 #include "EvtGenBase/EvtReport.hh"
  35 #include "EvtGenModels/EvtItgAbsFunction.hh"
  36 using std::endl;
  37
  38
  39 EvtItgAbsIntegrator::EvtItgAbsIntegrator(const EvtItgAbsFunction &theFunction):
  40   _myFunction(theFunction)
  41 {}
  42
  43 EvtItgAbsIntegrator::~EvtItgAbsIntegrator()
  44 {}
  45
  46 double
  47 EvtItgAbsIntegrator::normalisation() const {
  48   return evaluateIt(_myFunction.lowerRange(), _myFunction.upperRange());
  49 }
  50
  51 double
  52 EvtItgAbsIntegrator::evaluate(double lower, double upper) const{
  53
  54   double newLower(lower), newUpper(upper);
  55
  56   boundsCheck(newLower, newUpper);
  57
  58   return evaluateIt(newLower, newUpper);
  59 }
  60
  61 double
  62 EvtItgAbsIntegrator::trapezoid(double lower, double higher, int n, double &result) const {
  63
  64   if (n==1) return 0.5*(higher-lower)*(_myFunction(lower) + _myFunction(higher));
  65
  66   int it, j;
  67
  68   for (it=1, j=1;j<n-1;j++) it <<=1;
  69
  70   double itDouble(it);
  71
  72   double sum(0.0);
  73
  74   double deltaX((higher - lower)/itDouble);
  75
  76   double x(lower + 0.5* deltaX);
  77
  78   for (j=1;j<=it;j++){
  79     sum+=_myFunction(x);
  80     x+=deltaX;
  81   }
  82
  83   result = 0.5*(result+(higher - lower)*sum/itDouble);
  84
  85   return result;
  86 }
  87
  88 void
  89 EvtItgAbsIntegrator::boundsCheck(double &lower, double &upper) const{
  90
  91   if (lower < _myFunction.lowerRange() ) {
  92     report(WARNING,"EvtGen") << "Warning in EvtItgAbsIntegrator::evaluate.  Lower bound " << lower << " of integral "
  93                     << " is less than lower bound " << _myFunction.lowerRange()
  94                     << " of function.  No contribution from this range will be counted." << endl;
  95     lower = _myFunction.lowerRange();
  96   }
  97
  98   if (upper > _myFunction.upperRange() ) {
  99     report(WARNING,"EvtGen") << "Warning in EvtItgAbsIntegrator::evaluate.  Upper bound " << upper << " of integral "
 100                     << " is greater than upper bound " << _myFunction.upperRange()
 101                     << " of function.  No contribution from this range will be counted." << endl;
 102     upper = _myFunction.upperRange();
 103   }
 104
 105 }