Herwig/Utilities/GSLBisection.h

0001 // -*- C++ -*-
0002 //
0003 // GSLBisection.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_GSLBisection_H
0010 #define HERWIG_GSLBisection_H
0011 //
0012 // This is the declaration of the GSLBisection class.
0013 //
0014
0015 #include "ThePEG/Pointer/ReferenceCounted.h"
0016 #include "Herwig/Utilities/GSLHelper.h"
0017 #include <gsl/gsl_errno.h>
0018 #include <gsl/gsl_math.h>
0019 #include <gsl/gsl_roots.h>
0020
0021 namespace Herwig {
0022
0023 using namespace ThePEG;
0024
0025 /** \ingroup Utilities
0026  * This class is designed to find the root of a given function between
0027  * 2 limits using bisection methods.
0028  *
0029  * \author Manuel B\"ahr
0030  *
0031  * The function is supplied using a templated class that must define
0032  * operator(argument). The units of the argument ArgType and return type
0033  * ValType must be supplied in the integrand class using a typedef. In
0034  * addition the baseunit should be supplied by static methods vUnit()
0035  * and aUnit() to avoid numerical problems that arise when the centrally
0036  * defined baseunit is several orders of magnitude off the one you
0037  * need. As an example see: <br>
0038  * <code> struct integrand { </code><br>
0039  * <code> ... </code> <BR>
0040  * <code>Energy operator(CrossSection arg) const;</code><BR>
0041  * <code>typedef CrossSection ArgType</code><BR>
0042  * <code>typedef Energy ValType</code><BR>
0043  * <code>static ArgType aUnit(){return 1.*millibarn;} </code> <BR>
0044  * <code>static ValType vUnit(){return 1.*MeV;} </code> <BR>
0045  * <code> ... </code> <BR>
0046  * <code>}</code> <BR>
0047  * This can be facilitated by deriving from the GSLHelper struct. Which
0048  * implents the vUnit() and aUnit() methods using the baseunit static
0049  * method. Also the typedefs are written there.
0050  */
0051 class GSLBisection : public Pointer::ReferenceCounted {
0052
0053 public:
0054
0055   /**
0056    * Struct that is used to throw and catch GSL errors
0057    */
0058   struct GSLerror {};
0059
0060   /**
0061    * Struct that is used to throw and catch GSL errors
0062    */
0063   struct IntervalError {};
0064
0065   /** @name Standard constructors and destructors. */
0066   //@{
0067   /**
0068    * Default Constructor
0069    */
0070   GSLBisection() : abserr_(0), relerr_(1.E-8), maxPoints_(100) {}
0071
0072   /**
0073    * Specify all the parameters.
0074    * @param abserr Absolute error.
0075    * @param relerr Relative error.
0076    * @param max Maximum number of intervals
0077    */
0078   inline GSLBisection(double abserr, double relerr, int max) :
0079     abserr_(abserr), relerr_(relerr), maxPoints_(max) {}
0080
0081   //@}
0082
0083   /**
0084    * Function to overwrite the default GSL error handling
0085    */
0086   static void GSLsubstHandler(const char *, const char *,
0087                   int, int){
0088     throw GSLerror();
0089   }
0090
0091   /**
0092    * The result of the root finding.
0093    * @param function The integrand class that defines operator()
0094    * @param lower The lower limit of integration.
0095    * @param upper The upper limit of integration.
0096    */
0097   template <class T>
0098   inline typename T::ArgType value(const T & function,
0099                    const typename T::ArgType lower,
0100                    const typename T::ArgType upper) const;
0101
0102 private:
0103
0104   /**
0105    * The assignment operator is private and must never be called.
0106    * In fact, it should not even be implemented.
0107    */
0108   GSLBisection & operator=(const GSLBisection &) = delete;
0109
0110 private:
0111
0112   /**
0113    * The parameters controlling the absolute error.
0114    */
0115   double abserr_;
0116
0117   /**
0118    * The parameters controlling the relatve error.
0119    */
0120   double relerr_;
0121
0122   /**
0123    * The maximum number of evaluations to use.
0124    */
0125   int maxPoints_;
0126 };
0127
0128 }
0129
0130 #include "GSLBisection.tcc"
0131
0132 #endif /* HERWIG_GSLBisection_H */