Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBtoXsgammaRootFinder.cpp
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 EvtBtoXsgammaRootFinder::EvtBtoXsgammaRootFinder() {}
44
45 EvtBtoXsgammaRootFinder::~EvtBtoXsgammaRootFinder( )
46 {}
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   double f1 = theFunc->value(lowerValue) - functionValue;
58   double f2 = theFunc->value(upperValue) - functionValue;
59
60   if ( f1*f2 > 0.0 ) {
61     report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;  
62     return 0;
63   }
64
65   // Already have root
66   if (fabs(f1) < precision) {
67     root = lowerValue;
68     return root;
69   }
70   if (fabs(f2) < precision) {
71     root = upperValue;
72     return root;
73   }
74   
75   // Orient search so that f(xLower) < 0
76   if (f1 < 0.0) {
77     xLower = lowerValue;
78     xUpper = upperValue;
79   } else {
80     xLower = upperValue;
81     xUpper = lowerValue;
82   }
83   
84   double rootGuess = 0.5*(lowerValue + upperValue);
85   double dxold = fabs(upperValue - lowerValue);
86   double dx = dxold;
87   
88   double f = theFunc->value(rootGuess) - functionValue;
89   
90   for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
91     
92       dxold = dx;
93       dx = 0.5*(xUpper-xLower);
94       rootGuess = xLower+dx;
95
96       // If change in root is negligible, take it as solution.
97       if (fabs(xLower - rootGuess) < precision) {
98         root = rootGuess;
99         return root;
100       }
101       
102       f = theFunc->value(rootGuess) - functionValue;
103  
104       if (f < 0.0) {
105         xLower = rootGuess;
106       } else {
107         xUpper = rootGuess;
108       }
109       
110   }
111   
112   report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
113                            <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!" 
114                            <<" Returning false."<<endl;
115   return 0;
116   
117 }
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   EvtItgAbsIntegrator *func1Integ = new EvtItgSimpsonIntegrator(*theFunc1, integ1Precision, maxLoop1);
127   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   theFunc1->setCoeff(1,0,lowerValue);
133   theFunc2->setCoeff(1,0,lowerValue);
134   
135   double f1 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
136   theFunc1->setCoeff(1,0,upperValue);
137   theFunc2->setCoeff(1,0,upperValue);
138   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   if ( f1*f2 > 0.0 ) {
144     report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;  
145     return false;
146   }
147
148   // Already have root
149   if (fabs(f1) < precision) {
150     root = lowerValue;
151     return root;
152   }
153   if (fabs(f2) < precision) {
154     root = upperValue;
155     return root;
156   }
157   
158   // Orient search so that f(xLower) < 0
159   if (f1 < 0.0) {
160     xLower = lowerValue;
161     xUpper = upperValue;
162   } else {
163     xLower = upperValue;
164     xUpper = lowerValue;
165   }
166   
167   double rootGuess = 0.5*(lowerValue + upperValue);
168   double dxold = fabs(upperValue - lowerValue);
169   double dx = dxold;
170   
171   theFunc1->setCoeff(1,0,rootGuess);
172   theFunc2->setCoeff(1,0,rootGuess);
173   double f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
174   
175   for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
176     
177     dxold = dx;
178     dx = 0.5*(xUpper-xLower);
179     rootGuess = xLower+dx;
180     
181     // If change in root is negligible, take it as solution.
182     if (fabs(xLower - rootGuess) < precision) {
183       root = rootGuess;
184       return root;
185     }
186     
187     theFunc1->setCoeff(1,0,rootGuess);
188     theFunc2->setCoeff(1,0,rootGuess);
189     f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
190     
191     if (f < 0.0) {
192       xLower = rootGuess;
193     } else {
194       xUpper = rootGuess;
195     }
196     
197   }
198   
199   report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
200                            <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!" 
201                            <<" Returning false."<<endl;
202   return 0;
203   
204 }
205
206