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