diff --git a/include/LHAPDF/AlphaS.h b/include/LHAPDF/AlphaS.h
--- a/include/LHAPDF/AlphaS.h
+++ b/include/LHAPDF/AlphaS.h
@@ -1,336 +1,337 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_AlphaS_H
 #define LHAPDF_AlphaS_H
 
 #include "LHAPDF/Utils.h"
 #include "LHAPDF/Exceptions.h"
 #include "LHAPDF/KnotArray.h"
 
 namespace LHAPDF {
 
+
   /// @brief Calculator interface for computing alpha_s(Q2) in various ways
   ///
   /// The design of the AlphaS classes is that they are substitutible
   /// (cf. polymorphism) and are entirely non-dependent on the PDF and Info
   /// objects: hence they can be used by external code that actually doesn't
   /// want to do anything at all with PDFs, but which just wants to do some
   /// alpha_s interpolation.
   class AlphaS {
   public:
 
     /// Base class constructor for default param setup
     AlphaS();
 
     /// Destructor
     virtual ~AlphaS() {};
 
     /// @name alpha_s values
-    //@{
+    ///@{
 
     /// Calculate alphaS(Q)
     double alphasQ(double q) const { return alphasQ2(q*q); }
 
     /// Calculate alphaS(Q2)
     /// @todo Throw error in this base method if Q < Lambda?
     virtual double alphasQ2(double q2) const = 0;
 
-    //@}
+    ///@}
 
 
     /// @name alpha_s metadata
-    //@{
+    ///@{
 
     /// Calculate the number of active flavours at energy scale Q
     int numFlavorsQ(double q) const { return numFlavorsQ2(q*q); }
 
     /// Calculate the number of active flavours at energy scale Q2
     virtual int numFlavorsQ2(double q2) const;
 
     /// Get a quark mass by PDG code
     double quarkMass(int id) const;
 
     /// @brief Set quark masses by PDG code
     ///
     /// Used in the analytic and ODE solvers.
     void setQuarkMass(int id, double value);
 
     /// @brief Get a flavor scale threshold by PDG code
     ///
     /// Used in the analytic and ODE solvers.
     double quarkThreshold(int id) const;
 
     /// @brief Set a flavor threshold by PDG code (= quark masses by default)
     ///
     /// Used in the analytic and ODE solvers.
     void setQuarkThreshold(int id, double value);
 
     /// Get the order of QCD (expressed as number of loops)
     ///
     /// Used explicitly in the analytic and ODE solvers.
     int orderQCD() { return _qcdorder; }
 
     /// @brief Set the order of QCD (expressed as number of loops)
     ///
     /// Used in the analytic and ODE solvers.
     void setOrderQCD(int order) { _qcdorder = order; }
 
     /// @brief Set the Z mass used in this alpha_s
     ///
     /// Used in the ODE solver.
     void setMZ(double mz) { _mz = mz; }
 
     /// @brief Set the alpha_s(MZ) used in this alpha_s
     ///
     /// Used in the ODE solver.
     void setAlphaSMZ(double alphas) { _alphas_mz = alphas; }
 
     /// @brief Set the Z mass used in this alpha_s
     ///
     /// Used in the ODE solver.
     void setMassReference(double mref) { _mreference = mref; _customref = true; }
 
     /// @brief Set the alpha_s(MZ) used in this alpha_s
     ///
     /// Used in the ODE solver.
     void setAlphaSReference(double alphas) { _alphas_reference = alphas; _customref = true; }
 
     /// @brief Set the @a {i}th Lambda constant for @a i active flavors
     ///
     /// Used in the analytic solver.
     virtual void setLambda(unsigned int, double) {};
 
-    //@}
-
-    /// enum of flavor schemes
+    /// Enum of flavor schemes
     enum FlavorScheme { FIXED, VARIABLE };
 
     /// Get the implementation type of this AlphaS
     virtual std::string type() const = 0;
 
     /// Set flavor scheme of alpha_s solver
     void setFlavorScheme(FlavorScheme scheme, int nf = -1);
 
     /// Get flavor scheme
     FlavorScheme flavorScheme() const { return _flavorscheme; }
 
+    ///@}
+
 
   protected:
 
     /// @name Calculating beta function values
-    //@{
+    ///@{
 
     /// Calculate the i'th beta function given the number of active flavours
     /// Currently limited to 0 <= i <= 3
     /// Calculated using the MSbar scheme
     double _beta(int i, int nf) const;
 
     /// Calculate a vector of beta functions given the number of active flavours
     /// Currently returns a 4-element vector of beta0 -- beta3
     std::vector<double> _betas(int nf) const;
 
-    //@}
+    ///@}
 
 
   protected:
 
     /// Order of QCD evolution (expressed as number of loops)
     int _qcdorder;
 
     /// Mass of the Z boson in GeV
     double _mz;
 
     /// Value of alpha_s(MZ)
     double _alphas_mz;
 
     /// Reference mass in GeV
     double _mreference;
 
     /// Value of alpha_s(reference mass)
     double _alphas_reference;
 
     /// Decides whether to use custom reference values or fall back on MZ/AlphaS_MZ
     bool _customref;
 
     /// Masses of quarks in GeV
     /// Used for working out flavor thresholds and the number of quarks that are
     /// active at energy scale Q.
     std::map<int, double> _quarkmasses, _flavorthresholds;
 
     /// The flavor scheme in use
     FlavorScheme _flavorscheme;
 
     /// The allowed numbers of flavours in a fixed scheme
     int _fixflav;
 
   };
 
 
 
   /// Calculate alpha_s(Q2) by an analytic approximation
   class AlphaS_Analytic : public AlphaS {
   public:
 
     /// Implementation type of this solver
     std::string type() const { return "analytic"; }
 
     /// Calculate alphaS(Q2)
     double alphasQ2(double q2) const;
 
     /// Analytic has its own numFlavorsQ2 which respects the min/max nf set by the Lambdas
     int numFlavorsQ2(double q2) const;
 
     /// Set lambda_i (for i = flavour number)
     void setLambda(unsigned int i, double lambda);
 
 
   private:
 
     /// Get lambdaQCD for nf
     double _lambdaQCD(int nf) const;
 
     /// Recalculate min/max flavors in case lambdas have changed
     void _setFlavors();
 
 
     /// LambdaQCD values.
     std::map<int, double> _lambdas;
 
     /// Max number of flavors
     int _nfmax;
     /// Min number of flavors
     int _nfmin;
 
   };
 
 
 
   /// Interpolate alpha_s from tabulated points in Q2 via metadata
   ///
   /// @todo Extrapolation: log-gradient xpol at low Q, const at high Q?
   class AlphaS_Ipol : public AlphaS {
   public:
 
     /// Implementation type of this solver
     std::string type() const { return "ipol"; }
 
     /// Calculate alphaS(Q2)
     double alphasQ2(double q2) const;
 
     /// Set the array of Q values for interpolation
     ///
     /// Writes to the same internal arrays as setQ2Values, appropriately transformed.
     void setQValues(const std::vector<double>& qs);
 
     /// Set the array of Q2 values for interpolation
     ///
     /// Subgrids are represented by repeating the values which are the end of
     /// one subgrid and the start of the next. The supplied vector must match
     /// the layout of alpha_s values.
     void setQ2Values(const std::vector<double>& q2s) { _q2s = q2s; }
 
     /// Set the array of alpha_s(Q2) values for interpolation
     ///
     /// The supplied vector must match the layout of Q2 knots.  Subgrids may
     /// have discontinuities, i.e. different alpha_s values on either side of a
     /// subgrid boundary (for the same Q values).
     void setAlphaSValues(const std::vector<double>& as) { _as = as; }
 
 
   private:
 
     /// Standard cubic interpolation formula
     double _interpolateCubic(double T, double VL, double VDL, double VH, double VDH) const;
     /// Get the gradient for a patch in the middle of the grid
     double _ddq_central( size_t i ) const;
     /// Get the gradient for a patch at the low end of the grid
     double _ddq_forward( size_t i ) const;
     /// Get the gradient for a patch at the high end of the grid
     double _ddq_backward( size_t i ) const;
 
     /// Synchronise the contents of the single Q2 / alpha_s vectors into subgrid objects
     /// @note This is const so it can be called silently from a const method
     void _setup_grids() const;
 
 
     /// Map of AlphaSArrays "binned" for lookup by low edge in (log)Q2
     /// @note This is mutable so it can be initialized silently from a const method
     mutable std::map<double, AlphaSArray> _knotarrays;
 
     /// Array of ipol knots in Q2
     std::vector<double> _q2s;
     /// Array of alpha_s values for the Q2 knots
     std::vector<double> _as;
 
   };
 
 
 
   /// Solve the differential equation in alphaS using an implementation of RK4
   class AlphaS_ODE : public AlphaS {
   public:
 
     /// Implementation type of this solver
     std::string type() const { return "ode"; }
 
     /// Calculate alphaS(Q2)
     double alphasQ2( double q2 ) const;
 
     /// Set MZ, and also the caching flag
     void setMZ( double mz ) { _mz = mz; _calculated = false; }
 
     /// Set alpha_s(MZ), and also the caching flag
     void setAlphaSMZ( double alphas ) { _alphas_mz = alphas; _calculated = false; }
 
     /// Set reference mass, and also the caching flag
     void setMassReference( double mref ) { _mreference = mref; _calculated = false; _customref = true; }
 
     /// Set alpha_s(MZ), and also the caching flag
     void setAlphaSReference( double alphas ) { _alphas_reference = alphas; _calculated = false; _customref = true; }
 
     /// Set the array of Q values for interpolation, and also the caching flag
     void setQValues(const std::vector<double>& qs);
 
     /// @brief Set the array of Q2 values for interpolation, and also the caching flag
     ///
     /// Writes to the same internal array as setQValues, appropriately transformed.
     void setQ2Values( std::vector<double> q2s ) { _q2s = q2s; _calculated = false; }
 
 
   private:
 
     /// Calculate the derivative at Q2 = t, alpha_S = y
     double _derivative(double t, double y, const std::vector<double>& beta) const;
 
     /// Calculate the decoupling relation when going from num. flav. = ni -> nf
     /// abs(ni - nf) must be = 1
     double _decouple(double y, double t, unsigned int ni, unsigned int nf) const;
 
     /// Calculate the next step using RK4 with adaptive step size
     void _rk4(double& t, double& y, double h, const double allowed_change, const vector<double>& bs) const;
 
     /// Solve alpha_s for q2 using RK4
     void _solve(double q2, double& t, double& y, const double& allowed_relative, double h, double accuracy) const;
 
     /// Create interpolation grid
     void _interpolate() const;
 
 
     /// Vector of Q2s in case specific anchor points are used
     mutable std::vector<double> _q2s;
 
     /// Whether or not the ODE has been solved yet
     mutable bool _calculated;
 
     /// The interpolation used to get Alpha_s after the ODE has been solved
     mutable AlphaS_Ipol _ipol;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/Config.h b/include/LHAPDF/Config.h
--- a/include/LHAPDF/Config.h
+++ b/include/LHAPDF/Config.h
@@ -1,71 +1,69 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Config_H
 #define LHAPDF_Config_H
 
 #include "LHAPDF/Info.h"
 
 namespace LHAPDF {
 
 
   /// Class for PDF set metadata and manipulation
   class Config : public Info {
   public:
 
     /// @name Fetching/creation
-    //@{
+    ///@{
 
     /// Get the global configuration object
     ///
     /// The global config is populated by reading from lhapdf.conf if it is
     /// found in the search paths. It is a singleton, hence the 'get' accessor
     /// rather than a constructor.
     ///
     /// @note The LHAPDF system is responsible for deletion of the returned
     /// object. Do NOT delete it yourself!
     static Config& get();
-    //@}
+    ///@}
 
 
     /// Config destructor, used for end-of-run banner printing
     ~Config();
 
 
   private:
 
     /// Hide the default constructor
     Config() {
       // std::cout << "CONFIG CONSTRUCTION" << std::endl;
     }
 
-    //@}
-
   };
 
 
-  /// @name Convenient verbosity control
-  //@{
+  /// @defgroup verb Verbosity control
+  ///@{
 
   /// Convenient way to get the current verbosity level
   ///
   /// @note Verbosity, like any other flag, can also be set at lower levels. But who does that, really?!?
   inline int verbosity() {
     return Config::get().get_entry_as<int>("Verbosity", 1);
   }
 
   /// Convenient way to set the verbosity level
   ///
   /// @note Verbosity, like any other flag, can also be set at lower levels. But who does that, really?!?
   inline void setVerbosity(int v) {
     Config::get().set_entry("Verbosity", v);
   }
 
-  //@}
+  ///@}
 
 
 }
 #endif
diff --git a/include/LHAPDF/Exceptions.h b/include/LHAPDF/Exceptions.h
--- a/include/LHAPDF/Exceptions.h
+++ b/include/LHAPDF/Exceptions.h
@@ -1,128 +1,128 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Exceptions_H
 #define LHAPDF_Exceptions_H
 
 #include <exception>
 #include <stdexcept>
 
 namespace LHAPDF {
 
 
-  /// @name Exception classes for error handling
+  /// @defgroup exceptions Exception classes for error handling
   //@{
 
   /// @brief Generic unspecialised LHAPDF runtime error.
   ///
   /// NB. We don't use "Error" because that has a physics meaning!
   class Exception : public std::runtime_error {
   public:
     /// Constructor with error description string
     Exception(const std::string& what) : std::runtime_error(what) {}
   };
 
 
   /// Error for general PDF grid problems.
   class GridError : public Exception {
   public:
     /// Constructor with error description string
     GridError(const std::string& what) : Exception(what) {}
   };
 
 
   /// Error to be thrown when out of the valid range of a PDF.
   class RangeError : public Exception {
   public:
     /// Constructor with error description string
     RangeError(const std::string& what) : Exception(what) {}
   };
 
 
   /// Error for places where it should not have been possible to get to!
   class LogicError : public Exception {
   public:
     /// Constructor with error description string
     LogicError(const std::string& what) : Exception(what) {}
   };
 
 
   /// @brief Error for unfound or broken metadata entries.
   class MetadataError : public Exception {
   public:
     /// Constructor with error description string
     MetadataError(const std::string& what) : Exception(what) {}
   };
 
 
   /// @brief Error for file reading errors.
   class ReadError : public Exception {
   public:
     /// Constructor with error description string
     ReadError(const std::string& what) : Exception(what) {}
   };
 
 
   /// @brief Error for requests for unsupported/invalid flavour PIDs.
   class FlavorError : public Exception {
   public:
     /// Constructor with error description string
     FlavorError(const std::string& what) : Exception(what) {}
   };
 
 
   /// @brief Error to be raised by object factories given invalid requests.
   class FactoryError : public Exception {
   public:
     /// Constructor with error description string
     FactoryError(const std::string& what) : Exception(what) {}
   };
 
 
   /// @brief Error to be raised when a LHAPDF ID indexing fails
   class IndexError : public Exception {
   public:
     /// Constructor with error description string
     IndexError(const std::string& what) : Exception(what) {}
   };
 
 
   /// Error for general AlphaS computation problems.
   class AlphaSError : public Exception {
   public:
     /// Constructor with error description string
     AlphaSError(const std::string& what) : Exception(what) {}
   };
 
 
   /// @brief Error to be raised when a newer LHAPDF version is needed
   class VersionError : public Exception {
   public:
     /// Constructor with error description string
     VersionError(const std::string& what) : Exception(what) {}
   };
 
 
-  /// Problem exists between keyboard and chair.
+  /// Problem exists between keyboard and chair
   class UserError : public Exception {
   public:
     /// Constructor with error description string
     UserError(const std::string& what) : Exception(what) {}
   };
 
 
   /// This feature doesn't exist yet
   class NotImplementedError : public Exception {
   public:
     /// Constructor with error description string
     NotImplementedError(const std::string& what) : Exception(what) {}
   };
 
-  //@}
+  ///@}
 
 
 }
 #endif
diff --git a/include/LHAPDF/Extrapolator.h b/include/LHAPDF/Extrapolator.h
--- a/include/LHAPDF/Extrapolator.h
+++ b/include/LHAPDF/Extrapolator.h
@@ -1,80 +1,80 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Extrapolator_H
 #define LHAPDF_Extrapolator_H
 
 #include "LHAPDF/Utils.h"
 
 namespace LHAPDF {
 
 
   // Forward declaration
   class GridPDF;
 
 
   /// The general interface for extrapolating beyond grid boundaries
   class Extrapolator {
   public:
 
     /// Destructor to allow inheritance
     virtual ~Extrapolator() { }
 
 
     /// @name Binding to a PDF object
-    //@{
+    ///@{
 
     /// Bind to a GridPDF
     void bind(const GridPDF* pdf) { _pdf = pdf; }
 
     /// Unbind from GridPDF
     void unbind() { _pdf = 0; }
 
     /// Identify whether this Extrapolator has an associated PDF
     bool hasPDF() { return _pdf != 0; }
 
     /// Get the associated GridPDF
     const GridPDF& pdf() const { return *_pdf; }
 
-    //@}
+    ///@}
 
 
     /// @name Extrapolation methods
-    //@{
+    ///@{
 
     /// Extrapolate a single-point in (x,Q)
     ///
     /// @param id PDG parton ID
     /// @param x Momentum fraction
     /// @param q Energy scale
     /// @return The xf value at (x,q2)
     double extrapolateXQ(int id, double x, double q) const {
       return extrapolateXQ2(id, x, q*q );
     }
 
     /// Extrapolate a single-point in (x,Q2)
     ///
     /// @param id PDG parton ID
     /// @param x Momentum fraction
     /// @param q2 Squared energy scale
     /// @return The xf value at (x,q2)
     virtual double extrapolateXQ2(int id, double x, double q2) const = 0;
 
 
     /// @todo Make an all-PID version of extrapolateQ and Q2?
 
-    //@}
+    ///@}
 
 
   private:
 
     const GridPDF* _pdf;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/GridPDF.h b/include/LHAPDF/GridPDF.h
--- a/include/LHAPDF/GridPDF.h
+++ b/include/LHAPDF/GridPDF.h
@@ -1,242 +1,237 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_GridPDF_H
 #define LHAPDF_GridPDF_H
 
 #include "LHAPDF/PDF.h"
 #include "LHAPDF/Interpolator.h"
 #include "LHAPDF/Extrapolator.h"
 #include "LHAPDF/KnotArray.h"
 
 namespace LHAPDF {
 
 
   /// @brief A PDF defined via an interpolation grid
   class GridPDF : public PDF {
   public:
 
-    /// @name Creation and deletion
-    //@{
-
     /// Default constructor, making an empty PDF to be populated by hand.
     GridPDF() {
       _mempath = "";
       _info = PDFInfo();
       _forcePos = -1;
     }
 
     /// @brief Constructor from a file path
     ///
     /// We allow this to exist and be user-callable for testing and other
     /// special case uses, since if you are explicitly instantiating a GridPDF
     /// rather than acquiring it via a pointer/reference of PDF type, then you
     /// probably (hopefully) know what you're doing and aren't putting it into
     /// public production code!
     GridPDF(const std::string& path) {
       _loadInfo(path); // Sets _mempath
       _loadPlugins();
       _loadData(_mempath);
       _forcePos = -1;
     }
 
     /// Constructor from a set name and member ID
     GridPDF(const std::string& setname, int member) {
       _loadInfo(setname, member); // Sets _mempath
       _loadPlugins();
       _loadData(_mempath);
       _forcePos = -1;
     }
 
     /// Constructor from an LHAPDF ID
     GridPDF(int lhaid) {
       _loadInfo(lhaid); // Sets _mempath
       _loadPlugins();
       _loadData(_mempath);
       _forcePos = -1;
     }
 
     /// Virtual destructor to allow inheritance
     virtual ~GridPDF() { }
 
-    //@}
-
 
   protected:
 
     /// Load the interpolator, based on current metadata
     void _loadInterpolator();
 
     /// Load the PDF grid data block, based on current metadata
     void _loadExtrapolator();
 
     /// Load the alphaS, interpolator, and extrapolator based on current metadata
     void _loadPlugins() {
       _loadAlphaS();
       _loadInterpolator();
       _loadExtrapolator();
     }
 
     /// Load the PDF grid data block (not the metadata) from the given PDF member file
     void _loadData(const std::string& mempath);
 
 
   public:
 
     /// @name Interpolators and extrapolators
-    //@{
+    ///@{
 
     /// @brief Set the interpolator by pointer
     ///
     /// The provided Interpolator must have been new'd, as it will not be copied
     /// and ownership passes to this GridPDF: delete will be called on this ptr
     /// when this GridPDF goes out of scope or another setInterpolator call is made.
     void setInterpolator(Interpolator* ipol);
 
     /// @brief Set the interpolator by value
     ///
     /// The passed value must be a concrete instantiation of the Interpolator
     /// interface. It will be copied and heap-assigned for use inside this GridPDF.
     ///
     /// @todo Use SFINAE magic to restrict INTERPOLATOR to subclasses of Interpolator?
     template <typename INTERPOLATOR>
     void setInterpolator(INTERPOLATOR ipol) {
       setInterpolator(new INTERPOLATOR(ipol));
     }
 
     /// @brief Set the interpolator by name
     ///
     /// Use the interpolator specified by the given name, as passed to the
     /// createInterpolator factory function.
     void setInterpolator(const std::string& ipolname);
 
     /// Find whether an extrapolator has been set on this PDF
     bool hasInterpolator() const { return bool(_interpolator); }
 
     /// Get the current interpolator
     const Interpolator& interpolator() const;
 
 
     /// @brief Set the extrapolator by pointer
     ///
     /// The provided Extrapolator must have been new'd, as it will not be copied
     /// and ownership passes to this GridPDF: delete will be called on this ptr
     /// when this GridPDF goes out of scope or another setExtrapolator call is made.
     void setExtrapolator(Extrapolator* xpol);
 
     /// @brief Set the extrapolator by value
     ///
     /// The passed value must be a concrete instantiation of the Extrapolator
     /// interface. It will be copied and heap-assigned for use inside this GridPDF.
     ///
     /// @todo Use SFINAE magic to restrict EXTRAPOLATOR to subclasses of Extrapolator?
     template <typename EXTRAPOLATOR>
     void setExtrapolator(EXTRAPOLATOR xpol) {
       setExtrapolator(new EXTRAPOLATOR(xpol));
     }
 
     /// @brief Set the extrapolator by name
     ///
     /// Use the extrapolator specified by the given name, as passed to the
     /// createExtrapolator factory function.
     void setExtrapolator(const std::string& xpolname);
 
     /// Find whether an extrapolator has been set on this PDF
     bool hasExtrapolator() const { return bool(_extrapolator); }
 
     /// Get the current extrapolator
     const Extrapolator& extrapolator() const;
 
-    //@}
+    ///@}
 
 
   protected:
 
     /// @brief Get PDF xf(x,Q2) value (via grid inter/extrapolators)
     double _xfxQ2(int id, double x, double q2) const;
 
 
   public:
 
     /// @name Info about the grid, and access to the raw data points
-    //@{
+    ///@{
 
     /// Directly access the knot arrays in non-const mode, for programmatic filling
     std::map<double, KnotArrayNF>& knotarrays() {
       return _knotarrays;
     }
 
     /// Get the N-flavour subgrid containing Q2 = q2
     const KnotArrayNF& subgrid(double q2) const;
 
     /// Get the 1-flavour subgrid for PID=id containing Q2 = q2
     const KnotArray1F& subgrid(int id, double q2) const {
       return subgrid(q2).get_pid(id);
     }
 
     /// @brief Return a representative list of interpolation knots in x
     ///
     /// The x knot array for the first flavor grid of the lowest-Q2 subgrid is returned.
     const vector<double>& xKnots() const {
       const KnotArrayNF& subgrid1 = _knotarrays.begin()->second;
       const KnotArray1F& grid1 = subgrid1.get_first();
       return grid1.xs();
     }
 
     /// @brief Return a representative list of interpolation knots in Q2
     ///
     /// Constructed and cached by walking over all subgrids and concatenating their Q2 lists: expensive!
     const vector<double>& q2Knots() const;
 
   public:
 
     /// Check if x is in the grid range
     bool inRangeX(double x) const {
       assert(!xKnots().empty());
       if (x < xKnots().front()) return false;
       if (x > xKnots().back()) return false;
       return true;
     }
 
     /// Check if q2 is in the grid range
     bool inRangeQ2(double q2) const {
       assert(!q2Knots().empty());
       if (q2 < q2Knots().front()) return false;
       if (q2 > q2Knots().back()) return false;
       return true;
     }
 
-    //@}
+    ///@}
 
 
   private:
 
     /// Map of multi-flavour KnotArrays "binned" for lookup by low edge in Q2
     std::map<double, KnotArrayNF> _knotarrays;
 
     // /// Caching vector of x knot values
     // mutable std::vector<double> _xknots;
 
     /// Caching vector of Q2 knot values
     mutable std::vector<double> _q2knots;
 
     /// Typedef of smart pointer for ipol memory handling
     typedef unique_ptr<Interpolator> InterpolatorPtr;
 
     /// Typedef of smart pointer for xpol memory handling
     typedef unique_ptr<Extrapolator> ExtrapolatorPtr;
 
     /// Associated interpolator (mutable to allow laziness)
     mutable InterpolatorPtr _interpolator;
 
     /// Associated extrapolator (mutable to allow laziness)
     mutable ExtrapolatorPtr _extrapolator;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/Info.h b/include/LHAPDF/Info.h
--- a/include/LHAPDF/Info.h
+++ b/include/LHAPDF/Info.h
@@ -1,212 +1,208 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Info_H
 #define LHAPDF_Info_H
 
 #include "LHAPDF/Utils.h"
 #include "LHAPDF/Paths.h"
 #include "LHAPDF/Exceptions.h"
 
 namespace LHAPDF {
 
 
   /// Get the singleton global configuration object
   ///
   /// @todo Move this out of Info. To Factories.h or SystemConfig.h?
   ///
   /// The global config is populated by reading from lhapdf.conf if it is
   /// found in the search paths.
   // class Info;
   // Info& config();
 
 
 
   /// Metadata base class for PDFs, PDF sets, or global configuration
   class Info {
   public:
 
-    /// @name Creation and deletion
-    //@{
-
     /// Default constructor
     Info() { }
 
     /// Constructor
     Info(const std::string& path) {
       load(path);
     }
 
     /// Virtual destructor to allow inheritance
     virtual ~Info() { }
 
-    //@}
-
 
     /// @name Loading info from YAML files
-    //@{
+    ///@{
 
     /// Populate this info object from the specified YAML file path.
     ///
     /// This function may be called several times to read metadata from several
     /// YAML source files. Values for existing keys will be overwritten.
     void load(const std::string& filepath);
 
-    //@}
+    ///@}
 
 
     /// @name General metadata accessors
-    //@{
+    ///@{
 
     // /// Get all metadata as a map
     // const std::map<std::string, std::string>& metadata() const {
     //   return _metadict;
     // }
 
     // /// Get all metadata as a map (non-const)
     // std::map<std::string, std::string>& metadata() {
     //   return _metadict;
     // }
 
 
     /// Is a value defined for the given key on this specific object?
     bool has_key_local(const std::string& key) const {
       return _metadict.find(key) != _metadict.end();
     }
 
     /// Can this object return a value for the given key?
     ///
     /// The given key may be defined non-locally, in which case the cascading
     /// member -> set -> config info lookup is needed. These are implemented
     /// using has_key_local() and metadata_local().
     ///
     /// The default implementation is equivalent to has_key_local(). This is
     /// appropriate for Config.
     virtual bool has_key(const std::string& key) const {
       return has_key_local(key);
     }
 
 
     /// Retrieve a metadata string by key name, as defined on this specific object
     const std::string& get_entry_local(const std::string& key) const {
       if (has_key_local(key)) return _metadict.find(key)->second;
       throw MetadataError("Metadata for key: " + key + " not found.");
     }
 
     /// Retrieve a metadata string by key name
     ///
     /// The given key may be defined non-locally, in which case the cascading
     /// member -> set -> config info lookup is needed. These are implemented
     /// using has_key_local() and get_entry_local().
     ///
     /// The default implementation is equivalent to get_entry_local(). This is
     /// appropriate for Config.
     virtual const std::string& get_entry(const std::string& key) const {
       return get_entry_local(key);
     }
 
 
     /// Retrieve a metadata string by key name, with a default fallback
     virtual const std::string& get_entry(const std::string& key, const std::string& fallback) const {
       try {
         return get_entry(key);
       } catch (...) {
         return fallback;
       }
     }
 
 
     /// Retrieve a metadata entry by key name, with an inline type cast
     ///
     /// Specialisations are defined below for unpacking of comma-separated lists
     /// of strings, ints, and doubles.
     template <typename T>
     T get_entry_as(const std::string& key) const {
       const string& s = get_entry(key);
       return lexical_cast<T>(s);
     }
 
 
     /// Retrieve a metadata entry by key name, with an inline type cast and default fallback
     template <typename T>
     T get_entry_as(const std::string& key, const T& fallback) const {
       try {
         return get_entry_as<T>(key);
       } catch (...) {
         return fallback;
       }
     }
 
 
     /// Set a keyed value entry
     template <typename T>
     void set_entry(const std::string& key, const T& val) {
       _metadict[key] = to_str(val);
     }
 
-    //@}
+    ///@}
 
 
   protected:
 
     /// The string -> string native metadata storage container
     std::map<std::string, std::string> _metadict;
 
   };
 
 
+
   /// @name Info metadata function template specialisations
-  //@{
+  ///@{
 
   template <>
   inline bool Info::get_entry_as(const std::string& key) const {
     const string& s = get_entry(key);
     try {
       bool rtn = lexical_cast<bool>(s);
       return rtn;
     } catch (...) {
       if (s == "true" || s == "on" || s == "yes") return true;
       if (s == "false" || s == "off" || s == "no") return false;
     }
     throw MetadataError("'" + s + "' is not a valid string for conversion to bool type");
   }
 
   template <>
   inline std::vector<std::string> Info::get_entry_as(const std::string& key) const {
     static const string delim = ",";
     string strval = trim(get_entry(key));
     // cout << "@@ " << strval << endl;
     if (startswith(strval, "[")) strval = strval.substr(1, strval.size()-1);
     if (endswith(strval, "]")) strval = strval.substr(0, strval.size()-1);
     // cout << "## " << strval << endl;
     return split(strval, delim);
   }
 
   template <>
   inline std::vector<int> Info::get_entry_as(const std::string& key) const {
     const vector<string> strs = get_entry_as< vector<string> >(key);
     vector<int> rtn;
     rtn.reserve(strs.size());
     for (const string& s : strs) rtn.push_back( lexical_cast<int>(s) );
     assert(rtn.size() == strs.size());
     return rtn;
   }
 
   template <>
   inline std::vector<double> Info::get_entry_as(const std::string& key) const {
     const vector<string> strs = get_entry_as< vector<string> >(key);
     vector<double> rtn;
     rtn.reserve(strs.size());
     for (const string& s : strs) rtn.push_back( lexical_cast<double>(s) );
     assert(rtn.size() == strs.size());
     return rtn;
   }
 
-  //@}
+  ///@}
 
 
 }
 #endif
diff --git a/include/LHAPDF/Interpolator.h b/include/LHAPDF/Interpolator.h
--- a/include/LHAPDF/Interpolator.h
+++ b/include/LHAPDF/Interpolator.h
@@ -1,86 +1,86 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Interpolator_H
 #define LHAPDF_Interpolator_H
 
 #include "LHAPDF/Utils.h"
 #include "LHAPDF/KnotArray.h"
 
 namespace LHAPDF {
 
 
   // Forward declaration
   class GridPDF;
 
 
   /// The general interface for interpolating between grid points
   class Interpolator {
   public:
 
     /// Destructor to allow inheritance
     virtual ~Interpolator() { }
 
 
     /// @name Binding to a PDF object
-    //@{
+    ///@{
 
     /// Bind to a GridPDF
     void bind(const GridPDF* pdf) { _pdf = pdf; }
 
     /// Unbind from GridPDF
     void unbind() { _pdf = 0; }
 
     /// Identify whether this Interpolator has an associated PDF
     bool hasPDF() { return _pdf != 0; }
 
     /// Get the associated GridPDF
     const GridPDF& pdf() const { return *_pdf; }
 
-    //@}
+    ///@}
 
 
     /// @name Interpolation methods
-    //@{
+    ///@{
 
     /// Interpolate a single-point in (x,Q)
     double interpolateXQ(int id, double x, double q) const {
       return interpolateXQ2(id, x, q*q);
     }
 
     /// Interpolate a single-point in (x,Q2)
     double interpolateXQ2(int id, double x, double q2) const;
 
 
     /// @todo Make an all-PID version of interpolateQ and Q2?
 
+    ///@}
+
 
   protected:
 
     /// @brief Interpolate a single-point in (x,Q2), given x/Q2 values and subgrid indices.
     ///
     /// The key function to be overridden in derived classes: the subgrid and
     /// x/Q2 index lookup (and their caching) are done centrally in the
     /// Interpolator base class so do not need to be re-implemented in each
     /// flavour of interpolator.
     virtual double _interpolateXQ2(const KnotArray1F& subgrid, double x, size_t ix, double q2, size_t iq2) const = 0;
 
     /// @todo Implement this NF version, with a cached KnotArrayNF?
     // virtual double _interpolateXQ2(const KnotArrayNF& subgrid, int id, double x, size_t ix, double q2, size_t iq2) const;
 
-    //@}
-
 
   private:
 
     const GridPDF* _pdf;
 
   };
 
 
 }
 
 #endif
diff --git a/include/LHAPDF/KnotArray.h b/include/LHAPDF/KnotArray.h
--- a/include/LHAPDF/KnotArray.h
+++ b/include/LHAPDF/KnotArray.h
@@ -1,338 +1,333 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_KnotArray_H
 #define LHAPDF_KnotArray_H
 
 #include "LHAPDF/Exceptions.h"
 
 namespace LHAPDF {
 
 
   /// @brief Internal storage class for PDF data point grids
   ///
   /// We use "array" to refer to the "raw" knot grid, while "grid" means a grid-based PDF.
   /// The "1F" means that this is a single-flavour array
   class KnotArray1F {
   public:
 
-    /// @name Construction etc.
-    //@{
-
     /// Default constructor just for std::map insertability
     KnotArray1F() {}
 
     /// Constructor from x and Q2 knot values, and an xf value grid as strided list
     KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots, const std::vector<double>& xfs)
       : _xs(xknots), _q2s(q2knots), _xfs(xfs)
     {
       assert(_xfs.size() == size());
       _synclogs();
     }
 
     /// Constructor of a zero-valued array from x and Q2 knot values
     KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots)
       : _xs(xknots), _q2s(q2knots),
         _xfs(size(), 0.0)
     {
       assert(_xfs.size() == size());
       _synclogs();
     }
 
-    //@}
-
 
     /// @name x stuff
-    //@{
+    ///@{
 
     /// x knot setter
     /// @note Also zeros the xfs array, which is invalidated by resetting the x knots
     void setxs(const std::vector<double>& xs) {
       _xs = xs;
       _synclogs();
       _xfs = std::vector<double>(size(), 0.0);
     }
 
     /// Number of x knots
     size_t xsize() const { return _xs.size(); }
 
     /// x knot accessor
     const std::vector<double>& xs() const { return _xs; }
 
     /// log(x) knot accessor
     const std::vector<double>& logxs() const { return _logxs; }
 
     /// @brief Get the index of the closest x knot row <= x
     ///
     /// If the value is >= x_max, return i_max-1 (for polynomial spine construction)
     size_t ixbelow(double x) const {
       // Test that x is in the grid range
       if (x < xs().front()) throw GridError("x value " + to_str(x) + " is lower than lowest-x grid point at " + to_str(xs().front()));
       if (x > xs().back()) throw GridError("x value " + to_str(x) + " is higher than highest-x grid point at " + to_str(xs().back()));
       // Find the closest knot below the requested value
       size_t i = upper_bound(xs().begin(), xs().end(), x) - xs().begin();
       if (i == xs().size()) i -= 1; // can't return the last knot index
       i -= 1; // have to step back to get the knot <= x behaviour
       return i;
     }
 
-    //@}
+    ///@}
 
 
     /// @name Q2 stuff
-    //@{
+    ///@{
 
     /// Q2 knot setter
     /// @note Also zeros the xfs array, which is invalidated by resetting the Q2 knots
     void setq2s(const std::vector<double>& q2s) {
       _q2s = q2s;
       _synclogs();
       _xfs = std::vector<double>(size(), 0.0);
     }
 
     /// Number of Q2 knots
     size_t q2size() const { return _q2s.size(); }
 
     /// Q2 knot accessor
     const std::vector<double>& q2s() const { return _q2s; }
 
     /// log(Q2) knot accessor
     const std::vector<double>& logq2s() const { return _logq2s; }
 
     /// Get the index of the closest Q2 knot row <= q2
     ///
     /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
     size_t iq2below(double q2) const {
       // Test that Q2 is in the grid range
       if (q2 < q2s().front()) throw GridError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
       if (q2 > q2s().back()) throw GridError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
       /// Find the closest knot below the requested value
       size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
       if (i == q2s().size()) i -= 1; // can't return the last knot index
       i -= 1; // have to step back to get the knot <= q2 behaviour
       return i;
     }
 
-    //@}
+    ///@}
 
 
     /// @name PDF values at (x, Q2) points
-    //@{
+    ///@{
 
     /// Number of x knots
     size_t size() const { return xsize()*q2size(); }
 
     /// xf value accessor (const)
     const std::vector<double>& xfs() const { return _xfs; }
     /// xf value accessor (non-const)
     std::vector<double>& xfs() { return _xfs; }
     /// xf value setter
     void setxfs(const std::vector<double>& xfs) { _xfs = xfs; }
 
     /// Get the xf value at a particular indexed x,Q2 knot
     const double& xf(size_t ix, size_t iq2) const { return _xfs[ix*q2size() + iq2]; }
 
-    //@}
+    ///@}
 
 
   private:
 
     /// Synchronise log(x) and log(Q2) arrays from the x and Q2 ones
     void _synclogs() {
       _logxs.resize(_xs.size());
       _logq2s.resize(_q2s.size());
       for (size_t i = 0; i < _xs.size(); ++i) _logxs[i] = log(_xs[i]);
       for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
     }
 
     /// List of x knots
     std::vector<double> _xs;
     /// List of Q2 knots
     std::vector<double> _q2s;
     /// List of log(x) knots
     std::vector<double> _logxs;
     /// List of log(Q2) knots
     std::vector<double> _logq2s;
     /// List of xf values across the 2D knot array, stored as a strided [ix][iQ2] 1D array
     std::vector<double> _xfs;
 
   };
 
 
   /// @brief A collection of {KnotArray1F}s accessed by PID code
   ///
   /// The "NF" means "> 1 flavour", cf. the KnotArray1F name for a single flavour data array.
   class KnotArrayNF {
   public:
 
     /// How many {KnotArray1F}s are stored in this container?
     size_t size() const { return _map.size(); }
 
     /// Is this container empty?
     bool empty() const { return _map.empty(); }
 
     /// Does this contain a KnotArray1F for PID code @a id?
     bool has_pid(int id) const {
       return _map.find(id) != _map.end();
     }
 
     /// Get the KnotArray1F for PID code @a id
     const KnotArray1F& get_pid(int id) const {
       if (!has_pid(id)) throw FlavorError("Undefined particle ID requested: " + to_str(id));
       return _map.find(id)->second;
     }
 
     /// Convenience accessor for any valid subgrid, to get access to the x/Q2/etc. arrays
     const KnotArray1F& get_first() const {
       if (empty()) throw GridError("Tried to access grid indices when no flavour grids were loaded");
       return _map.begin()->second;
     }
 
     /// Get the KnotArray1F for PID code @a id
     void set_pid(int id, const KnotArray1F& ka) {
       _map[id] = ka;
     }
 
     /// Indexing operator (non-const)
     KnotArray1F& operator[](int id) { return _map[id]; }
 
     /// Access the xs array
     const std::vector<double>& xs() const { return get_first().xs(); }
     /// Access the log(x)s array
     const std::vector<double>& logxs() const { return get_first().logxs(); }
     /// Get the index of the closest x knot column <= x (see KnotArray1F)
     size_t ixbelow(double x) const { return get_first().ixbelow(x); }
 
     /// Access the Q2s array
     const std::vector<double>& q2s() const { return get_first().q2s(); }
     /// Access the log(Q2)s array
     const std::vector<double>& logq2s() const { return get_first().logq2s(); }
     /// Get the index of the closest Q2 knot row <= q2 (see KnotArray1F)
     size_t iq2below(double q2) const { return get_first().iq2below(q2); }
 
   private:
 
     /// Storage
     std::map<int, KnotArray1F> _map;
 
   };
 
 
 
   /// Internal storage class for alpha_s interpolation grids
   class AlphaSArray {
   public:
 
     /// @name Construction etc.
-    //@{
+    ///@{
 
     /// Default constructor just for std::map insertability
     AlphaSArray() {}
 
     /// Constructor from Q2 knot values and alpha_s values
     AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
       : _q2s(q2knots), _as(as)
     {
       _synclogs();
     }
 
-    //@}
+    ///@}
 
 
     /// @name Q2 stuff
-    //@{
+    ///@{
 
     /// Q2 knot vector accessor
     const std::vector<double>& q2s() const { return _q2s; }
 
     /// log(Q2) knot vector accessor
     const std::vector<double>& logq2s() const { return _logq2s; }
 
     /// Get the index of the closest Q2 knot row <= q2
     ///
     /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
     size_t iq2below(double q2) const {
       // Test that Q2 is in the grid range
       if (q2 < q2s().front()) throw AlphaSError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
       if (q2 > q2s().back()) throw AlphaSError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
       /// Find the closest knot below the requested value
       size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
       if (i == q2s().size()) i -= 1; // can't return the last knot index
       i -= 1; // have to step back to get the knot <= q2 behaviour
       return i;
     }
 
     /// Get the index of the closest logQ2 knot row <= logq2
     ///
     /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
     size_t ilogq2below(double logq2) const {
       // Test that log(Q2) is in the grid range
       if (logq2 < logq2s().front()) throw GridError("logQ2 value " + to_str(logq2) + " is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
       if (logq2 > logq2s().back()) throw GridError("logQ2 value " + to_str(logq2) + " is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
       /// Find the closest knot below the requested value
       size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
       if (i == logq2s().size()) i -= 1; // can't return the last knot index
       i -= 1; // have to step back to get the knot <= q2 behaviour
       return i;
     }
 
-    //@}
+    ///@}
 
 
     /// @name alpha_s values at Q2 points
-    //@{
+    ///@{
 
     /// alpha_s value accessor (const)
     const std::vector<double>& alphas() const { return _as; }
     // /// alpha_s value accessor (non-const)
     // std::vector<double>& alphas() { return _as; }
     // /// alpha_s value setter
     // void setalphas(const valarray& xfs) { _as = as; }
 
-    //@}
+    ///@}
 
 
     /// @name alpha_s derivatives vs (log)Q2, useful for interpolation
-    //@{
+    ///@{
 
     /// Forward derivative w.r.t. logQ2
     double ddlogq_forward(size_t i) const {
       return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
     }
 
     /// Backward derivative w.r.t. logQ2
     double ddlogq_backward(size_t i) const {
       return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
     }
 
     /// Central (avg of forward and backward) derivative w.r.t. logQ2
     double ddlogq_central(size_t i) const {
       return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
     }
 
-    //@}
+    ///@}
 
 
   private:
 
     /// Synchronise the log(Q2) array from the Q2 one
     void _synclogs() {
       _logq2s.resize(_q2s.size());
       for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
     }
 
     /// List of Q2 knots
     std::vector<double> _q2s;
     /// List of log(Q2) knots
     std::vector<double> _logq2s;
     /// List of alpha_s values across the knot array
     std::vector<double> _as;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/PDF.h b/include/LHAPDF/PDF.h
--- a/include/LHAPDF/PDF.h
+++ b/include/LHAPDF/PDF.h
@@ -1,524 +1,522 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_PDF_H
 #define LHAPDF_PDF_H
 
 #include "LHAPDF/PDFInfo.h"
 #include "LHAPDF/PDFIndex.h"
 #include "LHAPDF/Factories.h"
 #include "LHAPDF/AlphaS.h"
 #include "LHAPDF/Utils.h"
 #include "LHAPDF/Paths.h"
 #include "LHAPDF/Exceptions.h"
 #include "LHAPDF/Version.h"
 #include "LHAPDF/Config.h"
 
 namespace LHAPDF {
 
 
   /// @brief PDF is the general interface for access to parton density information.
   ///
   /// The PDF interface declares the general form of all PDF types, such as Grid based or analytic.
   class PDF {
   protected: //< These constructors should only be called by subclasses
 
     /// Internal convenience typedef for the AlphaS object handle
     /// @todo Reinstate this unique_ptr when C++98 header compatibility is no longer an issue
     // typedef unique_ptr<AlphaS> AlphaSPtr;
     typedef AlphaS* AlphaSPtr;
 
     /// Force initialization of the only non-class member.
     /// @todo Remove _alphas initialisation when it can be a smart ptr again
     PDF() : _alphas(0), _forcePos(0) { }
 
 
   public:
 
     /// Virtual destructor, to allow unfettered inheritance
     virtual ~PDF() {
       /// @todo Remove this delete when C++98 is gone, and unique_ptr can be reinstated
       delete _alphas;
     }
 
-    //@}
-
 
   protected:
 
 
     /// @name Helper methods for info loading / path setting, used by derived types
-    //@{
+    ///@{
 
     void _loadInfo(const std::string& mempath);
 
     void _loadInfo(const std::string& setname, int member) {
       const string searchpath = findpdfmempath(setname, member);
       if (searchpath.empty())
         throw UserError("Can't find a valid PDF " + setname + "/" + to_str(member));
       _loadInfo(searchpath);
     }
 
     void _loadInfo(int lhaid) {
       const pair<string,int> setname_memid = lookupPDF(lhaid);
       if (setname_memid.second == -1)
         throw IndexError("Can't find a PDF with LHAPDF ID = " + to_str(lhaid));
       _loadInfo(setname_memid.first, setname_memid.second);
     }
 
-    //@}
+    ///@}
 
 
   public:
 
     /// @name PDF values
-    //@{
+    ///@{
 
     /// @brief Get the PDF xf(x) value at (x,q2) for the given PID.
     ///
     /// All grids are defined in Q2 rather than Q since the natural value
     /// in MC programs is squared, so we typically avoid an expensive sqrt() call.
     ///
     /// @param id PDG parton ID
     /// @param x Momentum fraction
     /// @param q2 Squared energy (renormalization) scale
     /// @return The value of xf(x,q2)
     double xfxQ2(int id, double x, double q2) const;
 
 
     /// @brief Get the PDF xf(x) value at (x,q) for the given PID.
     ///
     /// xfxQ will square the given q and return the value from xfxQ2.
     /// All grids are defined in q2 rather than q since the natural value
     /// in MC programs is squared, so we typically avoid an expensive sqrt() call.
     ///
     /// @param id PDG parton ID
     /// @param x Momentum fraction
     /// @param q Energy (renormalization) scale
     /// @return The value of xf(x,q2)
     double xfxQ(int id, double x, double q) const {
       return xfxQ2(id, x, q*q);
     }
 
 
     /// @brief Get the PDF xf(x) value at (x,q2) for all supported PIDs.
     ///
     /// This version fills a user-supplied map to avoid container construction
     /// costs on every call.
     ///
     /// @param x Momentum fraction
     /// @param q2 Squared energy (renormalization) scale
     /// @param rtn Map of PDF xf(x,q2) values, to be filled
     void xfxQ2(double x, double q2, std::map<int, double>& rtn) const;
 
 
     /// @brief Get the PDF xf(x) value at (x,q) for all supported PIDs.
     ///
     /// This version fills a user-supplied map to avoid container construction
     /// costs on every call.
     ///
     /// @param x Momentum fraction
     /// @param q Energy (renormalization) scale
     /// @param rtn Map of PDF xf(x,q) values, to be filled
     void xfxQ(double x, double q, std::map<int, double>& rtn) const {
       xfxQ2(x, q*q, rtn);
     }
 
 
     /// @brief Get the PDF xf(x) value at (x,q2) for "standard" PIDs.
     ///
     /// This version fills a user-supplied vector to avoid container
     /// construction costs on every call.
     ///
     /// The filled vector follows the LHAPDF5 convention, with 13 entries
     /// running in the PDF ID order [-6, -5, ..., -1, 21, 1, ... 5, 6], i.e.
     /// quark PDF values will be at vector index pid+6 and the gluon at index 6.
     ///
     /// @param x Momentum fraction
     /// @param q2 Squared energy (renormalization) scale
     /// @param rtn Vector of PDF xf(x,q2) values, to be filled
     void xfxQ2(double x, double q2, std::vector<double>& rtn) const;
 
     /// @brief Get the PDF xf(x) value at (x,q) for "standard" PIDs.
     ///
     /// This version fills a user-supplied vector to avoid container
     /// construction costs on every call.
     ///
     /// The filled vector follows the LHAPDF5 convention, with 13 entries
     /// running in the PDF ID order [-6, -5, ..., -1, 21, 1, ... 5, 6], i.e.
     /// quark PDF values will be at vector index pid+6 and the gluon at index 6.
     ///
     /// @param x Momentum fraction
     /// @param q Energy (renormalization) scale
     /// @param rtn Vector of PDF xf(x,q) values, to be filled
     void xfxQ(double x, double q, std::vector<double>& rtn) const {
       xfxQ2(x, q*q, rtn);
     }
 
 
     /// @brief Get the PDF xf(x) value at (x,q2) for all supported PIDs.
     ///
     /// This version creates a new map on every call: prefer to use the
     /// fill-in-place version with a user-supplied map for many calls.
     ///
     /// @param x Momentum fraction
     /// @param q2 Squared energy (renormalization) scale
     /// @return A map of PDF xf(x,q2) values
     std::map<int, double> xfxQ2(double x, double q2) const;
 
     /// @brief Get the PDF xf(x) value at (x,q) for all supported PIDs.
     ///
     /// This version creates a new map on every call: prefer to use the
     /// fill-in-place version with a user-supplied map for many calls.
     ///
     /// xfxQ will square the given q and return the value from xfxQ2.
     /// All grids are defined in q2 rather than q since the natural value
     /// in MC programs is squared, so we typically avoid an expensive sqrt() call.
     ///
     /// @param x Momentum fraction
     /// @param q Energy (renormalization) scale
     /// @return A map of PDF xf(x,q) values
     std::map<int, double> xfxQ(double x, double q) const {
       return xfxQ2(x, q*q);
     }
 
 
   protected:
 
     /// @brief Calculate the PDF xf(x) value at (x,q2) for the given PID.
     ///
     /// This is the key function to be overridden in concrete PDF types, since
     /// it actually does the calculation of xf(x,Q2) by analytic, interpolation,
     /// or other means. The user-called xfxQ2 method exists so that the physical
     /// range and PID checks need only be done in one place, rather than need to
     /// be re-implemented in each concrete implementation.
     ///
     /// @param id Parton ID in the PDG scheme
     /// @param x Momentum fraction
     /// @param q2 Squared energy (renormalization) scale
     /// @return the value of xf(x,q2)
     virtual double _xfxQ2(int id, double x, double q2) const = 0;
 
-    //@}
+    ///@}
 
 
   public:
 
     /// @name Ranges of validity
-    //@{
+    ///@{
 
     /// Minimum valid x value for this PDF.
     virtual double xMin() {
       if (info().has_key("XMin"))
         return info().get_entry_as<double>("XMin");
       return numeric_limits<double>::epsilon();
     }
 
     /// Maximum valid x value for this PDF.
     virtual double xMax() {
       if (info().has_key("XMax"))
         return info().get_entry_as<double>("XMax");
       return 1.0;
     }
 
     /// Minimum valid Q value for this PDF (in GeV).
     /// @note This function calls sqrt(q2Min()). For better CPU efficiency and accuracy use q2Min() directly.
     virtual double qMin() {
       return info().get_entry_as<double>("QMin", 0);
     }
 
     /// @brief Maximum valid Q value for this PDF (in GeV).
     /// @note This function calls sqrt(q2Max()). For better CPU efficiency and accuracy use q2Max() directly.
     virtual double qMax() {
       return info().get_entry_as<double>("QMax", numeric_limits<double>::max());
     }
 
     /// Minimum valid Q2 value for this PDF (in GeV2).
     virtual double q2Min() {
       return sqr(this->qMin());
     }
 
     /// Maximum valid Q2 value for this PDF (in GeV2).
     virtual double q2Max() {
       // Explicitly re-access this from the info, to avoid an overflow from squaring double_max
       return (info().has_key("QMax")) ? sqr(info().get_entry_as<double>("QMax")) : numeric_limits<double>::max();
     }
 
     /// @brief Check whether PDF is set to only return positive (definite) values or not.
     ///
     /// This is to avoid overshooting in to negative values when
     /// interpolating/extrapolating PDFs that sharply decrease towards zero.
     /// 0 = unforced, 1 = forced positive, 2 = forced positive definite (>= 1e-10)
     int forcePositive() const {
       if (_forcePos < 0) //< Caching
         _forcePos = info().get_entry_as<unsigned int>("ForcePositive", 0);
       return _forcePos;
     }
 
     /// @brief Check whether the given x is physically valid
     ///
     /// Returns false for x less than 0 or greater than 1, since it
     /// is a momentum fraction and not valid outside those values.
     bool inPhysicalRangeX(double x) const {
       return x >= 0.0 && x <= 1.0;
     }
 
     /// @brief Check whether the given Q2 is physically valid
     ///
     /// Returns false for Q2 less than 0 (Q must be real-valued).
     bool inPhysicalRangeQ2(double q2) const {
       return q2 >= 0.0;
     }
 
     /// @brief Check whether the given Q is physically valid
     ///
     /// Returns false for Q less than 0 (Q must be positive).
     bool inPhysicalRangeQ(double q) const {
       return inPhysicalRangeQ2(q*q);
     }
 
     /// Check whether the given (x,Q2) is physically valid
     bool inPhysicalRangeXQ2(double x, double q2) const {
       return inPhysicalRangeX(x) && inPhysicalRangeQ2(q2);
     }
 
     /// Check whether the given (x,Q) is physically valid
     bool inPhysicalRangeXQ(double x, double q) const {
       return inPhysicalRangeX(x) && inPhysicalRangeQ(q);
     }
 
     /// @brief Grid range check for Q
     ///
     /// Return true when given Q is in the coverage range of this PDF.
     /// It actually squares the given Q and returns value from inRangeQ2.
     ///
     /// @param q Energy scale
     /// @return Whether q is in range
     virtual bool inRangeQ(double q) const {
       return inRangeQ2(q*q);
     }
 
     /// @brief Grid range check for Q2
     ///
     /// Return true when given Q2 is in the coverage range of this PDF.
     ///
     /// @param q2 Squared energy scale
     /// @return Whether q2 is in range
     virtual bool inRangeQ2(double q2) const = 0;
 
     /// @brief Grid range check for x
     ///
     /// Return true when given x is in the coverage range of this PDF.
     ///
     /// @param x Momentum fraction
     /// @return Whether x is in range
     virtual bool inRangeX(double x) const = 0;
 
     /// Combined range check for x and Q
     virtual bool inRangeXQ(double x, double q) const {
       return inRangeX(x) && inRangeQ(q);
     }
 
     /// Combined range check for x and Q2
     bool inRangeXQ2(double x, double q2) const {
       return inRangeX(x) && inRangeQ2(q2);
     }
 
-    //@}
+    ///@}
 
 
     /// @name Generic member-level metadata (including cascaded metadata from set & config level)
-    //@{
+    ///@{
 
     /// Get the info class that actually stores and handles the metadata
     PDFInfo& info() { return _info; }
 
     /// Get the info class that actually stores and handles the metadata (const version)
     const PDFInfo& info() const { return _info; }
 
     /// @brief Get the PDF set of which this is a member
     ///
     /// Obtained from the member file path, not Info-based metadata.
     PDFSet& set() const {
       return getPDFSet(_setname());
     }
 
-    //@}
+    ///@}
 
 
     /// @name Member-level metadata
-    //@{
+    ///@{
 
     /// @brief PDF member local ID number
     ///
     /// Obtained from the member file path, not Info-based metadata.
     int memberID() const {
       const string memname = file_stem(_mempath);
       assert(memname.length() > 5); // There must be more to the filename stem than just the _nnnn suffix
       const int memid = lexical_cast<int>(memname.substr(memname.length()-4)); //< Last 4 chars should be the member number
       return memid;
     }
 
     /// @brief PDF member global LHAPDF ID number
     ///
     /// Obtained from the member ID and the set's LHAPDF ID index
     int lhapdfID() const;
 
     /// Description of this PDF member
     std::string description() const {
       return info().get_entry("PdfDesc", "");
     }
 
     /// Version of this PDF's data file
     int dataversion() const {
       return info().get_entry_as<int>("DataVersion", -1);
     }
 
     /// Get the type of PDF member that this object represents (central, error)
     std::string type() const {
       return to_lower(info().get_entry("PdfType"));
     }
 
-    //@}
+    ///@}
 
 
     /// Summary printout
     void print(std::ostream& os=std::cout, int verbosity=1) const;
 
 
     /// @name Parton content and QCD parameters
-    //@{
+    ///@{
 
     /// @brief List of flavours defined by this PDF set.
     ///
     /// This list is stored locally after its initial read from the Info object
     /// to avoid unnecessary lookups and string decoding, since e.g. it is
     /// looked at by every call to the GridPDF's Interpolator and Extrapolator
     /// classes.
     ///
     /// @todo Make virtual for AnalyticPDF? Or allow manual setting of the Info?
     virtual const std::vector<int>& flavors() const {
       if (_flavors.empty()) {
         _flavors = info().get_entry_as< vector<int> >("Flavors");
         sort(_flavors.begin(), _flavors.end());
       }
       return _flavors;
     }
 
     /// Checks whether @a id is a valid parton for this PDF.
     bool hasFlavor(int id) const;
 
     /// @brief Order of QCD at which this PDF has been constructed
     ///
     /// "Order" is defined here and throughout LHAPDF as the maximum number of
     /// loops included in the matrix elements, in order to have an integer value
     /// for easy use in comparisons, as opposed to "LO", "NLO", etc. strings.
     int orderQCD() const {
       return info().get_entry_as<int>("OrderQCD");
     }
     /// @deprecated Use orderQCD instead
     int qcdOrder() const { return orderQCD(); }
 
     /// @brief Get a quark mass in GeV by PDG code (|PID| = 1-6 only)
     ///
     /// Convenience interface to the Mass* info keywords.
     /// Returns -1 for an undefined PID.
     double quarkMass(int id) const;
 
     /// @brief Get a flavor scale threshold in GeV by PDG code (|PID| = 1-6 only)
     /// Convenience interface to the Mass* and Threshold* info keywords.
     /// Returns -1 for an undefined PID.
     double quarkThreshold(int id) const;
 
-    //@}
+    ///@}
 
 
     /// @name QCD running coupling calculation
-    //@{
+    ///@{
 
     /// @brief Set the AlphaS calculator by pointer
     ///
     /// The provided AlphaS must have been new'd, as it will not be copied
     /// and ownership passes to this GridPDF: delete will be called on this ptr
     /// when this PDF goes out of scope or another setAlphaS call is made.
     void setAlphaS(AlphaS* alphas) {
       // _alphas.reset(alphas);
       if (hasAlphaS()) delete _alphas;
       _alphas = alphas;
     }
 
     /// @brief Check if an AlphaS calculator is set
     bool hasAlphaS() const {
       return _alphas;
     }
 
     /// @brief Retrieve the AlphaS object for this PDF
     AlphaS& alphaS() {
       return *_alphas;
     }
 
     /// @brief Retrieve the AlphaS object for this PDF (const)
     const AlphaS& alphaS() const {
       return *_alphas;
     }
 
     /// @brief Value of alpha_s(Q2) used by this PDF
     ///
     /// Calculated numerically, analytically, or interpolated according to
     /// metadata, using the AlphaS classes.
     double alphasQ(double q) const {
       return alphasQ2(q*q);
     }
 
     /// @brief Value of alpha_s(Q2) used by this PDF
     ///
     /// Calculated numerically, analytically, or interpolated according to
     /// metadata, using the AlphaS classes.
     double alphasQ2(double q2) const {
       if (!hasAlphaS()) throw Exception("No AlphaS pointer has been set");
       return _alphas->alphasQ2(q2);
     }
 
-    //@}
+    ///@}
 
 
   protected:
 
     void _loadAlphaS() {
       // _alphas.reset( mkAlphaS(info()) );
       if (hasAlphaS()) delete _alphas;
       _alphas = mkAlphaS(info());
     }
 
     /// Get the set name from the member data file path (for internal use only)
     std::string _setname() const {
       return basename(dirname(_mempath));
     }
 
     /// Member data file path
     std::string _mempath;
 
     /// Metadata container
     PDFInfo _info;
 
     /// Locally cached list of supported PIDs
     mutable vector<int> _flavors;
 
     /// Optionally loaded AlphaS object
     mutable AlphaSPtr _alphas;
 
     /// @brief Cached flag for whether to return only positive (or postive definite) PDF values
     ///
     /// A negative value indicates that the flag has not been set. 0 = no
     /// forcing, 1 = force positive (i.e. 0 is permitted, negative values are
     /// not), 2 = force positive definite (i.e. no values less than 1e-10).
     mutable int _forcePos;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/PDFIndex.h b/include/LHAPDF/PDFIndex.h
--- a/include/LHAPDF/PDFIndex.h
+++ b/include/LHAPDF/PDFIndex.h
@@ -1,48 +1,48 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_PDFIndex_H
 #define LHAPDF_PDFIndex_H
 
 #include "LHAPDF/Utils.h"
 
 namespace LHAPDF {
 
 
-  /// @name Functions for PDF lookup by LHAPDF ID index file
-  //@{
+  /// @defgroup index PDF lookup in the LHAPDF ID index
+  ///@{
 
   /// Get the singleton LHAPDF set ID -> PDF index map
   std::map<int, std::string>& getPDFIndex();
 
   /// Look up a PDF set name and member ID by the LHAPDF ID code
   ///
   /// The set name and member ID are returned as an std::pair.
   /// If lookup fails, a pair ("", -1) is returned.
   std::pair<std::string, int> lookupPDF(int lhaid);
 
   /// @brief Decode a single PDF member ID string into a setname,memid pair
   ///
   /// @note A trivial <SET,MEM> decoding rather than a "real lookup", for convenience & uniformity.
   std::pair<std::string, int> lookupPDF(const std::string& pdfstr);
 
   /// Look up the member's LHAPDF index from the set name and member ID.
   ///
   /// If lookup fails, -1 is returned, otherwise the LHAPDF ID code.
   /// NB. This function is relatively slow, since it requires std::map reverse lookup.
   int lookupLHAPDFID(const std::string& setname, int nmem);
 
   /// Look up the member's LHAPDF index from a setname/member string.
   inline int lookupLHAPDFID(const std::string& setname_nmem) {
     const std::pair<string,int> idpair = lookupPDF(setname_nmem);
     return lookupLHAPDFID(idpair.first, idpair.second);
   }
 
-  //@}
+  ///@}
 
 
 }
 #endif
diff --git a/include/LHAPDF/PDFInfo.h b/include/LHAPDF/PDFInfo.h
--- a/include/LHAPDF/PDFInfo.h
+++ b/include/LHAPDF/PDFInfo.h
@@ -1,76 +1,71 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_PDFInfo_H
 #define LHAPDF_PDFInfo_H
 
 #include "LHAPDF/Info.h"
 #include "LHAPDF/Factories.h"
 #include "LHAPDF/PDFIndex.h"
 
 namespace LHAPDF {
 
 
   /// Metadata class for PDF members
   class PDFInfo : public Info {
   public:
 
-    /// @name Creation and deletion
-    //@{
-
     /// Default constructor (for container compatibility)
     ///
     /// @note Don't use explicitly!
     ///
     /// @todo Remove?
     PDFInfo() { }
 
     /// Constructor from a PDF member's data path.
     ///
     /// @todo Bypasses standard path searching hence used by the path-based
     /// GridPDF constructor, for example.
     PDFInfo(const std::string& mempath);
 
     /// Constructor from a set name and member ID.
     PDFInfo(const std::string& setname, int member);
 
     /// Constructor from an LHAPDF ID code.
     PDFInfo(int lhaid);
 
-    //@}
-
 
     /// @name Metadata accessors
-    //@{
+    ///@{
 
     /// Can this Info object return a value for the given key? (it may be defined non-locally)
     bool has_key(const std::string& key) const;
 
     /// Retrieve a metadata string by key name
     const std::string& get_entry(const std::string& key) const;
 
     /// Retrieve a metadata string by key name, with a fallback
     const std::string& get_entry(const std::string& key, const std::string& fallback) const {
       return Info::get_entry(key, fallback);
     }
 
-    //@}
+    ///@}
 
 
   private:
 
     /// Name of the set in which this PDF is contained (for PDFSet lookup)
     std::string _setname;
 
     /// Member ID in PDF set
     /// @note Not currently used, but could be useful if a memberID method is exposed.
     int _member;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/PDFSet.h b/include/LHAPDF/PDFSet.h
--- a/include/LHAPDF/PDFSet.h
+++ b/include/LHAPDF/PDFSet.h
@@ -1,298 +1,293 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_PDFSet_H
 #define LHAPDF_PDFSet_H
 
 #include "LHAPDF/Info.h"
 #include "LHAPDF/Factories.h"
 #include "LHAPDF/Version.h"
 #include "LHAPDF/Config.h"
 #include "LHAPDF/Utils.h"
 
 namespace LHAPDF {
 
 
   // Forward declaration
   class PDF;
 
 
   /// Structure for storage of uncertainty info calculated over a PDF error set
   struct PDFUncertainty {
     /// Constructor
     PDFUncertainty(double cent=0, double eplus=0, double eminus=0, double esymm=0, double scalefactor=1,
 		   double eplus_pdf=0, double eminus_pdf=0, double esymm_pdf=0, double e_par=0)
       : central(cent), errplus(eplus), errminus(eminus), errsymm(esymm), scale(scalefactor),
         errplus_pdf(eplus_pdf), errminus_pdf(eminus_pdf), errsymm_pdf(esymm_pdf), err_par(e_par)
     {    }
     /// Variables for the central value, +ve, -ve & symmetrised errors, and a CL scalefactor
     double central, errplus, errminus, errsymm, scale;
     /// Add extra variables for separate PDF and parameter variation errors with combined sets
     double errplus_pdf, errminus_pdf, errsymm_pdf, err_par;
   };
 
 
   /// Class for PDF set metadata and manipulation
   class PDFSet : public Info {
   public:
 
-    /// @name Creation and deletion
-    //@{
-
     /// Default constructor (for container compatibility)
     /// @todo Remove?
     PDFSet() { }
 
     /// Constructor from a set name
     /// @todo Remove?
     PDFSet(const std::string& setname);
 
-    //@}
-
 
     /// @name PDF set metadata specialisations
-    //@{
+    ///@{
 
     /// @brief PDF set name
     ///
     /// @note _Not_ taken from the .info metadata file.
     std::string name() const {
       return _setname;
     }
 
     /// Description of the set
     std::string description() const {
       return get_entry("SetDesc");
     }
 
     /// First LHAPDF global index in this PDF set
     int lhapdfID() const {
       return get_entry_as<int>("SetIndex", -1);
     }
 
     /// Version of this PDF set's data files
     int dataversion() const {
       return get_entry_as<int>("DataVersion", -1);
     }
 
     /// Get the type of PDF errors in this set (replicas, symmhessian, hessian, custom, etc.)
     std::string errorType() const {
       return to_lower(get_entry("ErrorType", "UNKNOWN"));
     }
 
     /// @brief Get the confidence level of the Hessian eigenvectors, in percent.
     ///
     /// If not defined, assume 1-sigma = erf(1/sqrt(2)) =~ 68.268949% by default,
     /// unless this is a replica set for which return -1.
     double errorConfLevel() const;
 
     /// Number of members in this set
     // int numMembers() const {
     //   return get_entry_as<int>("NumMembers");
     // }
     size_t size() const {
       return get_entry_as<unsigned int>("NumMembers");
     }
 
-    //@}
+    ///@}
 
 
     /// Summary printout
     void print(std::ostream& os=std::cout, int verbosity=1) const;
 
 
     /// @name Creating PDF members
-    //@{
+    ///@{
 
     /// Make the nth PDF member in this set, returned by pointer
     ///
     /// @note As with the mkPDF factory method, the PDF pointer returned by this
     /// method is heap allocated and its memory management is now the
     /// responsibility of the caller.
     PDF* mkPDF(int member) const {
       return LHAPDF::mkPDF(name(), member);
     }
 
 
     /// Make all the PDFs in this set, filling a supplied vector with PDF pointers
     ///
     /// This version may be preferred in many circumstances, since it can avoid
     /// the overhead of creating a new temporary vector.
     ///
     /// A vector of *smart* pointers can be used, for any smart pointer type which
     /// supports construction from a raw pointer argument, e.g. unique_ptr<PDF>(PDF*).
     ///
     /// @note The supplied vector will be cleared before filling!
     ///
     /// @note As with the mkPDF method and factory function, the PDF pointers
     /// returned by this method are heap allocated and their memory management
     /// is now the responsibility of the caller, either manually for raw pointers
     /// or automatically if smart pointers are used.
     ///
     /// @note Use an *appropriate* smart pointer, of course! This depends in
     /// detail on how you will use the PDF objects (do you want shared or unique
     /// pointers?), but they also need to be compatible with storage in STL
     /// containers, e.g. std::unique_ptr or std::shared_ptr but *not* the
     /// deprecated std::auto_ptr.
     //
     /// @todo Needs to be implemented in the header since the arg type is templated.
     template <typename PTR>
     void mkPDFs(std::vector<PTR>& pdfs) const {
       const int v = verbosity();
       if (v > 0) {
         std::cout << "LHAPDF " << version() << " loading all " << size() << " PDFs in set " << name() << std::endl;
         this->print(std::cout, v);
         if (this->has_key("Note")) std::cout << get_entry("Note") << std::endl;
       }
       pdfs.clear();
       pdfs.reserve(size());
       if (v < 2) setVerbosity(0); //< Disable every-member printout unless verbosity level is high
       for (size_t i = 0; i < size(); ++i) {
         /// @todo Need to use an std::move here, or write differently, for unique_ptr to work?
         pdfs.push_back( PTR(mkPDF(i)) );
       }
       setVerbosity(v);
     }
 
     /// Make all the PDFs in this set, returning as a vector of PDF pointers
     ///
     /// @note As with the mkPDF method and factory function, the PDF pointers
     /// returned by this method are heap allocated and their memory management
     /// is now the responsibility of the caller.
     std::vector<PDF*> mkPDFs() const {
       std::vector<PDF*> rtn;
       mkPDFs(rtn);
       return rtn;
     }
 
     /// @todo Use the following with default function template args if C++11 is being used
     // template <typename PTR=PDF*>
     template <typename PTR>
     std::vector<PTR> mkPDFs() const {
       std::vector<PTR> rtn;
       mkPDFs(rtn);
       return rtn;
     }
 
-    //@}
+    ///@}
 
 
     /// @todo Add AlphaS getter for set-level alphaS?
 
 
     /// @name Generic metadata cascading mechanism
-    //@{
+    ///@{
 
     /// Can this Info object return a value for the given key? (it may be defined non-locally)
     bool has_key(const std::string& key) const {
       return has_key_local(key) || getConfig().has_key(key);
     }
 
     /// Retrieve a metadata string by key name
     const std::string& get_entry(const std::string& key) const {
       if (has_key_local(key)) return get_entry_local(key); //< value is defined locally
       return getConfig().get_entry(key); //< fall back to the global config
     }
 
     /// Retrieve a metadata string by key name, with a fallback
     const std::string& get_entry(const std::string& key, const std::string& fallback) const {
       return Info::get_entry(key, fallback);
     }
 
-    //@}
+    ///@}
 
 
     /// @name PDF set uncertainty functions
-    //@{
+    ///@{
 
     /// @brief Calculate central value and error from vector @c values with appropriate formulae for this set
     ///
     /// In the Hessian approach, the central value is the best-fit
     /// "values[0]" and the uncertainty is given by either the symmetric or
     /// asymmetric formula using eigenvector PDF sets.
     ///
     /// If the PDF set is given in the form of replicas, by default, the central value is
     /// given by the mean and is not necessarily "values[0]" for quantities with a non-linear
     /// dependence on PDFs, while the uncertainty is given by the standard deviation.
     ///
     /// Optional argument @c clpct is used to rescale uncertainties to a
     /// particular confidence level (in percent); a negative number will rescale to the
     /// default CL for this set.
     ///
     /// @note If @c cl is omitted, automatically rescale to normal 1-sigma ~ 68.268949% uncertainties.
     ///
     /// If the PDF set is given in the form of replicas, then optional argument
     /// @c alternative equal to true (default: false) will construct a confidence
     /// interval from the probability distribution of replicas, with the central
     /// value given by the median.
     ///
     /// For a combined set, a breakdown of the separate PDF and parameter
     /// variation uncertainties is available.  The parameter variation
     /// uncertainties are computed from the last 2*n members of the set, with n
     /// the number of parameters.
     PDFUncertainty uncertainty(const std::vector<double>& values,
                                double cl=100*erf(1/sqrt(2)), bool alternative=false) const;
 
     /// Calculate PDF uncertainties (as above), with efficient no-copy return to the @c rtn argument.
     /// @todo For real efficiency, the chaining of these functions should be the other way around
     void uncertainty(PDFUncertainty& rtn,
                      const std::vector<double>& values,
                      double cl=100*erf(1/sqrt(2)), bool alternative=false) const {
       rtn = uncertainty(values, cl, alternative);
     }
 
     /// @brief Calculate the PDF correlation between @c valuesA and @c valuesB using appropriate formulae for this set.
     ///
     /// The correlation can vary between -1 and +1 where values close to {-1,0,+1} mean that the two
     /// quantities A and B are {anticorrelated,uncorrelated,correlated}, respectively.
     ///
     /// For a combined set, the parameter variations are not included in the calculation of the correlation.
     double correlation(const std::vector<double>& valuesA, const std::vector<double>& valuesB) const;
 
     /// @brief Generate a random value from Hessian @c values and Gaussian random numbers.
     ///
     /// @note This routine is intended for advanced users!
     ///
     /// See Section 6 of G. Watt and R.S. Thorne, JHEP 1208 (2012) 052 [arXiv:1205.4024 [hep-ph]].
     ///
     /// Pass a vector @c values containing a value for each member of the Hessian PDF set.
     /// Pass a vector @c randoms containing neigen random numbers, where neigen is the number of distinct eigenvectors.
     ///
     /// Option @c symmetrise equal to true will symmetrise the random values (in the case of an asymmetric Hessian set)
     /// using a corrected Eq. (6.5) of arXiv:1205.4024v2, so that the average tends to the best-fit for a large number of replicas.
     ///
     /// Option @c symmetrise equal to false will use Eq. (6.4) of arXiv:1205.4024v2 (for an asymmetric Hessian set),
     /// then the average differs from the best-fit.  Option @c symmetrise has no effect for a symmetric Hessian set.
     ///
     /// Random values generated in this way can subsequently be used for applications such as Bayesian reweighting
     /// or combining predictions from different groups (as an alternative to taking the envelope).
     /// See, for example, supplementary material at http://mstwpdf.hepforge.org/random/.
     ///
     /// Use of this routine with a non-Hessian PDF set will throw a UserError.
     ///
     /// For a combined set, the parameter variations are not included in the generation of the random value.
     double randomValueFromHessian(const std::vector<double>& values, const std::vector<double>& randoms, bool symmetrise=true) const;
 
 
     /// Check that the PdfType of each member matches the ErrorType of the set.
     /// @todo We need to make the signature clearer -- what is the arg? Why not
     ///   automatically check the members? Why not a plural name? Why not on PDF?
     ///   "Hiding" the name for now with the leading underscore.
     void _checkPdfType(const std::vector<string>& pdftypes) const;
 
-    //@}
+    ///@}
 
 
   private:
 
     /// Name of this set
     std::string _setname;
 
   };
 
 
 }
 #endif
diff --git a/include/LHAPDF/Paths.h b/include/LHAPDF/Paths.h
--- a/include/LHAPDF/Paths.h
+++ b/include/LHAPDF/Paths.h
@@ -1,96 +1,100 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Paths_H
 #define LHAPDF_Paths_H
 
 #include "LHAPDF/Utils.h"
 
 namespace LHAPDF {
 
 
-  /// @name File searching and search path handling functions
-  //@{
+  /// @defgroup search Searching for PDF data
+  ///@{
+
+
+  /// @name Search path handling functions
+  ///@{
 
   /// @brief Get the ordered list of search paths, from $LHAPDF_DATA_PATH and the install location
   /// @note The install prefix will be appended *unless* $LHAPDF_DATA_PATH ends with a double colon, i.e. '::'
   std::vector<std::string> paths();
 
   /// Set the search paths list as a colon-separated string
   void setPaths(const std::string& pathstr);
   /// Set the search paths list
   inline void setPaths(std::vector<string> paths) {
     setPaths(join(paths, ":"));
   }
 
 
   /// Prepend to the search paths list
   inline void pathsPrepend(const std::string& p) {
     vector<string> ps = paths();
     ps.insert(ps.begin(), p);
     //ps.pop_back(); //< Discard the auto-added fallback path to the installed data prefix
     setPaths(ps);
   }
 
 
   /// Append to the search paths list
   inline void pathsAppend(const std::string& p) {
     vector<string> ps = paths();
     //ps.pop_back(); //< Discard the auto-added fallback path to the installed data prefix
     ps.push_back(p);
     setPaths(ps);
   }
 
 
   /// Return the first location in which a file is found
   ///
   /// If no matching file is found, return an empty path.
   std::string findFile(const std::string& target);
-  //@}
 
+  ///@}
 
 
-  /// @name Functions for handling standard LHAPDF filename structures
-  //@{
+  /// @name PDF-specific path functions
+  ///@{
 
   inline std::string pdfmempath(const std::string& setname, int member) {
     const string memname = setname + "_" + to_str_zeropad(member) + ".dat";
     const string mempath = setname / memname;
     return mempath;
   }
 
   inline std::string findpdfmempath(const std::string& setname, int member) {
     return findFile(pdfmempath(setname, member));
   }
 
   inline std::string pdfsetinfopath(const std::string& setname) {
     const string infoname = setname + ".info";
     const string setinfo = setname / infoname;
     return setinfo;
   }
 
   inline std::string findpdfsetinfopath(const std::string& setname) {
     /// @todo Check that set info and mem=0 file are in same dir?
     return findFile(pdfsetinfopath(setname));
   }
 
-  //@}
-
 
   /// @brief Get the names of all available PDF sets in the search path
   ///
   /// @note Taken from scanning the directories in the search path
   /// (i.e. LHAPDF_DATA_PATH) for viable PDF sets.
   ///
   /// @note The result is cached when first called, to avoid repeated filesystem
   /// walking. It's assumed that new PDFs will not appear on the filesystem
   /// during a run: please let the authors know if that's not a good assumption!
   const std::vector<std::string>& availablePDFSets();
 
+  ///@}
+
 
 }
 
 #endif
diff --git a/include/LHAPDF/Reweighting.h b/include/LHAPDF/Reweighting.h
--- a/include/LHAPDF/Reweighting.h
+++ b/include/LHAPDF/Reweighting.h
@@ -1,103 +1,107 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Reweighting_H
 #define LHAPDF_Reweighting_H
 
 #include "LHAPDF/PDF.h"
 #include "LHAPDF/PDFSet.h"
 
 namespace LHAPDF {
 
 
+  /// @name reweight PDF reweighting
+  ///@{
+
   namespace {
     inline bool _checkAlphasQ2(double Q2, const PDF& pdfa, const PDF& pdfb, double aschk) {
       if (aschk < 0) return true;
       const double as_a = pdfa.alphasQ2(Q2);
       const double as_b = pdfb.alphasQ2(Q2);
       if (2 * std::abs(as_a - as_b) / (std::abs(as_a) + std::abs(as_b)) < aschk) return true;
       std::cerr << "WARNING: alpha_s(Q2) mismatch in PDF reweighting "
                 << "at Q2 = " << Q2 << " GeV2:\n  "
                 << as_a << " for " << pdfa.set().name() << "/" << pdfa.memberID() << " vs. "
                 << as_b << " for " << pdfb.set().name() << "/" << pdfb.memberID()
                 << std::endl;
       return false;
     }
   }
 
 
-  /// @name Single beam reweighting functions
-  //@{
+  /// @name Single-beam reweighting functions
+  ///@{
 
   /// Get the PDF reweighting factor for a beam with id,x,Q parameters, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   inline double weightxQ2(int id, double x, double Q2, const PDF& basepdf, const PDF& newpdf, double aschk=5e-2) {
     if (aschk >= 0) _checkAlphasQ2(Q2, basepdf, newpdf, aschk);
     const double xf_base = basepdf.xfxQ2(id, x, Q2);
     const double xf_new = newpdf.xfxQ2(id, x, Q2);
     return xf_new / xf_base;
   }
 
   /// Get the PDF reweighting factor for a beam with id,x,Q parameters, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   template <typename PDFPTR>
   inline double weightxQ2(int id, double x, double Q2, const PDFPTR basepdf, const PDFPTR newpdf, double aschk=5e-2) {
     return weightxQ2(id, x, Q2, *basepdf, *newpdf, aschk);
   }
 
   /// Get the PDF reweighting factor for a beam with id,x,Q parameters, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   inline double weightxQ(int id, double x, double Q, const PDF& basepdf, const PDF& newpdf, double aschk=5e-2) {
     return weightxQ2(id, x, sqr(Q), basepdf, newpdf, aschk);
   }
 
   /// Get the PDF reweighting factor for a beam with id,x,Q parameters, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   template <typename PDFPTR>
   inline double weightxQ(int id, double x, double Q, const PDFPTR basepdf, const PDFPTR newpdf, double aschk=5e-2) {
     return weightxQ(id, x, Q, *basepdf, *newpdf, aschk);
   }
 
-  //@}
+  ///@}
 
 
   /// @name Two-beam reweighting functions
-  //@{
+  ///@{
 
   /// Get the PDF reweighting factor for two beams, one with id1,x1 and the other with id2,x2, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   inline double weightxxQ2(int id1, int id2, double x1, double x2, double Q2, const PDF& basepdf, const PDF& newpdf, double aschk=5e-2) {
     if (aschk >= 0) _checkAlphasQ2(Q2, basepdf, newpdf, aschk);
     const double w1 = weightxQ2(id1, x1, Q2, basepdf, newpdf, -1);
     const double w2 = weightxQ2(id2, x2, Q2, basepdf, newpdf, -1);
     return w1 * w2;
   }
 
   /// Get the PDF reweighting factor for two beams, one with id1,x1 and the other with id2,x2, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   template <typename PDFPTR>
   inline double weightxxQ2(int id1, int id2, double x1, double x2, double Q2, const PDFPTR basepdf, const PDFPTR newpdf, double aschk=5e-2) {
     return weightxxQ2(id1, id2, x1, x2, Q2, *basepdf, *newpdf, aschk);
   }
 
   /// Get the PDF reweighting factor for two beams, one with id1,x1 and the other with id2,x2, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   inline double weightxxQ(int id1, int id2, double x1, double x2, double Q, const PDF& basepdf, const PDF& newpdf, double aschk=5e-2) {
     return weightxxQ2(id1, id2, x1, x2, sqr(Q), basepdf, newpdf, aschk);
   }
 
   /// Get the PDF reweighting factor for two beams, one with id1,x1 and the other with id2,x2, from basepdf to newpdf
   /// @note For NLO calculations, in general different PDF values enter for each counterterm: be careful.
   template <typename PDFPTR>
   inline double weightxxQ(int id1, int id2, double x1, double x2, double Q, const PDFPTR basepdf, const PDFPTR newpdf, double aschk=5e-2) {
     return weightxxQ(id1, id2, x1, x2, Q, *basepdf, *newpdf, aschk);
   }
 
-  //@}
+  ///@}
 
+  ///@}
 
 }
 #endif
diff --git a/include/LHAPDF/Utils.h b/include/LHAPDF/Utils.h
--- a/include/LHAPDF/Utils.h
+++ b/include/LHAPDF/Utils.h
@@ -1,297 +1,301 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #pragma once
 #ifndef LHAPDF_Utils_H
 #define LHAPDF_Utils_H
 
 // STL includes
 #include <cassert>
 #include <stdexcept>
 #include <vector>
 #include <map>
 #include <string>
 #include <algorithm>
 #include <memory>
 #include <iostream>
 #include <iomanip>
 #include <sstream>
 #include <limits>
 #include <cmath>
 
 
 /// Namespace for all LHAPDF functions and classes
 namespace LHAPDF {
 
 
   // Allow implicit use of the std namespace within namespace LHAPDF
   using namespace std;
 
 
+  /// @defgroup utils Internal utility functions
+  ///@{
+
   /// @name String handling utility functions
-  //@{
+  ///@{
 
   /// When lexical_cast goes bad
   struct bad_lexical_cast : public std::runtime_error {
+    /// Constructor
     bad_lexical_cast(const std::string& what) : std::runtime_error(what) {}
   };
 
-  /// @brief Convert between any types via stringstream
+  /// @brief Convert between types via stringstream
   template<typename T, typename U>
   T lexical_cast(const U& in) {
     try {
       std::stringstream ss;
       ss << in;
       T out;
       ss >> out;
       return out;
     } catch (const std::exception& e) {
       throw bad_lexical_cast(e.what());
     }
   }
 
   /// Make a string representation of @a val
   template <typename T>
   inline std::string to_str(const T& val) {
     return lexical_cast<string>(val);
   }
 
   /// Make a string representation of a vector @a vec
   template <typename T>
   inline std::string to_str(const std::vector<T>& vec) {
     string rtn = "[";
     for (size_t i = 0; i < vec.size(); ++i) {
       rtn += to_str(vec[i]);
       if (i < vec.size()-1) rtn += ", ";
     }
     rtn += "]";
     return rtn;
   }
 
   /// Format an integer @a val as a zero-padded string of length @a nchars
   inline std::string to_str_zeropad(int val, size_t nchars=4) {
     stringstream ss;
     ss << setfill('0') << setw(nchars) << val;
     return ss.str();
   }
 
   /// Concatenate strings with separator strings between each element
   inline std::string join(const std::vector<std::string>& svec, const std::string& sep) {
     string rtn;
     for (size_t i = 0; i < svec.size(); ++i) {
       rtn += svec[i];
       if (i < svec.size()-1) rtn += sep;
     }
     return rtn;
   }
 
   /// Split a string by a given separator
   inline std::vector<std::string> split(const std::string& s, const std::string& sep) {
     vector<string> rtn;
     string tmp = s; // non-const working copy, to be incrementally truncated
     while (true) {
       const size_t delim_pos = tmp.find(sep);
       if (delim_pos == string::npos) break;
       const string s = tmp.substr(0, delim_pos);
       if (!s.empty()) rtn.push_back(s); // Don't insert "empties"
       tmp.replace(0, delim_pos+1, ""); // Remove already-processed part
     }
     if (!tmp.empty()) rtn.push_back(tmp); // Don't forget the trailing component!
     return rtn;
   }
 
   /// Does a string @a s contain the @a sub substring?
   inline bool contains(const std::string& s, const std::string& sub) {
     return s.find(sub) != string::npos;
   }
 
   /// Does a string @a s start with the @a sub substring?
   inline bool startswith(const std::string& s, const std::string& sub) {
     return s.find(sub) == 0;
   }
 
   /// Does a string @a s end with the @a sub substring?
   inline bool endswith(const std::string& s, const std::string& sub) {
     return s.find(sub) == s.length()-sub.length();
   }
 
   /// How many times does a string @a s contain the character @a c?
   inline size_t countchar(const std::string& s, const char c) {
     return std::count(s.begin(), s.end(), c);
   }
 
   /// Strip leading and trailing spaces (not in-place)
   inline std::string trim(const std::string& s) {
     const size_t firstnonspacepos = s.find_first_not_of(" ");
     const size_t lastnonspacepos = s.find_last_not_of(" ");
     if (firstnonspacepos == std::string::npos) return "";
     return s.substr(firstnonspacepos, lastnonspacepos-firstnonspacepos+1);
   }
 
   /// Convert a string to lower-case (not in-place)
   inline std::string to_lower(const std::string& s) {
     string rtn = s;
     transform(rtn.begin(), rtn.end(), rtn.begin(), (int(*)(int)) tolower);
     return rtn;
   }
 
   /// Convert a string to upper-case (not in-place)
   inline std::string to_upper(const std::string& s) {
     string rtn = s;
     transform(rtn.begin(), rtn.end(), rtn.begin(), (int(*)(int)) toupper);
     return rtn;
   }
 
-  //@}
+  ///@}
 
 
-  /// @name Generic path functions in the LHAPDF namespace
-  //@{
+  /// @name Filesystem utils
+  ///@{
 
   /// Check if a path @a p (either file or dir) exists
   bool path_exists(const std::string& p,int mode=0);
 
   /// Check if a file @a p exists
   bool file_exists(const std::string& p,int mode=0);
 
   /// Check if a dir @a p exists
   bool dir_exists(const std::string& p,int mode=0);
 
   /// Operator for joining strings @a a and @a b with filesystem separators
   inline std::string operator / (const std::string& a, const std::string& b) {
     // Ensure that a doesn't end with a slash, and b doesn't start with one, to avoid "//"
     const string anorm = (a.find("/") != std::string::npos) ? a.substr(0, a.find_last_not_of("/")+1) : a;
     const string bnorm = (b.find("/") != std::string::npos) ? b.substr(b.find_first_not_of("/")) : b;
     return anorm + "/" + bnorm;
   }
 
   /// Get the basename (i.e. terminal file name) from a path @a p
   inline std::string basename(const std::string& p) {
     if (!contains(p, "/")) return p;
     return p.substr(p.rfind("/")+1);
   }
 
   /// Get the dirname (i.e. path to the penultimate directory) from a path @a p
   inline std::string dirname(const std::string& p) {
     if (!contains(p, "/")) return "";
     return p.substr(0, p.rfind("/"));
   }
 
   /// Get the stem (i.e. part without a file extension) from a filename @a f
   inline std::string file_stem(const std::string& f) {
     if (!contains(f, ".")) return f;
     return f.substr(0, f.rfind("."));
   }
 
   /// Get the file extension from a filename @a f
   inline std::string file_extn(const std::string& f) {
     if (!contains(f, ".")) return "";
     return f.substr(f.rfind(".")+1);
   }
 
   /// @todo Add an abspath(p) function
 
-  //@}
+  ///@}
 
 
-  /// @name Math functions in the LHAPDF namespace
-  //@{
+  /// @name Math functions
+  ///@{
 
   /// Convenience function for squaring (of any type)
   template <typename N>
   inline N sqr(const N& x) { return x*x; }
 
   /// Get the sign of a number
   template <typename N>
   inline int sgn(N val) { return (N(0) < val) - (val < N(0)); }
 
   /// Check if a number is in a range (closed-open)
   inline int in_range(double x, double low, double high) { return x >= low && x < high; }
 
   /// Check if a number is in a range (closed-closed)
   inline int in_closed_range(double x, double low, double high) { return x >= low && x <= high; }
 
   /// Check if a number is in a range (open-open)
   inline int in_open_range(double x, double low, double high) { return x > low && x < high; }
 
   /// @todo Add iszero() & equals(,) functions?
 
   /// Quantiles of the standard normal probability distribution function
   double norm_quantile(double p);
 
   /// Quantiles of the chi-squared probability distribution function
   double chisquared_quantile(double p, double ndf);
 
-  //@}
+  ///@}
 
 
-  /// @name Container handling helpers
-  //@{
+  /// @name Container utils
+  ///@{
 
   /// Does the vector<T> @a container contain @a item?
   template <typename T>
   inline bool contains(const std::vector<T>& container, const T& item) {
     return find(container.begin(), container.end(), item) != container.end();
   }
 
   // /// Does the set<T> @a container contain @a item?
   // template <typename T>
   // inline bool contains(const std::set<T>& container, const T& item) {
   //   return container.find(item) != container.end();
   // }
 
   /// Does the map<K,T> @a container have a key K @a key?
   template <typename K, typename T>
   inline bool has_key(const std::map<K,T>& container, const K& key) {
     return container.find(key) != container.end();
   }
 
   // /// @name Implementation of generic begin/end container identification by traits
   // /// taken from http://stackoverflow.com/a/9407420/91808 . Needs C++11 (or maybe just C++0x).
   // //@{
 
   // #include <type_traits>
 
   // template<typename T>
   // struct has_const_iterator {
   // private:
   //   typedef char                      yes;
   //   typedef struct { char array[2]; } no;
   //   template<typename C> static yes test(typename C::const_iterator*);
   //   template<typename C> static no  test(...);
   // public:
   //   static const bool value = sizeof(test<T>(0)) == sizeof(yes);
   //   typedef T type;
   // };
 
   // template <typename T>
   // struct has_begin_end {
   //   template<typename C> static char (&f(typename std::enable_if<
   //     std::is_same<decltype(static_cast<typename C::const_iterator (C::*)() const>(&C::begin)),
   //     typename C::const_iterator(C::*)() const>::value, void>::type*))[1];
 
   //   template<typename C> static char (&f(...))[2];
 
   //   template<typename C> static char (&g(typename std::enable_if<
   //     std::is_same<decltype(static_cast<typename C::const_iterator (C::*)() const>(&C::end)),
   //     typename C::const_iterator(C::*)() const>::value, void>::type*))[1];
 
   //   template<typename C> static char (&g(...))[2];
 
   //   static bool const beg_value = sizeof(f<T>(0)) == 1;
   //   static bool const end_value = sizeof(g<T>(0)) == 1;
   // };
 
   // template<typename T>
   // struct is_container
   //   : std::integral_constant<bool, has_const_iterator<T>::value &&
   //                            has_begin_end<T>::beg_value && has_begin_end<T>::end_value>
   // { };
 
-  //@}
+  ///@}
 
 
 }
 #endif
diff --git a/src/Utils.cc b/src/Utils.cc
--- a/src/Utils.cc
+++ b/src/Utils.cc
@@ -1,362 +1,359 @@
 // -*- C++ -*-
 //
 // This file is part of LHAPDF
 // Copyright (C) 2012-2019 The LHAPDF collaboration (see AUTHORS for details)
 //
 #include "LHAPDF/Utils.h"
 
 #include <sys/stat.h>
 #ifdef HAVE_MPI
 #include <mpi.h>
 #endif
 
 namespace LHAPDF {
 
 
   /// Check if a path @a p (either file or dir) exists
   bool path_exists(const std::string& p,int mode)
   {
     return file_exists(p,mode) || dir_exists(p,mode);
   }
 
   /// Check if a file @a p exists
   bool file_exists(const std::string& file,int mode)
   {
     int exists(false);
     #ifdef HAVE_MPI
     if (mode==1 || MPI::COMM_WORLD.Get_rank()==0) {
       struct stat fst;
       if (stat(file.c_str(), &fst) != -1) exists = (fst.st_mode & S_IFMT) == S_IFREG;
     }
     if (mode!=1) MPI::COMM_WORLD.Bcast(&exists,1,MPI_INT,0);
     #else
     struct stat fst;
     if (stat(file.c_str(), &fst) != -1) exists = (fst.st_mode & S_IFMT) == S_IFREG;
     #endif
     return exists;
   }
 
   /// Check if a dir @a p exists
   bool dir_exists(const std::string& dir,int mode)
   {
     int exists(false);
     #ifdef HAVE_MPI
     if (mode==1 || MPI::COMM_WORLD.Get_rank()==0) {
       struct stat fst;
       if (stat(dir.c_str(), &fst) != -1) exists = (fst.st_mode & S_IFMT) == S_IFDIR;
     }
     if (mode!=1) MPI::COMM_WORLD.Bcast(&exists,1,MPI_INT,0);
     #else
     struct stat fst;
     if (stat(dir.c_str(), &fst) != -1) exists = (fst.st_mode & S_IFMT) == S_IFDIR;
     #endif
     return exists;
   }
 
+
+
   namespace {
 
-    /// @todo Tidy up
+    /// gamma functions from Cephes library -- http://www.netlib.org/cephes
+    /// Copyright 1985, 1987, 2000 by Stephen L. Moshier
+
     static const double kMACHEP = 1.11022302462515654042363166809e-16;
     static const double kMAXLOG = 709.782712893383973096206318587;
     static const double kBig = 4.503599627370496e15;
     static const double kBiginv =  2.22044604925031308085e-16;
 
     double igamc(double a, double x);
     double igam(double a, double x);
 
-    /// @name gamma functions from Cephes library -- http://www.netlib.org/cephes
-    ///
-    /// Copyright 1985, 1987, 2000 by Stephen L. Moshier
-    //@{
 
     /// @brief Incomplete gamma function (complement integral)
     ///
     /// \f$ \gamma_c(a,x) = 1 - \gamma(a,x) \f$
     /// \f$ \gamma_c(a,x) = 1/\Gamma(a) \int_x^\inf e^-t t^(a-1) dt \f$
     ///
     /// In this implementation both arguments must be positive.
     /// The integral is evaluated by either a power series or
     /// continued fraction expansion, depending on the relative
     /// values of a and x.
     double igamc(double a, double x) {
       double ans, ax, c, yc, r, t, y, z;
       double pk, pkm1, pkm2, qk, qkm1, qkm2;
 
       // LM: for negative values returns 0.0
       // This is correct if a is a negative integer since Gamma(-n) = +/- inf
       if (a <= 0)  return 0.0;
 
       if (x <= 0) return 1.0;
 
       if( (x < 1.0) || (x < a) )
         return( 1.0 - igam(a,x) );
 
       ax = a * log(x) - x - lgamma(a);
       if( ax < -kMAXLOG )
         return( 0.0 );
 
       ax = exp(ax);
 
       // Continued fraction
       y = 1.0 - a;
       z = x + y + 1.0;
       c = 0.0;
       pkm2 = 1.0;
       qkm2 = x;
       pkm1 = x + 1.0;
       qkm1 = z * x;
       ans = pkm1/qkm1;
 
       do
         {
           c += 1.0;
           y += 1.0;
           z += 2.0;
           yc = y * c;
           pk = pkm1 * z  -  pkm2 * yc;
           qk = qkm1 * z  -  qkm2 * yc;
           if(qk)
             {
               r = pk/qk;
               t = fabs( (ans - r)/r );
               ans = r;
             }
           else
             t = 1.0;
           pkm2 = pkm1;
           pkm1 = pk;
           qkm2 = qkm1;
           qkm1 = qk;
           if( fabs(pk) > kBig )
             {
               pkm2 *= kBiginv;
               pkm1 *= kBiginv;
               qkm2 *= kBiginv;
               qkm1 *= kBiginv;
             }
         }
       while (t > kMACHEP);
 
       return ans*ax;
     }
 
 
     /// @brief Left tail of incomplete gamma function
     ///
     /// \f$ \gamma(a,x) = x^a e^-x \sum_k=0^\inf x^k / \Gamma(a+k+1) \f$
     double igam( double a, double x )
     {
       double ans, ax, c, r;
 
       // LM: for negative values returns 1.0 instead of zero
       // This is correct if a is a negative integer since Gamma(-n) = +/- inf
       if (a <= 0)  return 1.0;
 
       if (x <= 0)  return 0.0;
 
       if( (x > 1.0) && (x > a ) )
         return( 1.0 - igamc(a,x) );
 
       /* Compute  x**a * exp(-x) / gamma(a)  */
       ax = a * log(x) - x - lgamma(a);
       if( ax < -kMAXLOG )
         return( 0.0 );
 
       ax = exp(ax);
 
       /* power series */
       r = a;
       c = 1.0;
       ans = 1.0;
 
       do
         {
           r += 1.0;
           c *= x/r;
           ans += c;
         }
       while( c/ans > kMACHEP );
 
       return( ans * ax/a );
     }
 
-    //@}
-
   }
 
 
-
   /// @brief Compute quantiles for standard normal distribution N(0, 1) at probability p
   ///
   /// ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.
   double norm_quantile(double p) {
 
     /// @todo Return +-inf
     if (p <=0 || p >= 1) {
       cerr << "norm_quantile: probability outside (0, 1)" << endl;
       return 0;
     }
 
     const double  a0 = 3.3871328727963666080e0;
     const double  a1 = 1.3314166789178437745e+2;
     const double  a2 = 1.9715909503065514427e+3;
     const double  a3 = 1.3731693765509461125e+4;
     const double  a4 = 4.5921953931549871457e+4;
     const double  a5 = 6.7265770927008700853e+4;
     const double  a6 = 3.3430575583588128105e+4;
     const double  a7 = 2.5090809287301226727e+3;
     const double  b1 = 4.2313330701600911252e+1;
     const double  b2 = 6.8718700749205790830e+2;
     const double  b3 = 5.3941960214247511077e+3;
     const double  b4 = 2.1213794301586595867e+4;
     const double  b5 = 3.9307895800092710610e+4;
     const double  b6 = 2.8729085735721942674e+4;
     const double  b7 = 5.2264952788528545610e+3;
     const double  c0 = 1.42343711074968357734e0;
     const double  c1 = 4.63033784615654529590e0;
     const double  c2 = 5.76949722146069140550e0;
     const double  c3 = 3.64784832476320460504e0;
     const double  c4 = 1.27045825245236838258e0;
     const double  c5 = 2.41780725177450611770e-1;
     const double  c6 = 2.27238449892691845833e-2;
     const double  c7 = 7.74545014278341407640e-4;
     const double  d1 = 2.05319162663775882187e0;
     const double  d2 = 1.67638483018380384940e0;
     const double  d3 = 6.89767334985100004550e-1;
     const double  d4 = 1.48103976427480074590e-1;
     const double  d5 = 1.51986665636164571966e-2;
     const double  d6 = 5.47593808499534494600e-4;
     const double  d7 = 1.05075007164441684324e-9;
     const double  e0 = 6.65790464350110377720e0;
     const double  e1 = 5.46378491116411436990e0;
     const double  e2 = 1.78482653991729133580e0;
     const double  e3 = 2.96560571828504891230e-1;
     const double  e4 = 2.65321895265761230930e-2;
     const double  e5 = 1.24266094738807843860e-3;
     const double  e6 = 2.71155556874348757815e-5;
     const double  e7 = 2.01033439929228813265e-7;
     const double  f1 = 5.99832206555887937690e-1;
     const double  f2 = 1.36929880922735805310e-1;
     const double  f3 = 1.48753612908506148525e-2;
     const double  f4 = 7.86869131145613259100e-4;
     const double  f5 = 1.84631831751005468180e-5;
     const double  f6 = 1.42151175831644588870e-7;
     const double  f7 = 2.04426310338993978564e-15;
 
     const double split1 = 0.425;
     const double split2=5.;
     const double konst1=0.180625;
     const double konst2=1.6;
 
     double q, r, quantile;
     q = p - 0.5;
     if (fabs(q) < split1) {
       r = konst1 - q*q;
       quantile = q* (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
                        * r + a2) * r + a1) * r + a0) /
         (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
            * r + b2) * r + b1) * r + 1.);
     } else {
       r = (q < 0) ? p : 1-p;
       //error case
       if (r<=0)
         quantile=0;
       else {
         r=sqrt(-log(r));
         if (r<=split2) {
           r=r-konst2;
           quantile=(((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
                       * r + c2) * r + c1) * r + c0) /
             (((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
                * r + d2) * r + d1) * r + 1);
         } else{
           r=r-split2;
           quantile=(((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
                       * r + e2) * r + e1) * r + e0) /
             (((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
                * r + f2) * r + f1) * r + 1);
         }
         if (q<0) quantile=-quantile;
       }
     }
     return quantile;
   }
 
 
   /// @brief Compute quantiles of the chi-squared probability distribution function
   ///
   /// Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35
   /// implemented by Anna Kreshuk.
   /// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
   /// Parameters:
   ///   @arg p   - the probability value, at which the quantile is computed
   ///   @arg ndf - number of degrees of freedom
   double chisquared_quantile(double p, double ndf) {
     static const double c[] = {0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
                                4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
                                84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
                                210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
                                462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
                                932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
                                2520.0, 5040.0};
     static const double e = 5e-7;
     static const double aa = 0.6931471806;
     static const int maxit = 20;
 
     /// @todo Tidy
     double ch, p1, p2, q, t, a, b, x;
     double s1, s2, s3, s4, s5, s6;
 
     if (ndf <= 0) return 0;
 
     const double g = lgamma(0.5*ndf);
     const double xx = 0.5 * ndf;
     const double cp = xx - 1;
 
     if (ndf >= log(p)*(-c[5])){
       // Starting approximation for ndf less than or equal to 0.32
       if (ndf > c[3]) {
         x = norm_quantile(p);
         // Starting approximation using Wilson and Hilferty estimate
         p1 = c[2]/ndf;
         ch = ndf*pow((x*sqrt(p1) + 1 - p1), 3);
         if (ch > c[6]*ndf + 6)
           ch = -2 * (log(1-p) - cp * log(0.5 * ch) + g);
       } else {
         ch = c[4];
         a = log(1-p);
         do {
           q = ch;
           p1 = 1 + ch * (c[7]+ch);
           p2 = ch * (c[9] + ch * (c[8] + ch));
           t = -0.5 + (c[7] + 2 * ch) / p1 - (c[9] + ch * (c[10] + 3 * ch)) / p2;
           ch = ch - (1 - exp(a + g + 0.5 * ch + cp * aa) *p2 / p1) / t;
         } while (fabs(q/ch - 1) > c[1]);
       }
     } else {
       ch = pow((p * xx * exp(g + xx * aa)),(1./xx));
       if (ch < e) return ch;
     }
     // Call to algorithm AS 239 and calculation of seven term Taylor series
     for (int i = 0; i < maxit; ++i) {
       q = ch;
       p1 = 0.5 * ch;
       p2 = p - igam(xx, p1);
 
       t = p2 * exp(xx * aa + g + p1 - cp * log(ch));
       b = t / ch;
       a = 0.5 * t - b * cp;
       s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] +c[11] * a))))) / c[24];
       s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
       s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
       s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
       s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
       s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
       ch = ch + t * (1 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
       if (fabs(q / ch - 1) > e) break;
     }
     return ch;
   }
 
 
 }