LHAPDF  6.5.5
AlphaS.h
1 // -*- C++ -*-
2 //
3 // This file is part of LHAPDF
4 // Copyright (C) 2012-2024 The LHAPDF collaboration (see AUTHORS for details)
5 //
6 #pragma once
7 #ifndef LHAPDF_AlphaS_H
8 #define LHAPDF_AlphaS_H
9 
10 #include "LHAPDF/Utils.h"
11 #include "LHAPDF/Exceptions.h"
12 #include "LHAPDF/KnotArray.h"
13 
14 namespace LHAPDF {
15 
16 
17  /// @brief Calculator interface for computing alpha_s(Q2) in various ways
18  ///
19  /// The design of the AlphaS classes is that they are substitutible
20  /// (cf. polymorphism) and are entirely non-dependent on the PDF and Info
21  /// objects: hence they can be used by external code that actually doesn't
22  /// want to do anything at all with PDFs, but which just wants to do some
23  /// alpha_s interpolation.
24  class AlphaS {
25  public:
26 
27  /// Base class constructor for default param setup
28  AlphaS();
29 
30  /// Destructor
31  virtual ~AlphaS() {};
32 
33  /// @name alpha_s values
34  ///@{
35 
36  /// Calculate alphaS(Q)
37  double alphasQ(double q) const { return alphasQ2(q*q); }
38 
39  /// Calculate alphaS(Q2)
40  /// @todo Throw error in this base method if Q < Lambda?
41  virtual double alphasQ2(double q2) const = 0;
42 
43  ///@}
44 
45 
46  /// @name alpha_s metadata
47  ///@{
48 
49  /// Calculate the number of active flavours at energy scale Q
50  int numFlavorsQ(double q) const { return numFlavorsQ2(q*q); }
51 
52  /// Calculate the number of active flavours at energy scale Q2
53  virtual int numFlavorsQ2(double q2) const;
54 
55  /// Get a quark mass by PDG code
56  double quarkMass(int id) const;
57 
58  /// @brief Set quark masses by PDG code
59  ///
60  /// Used in the analytic and ODE solvers.
61  void setQuarkMass(int id, double value);
62 
63  /// @brief Get a flavor scale threshold by PDG code
64  ///
65  /// Used in the analytic and ODE solvers.
66  double quarkThreshold(int id) const;
67 
68  /// @brief Set a flavor threshold by PDG code (= quark masses by default)
69  ///
70  /// Used in the analytic and ODE solvers.
71  void setQuarkThreshold(int id, double value);
72 
73  /// Get the order of QCD (expressed as number of loops)
74  ///
75  /// Used explicitly in the analytic and ODE solvers.
76  int orderQCD() { return _qcdorder; }
77 
78  /// @brief Set the order of QCD (expressed as number of loops)
79  ///
80  /// Used in the analytic and ODE solvers.
81  void setOrderQCD(int order) { _qcdorder = order; }
82 
83  /// @brief Set the Z mass used in this alpha_s
84  ///
85  /// Used in the ODE solver.
86  void setMZ(double mz) { _mz = mz; }
87 
88  /// @brief Set the alpha_s(MZ) used in this alpha_s
89  ///
90  /// Used in the ODE solver.
91  void setAlphaSMZ(double alphas) { _alphas_mz = alphas; }
92 
93  /// @brief Set the Z mass used in this alpha_s
94  ///
95  /// Used in the ODE solver.
96  void setMassReference(double mref) { _mreference = mref; _customref = true; }
97 
98  /// @brief Set the alpha_s(MZ) used in this alpha_s
99  ///
100  /// Used in the ODE solver.
101  void setAlphaSReference(double alphas) { _alphas_reference = alphas; _customref = true; }
102 
103  /// @brief Set the @a {i}th Lambda constant for @a i active flavors
104  ///
105  /// Used in the analytic solver.
106  virtual void setLambda(unsigned int, double) {};
107 
108  /// Enum of flavor schemes
109  enum FlavorScheme { FIXED, VARIABLE };
110 
111  /// Get the implementation type of this AlphaS
112  virtual std::string type() const = 0;
113 
114  /// Set flavor scheme of alpha_s solver
115  void setFlavorScheme(FlavorScheme scheme, int nf = -1);
116 
117  /// Get flavor scheme
119 
120  ///@}
121 
122 
123  protected:
124 
125  /// @name Calculating beta function values
126  ///@{
127 
128  /// Calculate the i'th beta function given the number of active flavours
129  /// Currently limited to 0 <= i <= 3
130  /// Calculated using the MSbar scheme
131  double _beta(int i, int nf) const;
132 
133  /// Calculate a vector of beta functions given the number of active flavours
134  /// Currently returns a 4-element vector of beta0 -- beta3
135  std::vector<double> _betas(int nf) const;
136 
137  ///@}
138 
139 
140  protected:
141 
142  /// Order of QCD evolution (expressed as number of loops)
144 
145  /// Mass of the Z boson in GeV
146  double _mz;
147 
148  /// Value of alpha_s(MZ)
149  double _alphas_mz;
150 
151  /// Reference mass in GeV
152  double _mreference;
153 
154  /// Value of alpha_s(reference mass)
156 
157  /// Decides whether to use custom reference values or fall back on MZ/AlphaS_MZ
159 
160  /// Masses of quarks in GeV
161  /// Used for working out flavor thresholds and the number of quarks that are
162  /// active at energy scale Q.
163  std::map<int, double> _quarkmasses, _flavorthresholds;
164 
165  /// The flavor scheme in use
167 
168  /// The allowed numbers of flavours in a fixed scheme
169  int _fixflav;
170 
171  };
172 
173 
174 
175  /// Calculate alpha_s(Q2) by an analytic approximation
176  class AlphaS_Analytic : public AlphaS {
177  public:
178 
179  /// Implementation type of this solver
180  std::string type() const { return "analytic"; }
181 
182  /// Calculate alphaS(Q2)
183  double alphasQ2(double q2) const;
184 
185  /// Analytic has its own numFlavorsQ2 which respects the min/max nf set by the Lambdas
186  int numFlavorsQ2(double q2) const;
187 
188  /// Set lambda_i (for i = flavour number)
189  void setLambda(unsigned int i, double lambda);
190 
191 
192  private:
193 
194  /// Get lambdaQCD for nf
195  double _lambdaQCD(int nf) const;
196 
197  /// Recalculate min/max flavors in case lambdas have changed
198  void _setFlavors();
199 
200 
201  /// LambdaQCD values.
202  std::map<int, double> _lambdas;
203 
204  /// Max number of flavors
205  int _nfmax;
206  /// Min number of flavors
207  int _nfmin;
208 
209  };
210 
211 
212 
213  /// Interpolate alpha_s from tabulated points in Q2 via metadata
214  ///
215  /// @todo Extrapolation: log-gradient xpol at low Q, const at high Q?
216  class AlphaS_Ipol : public AlphaS {
217  public:
218 
219  /// Implementation type of this solver
220  std::string type() const { return "ipol"; }
221 
222  /// Calculate alphaS(Q2)
223  double alphasQ2(double q2) const;
224 
225  /// Set the array of Q values for interpolation
226  ///
227  /// Writes to the same internal arrays as setQ2Values, appropriately transformed.
228  void setQValues(const std::vector<double>& qs);
229 
230  /// Set the array of Q2 values for interpolation
231  ///
232  /// Subgrids are represented by repeating the values which are the end of
233  /// one subgrid and the start of the next. The supplied vector must match
234  /// the layout of alpha_s values.
235  void setQ2Values(const std::vector<double>& q2s) { _q2s = q2s; }
236 
237  /// Set the array of alpha_s(Q2) values for interpolation
238  ///
239  /// The supplied vector must match the layout of Q2 knots. Subgrids may
240  /// have discontinuities, i.e. different alpha_s values on either side of a
241  /// subgrid boundary (for the same Q values).
242  void setAlphaSValues(const std::vector<double>& as) { _as = as; }
243 
244 
245  private:
246 
247  /// Standard cubic interpolation formula
248  double _interpolateCubic(double T, double VL, double VDL, double VH, double VDH) const;
249  /// Get the gradient for a patch in the middle of the grid
250  double _ddq_central( size_t i ) const;
251  /// Get the gradient for a patch at the low end of the grid
252  double _ddq_forward( size_t i ) const;
253  /// Get the gradient for a patch at the high end of the grid
254  double _ddq_backward( size_t i ) const;
255 
256  /// Synchronise the contents of the single Q2 / alpha_s vectors into subgrid objects
257  /// @note This is const so it can be called silently from a const method
258  void _setup_grids() const;
259 
260 
261  /// Map of AlphaSArrays "binned" for lookup by low edge in (log)Q2
262  /// @note This is mutable so it can be initialized silently from a const method
263  mutable std::map<double, AlphaSArray> _knotarrays;
264 
265  /// Array of ipol knots in Q2
266  std::vector<double> _q2s;
267  /// Array of alpha_s values for the Q2 knots
268  std::vector<double> _as;
269 
270  };
271 
272 
273 
274  /// Solve the differential equation in alphaS using an implementation of RK4
275  class AlphaS_ODE : public AlphaS {
276  public:
277 
278  /// Implementation type of this solver
279  std::string type() const { return "ode"; }
280 
281  /// Calculate alphaS(Q2)
282  double alphasQ2( double q2 ) const;
283 
284  /// Set MZ, and also the caching flag
285  void setMZ( double mz ) { _mz = mz; _calculated = false; }
286 
287  /// Set alpha_s(MZ), and also the caching flag
288  void setAlphaSMZ( double alphas ) { _alphas_mz = alphas; _calculated = false; }
289 
290  /// Set reference mass, and also the caching flag
291  void setMassReference( double mref ) { _mreference = mref; _calculated = false; _customref = true; }
292 
293  /// Set alpha_s(MZ), and also the caching flag
294  void setAlphaSReference( double alphas ) { _alphas_reference = alphas; _calculated = false; _customref = true; }
295 
296  /// Set the array of Q values for interpolation, and also the caching flag
297  void setQValues(const std::vector<double>& qs);
298 
299  /// @brief Set the array of Q2 values for interpolation, and also the caching flag
300  ///
301  /// Writes to the same internal array as setQValues, appropriately transformed.
302  void setQ2Values( std::vector<double> q2s ) { _q2s = q2s; _calculated = false; }
303 
304 
305  private:
306 
307  /// Calculate the derivative at Q2 = t, alpha_S = y
308  double _derivative(double t, double y, const std::vector<double>& beta) const;
309 
310  /// Calculate the decoupling relation when going from num. flav. = ni -> nf
311  /// abs(ni - nf) must be = 1
312  double _decouple(double y, double t, unsigned int ni, unsigned int nf) const;
313 
314  /// Calculate the next step using RK4 with adaptive step size
315  void _rk4(double& t, double& y, double h, const double allowed_change, const vector<double>& bs) const;
316 
317  /// Solve alpha_s for q2 using RK4
318  void _solve(double q2, double& t, double& y, const double& allowed_relative, double h, double accuracy) const;
319 
320  /// Create interpolation grid
321  void _interpolate() const;
322 
323 
324  /// Vector of Q2s in case specific anchor points are used
325  mutable std::vector<double> _q2s;
326 
327  /// Whether or not the ODE has been solved yet
328  mutable bool _calculated;
329 
330  /// The interpolation used to get Alpha_s after the ODE has been solved
332 
333  };
334 
335 
336 }
337 #endif