AlphaS.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2023 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 {
 
 
  class AlphaS {
  public:
 
    AlphaS();
 
    virtual ~AlphaS() {};
 
 
    double alphasQ(double q) const { return alphasQ2(q*q); }
 
    virtual double alphasQ2(double q2) const = 0;
 
 
 
 
    int numFlavorsQ(double q) const { return numFlavorsQ2(q*q); }
 
    virtual int numFlavorsQ2(double q2) const;
 
    double quarkMass(int id) const;
 
    void setQuarkMass(int id, double value);
 
    double quarkThreshold(int id) const;
 
    void setQuarkThreshold(int id, double value);
 
    int orderQCD() { return _qcdorder; }
 
    void setOrderQCD(int order) { _qcdorder = order; }
 
    void setMZ(double mz) { _mz = mz; }
 
    void setAlphaSMZ(double alphas) { _alphas_mz = alphas; }
 
    void setMassReference(double mref) { _mreference = mref; _customref = true; }
 
    void setAlphaSReference(double alphas) { _alphas_reference = alphas; _customref = true; }
 
    virtual void setLambda(unsigned int, double) {};
 
    enum FlavorScheme { FIXED, VARIABLE };
 
    virtual std::string type() const = 0;
 
    void setFlavorScheme(FlavorScheme scheme, int nf = -1);
 
    FlavorScheme flavorScheme() const { return _flavorscheme; }
 
 
 
  protected:
 
 
    double _beta(int i, int nf) const;
 
    std::vector<double> _betas(int nf) const;
 
 
 
  protected:
 
    int _qcdorder;
 
    double _mz;
 
    double _alphas_mz;
 
    double _mreference;
 
    double _alphas_reference;
 
    bool _customref;
 
    std::map<int, double> _quarkmasses, _flavorthresholds;
 
    FlavorScheme _flavorscheme;
 
    int _fixflav;
 
  };
 
 
 
  class AlphaS_Analytic : public AlphaS {
  public:
 
    std::string type() const { return "analytic"; }
 
    double alphasQ2(double q2) const;
 
    int numFlavorsQ2(double q2) const;
 
    void setLambda(unsigned int i, double lambda);
 
 
  private:
 
    double _lambdaQCD(int nf) const;
 
    void _setFlavors();
 
 
    std::map<int, double> _lambdas;
 
    int _nfmax;
    int _nfmin;
 
  };
 
 
 
  class AlphaS_Ipol : public AlphaS {
  public:
 
    std::string type() const { return "ipol"; }
 
    double alphasQ2(double q2) const;
 
    void setQValues(const std::vector<double>& qs);
 
    void setQ2Values(const std::vector<double>& q2s) { _q2s = q2s; }
 
    void setAlphaSValues(const std::vector<double>& as) { _as = as; }
 
 
  private:
 
    double _interpolateCubic(double T, double VL, double VDL, double VH, double VDH) const;
    double _ddq_central( size_t i ) const;
    double _ddq_forward( size_t i ) const;
    double _ddq_backward( size_t i ) const;
 
    void _setup_grids() const;
 
 
    mutable std::map<double, AlphaSArray> _knotarrays;
 
    std::vector<double> _q2s;
    std::vector<double> _as;
 
  };
 
 
 
  class AlphaS_ODE : public AlphaS {
  public:
 
    std::string type() const { return "ode"; }
 
    double alphasQ2( double q2 ) const;
 
    void setMZ( double mz ) { _mz = mz; _calculated = false; }
 
    void setAlphaSMZ( double alphas ) { _alphas_mz = alphas; _calculated = false; }
 
    void setMassReference( double mref ) { _mreference = mref; _calculated = false; _customref = true; }
 
    void setAlphaSReference( double alphas ) { _alphas_reference = alphas; _calculated = false; _customref = true; }
 
    void setQValues(const std::vector<double>& qs);
 
    void setQ2Values( std::vector<double> q2s ) { _q2s = q2s; _calculated = false; }
 
 
  private:
 
    double _derivative(double t, double y, const std::vector<double>& beta) const;
 
    double _decouple(double y, double t, unsigned int ni, unsigned int nf) const;
 
    void _rk4(double& t, double& y, double h, const double allowed_change, const vector<double>& bs) const;
 
    void _solve(double q2, double& t, double& y, const double& allowed_relative, double h, double accuracy) const;
 
    void _interpolate() const;
 
 
    mutable std::vector<double> _q2s;
 
    mutable bool _calculated;
 
    mutable AlphaS_Ipol _ipol;
 
  };
 
 
}
#endif