]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //-------------------------------------------------------------------------- |
2 | // | |
3 | // Copyright Information: See EvtGen/COPYRIGHT | |
4 | // | |
5 | // Environment: | |
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. | |
9 | // | |
10 | // Module: EvtItgIntegrator.cc | |
11 | // | |
12 | // Description: | |
13 | // Simpson integrator (Stolen and modified from | |
14 | // the BaBar IntegrationUtils package - author: Phil Strother). | |
15 | // | |
16 | // Modification history: | |
17 | // | |
18 | // Jane Tinslay March 21, 2001 Module adapted for use in | |
19 | // EvtGen | |
20 | // | |
21 | //------------------------------------------------------------------------ | |
22 | #include "EvtGenBase/EvtPatches.hh" | |
23 | #include "EvtGenModels/EvtItgAbsIntegrator.hh" | |
24 | ||
25 | //------------- | |
26 | // C Headers -- | |
27 | //------------- | |
28 | extern "C" { | |
29 | } | |
30 | ||
31 | #include <math.h> | |
32 | #include <iostream> | |
33 | ||
34 | #include "EvtGenBase/EvtReport.hh" | |
35 | #include "EvtGenModels/EvtItgAbsFunction.hh" | |
36 | using std::endl; | |
37 | ||
38 | ||
39 | EvtItgAbsIntegrator::EvtItgAbsIntegrator(const EvtItgAbsFunction &theFunction): | |
40 | _myFunction(theFunction) | |
41 | {} | |
42 | ||
43 | EvtItgAbsIntegrator::~EvtItgAbsIntegrator() | |
44 | {} | |
45 | ||
46 | double | |
47 | EvtItgAbsIntegrator::normalisation() const { | |
48 | return evaluateIt(_myFunction.lowerRange(), _myFunction.upperRange()); | |
49 | } | |
50 | ||
51 | double | |
52 | EvtItgAbsIntegrator::evaluate(double lower, double upper) const{ | |
53 | ||
54 | double newLower(lower), newUpper(upper); | |
55 | ||
56 | boundsCheck(newLower, newUpper); | |
57 | ||
58 | return evaluateIt(newLower, newUpper); | |
59 | } | |
60 | ||
61 | double | |
62 | EvtItgAbsIntegrator::trapezoid(double lower, double higher, int n, double &result) const { | |
63 | ||
64 | if (n==1) return 0.5*(higher-lower)*(_myFunction(lower) + _myFunction(higher)); | |
65 | ||
66 | int it, j; | |
67 | ||
68 | for (it=1, j=1;j<n-1;j++) it <<=1; | |
69 | ||
70 | double itDouble(it); | |
71 | ||
72 | double sum(0.0); | |
73 | ||
74 | double deltaX((higher - lower)/itDouble); | |
75 | ||
76 | double x(lower + 0.5* deltaX); | |
77 | ||
78 | for (j=1;j<=it;j++){ | |
79 | sum+=_myFunction(x); | |
80 | x+=deltaX; | |
81 | } | |
82 | ||
83 | result = 0.5*(result+(higher - lower)*sum/itDouble); | |
84 | ||
85 | return result; | |
86 | } | |
87 | ||
88 | void | |
89 | EvtItgAbsIntegrator::boundsCheck(double &lower, double &upper) const{ | |
90 | ||
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(); | |
96 | } | |
97 | ||
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(); | |
103 | } | |
104 | ||
105 | } |