correcting alphasDown with alphasDown in MatchPDFs.f + adding beta fu…

…nction terms + updating documentation

doc/sources/QCD_QED_common_basis.tex

@@ -57,7 +57,7 @@
 
 \newpage{}
 
-\section{The Structure of the DGLAP Equation}
+\section{The structure of the DGLAP equation}
 
 The DGLAP equation that governs the PDF evolution has a general
 structure that in QCD holds at any perturbative order. Suppose one
@@ -637,7 +637,7 @@ \section{The Structure of the DGLAP Equation}
 \end{equation}
 
 
-\section{Evolution Basis}
+\section{Evolution basis}
 
 In order to diagonlize as much as possible the evolution matrix in the
 presence of QED corrections avoiding unnecessery couplings between parton distributions, we propose
@@ -908,7 +908,7 @@ \section{QED corrections at LO}
 $2e_\Sigma^2=\theta^-$ and $2\delta_e^2=\theta^+$, is consistent with
 eq.~(9) of the APFEL paper.
 
-\section{Including the Lepton PDFs}
+\section{Including the lepton PDFs}
 
 In order to include the lepton PDFs in the coupled QCD$\times$QED
 DGLAP evolution we first need to make some preliminary
@@ -1408,7 +1408,7 @@ \section{Including the Lepton PDFs}
 \end{equation}
 
 
-\subsection{Evolution Equations at LO}
+\subsection{Evolution equations at LO}
 
 At LO in QED, considering that $\mathcal{P}_{q\ell}^S = 0$ and that $\mathcal{P}_{ij}=\widetilde{P}_{ij}$, the
 evolution equations above reduce to:
@@ -1529,8 +1529,79 @@ \subsection{Evolution Equations at LO}
 that we need to introduce a new threshold between the charm and the
 bottom thresholds and this will complicate the structure of the code.
 
-\bibliographystyle{unsrt}
-\phantomsection\addcontentsline{toc}{section}{\refname}\nocite{*}
-\bibliography{bibliography}
+\section{QED corrections at NLO}
+
+In this section we discuss the details of the implementation of the
+NLO QED corrections.  While the inclusion of the LO corrections
+presents many simplifications, $e.g.$ QED and QCD corrections do not
+mix and thus the DGLAP equations as well as the $\alpha_s$ and
+$\alpha$ evolution equantions are decoupled, when including NLO
+corrections QED and QCD corrections mix both in the DGLAP and in the
+coupling evolution equations. In addition, as far as DIS structure
+functions are concerned, such corrections induce photon initiated
+diagrams that have to be included in order to have the full set of
+corrections. In the following we will first discuss how to generalize
+the coupling evolution equations, we will then turn to the DGLAP, and
+finally we will consider the DIS structure functions.
+
+\subsection{Evolution equations for the couplings}
+
+As already mentioned, NLO QED corrections induce a mixing with QCD.
+At the level of the couplings, this essentially means that the QCD
+$\beta$-function will get corrections proportional to $\alpha$ and
+vice-versa, the QED $\beta$-function will get corrections proportional
+to $\alpha_s$, that is:
+\begin{equation}
+\begin{array}{rcl}
+\displaystyle \mu^2\frac{\partial \alpha_s}{\partial \mu^2} &=& \displaystyle
+                                                \beta^{\rm QCD}(\alpha_s,\alpha)\,.\\
+\\
+\displaystyle \mu^2\frac{\partial \alpha}{\partial \mu^2} &=& \displaystyle \beta^{\rm QED}(\alpha_s,\alpha)\,.
+\end{array}
+\end{equation}
+As a consequence, these evolution equations form a set of coupled
+differential equations.
+
+Before discussing the numerical solution of these equations, we first
+need to know explicitly the new contributions to the
+$\beta$-functions. In particular, up to NLO one has: 
+\begin{equation}
+\beta^{\rm QCD}(\alpha_s,\alpha) = -\alpha_s\left[\beta_0^{(\alpha_s)}\left(\frac{\alpha_s}{4\pi}\right)+\beta_1^{(\alpha_s\alpha)}\left(\frac{\alpha_s}{4\pi}\right) \left(\frac{\alpha}{4\pi}\right)+\beta_1^{(\alpha_s^2)}\left(\frac{\alpha_s}{4\pi}\right)^2+\dots\right]\,,
+\end{equation}
+and:
+\begin{equation}
+\beta^{\rm QED}(\alpha_s,\alpha) = -\alpha\left[\beta_0^{(\alpha)}\left(\frac{\alpha}{4\pi}\right)+\beta_1^{(\alpha\alpha_s)}\left(\frac{\alpha}{4\pi}\right) \left(\frac{\alpha_s}{4\pi}\right)+\beta_1^{(\alpha^2)}\left(\frac{\alpha}{4\pi}\right)^2+\dots\right]\,,
+\end{equation}
+where the new terms $\beta_1^{(\alpha_s\alpha)}$ and
+$\beta_1^{(\alpha\alpha_s)}$ can be taken from
+Ref.~\cite{Surguladze:1996hx} and, taking into account the additional
+factor four in the definition of the expansion parameter and writing
+explicitly the color factors one finds:
+\begin{equation}
+\beta_1^{(\alpha_s\alpha)} = -2\sum_{i=1}^{n_f}
+q_i^2\quad\mbox{and}\quad \beta_1^{(\alpha\alpha_s)} = -\frac{16}{3}N_c\sum_{i=1}^{n_f} q_i^2\,.
+\end{equation}
+where $n_f$ is the number of active quark flavors and $N_c=3$ is the
+number of colors.  We also need the term $\beta_1^{(\alpha^2)}$ which
+can again be taken from the same paper:
+\begin{equation}
+\beta_1^{(\alpha^2)} = -4\left(n_l+N_c\sum_{i=1}^{n_f} q_i^2\right)\,.
+\end{equation}
+where $n_l$ is the number of active lepton flavors.
+
+
+
+\begin{thebibliography}{alp}
+
+%\cite{Surguladze:1996hx}
+\bibitem{Surguladze:1996hx}
+  L.~R.~Surguladze,
+  %``O(alpha**n alpha-s**m) corrections in e+ e- annihilation and tau decay,''
+  hep-ph/9803211.
+  %%CITATION = HEP-PH/9803211;%%
+  %3 citations counted in INSPIRE as of 29 Apr 2016
+
+\end{thebibliography}
+
 
 \end{document}

src/Evolution/MatchPDFs.f

@@ -33,10 +33,10 @@ subroutine MatchPDFs(nf,sgn,fevQCD)
 *
       double precision fevQCD(0:13,0:nint_max)
 *
-*     Get alphas value at the heavy quark threshold (with nf active flavours)
+*     Get alphas value at the heavy quark threshold (with nf+1 active flavours)
 *
-c      coup = asthUp(nf)
-      coup = asthDown(nf)
+      coup = asthUp(nf)
+c      coup = asthDown(nf)
 *
 *     Singlet map
 *

src/Evolution/a_QCD.f

@@ -339,7 +339,7 @@ function fbeta(a,nf,ipt)
       elseif(ipt.eq.2)then
          fbeta = - a**2 * ( beta0apf(nf)
      1           + a * ( beta1apf(nf) + a * beta2apf(nf) ) )
-      elseif(ipt.eq.3)then
+      elseif(ipt.ge.3)then
          fbeta = - a**2 * ( beta0apf(nf)
      1           + a * ( beta1apf(nf)
      2           + a * ( beta2apf(nf) + a * beta3apf(nf) ) ) )

src/Evolution/a_QED.f

@@ -206,7 +206,7 @@ FUNCTION ALPHAQEDEV(NF,NL,Q2F,Q2I,ALPHAREF)
       BETA0 = BETA0QED(NF,NL)
       L = DLOG( Q2F / Q2I )
 *
-      ALPHAQEDEV = ALPHAREF / ( 1D0 - ALPHAREF * BETA0 * L )
+      ALPHAQEDEV = ALPHAREF / ( 1D0 + ALPHAREF * BETA0 * L )
 *
       RETURN
       END
@@ -240,7 +240,105 @@ function beta0qed(nf,nl)
       sumch2(5) = 11d0 / 9d0
       sumch2(6) = 5d0 / 3d0
 *
-      beta0qed = 8d0 / 3d0 * ( nc * sumch2(nf) + nl )
+      beta0qed = - 8d0 / 3d0 * ( nc * sumch2(nf) + nl )
+c      beta0qed = - 4d0 / 3d0 * ( nc * sumch2(nf) + nl )
+*
+      return
+      end
+*
+****************************************************************
+      function beta1qed(nf,nl)
+*
+      implicit none
+**
+*     Input Variables
+*
+      integer nf,nl
+**
+*     Internal Variables
+*
+      integer nc
+      double precision sumch4(3:6)
+**
+*     Output Variables
+*
+      double precision beta1qed
+*
+*     Number of colours
+*
+      nc = 3
+*
+*     Sum of the first 3, 4, 5 and 6 electric charges to the fouth
+*
+      sumch4(3) = 18d0 / 81d0
+      sumch4(4) = 34d0 / 81d0
+      sumch4(5) = 35d0 / 81d0
+      sumch4(6) = 51d0 / 81d0
+*
+      beta1qed = - 16d0 / 4d0 * ( nc * sumch4(nf) + nl )
+*
+      return
+      end
+*
+****************************************************************
+      function beta1qcdqed(nf)
+*
+      implicit none
+**
+*     Input Variables
+*
+      integer nf
+**
+*     Internal Variables
+*
+      double precision sumch2(3:6)
+**
+*     Output Variables
+*
+      double precision beta1qcdqed
+*
+*     Sum of the first 3, 4, 5 and 6 squared electric charges
+*
+      sumch2(3) = 2d0 / 3d0
+      sumch2(4) = 10d0 / 9d0
+      sumch2(5) = 11d0 / 9d0
+      sumch2(6) = 5d0 / 3d0
+*
+      beta1qcdqed = - 2d0 * sumch2(nf)
+*
+      return
+      end
+*
+****************************************************************
+      function beta1qedqcd(nf)
+*
+      implicit none
+**
+*     Input Variables
+*
+      integer nf
+**
+*     Internal Variables
+*
+      integer nc
+      double precision sumch2(3:6)
+**
+*     Output Variables
+*
+      double precision beta1qedqcd
+*
+*     Number of colours
+*
+      nc = 3
+*
+*     Sum of the first 3, 4, 5 and 6 squared electric charges
+*
+      sumch2(3) = 2d0 / 3d0
+      sumch2(4) = 10d0 / 9d0
+      sumch2(5) = 11d0 / 9d0
+      sumch2(6) = 5d0 / 3d0
+*
+      beta1qedqcd = - 16d0 / 3d0 * nc * sumch2(nf)
 *
       return
       end
@@ -276,10 +374,10 @@ function a_as_exact(nf,nl,as0,a0,as,ipt)
 *     Fourth-order runge-kutta beyond the leading order
 *
       do n=1,nstep
-         k1 = h * fbetaQED(a,nf,nl) / fbeta(as,nf,ipt)
-         k2 = h * fbetaQED(a+k1/2d0,nf,nl) / fbeta(as+h/2d0,nf,ipt)
-         k3 = h * fbetaQED(a+k2/2d0,nf,nl) / fbeta(as+h/2d0,nf,ipt)
-         k4 = h * fbetaQED(a+k3,nf,nl) / fbeta(as+h,nf,ipt)
+         k1 = h * fbetaQED(a,nf,nl,ipt) / fbeta(as,nf,ipt)
+         k2 = h * fbetaQED(a+k1/2d0,nf,nl,ipt) / fbeta(as+h/2d0,nf,ipt)
+         k3 = h * fbetaQED(a+k2/2d0,nf,nl,ipt) / fbeta(as+h/2d0,nf,ipt)
+         k4 = h * fbetaQED(a+k3,nf,nl,ipt) / fbeta(as+h,nf,ipt)
 *
          a  = a + ( k1 + 2d0 * k2 + 2d0 * k3 + k4 ) / 6d0
          as = as + h
@@ -295,24 +393,29 @@ function a_as_exact(nf,nl,as0,a0,as,ipt)
 *     QED beta function.
 *
 ****************************************************************************
-      function fbetaQED(a,nf,nl)
+      function fbetaQED(a,nf,nl,ipt)
 *
       implicit none
 **
 *     Input Variables
 *
       double precision a
-      integer nf,nl
+      integer nf,nl,ipt
 **
 *     Internal Variables
 *
-      double precision beta0qed
+      double precision beta0qed!,beta1qed
 **
 *     Output Variables
 *
       double precision fbetaQED
 *
-      fbetaQED = a**2 * beta0qed(nf,nl)
+      fbetaQED = - a**2 * beta0qed(nf,nl)
+c      if(ipt.eq.0)then
+c         fbetaQED = - a**2 * beta0qed(nf,nl)
+c      elseif(ipt.ge.1)then
+c         fbetaQED = - a**2 * ( beta0qed(nf,nl) + a * beta1qed(nf,nl) )
+c      endif
 *
       return
       end

src/Evolution/integralsQED.f

@@ -51,7 +51,7 @@ function integralsQED(alpha,beta,coup,kk)
       if(PDFEvol.eq."exactmu".or.Th.eq."QUniD")then
          integralsQED = coup * SQ(igrid,pnf,wnl,kk,alpha,beta)
       else
-         integralsQED = SQ(igrid,pnf,wnl,kk,alpha,beta) 
+         integralsQED = - SQ(igrid,pnf,wnl,kk,alpha,beta)
      1                / beta0qed(bnf,wnl) / coup
       endif
 *