Converting TEvtGen to native cmake

[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtBtoXsgammaRootFinder.cpp

diff --git a/TEvtGen/EvtGen/EvtGenModels/EvtBtoXsgammaRootFinder.cpp b/TEvtGen/EvtGen/EvtGenModels/EvtBtoXsgammaRootFinder.cpp
new file mode 100644 (file)
index 0000000..6048b77
--- /dev/null
+++ b/TEvtGen/EvtGen/EvtGenModels/EvtBtoXsgammaRootFinder.cpp
@@ -0,0 +1,206 @@
+//--------------------------------------------------------------------------
+//
+//
+// Copyright Information: See EvtGen/COPYRIGHT
+//
+// Environment:
+//      This software is part of the EvtGen package developed jointly
+//      for the BaBar and CLEO collaborations.  If you use all or part
+//      of it, please give an appropriate acknowledgement.
+//
+// Module: EvtBtoXsgammaRootFinder.cc
+//
+// Description:
+//      Root finders for EvtBtoXsgammaKagan module.
+//
+// Modification history:
+//
+//      Jane Tinslay       March 21, 2001       Module created
+//------------------------------------------------------------------------
+#include "EvtGenBase/EvtPatches.hh"
+
+#include "EvtGenModels/EvtBtoXsgammaRootFinder.hh"
+#include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
+#include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
+#include "EvtGenBase/EvtReport.hh"
+#include <math.h>
+using std::endl;
+
+//-------------
+// C Headers --
+//-------------
+extern "C" {
+}
+
+//-----------------------------------------------------------------------
+// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
+//-----------------------------------------------------------------------
+
+#define EVTITGROOTFINDER_MAXIT 100
+#define EVTITGROOTFINDER_RELATIVEPRECISION 1.0e-16
+
+
+EvtBtoXsgammaRootFinder::EvtBtoXsgammaRootFinder() {}
+
+EvtBtoXsgammaRootFinder::~EvtBtoXsgammaRootFinder( )
+{}
+
+double
+EvtBtoXsgammaRootFinder::GetRootSingleFunc(const EvtItgAbsFunction* theFunc, double functionValue, double lowerValue, double upperValue, double precision) {
+  
+  // Use the bisection to find the root.
+  // Iterates until find root to the accuracy of precision
+
+  double xLower = 0.0, xUpper = 0.0;
+  double root=0;
+
+  double f1 = theFunc->value(lowerValue) - functionValue;
+  double f2 = theFunc->value(upperValue) - functionValue;
+
+  if ( f1*f2 > 0.0 ) {
+    report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;  
+    return 0;
+  }
+
+  // Already have root
+  if (fabs(f1) < precision) {
+    root = lowerValue;
+    return root;
+  }
+  if (fabs(f2) < precision) {
+    root = upperValue;
+    return root;
+  }
+  
+  // Orient search so that f(xLower) < 0
+  if (f1 < 0.0) {
+    xLower = lowerValue;
+    xUpper = upperValue;
+  } else {
+    xLower = upperValue;
+    xUpper = lowerValue;
+  }
+  
+  double rootGuess = 0.5*(lowerValue + upperValue);
+  double dxold = fabs(upperValue - lowerValue);
+  double dx = dxold;
+  
+  double f = theFunc->value(rootGuess) - functionValue;
+  
+  for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
+    
+      dxold = dx;
+      dx = 0.5*(xUpper-xLower);
+      rootGuess = xLower+dx;
+
+      // If change in root is negligible, take it as solution.
+      if (fabs(xLower - rootGuess) < precision) {
+       root = rootGuess;
+       return root;
+      }
+      
+      f = theFunc->value(rootGuess) - functionValue;
+ 
+      if (f < 0.0) {
+       xLower = rootGuess;
+      } else {
+       xUpper = rootGuess;
+      }
+      
+  }
+  
+  report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
+                          <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!" 
+                          <<" Returning false."<<endl;
+  return 0;
+  
+}
+
+double
+EvtBtoXsgammaRootFinder::GetGaussIntegFcnRoot(EvtItgAbsFunction *theFunc1, EvtItgAbsFunction *theFunc2, double integ1Precision, double integ2Precision, int maxLoop1, int maxLoop2, double integLower, double integUpper, double lowerValue, double upperValue, double precision) {
+ 
+  // Use the bisection to find the root.
+  // Iterates until find root to the accuracy of precision
+  
+  //Need to work with integrators
+  EvtItgAbsIntegrator *func1Integ = new EvtItgSimpsonIntegrator(*theFunc1, integ1Precision, maxLoop1);
+  EvtItgAbsIntegrator *func2Integ = new EvtItgSimpsonIntegrator(*theFunc2, integ2Precision, maxLoop2);
+  
+  
+  //coefficient 1 of the integrators is the root to be found
+  //need to set this to lower value to start off with
+  theFunc1->setCoeff(1,0,lowerValue);
+  theFunc2->setCoeff(1,0,lowerValue);
+  
+  double f1 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
+  theFunc1->setCoeff(1,0,upperValue);
+  theFunc2->setCoeff(1,0,upperValue);
+  double f2 = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
+  
+  double xLower = 0.0, xUpper = 0.0;
+  double root=0;
+
+  if ( f1*f2 > 0.0 ) {
+    report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: No root in specified range !"<<endl;  
+    return false;
+  }
+
+  // Already have root
+  if (fabs(f1) < precision) {
+    root = lowerValue;
+    return root;
+  }
+  if (fabs(f2) < precision) {
+    root = upperValue;
+    return root;
+  }
+  
+  // Orient search so that f(xLower) < 0
+  if (f1 < 0.0) {
+    xLower = lowerValue;
+    xUpper = upperValue;
+  } else {
+    xLower = upperValue;
+    xUpper = lowerValue;
+  }
+  
+  double rootGuess = 0.5*(lowerValue + upperValue);
+  double dxold = fabs(upperValue - lowerValue);
+  double dx = dxold;
+  
+  theFunc1->setCoeff(1,0,rootGuess);
+  theFunc2->setCoeff(1,0,rootGuess);
+  double f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
+  
+  for (int j = 0; j< EVTITGROOTFINDER_MAXIT; j++) {
+    
+    dxold = dx;
+    dx = 0.5*(xUpper-xLower);
+    rootGuess = xLower+dx;
+    
+    // If change in root is negligible, take it as solution.
+    if (fabs(xLower - rootGuess) < precision) {
+      root = rootGuess;
+      return root;
+    }
+    
+    theFunc1->setCoeff(1,0,rootGuess);
+    theFunc2->setCoeff(1,0,rootGuess);
+    f = func1Integ->evaluate(integLower,integUpper) - theFunc2->getCoeff(1,2)*func2Integ->evaluate(integLower,integUpper);
+    
+    if (f < 0.0) {
+      xLower = rootGuess;
+    } else {
+      xUpper = rootGuess;
+    }
+    
+  }
+  
+  report(WARNING,"EvtGen") << "EvtBtoXsgammaRootFinder: Maximum number of iterations "
+                          <<"in EvtBtoXsgammaRootFinder::foundRoot exceeded!" 
+                          <<" Returning false."<<endl;
+  return 0;
+  
+}
+
+