[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtItgSimpsonIntegrator.cxx

//--------------------------------------------------------------------------
//
//
// 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: EvtItgSimpsonIntegrator.hh
//
// Description:
//      Abstraction of a generic function for use in integration methods elsewhere
//      in this package. (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 "EvtGenBase/EvtPatches.hh"

#include "EvtGenModels/EvtItgSimpsonIntegrator.hh"

//-------------
// C Headers --
//-------------
extern "C" {
}

//---------------
// C++ Headers --
//---------------

#include <math.h>

//-------------------------------
// Collaborating Class Headers --
//-------------------------------

#include "EvtGenModels/EvtItgAbsFunction.hh"
#include "EvtGenBase/EvtReport.hh"
using std::endl;


EvtItgSimpsonIntegrator::EvtItgSimpsonIntegrator(const EvtItgAbsFunction &theFunction, double precision, int maxLoop):
  EvtItgAbsIntegrator(theFunction),
  _precision(precision),
  _maxLoop(maxLoop)
{}


//--------------
// Destructor --
//--------------

EvtItgSimpsonIntegrator::~EvtItgSimpsonIntegrator()
{}

double
EvtItgSimpsonIntegrator::evaluateIt(double lower, double higher) const{
  
  // report(INFO,"EvtGen")<<"in evaluate"<<endl;
  int j;
  double result(0.0);
  double s, st, ost(0.0);
  for (j=1;j<4;j++) {
    st = trapezoid(lower, higher, j, result);
    s = (4.0 * st - ost)/3.0;
    ost=st;
  }

  double olds(s);
  st = trapezoid(lower, higher, j, result);
  s = (4.0 * st - ost)/3.0;

  if (fabs(s - olds) < _precision*fabs(olds) || (s==0.0 && olds==0.0))     return s;
  
  ost=st;

  for (j=5;j<_maxLoop;j++){

    st = trapezoid(lower, higher, j, result);
    s = (4.0 * st - ost)/3.0;
    
    if (fabs(s - olds) < _precision*fabs(olds) || (s==0.0 && olds==0.0))    return s;
    olds=s;
    ost=st;
  }
  
  report(ERROR,"EvtGen") << "Severe error in EvtItgSimpsonIntegrator.  Failed to converge after loop with 2**"
		 << _maxLoop << " calls to the integrand in." << endl;
  
  return 0.0;
    
}
Commit	Line	Data
da0e9ce3	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	}