Line data Source code
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: EvtBtoXsgammaRootFinder.cc
12 : //
13 : // Description:
14 : // Root finders for EvtBtoXsgammaKagan module.
15 : //
16 : // Modification history:
17 : //
18 : // Jane Tinslay March 21, 2001 Module created
19 : //------------------------------------------------------------------------
20 : #include "EvtGenBase/EvtPatches.hh"
21 :
22 : #include "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
23 : #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
24 : #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
25 : #include "EvtGenBase/EvtReport.hh"
26 : #include <math.h>
27 : using std::endl;
28 :
29 : //-------------
30 : // C Headers --
31 : //-------------
32 : extern "C" {
33 : }
34 :
35 : //-----------------------------------------------------------------------
36 : // Local Macros, Typedefs, Structures, Unions and Forward Declarations --
37 : //-----------------------------------------------------------------------
38 :
39 : #define EVTITGROOTFINDER_MAXIT 100
40 : #define EVTITGROOTFINDER_RELATIVEPRECISION 1.0e-16
41 :
42 :
43 0 : EvtBtoXsgammaRootFinder::EvtBtoXsgammaRootFinder() {}
44 :
45 : EvtBtoXsgammaRootFinder::~EvtBtoXsgammaRootFinder( )
46 0 : {}
47 :
48 : double
49 : EvtBtoXsgammaRootFinder::GetRootSingleFunc(const EvtItgAbsFunction* theFunc, double functionValue, double lowerValue, double upperValue, double precision) {
50 :
51 : // Use the bisection to find the root.
52 : // Iterates until find root to the accuracy of precision
53 :
54 : double xLower = 0.0, xUpper = 0.0;
55 : double root=0;
56 :
57 0 : double f1 = theFunc->value(lowerValue) - functionValue;
58 0 : double f2 = theFunc->value(upperValue) - functionValue;
59 :
60 0 : if ( f1*f2 > 0.0 ) {
61 0 : report(Severity::Warning,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;
62 0 : return 0;
63 : }
64 :
65 : // Already have root
66 0 : if (fabs(f1) < precision) {
67 : root = lowerValue;
68 0 : return root;
69 : }
70 0 : if (fabs(f2) < precision) {
71 : root = upperValue;
72 0 : return root;
73 : }
74 :
75 : // Orient search so that f(xLower) < 0
76 0 : if (f1 < 0.0) {
77 : xLower = lowerValue;
78 : xUpper = upperValue;
79 0 : } else {
80 : xLower = upperValue;
81 : xUpper = lowerValue;
82 : }
83 :
84 0 : double rootGuess = 0.5*(lowerValue + upperValue);
85 0 : double dxold = fabs(upperValue - lowerValue);
86 : double dx = dxold;
87 :
88 0 : double f = theFunc->value(rootGuess) - functionValue;
89 :
90 0 : for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
91 :
92 : dxold = dx;
93 0 : dx = 0.5*(xUpper-xLower);
94 0 : rootGuess = xLower+dx;
95 :
96 : // If change in root is negligible, take it as solution.
97 0 : if (fabs(xLower - rootGuess) < precision) {
98 : root = rootGuess;
99 0 : return root;
100 : }
101 :
102 0 : f = theFunc->value(rootGuess) - functionValue;
103 :
104 0 : if (f < 0.0) {
105 : xLower = rootGuess;
106 0 : } else {
107 : xUpper = rootGuess;
108 : }
109 :
110 : }
111 :
112 0 : report(Severity::Warning,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
113 0 : <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!"
114 0 : <<" Returning false."<<endl;
115 0 : return 0;
116 :
117 0 : }
118 :
119 : double
120 : EvtBtoXsgammaRootFinder::GetGaussIntegFcnRoot(EvtItgAbsFunction *theFunc1, EvtItgAbsFunction *theFunc2, double integ1Precision, double integ2Precision, int maxLoop1, int maxLoop2, double integLower, double integUpper, double lowerValue, double upperValue, double precision) {
121 :
122 : // Use the bisection to find the root.
123 : // Iterates until find root to the accuracy of precision
124 :
125 : //Need to work with integrators
126 0 : EvtItgAbsIntegrator *func1Integ = new EvtItgSimpsonIntegrator(*theFunc1, integ1Precision, maxLoop1);
127 0 : EvtItgAbsIntegrator *func2Integ = new EvtItgSimpsonIntegrator(*theFunc2, integ2Precision, maxLoop2);
128 :
129 :
130 : //coefficient 1 of the integrators is the root to be found
131 : //need to set this to lower value to start off with
132 0 : theFunc1->setCoeff(1,0,lowerValue);
133 0 : theFunc2->setCoeff(1,0,lowerValue);
134 :
135 0 : double f1 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
136 0 : theFunc1->setCoeff(1,0,upperValue);
137 0 : theFunc2->setCoeff(1,0,upperValue);
138 0 : double f2 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
139 :
140 : double xLower = 0.0, xUpper = 0.0;
141 : double root=0;
142 :
143 0 : if ( f1*f2 > 0.0 ) {
144 0 : report(Severity::Warning,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;
145 0 : return false;
146 : }
147 :
148 : // Already have root
149 0 : if (fabs(f1) < precision) {
150 : root = lowerValue;
151 0 : return root;
152 : }
153 0 : if (fabs(f2) < precision) {
154 : root = upperValue;
155 0 : return root;
156 : }
157 :
158 : // Orient search so that f(xLower) < 0
159 0 : if (f1 < 0.0) {
160 : xLower = lowerValue;
161 : xUpper = upperValue;
162 0 : } else {
163 : xLower = upperValue;
164 : xUpper = lowerValue;
165 : }
166 :
167 0 : double rootGuess = 0.5*(lowerValue + upperValue);
168 0 : double dxold = fabs(upperValue - lowerValue);
169 : double dx = dxold;
170 :
171 0 : theFunc1->setCoeff(1,0,rootGuess);
172 0 : theFunc2->setCoeff(1,0,rootGuess);
173 0 : double f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
174 :
175 0 : for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
176 :
177 : dxold = dx;
178 0 : dx = 0.5*(xUpper-xLower);
179 0 : rootGuess = xLower+dx;
180 :
181 : // If change in root is negligible, take it as solution.
182 0 : if (fabs(xLower - rootGuess) < precision) {
183 : root = rootGuess;
184 0 : return root;
185 : }
186 :
187 0 : theFunc1->setCoeff(1,0,rootGuess);
188 0 : theFunc2->setCoeff(1,0,rootGuess);
189 0 : f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
190 :
191 0 : if (f < 0.0) {
192 : xLower = rootGuess;
193 0 : } else {
194 : xUpper = rootGuess;
195 : }
196 :
197 : }
198 :
199 0 : report(Severity::Warning,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
200 0 : <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!"
201 0 : <<" Returning false."<<endl;
202 0 : return 0;
203 :
204 0 : }
205 :
206 :
|
