EvtItgAbsIntegrator.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Copyright Information: See EvtGen/COPYRIGHT
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Module: EvtItgIntegrator.cc
//
// Description:
//      Simpson integrator (Stolen and modified from
//      the BaBar IntegrationUtils package - author: Phil Strother).
//
// Modification history:
//
//    Jane Tinslay                March 21, 2001       Module adapted for use in
//                                                     EvtGen
//
//------------------------------------------------------------------------
#include "EvtItgAbsIntegrator.hh"
#include "../EvtGenBase/EvtPatches.hh"
 
//-------------
// C Headers --
//-------------
extern "C" {}
 
#include <iostream>
#include <math.h>
 
#include "../EvtGenBase/EvtReport.hh"
#include "EvtItgAbsFunction.hh"
using std::endl;
 
EvtItgAbsIntegrator::EvtItgAbsIntegrator( const EvtItgAbsFunction& theFunction )
    : _myFunction( theFunction ) {}
 
EvtItgAbsIntegrator::~EvtItgAbsIntegrator() {}
 
double EvtItgAbsIntegrator::normalisation() const {
  return evaluateIt( _myFunction.lowerRange(), _myFunction.upperRange() );
}
 
double EvtItgAbsIntegrator::evaluate( double lower, double upper ) const {
 
  double newLower( lower ), newUpper( upper );
 
  boundsCheck( newLower, newUpper );
 
  return evaluateIt( newLower, newUpper );
}
 
double EvtItgAbsIntegrator::trapezoid( double lower, double higher, int n,
                                       double& result ) const {
 
  if ( n == 1 )
    return 0.5 * ( higher - lower ) * ( _myFunction( lower ) + _myFunction( higher ) );
 
  int it, j;
 
  for ( it = 1, j = 1; j < n - 1; j++ ) it <<= 1;
 
  double itDouble( it );
 
  double sum( 0.0 );
 
  double deltaX( ( higher - lower ) / itDouble );
 
  double x( lower + 0.5 * deltaX );
 
  for ( j = 1; j <= it; j++ )
  {
    sum += _myFunction( x );
    x += deltaX;
  }
 
  result = 0.5 * ( result + ( higher - lower ) * sum / itDouble );
 
  return result;
}
 
void EvtItgAbsIntegrator::boundsCheck( double& lower, double& upper ) const {
 
  if ( lower < _myFunction.lowerRange() )
  {
    report( WARNING, "EvtGen" )
        << "Warning in EvtItgAbsIntegrator::evaluate.  Lower bound " << lower
        << " of integral "
        << " is less than lower bound " << _myFunction.lowerRange()
        << " of function.  No contribution from this range will be counted." << endl;
    lower = _myFunction.lowerRange();
  }
 
  if ( upper > _myFunction.upperRange() )
  {
    report( WARNING, "EvtGen" )
        << "Warning in EvtItgAbsIntegrator::evaluate.  Upper bound " << upper
        << " of integral "
        << " is greater than upper bound " << _myFunction.upperRange()
        << " of function.  No contribution from this range will be counted." << endl;
    upper = _myFunction.upperRange();
  }
}