TEvtGen/EvtGenModels/EvtBtoXsgammaRootFinder.cxx

   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