]>
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 | } |