Line data    Source code

       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           0 :   _myFunction(theFunction)
      41           0 : {}
      42             : 
      43             : EvtItgAbsIntegrator::~EvtItgAbsIntegrator()
      44           0 : {}
      45             :   
      46             : double 
      47             : EvtItgAbsIntegrator::normalisation() const {
      48           0 :   return evaluateIt(_myFunction.lowerRange(), _myFunction.upperRange());
      49             : }
      50             : 
      51             : double
      52             : EvtItgAbsIntegrator::evaluate(double lower, double upper) const{
      53             : 
      54           0 :   double newLower(lower), newUpper(upper);
      55             : 
      56           0 :   boundsCheck(newLower, newUpper);
      57             : 
      58           0 :   return evaluateIt(newLower, newUpper);
      59           0 : }
      60             : 
      61             : double 
      62             : EvtItgAbsIntegrator::trapezoid(double lower, double higher, int n, double &result) const {
      63             : 
      64           0 :   if (n==1) return 0.5*(higher-lower)*(_myFunction(lower) + _myFunction(higher));
      65             :   
      66             :   int it, j;
      67             :   
      68           0 :   for (it=1, j=1;j<n-1;j++) it <<=1;
      69             :   
      70           0 :   double itDouble(it);
      71             :   
      72             :   double sum(0.0);
      73             : 
      74           0 :   double deltaX((higher - lower)/itDouble);
      75             :   
      76           0 :   double x(lower + 0.5* deltaX);
      77             :     
      78           0 :   for (j=1;j<=it;j++){
      79           0 :     sum+=_myFunction(x);
      80           0 :     x+=deltaX;
      81             :   }
      82             :   
      83           0 :   result = 0.5*(result+(higher - lower)*sum/itDouble);
      84             : 
      85             :   return result;
      86           0 : }
      87             : 
      88             : void 
      89             : EvtItgAbsIntegrator::boundsCheck(double &lower, double &upper) const{
      90             : 
      91           0 :   if (lower < _myFunction.lowerRange() ) {
      92           0 :     report(Severity::Warning,"EvtGen") << "Warning in EvtItgAbsIntegrator::evaluate.  Lower bound " << lower << " of integral " 
      93           0 :                     << " is less than lower bound " << _myFunction.lowerRange() 
      94           0 :                     << " of function.  No contribution from this range will be counted." << endl;
      95           0 :     lower = _myFunction.lowerRange();
      96           0 :   }
      97             : 
      98           0 :   if (upper > _myFunction.upperRange() ) {
      99           0 :     report(Severity::Warning,"EvtGen") << "Warning in EvtItgAbsIntegrator::evaluate.  Upper bound " << upper << " of integral "
     100           0 :                     << " is greater than upper bound " << _myFunction.upperRange() 
     101           0 :                     << " of function.  No contribution from this range will be counted." << endl;  
     102           0 :     upper = _myFunction.upperRange();
     103           0 :   }
     104             : 
     105           0 : }