Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtItgSimpsonIntegrator.cpp
CommitLineData
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//-------------
31extern "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"
46using std::endl;
47
48
49EvtItgSimpsonIntegrator::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
60EvtItgSimpsonIntegrator::~EvtItgSimpsonIntegrator()
61{}
62
63double
64EvtItgSimpsonIntegrator::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}