1 //--------------------------------------------------------------------------
3 // Copyright Information: See EvtGen/COPYRIGHT
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.
10 // Module: EvtItgIntegrator.cc
13 // Simpson integrator (Stolen and modified from
14 // the BaBar IntegrationUtils package - author: Phil Strother).
16 // Modification history:
18 // Jane Tinslay March 21, 2001 Module adapted for use in
21 //------------------------------------------------------------------------
22 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenModels/EvtItgAbsIntegrator.hh"
34 #include "EvtGenBase/EvtReport.hh"
35 #include "EvtGenModels/EvtItgAbsFunction.hh"
39 EvtItgAbsIntegrator::EvtItgAbsIntegrator(const EvtItgAbsFunction &theFunction):
40 _myFunction(theFunction)
43 EvtItgAbsIntegrator::~EvtItgAbsIntegrator()
47 EvtItgAbsIntegrator::normalisation() const {
48 return evaluateIt(_myFunction.lowerRange(), _myFunction.upperRange());
52 EvtItgAbsIntegrator::evaluate(double lower, double upper) const{
54 double newLower(lower), newUpper(upper);
56 boundsCheck(newLower, newUpper);
58 return evaluateIt(newLower, newUpper);
62 EvtItgAbsIntegrator::trapezoid(double lower, double higher, int n, double &result) const {
64 if (n==1) return 0.5*(higher-lower)*(_myFunction(lower) + _myFunction(higher));
68 for (it=1, j=1;j<n-1;j++) it <<=1;
74 double deltaX((higher - lower)/itDouble);
76 double x(lower + 0.5* deltaX);
83 result = 0.5*(result+(higher - lower)*sum/itDouble);
89 EvtItgAbsIntegrator::boundsCheck(double &lower, double &upper) const{
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();
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();