Herwig/Utilities/GSLIntegrator.tcc

0001 // -*- C++ -*-
0002 //
0003 // GSLIntegrator.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 //
0010 // This is the implementation of the non-inlined templated member
0011 // functions of the GSLIntegrator class.
0012 //
0013 using namespace Herwig;
0014 using namespace ThePEG;
0015
0016 namespace {
0017
0018   template <class T> struct param {
0019
0020     //The integrand function
0021     const T & function;
0022
0023   };
0024
0025   template<class T> double integrand(double x , void * p) {
0026     //Units of the argument and return type
0027     const typename T::ValType ValUnit =
0028       TypeTraits<typename T::ValType>::baseunit();
0029     const typename T::ArgType ArgUnit =
0030       TypeTraits<typename T::ArgType>::baseunit();
0031
0032     const T & f = ((struct param<T> *)p)->function;
0033     return f(x * ArgUnit ) / ValUnit;
0034   }
0035
0036 }
0037
0038 namespace Herwig {
0039 using namespace ThePEG;
0040
0041
0042 template <class T>
0043 inline GSLIntegrator::ValT<T>
0044 GSLIntegrator::value(const T & fn,
0045                      const typename T::ArgType lower,
0046                      const typename T::ArgType upper) const
0047 {
0048   GSLIntegrator::ValT<T> error;
0049   return value(fn,lower,upper,error);
0050 }
0051
0052
0053 template <class T>
0054 inline GSLIntegrator::ValT<T>
0055 GSLIntegrator::value(const T & fn,
0056                      const typename T::ArgType lower,
0057                      const typename T::ArgType upper,
0058                      GSLIntegrator::ValT<T> & error) const
0059 {
0060   typedef typename T::ValType ValType;
0061   typedef typename T::ArgType ArgType;
0062   const ValType ValUnit = TypeTraits<ValType>::baseunit();
0063   const ArgType ArgUnit = TypeTraits<ArgType>::baseunit();
0064
0065   double result(0.), error2(0.);
0066
0067   param<T> parameters = { fn };
0068   gsl_function integrationFunction;
0069   integrationFunction.function = &integrand<T>;
0070   integrationFunction.params = &parameters;
0071
0072   gsl_integration_workspace * workspace =
0073     gsl_integration_workspace_alloc(_nbins);
0074   //do integration
0075   //Want to check error messages ourselves
0076   gsl_error_handler_t * oldhandler = gsl_set_error_handler_off();
0077   int status = gsl_integration_qags(&integrationFunction, lower/ArgUnit,
0078                                     upper/ArgUnit, _abserr, _relerr, _nbins,
0079                                     workspace, &result, &error2);
0080   if( status > 0 ) {
0081     CurrentGenerator::log() << "An error occurred in the GSL "
0082       "integration subroutine:\n";
0083     switch( status ) {
0084     case GSL_EMAXITER:
0085       CurrentGenerator::log() << "The maximum number of subdivisions "
0086         "was exceeded.\n";
0087       break;
0088     case GSL_EROUND:
0089       CurrentGenerator::log() << "Cannot reach tolerance because of "
0090         "roundoff error, or roundoff error was detected in the "
0091         "extrapolation table.\n";
0092       break;
0093     case GSL_ESING:
0094       CurrentGenerator::log() << "A non-integrable singularity or "
0095         "other bad integrand behavior was found in the integration "
0096         "interval.\n";
0097       break;
0098     case GSL_EDIVERGE:
0099       CurrentGenerator::log() << "The integral is divergent, "
0100         "or too slowly convergent to be integrated numerically.\n";
0101       break;
0102     default:
0103       CurrentGenerator::log() << "A general error occurred with code "
0104                               << status << '\n';
0105     }
0106     result = 0.;
0107   }
0108   gsl_set_error_handler(oldhandler);
0109   gsl_integration_workspace_free(workspace);
0110
0111   //fix units and return
0112   error = error2* ValUnit * ArgUnit;
0113   return result * ValUnit * ArgUnit;
0114 }
0115
0116 }