00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "EvtGenBase/EvtPatches.hh"
00021
00022 #include "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
00023 #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
00024 #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
00025 #include "EvtGenBase/EvtReport.hh"
00026 #include <math.h>
00027 using std::endl;
00028
00029
00030
00031
00032 extern "C" {
00033 }
00034
00035
00036
00037
00038
00039 #define EVTITGROOTFINDER_MAXIT 100
00040 #define EVTITGROOTFINDER_RELATIVEPRECISION 1.0e-16
00041
00042
00043 EvtBtoXsgammaRootFinder::EvtBtoXsgammaRootFinder() {}
00044
00045 EvtBtoXsgammaRootFinder::~EvtBtoXsgammaRootFinder( )
00046 {}
00047
00048 double
00049 EvtBtoXsgammaRootFinder::GetRootSingleFunc(const EvtItgAbsFunction* theFunc, double functionValue, double lowerValue, double upperValue, double precision) {
00050
00051
00052
00053
00054 double xLower = 0.0, xUpper = 0.0;
00055 double root=0;
00056
00057 double f1 = theFunc->value(lowerValue) - functionValue;
00058 double f2 = theFunc->value(upperValue) - functionValue;
00059
00060 if ( f1*f2 > 0.0 ) {
00061 report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;
00062 return 0;
00063 }
00064
00065
00066 if (fabs(f1) < precision) {
00067 root = lowerValue;
00068 return root;
00069 }
00070 if (fabs(f2) < precision) {
00071 root = upperValue;
00072 return root;
00073 }
00074
00075
00076 if (f1 < 0.0) {
00077 xLower = lowerValue;
00078 xUpper = upperValue;
00079 } else {
00080 xLower = upperValue;
00081 xUpper = lowerValue;
00082 }
00083
00084 double rootGuess = 0.5*(lowerValue + upperValue);
00085 double dxold = fabs(upperValue - lowerValue);
00086 double dx = dxold;
00087
00088 double f = theFunc->value(rootGuess) - functionValue;
00089
00090 for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
00091
00092 dxold = dx;
00093 dx = 0.5*(xUpper-xLower);
00094 rootGuess = xLower+dx;
00095
00096
00097 if (fabs(xLower - rootGuess) < precision) {
00098 root = rootGuess;
00099 return root;
00100 }
00101
00102 f = theFunc->value(rootGuess) - functionValue;
00103
00104 if (f < 0.0) {
00105 xLower = rootGuess;
00106 } else {
00107 xUpper = rootGuess;
00108 }
00109
00110 }
00111
00112 report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
00113 <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!"
00114 <<" Returning false."<<endl;
00115 return 0;
00116
00117 }
00118
00119 double
00120 EvtBtoXsgammaRootFinder::GetGaussIntegFcnRoot(EvtItgAbsFunction *theFunc1, EvtItgAbsFunction *theFunc2, double integ1Precision, double integ2Precision, int maxLoop1, int maxLoop2, double integLower, double integUpper, double lowerValue, double upperValue, double precision) {
00121
00122
00123
00124
00125
00126 EvtItgAbsIntegrator *func1Integ = new EvtItgSimpsonIntegrator(*theFunc1, integ1Precision, maxLoop1);
00127 EvtItgAbsIntegrator *func2Integ = new EvtItgSimpsonIntegrator(*theFunc2, integ2Precision, maxLoop2);
00128
00129
00130
00131
00132 theFunc1->setCoeff(1,0,lowerValue);
00133 theFunc2->setCoeff(1,0,lowerValue);
00134
00135 double f1 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
00136 theFunc1->setCoeff(1,0,upperValue);
00137 theFunc2->setCoeff(1,0,upperValue);
00138 double f2 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
00139
00140 double xLower = 0.0, xUpper = 0.0;
00141 double root=0;
00142
00143 if ( f1*f2 > 0.0 ) {
00144 report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;
00145 return false;
00146 }
00147
00148
00149 if (fabs(f1) < precision) {
00150 root = lowerValue;
00151 return root;
00152 }
00153 if (fabs(f2) < precision) {
00154 root = upperValue;
00155 return root;
00156 }
00157
00158
00159 if (f1 < 0.0) {
00160 xLower = lowerValue;
00161 xUpper = upperValue;
00162 } else {
00163 xLower = upperValue;
00164 xUpper = lowerValue;
00165 }
00166
00167 double rootGuess = 0.5*(lowerValue + upperValue);
00168 double dxold = fabs(upperValue - lowerValue);
00169 double dx = dxold;
00170
00171 theFunc1->setCoeff(1,0,rootGuess);
00172 theFunc2->setCoeff(1,0,rootGuess);
00173 double f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
00174
00175 for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
00176
00177 dxold = dx;
00178 dx = 0.5*(xUpper-xLower);
00179 rootGuess = xLower+dx;
00180
00181
00182 if (fabs(xLower - rootGuess) < precision) {
00183 root = rootGuess;
00184 return root;
00185 }
00186
00187 theFunc1->setCoeff(1,0,rootGuess);
00188 theFunc2->setCoeff(1,0,rootGuess);
00189 f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
00190
00191 if (f < 0.0) {
00192 xLower = rootGuess;
00193 } else {
00194 xUpper = rootGuess;
00195 }
00196
00197 }
00198
00199 report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
00200 <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!"
00201 <<" Returning false."<<endl;
00202 return 0;
00203
00204 }
00205
00206