Herwig/Utilities/GaussianIntegrator.h

0001 // -*- C++ -*-
0002 //
0003 // GaussianIntegrator.h 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 #ifndef HERWIG_GaussianIntegrator_H
0010 #define HERWIG_GaussianIntegrator_H
0011 //
0012 // This is the declaration of the GaussianIntegrator class.
0013 //
0014
0015 #include "ThePEG/Pointer/ReferenceCounted.h"
0016 #include "ThePEG/Repository/CurrentGenerator.h"
0017 #include <vector>
0018
0019 namespace Herwig {
0020
0021 using namespace ThePEG;
0022
0023 /** \ingroup Utilities
0024  *  \author  Peter Richardson
0025  *  This class is designed to perform the integral of a function
0026  *  using Gaussian quadrature.The method is adaptive based on using 6th,12th,
0027  *  24th,48th, or 96th order Gaussian quadrature combined with
0028  *  subdivision of the integral if this is insufficient.
0029  *
0030  *  The class is templated on a simple class which should provide a
0031  *  T::operator () (double) const which provides the integrand for the function.
0032  */
0033 class GaussianIntegrator : public Pointer::ReferenceCounted {
0034
0035 public:
0036
0037   /** @name Standard constructors and destructors. */
0038   //@{
0039   /**
0040    * Default Constructor
0041    */
0042   GaussianIntegrator()
0043     : _abserr(1.E-35), _relerr(5.E-5), _binwidth(1.E-5), _maxeval(100000) {
0044     // setup the weights and abscissae
0045     Init();
0046   }
0047
0048   /**
0049    * Specify all the parameters.
0050    * @param abserr Absolute error.
0051    * @param relerr Relative error.
0052    * @param binwidth Width of the bin as a fraction of the integration region.
0053    * @param maxeval Maximum number of function evaluations
0054    */
0055   GaussianIntegrator(double abserr, double relerr, double binwidth,
0056                      int maxeval)
0057     : _abserr(abserr), _relerr(relerr),
0058       _binwidth(binwidth),
0059       _maxeval(maxeval) {
0060     // setup the weights and abscissae
0061     Init();
0062   }
0063
0064   /// helper type for the integration result
0065   template <class T>
0066   using ValT = decltype(std::declval<typename T::ValType>()
0067                       * std::declval<typename T::ArgType>());
0068
0069   /**
0070    * The value of the integral
0071    * @param lower The lower limit of integration.
0072    * @param upper The upper limit of integration.
0073    */
0074   template <class T>
0075   inline ValT<T> value(const T &,
0076                    const typename T::ArgType lower,
0077                    const typename T::ArgType upper) const;
0078
0079 private:
0080
0081   /**
0082    * The assignment operator is private and must never be called.
0083    * In fact, it should not even be implemented.
0084    */
0085   GaussianIntegrator & operator=(const GaussianIntegrator &) = delete;
0086
0087   /**
0088    * Initialise the weights and abscissae.
0089    */
0090   void Init();
0091
0092 private:
0093
0094   /**
0095    * The weights for the gaussian quadrature.
0096    */
0097   std::vector< std::vector<double> > _weights;
0098
0099   /**
0100    * The abscissae.
0101    */
0102   std::vector< std::vector <double> > _abscissae;
0103
0104   /**
0105    * The parameters controlling the error.
0106    */
0107   double _abserr,_relerr;
0108
0109   /**
0110    * The minimum width of a bin as a fraction of the integration region.
0111    */
0112   double _binwidth;
0113
0114   /**
0115    * Maximum number of function evaluations.
0116    */
0117   int _maxeval;
0118
0119 };
0120
0121 }
0122
0123 #include "GaussianIntegrator.tcc"
0124
0125 #endif /* HERWIG_GaussianIntegrator_H */