diff --git a/gonsalves/Cteq61Pdf.f b/gonsalves/Cteq61Pdf.f
deleted file mode 100644
index 29a007d..0000000
--- a/gonsalves/Cteq61Pdf.f
+++ /dev/null
@@ -1,513 +0,0 @@
-C============================================================================
-C April 10, 2002, v6.01
-C February 23, 2003, v6.1
-C
-C Ref[1]: "New Generation of Parton Distributions with Uncertainties from Global QCD Analysis"
-C By: J. Pumplin, D.R. Stump, J.Huston, H.L. Lai, P. Nadolsky, W.K. Tung
-C JHEP 0207:012(2002), hep-ph/0201195
-C
-C Ref[2]: "Inclusive Jet Production, Parton Distributions, and the Search for New Physics"
-C By : D. Stump, J. Huston, J. Pumplin, W.K. Tung, H.L. Lai, S. Kuhlmann, J. Owens
-C hep-ph/0303013
-C
-C This package contains
-C (1) 4 standard sets of CTEQ6 PDF's (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) ;
-C (2) 40 up/down sets (with respect to CTEQ6M) for uncertainty studies from Ref[1];
-C (3) updated version of the above: CTEQ6.1M and its 40 up/down eigenvector sets from Ref[2].
-C
-C The CTEQ6.1M set provides a global fit that is almost equivalent in every respect
-C to the published CTEQ6M, Ref[1], although some parton distributions (e.g., the gluon)
-C may deviate from CTEQ6M in some kinematic ranges by amounts that are well within the
-C specified uncertainties.
-C The more significant improvements of the new version are associated with some of the
-C 40 eigenvector sets, which are made more symmetrical and reliable in (3), compared to (2).
-
-C Details about calling convention are:
-C ---------------------------------------------------------------------------
-C Iset PDF-set Description Alpha_s(Mz)**Lam4 Lam5 Table_File
-C ===========================================================================
-C Standard, "best-fit", sets:
-C --------------------------
-C 1 CTEQ6M Standard MSbar scheme 0.118 326 226 cteq6m.tbl
-C 2 CTEQ6D Standard DIS scheme 0.118 326 226 cteq6d.tbl
-C 3 CTEQ6L Leading Order 0.118** 326** 226 cteq6l.tbl
-C 4 CTEQ6L1 Leading Order 0.130** 215** 165 cteq6l1.tbl
-C ============================================================================
-C For uncertainty calculations using eigenvectors of the Hessian:
-C ---------------------------------------------------------------
-C central + 40 up/down sets along 20 eigenvector directions
-C -----------------------------
-C Original version, Ref[1]: central fit: CTEQ6M (=CTEQ6M.00)
-C -----------------------
-C 1xx CTEQ6M.xx +/- sets 0.118 326 226 cteq6m1xx.tbl
-C where xx = 01-40: 01/02 corresponds to +/- for the 1st eigenvector, ... etc.
-C e.g. 100 is CTEQ6M.00 (=CTEQ6M),
-C 101/102 are CTEQ6M.01/02, +/- sets of 1st eigenvector, ... etc.
-C -----------------------
-C Updated version, Ref[2]: central fit: CTEQ6.1M (=CTEQ61.00)
-C -----------------------
-C 2xx CTEQ61.xx +/- sets 0.118 326 226 ctq61.xx.tbl
-C where xx = 01-40: 01/02 corresponds to +/- for the 1st eigenvector, ... etc.
-C e.g. 200 is CTEQ61.00 (=CTEQ6.1M),
-C 201/202 are CTEQ61.01/02, +/- sets of 1st eigenvector, ... etc.
-C ===========================================================================
-C ** ALL fits are obtained by using the same coupling strength
-C \alpha_s(Mz)=0.118 and the NLO running \alpha_s formula, except CTEQ6L1
-C which uses the LO running \alpha_s and its value determined from the fit.
-C For the LO fits, the evolution of the PDF and the hard cross sections are
-C calculated at LO. More detailed discussions are given in the references.
-C
-C The table grids are generated for 10^-6 < x < 1 and 1.3 < Q < 10,000 (GeV).
-C PDF values outside of the above range are returned using extrapolation.
-C Lam5 (Lam4) represents Lambda value (in MeV) for 5 (4) flavors.
-C The matching alpha_s between 4 and 5 flavors takes place at Q=4.5 GeV,
-C which is defined as the bottom quark mass, whenever it can be applied.
-C
-C The Table_Files are assumed to be in the working directory.
-C
-C Before using the PDF, it is necessary to do the initialization by
-C Call SetCtq6(Iset)
-C where Iset is the desired PDF specified in the above table.
-C
-C The function Ctq6Pdf (Iparton, X, Q)
-C returns the parton distribution inside the proton for parton [Iparton]
-C at [X] Bjorken_X and scale [Q] (GeV) in PDF set [Iset].
-C Iparton is the parton label (5, 4, 3, 2, 1, 0, -1, ......, -5)
-C for (b, c, s, d, u, g, u_bar, ..., b_bar),
-C
-C For detailed information on the parameters used, e.q. quark masses,
-C QCD Lambda, ... etc., see info lines at the beginning of the
-C Table_Files.
-C
-C These programs, as provided, are in double precision. By removing the
-C "Implicit Double Precision" lines, they can also be run in single
-C precision.
-C
-C If you have detailed questions concerning these CTEQ6 distributions,
-C or if you find problems/bugs using this package, direct inquires to
-C Pumplin@pa.msu.edu or Tung@pa.msu.edu.
-C
-C===========================================================================
-
- Function Ctq6Pdf (Iparton, X, Q)
- Implicit Double Precision (A-H,O-Z)
- Logical Warn
- Common
- > / CtqPar2 / Nx, Nt, NfMx
- > / QCDtable / Alambda, Nfl, Iorder
-
- Data Warn /.true./
- save Warn
-
- If (X .lt. 0D0 .or. X .gt. 1D0) Then
- Print *, 'X out of range in Ctq6Pdf: ', X
- Stop
- Endif
- If (Q .lt. Alambda) Then
- Print *, 'Q out of range in Ctq6Pdf: ', Q
- Stop
- Endif
- If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
- If (Warn) Then
-C put a warning for calling extra flavor.
- Warn = .false.
- Print *, 'Warning: Iparton out of range in Ctq6Pdf! '
- Print *, 'Iparton, MxFlvN0: ', Iparton, NfMx
- Endif
- Ctq6Pdf = 0D0
- Return
- Endif
-
- Ctq6Pdf = PartonX6 (Iparton, X, Q)
- if(Ctq6Pdf.lt.0.D0) Ctq6Pdf = 0.D0
-
- Return
-
-C ********************
- End
-
- Subroutine SetCtq6 (Iset)
- Implicit Double Precision (A-H,O-Z)
- Parameter (Isetmax0=5)
- Character Flnm(Isetmax0)*6, nn*3, Tablefile*40
- Data (Flnm(I), I=1,Isetmax0)
- > / 'cteq6m', 'cteq6d', 'cteq6l', 'cteq6l','ctq61.'/
- Data Isetold, Isetmin0, Isetmin1, Isetmax1 /-987,1,100,140/
- Data Isetmin2,Isetmax2 /200,240/
- save
-
-C If data file not initialized, do so.
- If(Iset.ne.Isetold) then
- IU= NextUn()
- If (Iset.ge.Isetmin0 .and. Iset.le.3) Then
- Tablefile=Flnm(Iset)//'.tbl'
- Elseif (Iset.eq.4) Then
- Tablefile=Flnm(Iset)//'1.tbl'
- Elseif (Iset.ge.Isetmin1 .and. Iset.le.Isetmax1) Then
- write(nn,'(I3)') Iset
- Tablefile=Flnm(1)//nn//'.tbl'
- Elseif (Iset.ge.Isetmin2 .and. Iset.le.Isetmax2) Then
- write(nn,'(I3)') Iset
- Tablefile=Flnm(5)//nn(2:3)//'.tbl'
- Else
- Print *, 'Invalid Iset number in SetCtq6 :', Iset
- Stop
- Endif
- Open(IU, File=Tablefile, Status='OLD', Err=100)
- 21 Call ReadTbl (IU)
- Close (IU)
- Isetold=Iset
- Endif
- Return
-
- 100 Print *, ' Data file ', Tablefile, ' cannot be opened '
- >//'in SetCtq6!!'
- Stop
-C ********************
- End
-
- Subroutine ReadTbl (Nu)
- Implicit Double Precision (A-H,O-Z)
- Character Line*80
- PARAMETER (MXX = 96, MXQ = 20, MXF = 5)
- PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX)
- Common
- > / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX)
- > / CtqPar2 / Nx, Nt, NfMx
- > / XQrange / Qini, Qmax, Xmin
- > / QCDtable / Alambda, Nfl, Iorder
- > / Masstbl / Amass(6)
-
- Read (Nu, '(A)') Line
- Read (Nu, '(A)') Line
- Read (Nu, *) Dr, Fl, Al, (Amass(I),I=1,6)
- Iorder = Nint(Dr)
- Nfl = Nint(Fl)
- Alambda = Al
-
- Read (Nu, '(A)') Line
- Read (Nu, *) NX, NT, NfMx
-
- Read (Nu, '(A)') Line
- Read (Nu, *) QINI, QMAX, (TV(I), I =0, NT)
-
- Read (Nu, '(A)') Line
- Read (Nu, *) XMIN, (XV(I), I =0, NX)
-
- Do 11 Iq = 0, NT
- TV(Iq) = Log(Log (TV(Iq) /Al))
- 11 Continue
-C
-C Since quark = anti-quark for nfl>2 at this stage,
-C we Read out only the non-redundent data points
-C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
-
- Nblk = (NX+1) * (NT+1)
- Npts = Nblk * (NfMx+3)
- Read (Nu, '(A)') Line
- Read (Nu, *, IOSTAT=IRET) (UPD(I), I=1,Npts)
-
- Return
-C ****************************
- End
-
- Function NextUn()
-C Returns an unallocated FORTRAN i/o unit.
- Logical EX
-C
- Do 10 N = 10, 300
- INQUIRE (UNIT=N, OPENED=EX)
- If (.NOT. EX) then
- NextUn = N
- Return
- Endif
- 10 Continue
- Stop ' There is no available I/O unit. '
-C *************************
- End
-C
-
- SUBROUTINE POLINT (XA,YA,N,X,Y,DY)
-
- IMPLICIT DOUBLE PRECISION (A-H, O-Z)
-C Adapted from "Numerical Recipes"
- PARAMETER (NMAX=10)
- DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
- NS=1
- DIF=ABS(X-XA(1))
- DO 11 I=1,N
- DIFT=ABS(X-XA(I))
- IF (DIFT.LT.DIF) THEN
- NS=I
- DIF=DIFT
- ENDIF
- C(I)=YA(I)
- D(I)=YA(I)
-11 CONTINUE
- Y=YA(NS)
- NS=NS-1
- DO 13 M=1,N-1
- DO 12 I=1,N-M
- HO=XA(I)-X
- HP=XA(I+M)-X
- W=C(I+1)-D(I)
- DEN=HO-HP
- IF(DEN.EQ.0.)PAUSE
- DEN=W/DEN
- D(I)=HP*DEN
- C(I)=HO*DEN
-12 CONTINUE
- IF (2*NS.LT.N-M)THEN
- DY=C(NS+1)
- ELSE
- DY=D(NS)
- NS=NS-1
- ENDIF
- Y=Y+DY
-13 CONTINUE
- RETURN
- END
-
- Function PartonX6 (IPRTN, XX, QQ)
-
-c Given the parton distribution function in the array U in
-c COMMON / PEVLDT / , this routine interpolates to find
-c the parton distribution at an arbitray point in x and q.
-c
- Implicit Double Precision (A-H,O-Z)
-
- Parameter (MXX = 96, MXQ = 20, MXF = 5)
- Parameter (MXQX= MXQ * MXX, MXPQX = MXQX * (MXF+3))
-
- Common
- > / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX)
- > / CtqPar2 / Nx, Nt, NfMx
- > / XQrange / Qini, Qmax, Xmin
-
- Dimension fvec(4), fij(4)
- Dimension xvpow(0:mxx)
- Data OneP / 1.00001 /
- Data xpow / 0.3d0 / !**** choice of interpolation variable
- Data nqvec / 4 /
- Data ientry / 0 /
- Save ientry,xvpow
-
-c store the powers used for interpolation on first call...
- if(ientry .eq. 0) then
- ientry = 1
-
- xvpow(0) = 0D0
- do i = 1, nx
- xvpow(i) = xv(i)**xpow
- enddo
- endif
-
- X = XX
- Q = QQ
- tt = log(log(Q/Al))
-
-c ------------- find lower end of interval containing x, i.e.,
-c get jx such that xv(jx) .le. x .le. xv(jx+1)...
- JLx = -1
- JU = Nx+1
- 11 If (JU-JLx .GT. 1) Then
- JM = (JU+JLx) / 2
- If (X .Ge. XV(JM)) Then
- JLx = JM
- Else
- JU = JM
- Endif
- Goto 11
- Endif
-C Ix 0 1 2 Jx JLx Nx-2 Nx
-C |---|---|---|...|---|-x-|---|...|---|---|
-C x 0 Xmin x 1
-C
- If (JLx .LE. -1) Then
- Print '(A,1pE12.4)', 'Severe error: x <= 0 in PartonX6! x = ', x
- Stop
- ElseIf (JLx .Eq. 0) Then
- Jx = 0
- Elseif (JLx .LE. Nx-2) Then
-
-C For interrior points, keep x in the middle, as shown above
- Jx = JLx - 1
- Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
-
-C We tolerate a slight over-shoot of one (OneP=1.00001),
-C perhaps due to roundoff or whatever, but not more than that.
-C Keep at least 4 points >= Jx
- Jx = JLx - 2
- Else
- Print '(A,1pE12.4)', 'Severe error: x > 1 in PartonX6! x = ', x
- Stop
- Endif
-C ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
-
-C This is the variable to be interpolated in
- ss = x**xpow
-
- If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
-
-c initiation work for "interior bins": store the lattice points in s...
- svec1 = xvpow(jx)
- svec2 = xvpow(jx+1)
- svec3 = xvpow(jx+2)
- svec4 = xvpow(jx+3)
-
- s12 = svec1 - svec2
- s13 = svec1 - svec3
- s23 = svec2 - svec3
- s24 = svec2 - svec4
- s34 = svec3 - svec4
-
- sy2 = ss - svec2
- sy3 = ss - svec3
-
-c constants needed for interpolating in s at fixed t lattice points...
- const1 = s13/s23
- const2 = s12/s23
- const3 = s34/s23
- const4 = s24/s23
- s1213 = s12 + s13
- s2434 = s24 + s34
- sdet = s12*s34 - s1213*s2434
- tmp = sy2*sy3/sdet
- const5 = (s34*sy2-s2434*sy3)*tmp/s12
- const6 = (s1213*sy2-s12*sy3)*tmp/s34
-
- EndIf
-
-c --------------Now find lower end of interval containing Q, i.e.,
-c get jq such that qv(jq) .le. q .le. qv(jq+1)...
- JLq = -1
- JU = NT+1
- 12 If (JU-JLq .GT. 1) Then
- JM = (JU+JLq) / 2
- If (tt .GE. TV(JM)) Then
- JLq = JM
- Else
- JU = JM
- Endif
- Goto 12
- Endif
-
- If (JLq .LE. 0) Then
- Jq = 0
- Elseif (JLq .LE. Nt-2) Then
-C keep q in the middle, as shown above
- Jq = JLq - 1
- Else
-C JLq .GE. Nt-1 case: Keep at least 4 points >= Jq.
- Jq = Nt - 3
-
- Endif
-C This is the interpolation variable in Q
-
- If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
-c store the lattice points in t...
- tvec1 = Tv(jq)
- tvec2 = Tv(jq+1)
- tvec3 = Tv(jq+2)
- tvec4 = Tv(jq+3)
-
- t12 = tvec1 - tvec2
- t13 = tvec1 - tvec3
- t23 = tvec2 - tvec3
- t24 = tvec2 - tvec4
- t34 = tvec3 - tvec4
-
- ty2 = tt - tvec2
- ty3 = tt - tvec3
-
- tmp1 = t12 + t13
- tmp2 = t24 + t34
-
- tdet = t12*t34 - tmp1*tmp2
-
- EndIf
-
-
-c get the pdf function values at the lattice points...
-
- If (Iprtn .GE. 3) Then
- Ip = - Iprtn
- Else
- Ip = Iprtn
- EndIf
- jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
-
- Do it = 1, nqvec
-
- J1 = jtmp + it*(NX+1)
-
- If (Jx .Eq. 0) Then
-C For the first 4 x points, interpolate x^2*f(x,Q)
-C This applies to the two lowest bins JLx = 0, 1
-C We can not put the JLx.eq.1 bin into the "interrior" section
-C (as we do for q), since Upd(J1) is undefined.
- fij(1) = 0
- fij(2) = Upd(J1+1) * XV(1)**2
- fij(3) = Upd(J1+2) * XV(2)**2
- fij(4) = Upd(J1+3) * XV(3)**2
-C
-C Use Polint which allows x to be anywhere w.r.t. the grid
-
- Call Polint (XVpow(0), Fij(1), 4, ss, Fx, Dfx)
-
- If (x .GT. 0D0) Fvec(it) = Fx / x**2
-C Pdf is undefined for x.eq.0
- ElseIf (JLx .Eq. Nx-1) Then
-C This is the highest x bin:
-
- Call Polint (XVpow(Nx-3), Upd(J1), 4, ss, Fx, Dfx)
-
- Fvec(it) = Fx
-
- Else
-C for all interior points, use Jon's in-line function
-C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
- sf2 = Upd(J1+1)
- sf3 = Upd(J1+2)
-
- g1 = sf2*const1 - sf3*const2
- g4 = -sf2*const3 + sf3*const4
-
- Fvec(it) = (const5*(Upd(J1)-g1)
- & + const6*(Upd(J1+3)-g4)
- & + sf2*sy3 - sf3*sy2) / s23
-
- Endif
-
- enddo
-C We now have the four values Fvec(1:4)
-c interpolate in t...
-
- If (JLq .LE. 0) Then
-C 1st Q-bin, as well as extrapolation to lower Q
- Call Polint (TV(0), Fvec(1), 4, tt, ff, Dfq)
-
- ElseIf (JLq .GE. Nt-1) Then
-C Last Q-bin, as well as extrapolation to higher Q
- Call Polint (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
- Else
-C Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
-C which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
-C the full range QV(0:Nt) (in contrast to XV)
- tf2 = fvec(2)
- tf3 = fvec(3)
-
- g1 = ( tf2*t13 - tf3*t12) / t23
- g4 = (-tf2*t34 + tf3*t24) / t23
-
- h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
- & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
-
- ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
- EndIf
-
- PartonX6 = ff
-
- Return
-C ********************
- End
diff --git a/gonsalves/README.md b/gonsalves/README.md
deleted file mode 100644
index 5a84ad1..0000000
--- a/gonsalves/README.md
+++ /dev/null
@@ -1,4 +0,0 @@
-Gonsalves, Pawlowski, Wai calculation of dsigma/dqt
-Phys. Rev. D40, 2245 (1989)
-Code from
-http://www.physics.buffalo.edu/gonsalves/ewbqt/public.tar.gz
diff --git a/gonsalves/alphas.f b/gonsalves/alphas.f
deleted file mode 100644
index 4b8f004..0000000
--- a/gonsalves/alphas.f
+++ /dev/null
@@ -1,151 +0,0 @@
-* Extracted from the MRST demonstration program ggh.f
- FUNCTION ALPHAS(SCALE,LAMBDA,IORDER)
-C.
-C. The QCD coupling used in the MRS analysis.
-C. Automatically corrects for NF=3,4,5 with thresholds at m_Q.
-C. The following parameters must be supplied:
-C. SCALE = the QCD scale in GeV (real)
-C. LAMBDA = the 4 flavour MSbar Lambda parameter in GeV (real)
-C. IORDER = 0 for LO, 1 for NLO and 2 for NNLO
-C.
- IMPLICIT REAL*8 (A-H,O-Z)
- REAL*8 LAMBDA
- DATA TOL,QSCT,QSDT/5D-4,74D0,8.18D0/
- PI=4D0*DATAN(1D0)
- IORD=IORDER
- ITH=0
- FLAV=4D0
- AL=LAMBDA
- AL2=LAMBDA*LAMBDA
- QS=SCALE*SCALE
- T=DLOG(QS/AL2)
- TT=T
-C.
- qsctt=qsct/4.
- qsdtt=qsdt/4.
- if(qs.lt.0.5d0) then !! running stops below 0.5
- qs=0.5d0
- t=dlog(qs/al2)
- tt=t
- endif
-
- IF(QS.gt.QSCTT) GO TO 12
- IF(QS.lt.QSDTT) GO TO 312
- 11 CONTINUE
- B0=11-2.*FLAV/3.
- X1=4.*PI/B0
- 5 continue
- IF(IORD.eq.0) then
- ALPHAS=X1/T
- ELSEIF(IORD.eq.1) then
- B1=102.-38.*FLAV/3.
- X2=B1/B0**2
- AS2=X1/T*(1.-X2*dlog(T)/T)
- 95 AS=AS2
- F=-T+X1/AS-X2*dlog(X1/AS+X2)
- FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
- AS2=AS-F/FP
- DEL=ABS(F/FP/AS)
- IF((DEL-TOL).GT.0.) go to 95
- ALPHAS=AS2
- ELSEIF(IORD.eq.2) then
- ALPHAS=qwikalf(t,2,flav)
- ENDIF
- IF(ITH.EQ.0) RETURN
- GO TO (13,14,15) ITH
- print *, 'here'
- GO TO 5
- 12 ITH=1
- T=dlog(QSCTT/AL2)
- GO TO 11
- 13 ALFQC4=ALPHAS
- FLAV=5.
- ITH=2
- GO TO 11
- 14 ALFQC5=ALPHAS
- ITH=3
- T=TT
- GO TO 11
- 15 ALFQS5=ALPHAS
- ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
- ALPHAS=1./ALFINV
- RETURN
-
- 311 CONTINUE
- B0=11-2.*FLAV/3.
- X1=4.*PI/B0
- IF(IORD.eq.0) then
- ALPHAS=X1/T
- ELSEIF(IORD.eq.1) then
- B1=102.-38.*FLAV/3.
- X2=B1/B0**2
- AS2=X1/T*(1.-X2*dlog(T)/T)
- 35 AS=AS2
- F=-T+X1/AS-X2*dlog(X1/AS+X2)
- FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
- AS2=AS-F/FP
- DEL=ABS(F/FP/AS)
- IF((DEL-TOL).GT.0.) go to 35
- ELSEIF(IORD.eq.2) then
- ALPHAS=qwikalf(t,2,flav)
- ENDIF
- IF(ITH.EQ.0) RETURN
- GO TO (313,314,315) ITH
- 312 ITH=1
- T=dlog(QSDTT/AL2)
- GO TO 311
- 313 ALFQC4=ALPHAS
- FLAV=3.
- ITH=2
- GO TO 311
- 314 ALFQC3=ALPHAS
- ITH=3
- T=TT
- GO TO 311
- 315 ALFQS3=ALPHAS
- ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
- ALPHAS=1./ALFINV
- RETURN
- END
- function qwikalf(t,iord,flav)
- implicit real*8(a-h,o-z)
- dimension z3(6),z4(6),z5(6),zz3(6),zz4(6),zz5(6)
- data z3/ -.161667E+01,0.954244E+01,
- .0.768623E+01,0.101523E+00,-.360127E-02,0.457867E-04/
- data z4/ -.172239E+01,0.831185E+01,
- .0.721463E+01,0.835531E-01,-.285436E-02,0.349129E-04/
- data z5/ -.872190E+00,0.572816E+01,
- .0.716119E+01,0.195884E-01,-.300199E-03,0.151741E-05/
- data zz3/-.155611E+02,0.168406E+02,
- .0.603014E+01,0.257682E+00,-.970217E-02,0.127628E-03/
- data zz4/-.106762E+02,0.118497E+02,0.664964E+01,
- .0.112996E+00,-.317551E-02,0.302434E-04/
- data zz5/-.531860E+01,0.708503E+01,0.698352E+01,
- .0.274170E-01,-.426894E-03,0.217591E-05/
- data pi/3.14159/
- nfm2=flav-2.
- x=dsqrt(t)
- x2=x*x
- x3=x*x2
- x4=x*x3
- x5=x*x4
- go to (1,2) iord !!iord change!!
- 1 go to (3,4,5) nfm2
- 3 y=z3(1)+z3(2)*x+z3(3)*x2+z3(4)*x3+z3(5)*x4+z3(6)*x5
- go to 10
- 4 y=z4(1)+z4(2)*x+z4(3)*x2+z4(4)*x3+z4(5)*x4+z4(6)*x5
- go to 10
- 5 y=z5(1)+z5(2)*x+z5(3)*x2+z5(4)*x3+z5(5)*x4+z5(6)*x5
- go to 10
- 2 go to (6,7,8) nfm2
- 6 y=zz3(1)+zz3(2)*x+zz3(3)*x2+zz3(4)*x3+zz3(5)*x4+zz3(6)*x5
- go to 10
- 7 y=zz4(1)+zz4(2)*x+zz4(3)*x2+zz4(4)*x3+zz4(5)*x4+zz4(6)*x5
- go to 10
- 8 y=zz5(1)+zz5(2)*x+zz5(3)*x2+zz5(4)*x3+zz5(5)*x4+zz5(6)*x5
- go to 10
- 10 qwikalf=4.*pi/y
- return
- end
-
-
diff --git a/gonsalves/cteq.f b/gonsalves/cteq.f
deleted file mode 100644
index d25549d..0000000
--- a/gonsalves/cteq.f
+++ /dev/null
@@ -1,25 +0,0 @@
-* CTEQ Parton Densities
- SUBROUTINE cteq(nset,x,qsq,pden)
- IMPLICIT NONE
- INTEGER nset,mode,i
- DOUBLE PRECISION x,qsq,pden(-6:6)
- DOUBLE PRECISION q,upv,dnv,usea,dsea,str,chm,bot,glu
- DOUBLE PRECISION Ctq6Pdf
- EXTERNAL Ctq6Pdf
-
- q=dsqrt(qsq)
- IF (nset.EQ.1) THEN
- mode=1
- CALL SetCtq6(mode)
- DO i=-5,5
- pden(i) = Ctq6Pdf(i,x,q)
- ENDDO
- pden(6) = 0d0
- pden(-6) = 0d0
- ELSE
- PRINT *,' bad CTEQ mode number',nset
- STOP
- ENDIF
-
- RETURN
- END
diff --git a/gonsalves/dat.tar.gz b/gonsalves/dat.tar.gz
deleted file mode 100644
index 6b69f83..0000000
Binary files a/gonsalves/dat.tar.gz and /dev/null differ
diff --git a/gonsalves/dopden.f b/gonsalves/dopden.f
deleted file mode 100644
index c308c83..0000000
--- a/gonsalves/dopden.f
+++ /dev/null
@@ -1,206 +0,0 @@
-C-----------------------------------------------------------------------
-C Duke-Owens' Q**2 dependent parton densities
-C Phys. Rev. D 30, 49 (1984)
-C
- SUBROUTINE dopden (nset,x,qsq,pden)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION pden(-6:6)
-
-C--- Q for start of evolution
- q0=2.
-
-C--- The output distributions for the proton are
-C pden(0) = gluon
-C pden(1) = valence up quark + up sea
-C pden(2) = valence down quark + down sea
-C pden(3) = strange sea
-C pden(4) = charm sea
-C pden(5) = bottom sea
-C pden(6) = top sea
-C pden(-1) = anti-up sea
-C pden(-2) = anti-down sea
-C pden(-3) = anti-strange sea
-C pden(-4) = anti-charm sea
-C pden(-5) = anti-bottom sea
-C pden(-6) = anti-top sea
-
- IF (nset.LT.1.OR.nset.GT.2) THEN
- PRINT *,'Bad Duke Owens Set ',nset
- STOP
- ENDIF
- GOTO (1,2) nset
-1 CONTINUE
-C--- nset = 1 : set 1 Lambda = 0.2 GeV/c
-
- xlambd=0.2
- s=dlog(qsq/xlambd**2)
- s=s/(2.*dlog(q0/xlambd))
- s=dlog(s)
-
-C valence quark parameters
-
- eta1=0.419+0.004*s-0.007*s**2
- eta2=3.46+0.724*s-0.066*s**2
- gamud=4.40-4.86*s+1.33*s**2
- eta3=0.763-0.237*s+0.026*s**2
- eta4=4.00+0.627*s-0.019*s**2
- gamd=-0.421*s+0.033*s**2
-
-C up-down-strange sea quark parameters
-
- saa=1.265-1.132*s+0.293*s**2
- sa=-0.372*s-0.029*s**2
- sb=8.05+1.59*s-0.153*s**2
- salp=6.31*s-0.273*s**2
- sbet=-10.5*s-3.17*s**2
- sgam=14.7*s+9.80*s**2
-
-C charmed quark parameters
-
- chaa=0.135*s-0.075*s**2
- cha=-0.036-0.222*s-0.058*s**2
- chb=6.35+3.26*s-0.909*s**2
- chalp=-3.03*s+1.50*s**2
- chbet=17.4*s-11.3*s**2
- chgam=-17.9*s+15.6*s**2
-
-C gluon parameters
-
- glaa=1.56-1.71*s+0.638*s**2
- gla=-0.949*s+0.325*s**2
- glb=6.0+1.44*s-1.05*s**2
- glalp=9.0-7.19*s+0.255*s**2
- glbet=-16.5*s+10.9*s**2
- glgam=15.3*s-10.1*s**2
-
- GOTO 3
-
-2 CONTINUE
-C--- nset = 2 : set 2 Lambda = 0.4 GeV/c
-
- xlambd=0.4
- s=dlog(qsq/xlambd**2)
- s=s/(2.*dlog(q0/xlambd))
- s=dlog(s)
-
-C valence quark parameters
-
- eta1=0.374+0.014*s
- eta2=3.33+0.753*s-0.076*s**2
- gamud=6.03-6.22*s+1.56*s**2
- eta3=0.761-0.232*s+0.023*s**2
- eta4=3.83+0.627*s-0.019*s**2
- gamd=-0.418*s+0.036*s**2
-
-C up-down-strange sea quark parameters
-
- saa=1.67-1.92*s+0.582*s**2
- sa=-0.273*s-0.164*s**2
- sb=9.15+0.530*s-0.763*s**2
- salp=15.7*s-2.83*s**2
- sbet=-101.*s+44.7*s**2
- sgam=223.*s-117.*s**2
-
-C charmed quark parameters
-
- chaa=0.067*s-0.031*s**2
- cha=-0.120-0.233*s-0.023*s**2
- chb=3.51+3.66*s-0.453*s**2
- chalp=-0.474*s+0.358*s**2
- chbet=9.50*s-5.43*s**2
- chgam=-16.6*s+15.5*s**2
-
-C gluon parameters
-
- glaa=0.879-0.971*s+0.434*s**2
- gla=-1.16*s+0.476*s**2
- glb=4.0+1.23*s-0.254*s**2
- glalp=9.0-5.64*s-0.817*s**2
- glbet=-7.54*s+5.50*s**2
- glgam=-0.596*s+0.126*s**2
-
-3 CONTINUE
-
-C gluon density
-
- xglu=glaa*x**gla*(1.-x)**glb*(1.+glalp*x+glbet*x**2+glgam*x**3)
- pden(0)=xglu/x
-
-C valence quark densities
-
- xnud=3./beta(eta1,eta2+1.)
- xnud=xnud/(1.+gamud*eta1/(eta1+eta2+1.))
- xnd=1./beta(eta3,eta4+1.)
- xnd=xnd/(1.+gamd*eta3/(eta3+eta4+1.))
-
- xuvdv=xnud*x**eta1*(1.-x)**eta2*(1.+gamud*x)
- xdv=xnd*x**eta3*(1.-x)**eta4*(1.+gamd*x)
- pden(1)=(xuvdv-xdv)/x
- pden(2)=xdv/x
-
-C up-down-strange sea quark densities
-
- xsea=saa*x**sa*(1.-x)**sb*(1.+salp*x+sbet*x**2+sgam*x**3)
- pden(-1)=xsea/x/6.0
- pden(-2)=pden(3)
- pden(-3)=pden(3)
-
-C charmed quark density
-
- xchm=chaa*x**cha*(1.-x)**chb*(1.+chalp*x+chbet*x**2+chgam*x**3)
- pden(-4)=xchm/x
-
-C--- top and bottom are not implemented - return 0
- pden(-5)=0.
- pden(-6)=0.
-
-C--- rearrange
- pden(1)=pden(1)+pden(-1)
- pden(2)=pden(2)+pden(-2)
- do i=3,6
- pden(i)=pden(-i)
- enddo
-
- RETURN
- END
-C-----------------------------------------------------------------------
-C Logarithm of the gamma function - numerical recipes
-C
- DOUBLE PRECISION FUNCTION gammln (xx)
-C--- Routine supposed to be accurate for xx > 1.
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION cof(6)
- DATA cof,stp /76.18009173d0,-86.50532033d0,24.01409822d0,
- & -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
- DATA half,one,fpf /0.5d0,1.0d0,5.5d0/
- x=xx-one
- tmp=x+fpf
- tmp=(x+half)*dlog(tmp)-tmp
- ser=one
- DO 11 j=1,6
- x=x+one
- ser=ser+cof(j)/x
-11 CONTINUE
- gammln=tmp+dlog(stp*ser)
- RETURN
- END
-C-----------------------------------------------------------------------
-C The beta function : both arguments positive
-C
- DOUBLE PRECISION FUNCTION beta (xx,yy)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DATA one /1.0d0/
- x=xx
- y=yy
- IF (xx.LT.one) x=x+one
- IF (yy.LT.one) y=y+one
- beta=gammln(x)+gammln(y)-gammln(x+y)
- beta=dexp(beta)
- IF (xx.LT.one) THEN
- beta=beta*(one+yy/xx)
- IF (yy.LT.one) beta=beta*(one+x/yy)
- ELSE
- IF (yy.LT.one) beta=beta*(one+xx/yy)
- ENDIF
- RETURN
- END
diff --git a/gonsalves/ehlq.f b/gonsalves/ehlq.f
deleted file mode 100644
index 988ebb8..0000000
--- a/gonsalves/ehlq.f
+++ /dev/null
@@ -1,522 +0,0 @@
-C-----------------------------------------------------------------------
- SUBROUTINE ehlq (mode,x,qsq,pden)
- IMPLICIT NONE
- INTEGER mode
- DOUBLE PRECISION x,qsq,pden(-6:6)
- DOUBLE PRECISION q,upv,dnv,sea,str,chm,bot,top,glu
-
- q=dsqrt(qsq)
- IF (mode.EQ.1) THEN
- CALL sfehlq1(x,q,upv,dnv,sea,str,chm,bot,top,glu)
- ELSEIF (mode.EQ.2) THEN
- CALL sfehlq2(x,q,upv,dnv,sea,str,chm,bot,top,glu)
- ELSE
- PRINT *,' bad EHLQ mode number'
- STOP
- ENDIF
- pden(0) = glu/x
- pden(1) = (upv+sea)/x
- pden(-1) = sea/x
- pden(2) = (dnv+sea)/x
- pden(-2) = sea/x
- pden(3) = str/x
- pden(-3) = pden(3)
- pden(4) = chm/x
- pden(-4) = pden(4)
- pden(5) = bot/x
- pden(-5) = pden(5)
- pden(6) = top
- pden(-6) = pden(6)
-
- RETURN
- END
-C-----------------------------------------------------------------------
-* $Id: sfehlq1.F,v 1.1.1.2 1996/10/30 08:29:46 cernlib Exp $
-*
-* $Log: sfehlq1.F,v $
-* Revision 1.1.1.2 1996/10/30 08:29:46 cernlib
-* Version 7.04
-*
-* Revision 1.1.1.1 1996/04/12 15:29:35 plothow
-* Version 7.01
-*
-*
-* #include "pdf/pilot.h"
-C-----------------------------------------------------------------------
- SUBROUTINE SFEHLQ1(DX,DSCALE,DUP,DDN,DSEA,DSTR,DCHM,DBOT,DTOP,DGL)
-C ::::::: EHLQ SET 1 SUPERFAST DOUBLE PRECISION :::::::
-C
-* #include "pdf/expdp.inc"
- DOUBLE PRECISION
- + DX,DSCALE,DUP,DDN,DSEA,DSTR,DCHM,DBOT,DTOP,DGL
- PARAMETER (ALAM=0.2, Q02=5.0)
-* #include "pdf/w50511.inc"
- PARAMETER (BMAS=4.25,TMAS=178.1)
- DIMENSION XQ(9),TX(6),TT(6),TS(6),NEHLQ(8),CEHLQ(6,6,2,8)
-C
-C...The following data lines are coefficients needed in the
-C...Eichten, Hinchliffe, Lane, Quigg proton structure function
-C...parametrizations, see below.
-C...Powers of 1-x in different cases.
- DATA NEHLQ/3,4,7,5,7,7,7,7/
-C...Expansion coefficients for up valence quark distribution.
- DATA (((CEHLQ(IX,IT,NX,1),IX=1,6),IT=1,6),NX=1,2)/
- 1 7.677E-01,-2.087E-01,-3.303E-01,-2.517E-02,-1.570E-02,-1.000E-04,
- 2-5.326E-01,-2.661E-01, 3.201E-01, 1.192E-01, 2.434E-02, 7.620E-03,
- 3 2.162E-01, 1.881E-01,-8.375E-02,-6.515E-02,-1.743E-02,-5.040E-03,
- 4-9.211E-02,-9.952E-02, 1.373E-02, 2.506E-02, 8.770E-03, 2.550E-03,
- 5 3.670E-02, 4.409E-02, 9.600E-04,-7.960E-03,-3.420E-03,-1.050E-03,
- 6-1.549E-02,-2.026E-02,-3.060E-03, 2.220E-03, 1.240E-03, 4.100E-04,
- 1 2.395E-01, 2.905E-01, 9.778E-02, 2.149E-02, 3.440E-03, 5.000E-04,
- 2 1.751E-02,-6.090E-03,-2.687E-02,-1.916E-02,-7.970E-03,-2.750E-03,
- 3-5.760E-03,-5.040E-03, 1.080E-03, 2.490E-03, 1.530E-03, 7.500E-04,
- 4 1.740E-03, 1.960E-03, 3.000E-04,-3.400E-04,-2.900E-04,-1.800E-04,
- 5-5.300E-04,-6.400E-04,-1.700E-04, 4.000E-05, 6.000E-05, 4.000E-05,
- 6 1.700E-04, 2.200E-04, 8.000E-05, 1.000E-05,-1.000E-05,-1.000E-05/
-C...Expansion coefficients for down valence quark distribution.
- DATA (((CEHLQ(IX,IT,NX,2),IX=1,6),IT=1,6),NX=1,2)/
- 1 3.813E-01,-8.090E-02,-1.634E-01,-2.185E-02,-8.430E-03,-6.200E-04,
- 2-2.948E-01,-1.435E-01, 1.665E-01, 6.638E-02, 1.473E-02, 4.080E-03,
- 3 1.252E-01, 1.042E-01,-4.722E-02,-3.683E-02,-1.038E-02,-2.860E-03,
- 4-5.478E-02,-5.678E-02, 8.900E-03, 1.484E-02, 5.340E-03, 1.520E-03,
- 5 2.220E-02, 2.567E-02,-3.000E-05,-4.970E-03,-2.160E-03,-6.500E-04,
- 6-9.530E-03,-1.204E-02,-1.510E-03, 1.510E-03, 8.300E-04, 2.700E-04,
- 1 1.261E-01, 1.354E-01, 3.958E-02, 8.240E-03, 1.660E-03, 4.500E-04,
- 2 3.890E-03,-1.159E-02,-1.625E-02,-9.610E-03,-3.710E-03,-1.260E-03,
- 3-1.910E-03,-5.600E-04, 1.590E-03, 1.590E-03, 8.400E-04, 3.900E-04,
- 4 6.400E-04, 4.900E-04,-1.500E-04,-2.900E-04,-1.800E-04,-1.000E-04,
- 5-2.000E-04,-1.900E-04, 0.000E+00, 6.000E-05, 4.000E-05, 3.000E-05,
- 6 7.000E-05, 8.000E-05, 2.000E-05,-1.000E-05,-1.000E-05,-1.000E-05/
-C...Expansion coefficients for up and down sea quark distributions.
- DATA (((CEHLQ(IX,IT,NX,3),IX=1,6),IT=1,6),NX=1,2)/
- 1 6.870E-02,-6.861E-02, 2.973E-02,-5.400E-03, 3.780E-03,-9.700E-04,
- 2-1.802E-02, 1.400E-04, 6.490E-03,-8.540E-03, 1.220E-03,-1.750E-03,
- 3-4.650E-03, 1.480E-03,-5.930E-03, 6.000E-04,-1.030E-03,-8.000E-05,
- 4 6.440E-03, 2.570E-03, 2.830E-03, 1.150E-03, 7.100E-04, 3.300E-04,
- 5-3.930E-03,-2.540E-03,-1.160E-03,-7.700E-04,-3.600E-04,-1.900E-04,
- 6 2.340E-03, 1.930E-03, 5.300E-04, 3.700E-04, 1.600E-04, 9.000E-05,
- 1 1.014E+00,-1.106E+00, 3.374E-01,-7.444E-02, 8.850E-03,-8.700E-04,
- 2 9.233E-01,-1.285E+00, 4.475E-01,-9.786E-02, 1.419E-02,-1.120E-03,
- 3 4.888E-02,-1.271E-01, 8.606E-02,-2.608E-02, 4.780E-03,-6.000E-04,
- 4-2.691E-02, 4.887E-02,-1.771E-02, 1.620E-03, 2.500E-04,-6.000E-05,
- 5 7.040E-03,-1.113E-02, 1.590E-03, 7.000E-04,-2.000E-04, 0.000E+00,
- 6-1.710E-03, 2.290E-03, 3.800E-04,-3.500E-04, 4.000E-05, 1.000E-05/
-C...Expansion coefficients for gluon distribution.
- DATA (((CEHLQ(IX,IT,NX,4),IX=1,6),IT=1,6),NX=1,2)/
- 1 9.482E-01,-9.578E-01, 1.009E-01,-1.051E-01, 3.456E-02,-3.054E-02,
- 2-9.627E-01, 5.379E-01, 3.368E-01,-9.525E-02, 1.488E-02,-2.051E-02,
- 3 4.300E-01,-8.306E-02,-3.372E-01, 4.902E-02,-9.160E-03, 1.041E-02,
- 4-1.925E-01,-1.790E-02, 2.183E-01, 7.490E-03, 4.140E-03,-1.860E-03,
- 5 8.183E-02, 1.926E-02,-1.072E-01,-1.944E-02,-2.770E-03,-5.200E-04,
- 6-3.884E-02,-1.234E-02, 5.410E-02, 1.879E-02, 3.350E-03, 1.040E-03,
- 1 2.948E+01,-3.902E+01, 1.464E+01,-3.335E+00, 5.054E-01,-5.915E-02,
- 2 2.559E+01,-3.955E+01, 1.661E+01,-4.299E+00, 6.904E-01,-8.243E-02,
- 3-1.663E+00, 1.176E+00, 1.118E+00,-7.099E-01, 1.948E-01,-2.404E-02,
- 4-2.168E-01, 8.170E-01,-7.169E-01, 1.851E-01,-1.924E-02,-3.250E-03,
- 5 2.088E-01,-4.355E-01, 2.239E-01,-2.446E-02,-3.620E-03, 1.910E-03,
- 6-9.097E-02, 1.601E-01,-5.681E-02,-2.500E-03, 2.580E-03,-4.700E-04/
-C...Expansion coefficients for strange sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,5),IX=1,6),IT=1,6),NX=1,2)/
- 1 4.968E-02,-4.173E-02, 2.102E-02,-3.270E-03, 3.240E-03,-6.700E-04,
- 2-6.150E-03,-1.294E-02, 6.740E-03,-6.890E-03, 9.000E-04,-1.510E-03,
- 3-8.580E-03, 5.050E-03,-4.900E-03,-1.600E-04,-9.400E-04,-1.500E-04,
- 4 7.840E-03, 1.510E-03, 2.220E-03, 1.400E-03, 7.000E-04, 3.500E-04,
- 5-4.410E-03,-2.220E-03,-8.900E-04,-8.500E-04,-3.600E-04,-2.000E-04,
- 6 2.520E-03, 1.840E-03, 4.100E-04, 3.900E-04, 1.600E-04, 9.000E-05,
- 1 9.235E-01,-1.085E+00, 3.464E-01,-7.210E-02, 9.140E-03,-9.100E-04,
- 2 9.315E-01,-1.274E+00, 4.512E-01,-9.775E-02, 1.380E-02,-1.310E-03,
- 3 4.739E-02,-1.296E-01, 8.482E-02,-2.642E-02, 4.760E-03,-5.700E-04,
- 4-2.653E-02, 4.953E-02,-1.735E-02, 1.750E-03, 2.800E-04,-6.000E-05,
- 5 6.940E-03,-1.132E-02, 1.480E-03, 6.500E-04,-2.100E-04, 0.000E+00,
- 6-1.680E-03, 2.340E-03, 4.200E-04,-3.400E-04, 5.000E-05, 1.000E-05/
-C...Expansion coefficients for charm sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,6),IX=1,6),IT=1,6),NX=1,2)/
- 1 9.270E-03,-1.817E-02, 9.590E-03,-6.390E-03, 1.690E-03,-1.540E-03,
- 2 5.710E-03,-1.188E-02, 6.090E-03,-4.650E-03, 1.240E-03,-1.310E-03,
- 3-3.960E-03, 7.100E-03,-3.590E-03, 1.840E-03,-3.900E-04, 3.400E-04,
- 4 1.120E-03,-1.960E-03, 1.120E-03,-4.800E-04, 1.000E-04,-4.000E-05,
- 5 4.000E-05,-3.000E-05,-1.800E-04, 9.000E-05,-5.000E-05,-2.000E-05,
- 6-4.200E-04, 7.300E-04,-1.600E-04, 5.000E-05, 5.000E-05, 5.000E-05,
- 1 8.098E-01,-1.042E+00, 3.398E-01,-6.824E-02, 8.760E-03,-9.000E-04,
- 2 8.961E-01,-1.217E+00, 4.339E-01,-9.287E-02, 1.304E-02,-1.290E-03,
- 3 3.058E-02,-1.040E-01, 7.604E-02,-2.415E-02, 4.600E-03,-5.000E-04,
- 4-2.451E-02, 4.432E-02,-1.651E-02, 1.430E-03, 1.200E-04,-1.000E-04,
- 5 1.122E-02,-1.457E-02, 2.680E-03, 5.800E-04,-1.200E-04, 3.000E-05,
- 6-7.730E-03, 7.330E-03,-7.600E-04,-2.400E-04, 1.000E-05, 0.000E+00/
-C...Expansion coefficients for bottom sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,7),IX=1,6),IT=1,6),NX=1,2)/
- 1 9.010E-03,-1.401E-02, 7.150E-03,-4.130E-03, 1.260E-03,-1.040E-03,
- 2 6.280E-03,-9.320E-03, 4.780E-03,-2.890E-03, 9.100E-04,-8.200E-04,
- 3-2.930E-03, 4.090E-03,-1.890E-03, 7.600E-04,-2.300E-04, 1.400E-04,
- 4 3.900E-04,-1.200E-03, 4.400E-04,-2.500E-04, 2.000E-05,-2.000E-05,
- 5 2.600E-04, 1.400E-04,-8.000E-05, 1.000E-04, 1.000E-05, 1.000E-05,
- 6-2.600E-04, 3.200E-04, 1.000E-05,-1.000E-05, 1.000E-05,-1.000E-05,
- 1 8.029E-01,-1.075E+00, 3.792E-01,-7.843E-02, 1.007E-02,-1.090E-03,
- 2 7.903E-01,-1.099E+00, 4.153E-01,-9.301E-02, 1.317E-02,-1.410E-03,
- 3-1.704E-02,-1.130E-02, 2.882E-02,-1.341E-02, 3.040E-03,-3.600E-04,
- 4-7.200E-04, 7.230E-03,-5.160E-03, 1.080E-03,-5.000E-05,-4.000E-05,
- 5 3.050E-03,-4.610E-03, 1.660E-03,-1.300E-04,-1.000E-05, 1.000E-05,
- 6-4.360E-03, 5.230E-03,-1.610E-03, 2.000E-04,-2.000E-05, 0.000E+00/
-C...Expansion coefficients for top sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,8),IX=1,6),IT=1,6),NX=1,2)/
- 1 4.410E-03,-7.480E-03, 3.770E-03,-2.580E-03, 7.300E-04,-7.100E-04,
- 2 3.840E-03,-6.050E-03, 3.030E-03,-2.030E-03, 5.800E-04,-5.900E-04,
- 3-8.800E-04, 1.660E-03,-7.500E-04, 4.700E-04,-1.000E-04, 1.000E-04,
- 4-8.000E-05,-1.500E-04, 1.200E-04,-9.000E-05, 3.000E-05, 0.000E+00,
- 5 1.300E-04,-2.200E-04,-2.000E-05,-2.000E-05,-2.000E-05,-2.000E-05,
- 6-7.000E-05, 1.900E-04,-4.000E-05, 2.000E-05, 0.000E+00, 0.000E+00,
- 1 6.623E-01,-9.248E-01, 3.519E-01,-7.930E-02, 1.110E-02,-1.180E-03,
- 2 6.380E-01,-9.062E-01, 3.582E-01,-8.479E-02, 1.265E-02,-1.390E-03,
- 3-2.581E-02, 2.125E-02, 4.190E-03,-4.980E-03, 1.490E-03,-2.100E-04,
- 4 7.100E-04, 5.300E-04,-1.270E-03, 3.900E-04,-5.000E-05,-1.000E-05,
- 5 3.850E-03,-5.060E-03, 1.860E-03,-3.500E-04, 4.000E-05, 0.000E+00,
- 6-3.530E-03, 4.460E-03,-1.500E-03, 2.700E-04,-3.000E-05, 0.000E+00/
-C
-C
-C...Proton structure functions from Eichten, Hinchliffe, Lane, Quigg.
-C...Allowed variable range: 5 GeV^2 < Q^2 < 1E8 GeV^2; 1E-4 < x < 1
-C
-*...Set number of flavors - modify CERNLIB routine
- NFL=4
- IF (DSCALE.GT.BMAS) NFL=5
- IF (DSCALE.GT.TMAS) NFL=6
-C...Determine Lamdba and x and t expansion variables.
- X = DX
- SCALE = DSCALE
- KF = ABS(NFL)
- Q2 =SCALE**2
- TMIN=LOG(Q02/ALAM**2)
- TMAX=LOG(1E8/ALAM**2)
- T=LOG(MAX(1.,Q2/ALAM**2))
- VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN)))
- NX=1
- IF(X.LE.0.1) NX=2
- IF(NX.EQ.1) VX=(2.*X-1.1)/0.9
- IF(NX.EQ.2) VX=MAX(-1.,(2.*LOG(X)+11.51293)/6.90776)
- CXS=1.
-
-C...Chebyshev polynomials for x and t expansion.
- TX(1)=1.
- TX(2)=VX
- TX(3)=2.*VX**2-1.
- TX(4)=4.*VX**3-3.*VX
- TX(5)=8.*VX**4-8.*VX**2+1.
- TX(6)=16.*VX**5-20.*VX**3+5.*VX
- TT(1)=1.
- TT(2)=VT
- TT(3)=2.*VT**2-1.
- TT(4)=4.*VT**3-3.*VT
- TT(5)=8.*VT**4-8.*VT**2+1.
- TT(6)=16.*VT**5-20.*VT**3+5.*VT
-
-C...Calculate structure functions.
- DO 120 KFL=1,6
- XQSUM=0.
- DO 110 IT=1,6
- DO 110 IX=1,6
- 110 XQSUM=XQSUM+CEHLQ(IX,IT,NX,KFL)*TX(IX)*TT(IT)
- 120 XQ(KFL)=XQSUM*(1.-X)**NEHLQ(KFL)*CXS
-C
-C...Put into output variables
- DUP = XQ(1)
- DDN = XQ(2)
- DSEA = XQ(3)
- DGL = XQ(4)
- DSTR = XQ(5)
- DCHM = XQ(6)
- DBOT = 0.D0
- DTOP = 0.D0
-C
- IF (KF.GE.5) THEN
-C
-C...Special expansion for bottom (threshold effects).
- BOT = 0.
- TMIN=8.1905
- IF(T.LE.TMIN) GOTO 140
- VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN)))
- TT(1)=1.
- TT(2)=VT
- TT(3)=2.*VT**2-1.
- TT(4)=4.*VT**3-3.*VT
- TT(5)=8.*VT**4-8.*VT**2+1.
- TT(6)=16.*VT**5-20.*VT**3+5.*VT
- XQSUM=0.
- DO 130 IT=1,6
- DO 130 IX=1,6
- 130 XQSUM=XQSUM+CEHLQ(IX,IT,NX,7)*TX(IX)*TT(IT)
- BOT=XQSUM*(1.-X)**NEHLQ(7)*CXS
- 140 CONTINUE
- ENDIF
- DBOT = BOT
-C--
- IF (KF.GE.6) THEN
-C
-C...Special expansion for top (threshold effects).
- TOP = 0.
- TMIN=11.5528
- TMIN=TMIN+2.*LOG(TMAS/30.)
- TMAX=TMAX+2.*LOG(TMAS/30.)
- IF(T.LE.TMIN) GOTO 160
- VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN)))
- TT(1)=1.
- TT(2)=VT
- TT(3)=2.*VT**2-1.
- TT(4)=4.*VT**3-3.*VT
- TT(5)=8.*VT**4-8.*VT**2+1.
- TT(6)=16.*VT**5-20.*VT**3+5.*VT
- XQSUM=0.
- DO 150 IT=1,6
- DO 150 IX=1,6
- 150 XQSUM=XQSUM+CEHLQ(IX,IT,NX,8)*TX(IX)*TT(IT)
- TOP=XQSUM*(1.-X)**NEHLQ(8)*CXS
- 160 CONTINUE
- ENDIF
- DTOP = TOP
-C
- RETURN
- END
-C-----------------------------------------------------------------------
-* $Id: sfehlq2.F,v 1.1.1.2 1996/10/30 08:29:47 cernlib Exp $
-*
-* $Log: sfehlq2.F,v $
-* Revision 1.1.1.2 1996/10/30 08:29:47 cernlib
-* Version 7.04
-*
-* Revision 1.1.1.1 1996/04/12 15:29:35 plothow
-* Version 7.01
-*
-*
-* #include "pdf/pilot.h"
-C-----------------------------------------------------------------------
- SUBROUTINE SFEHLQ2(DX,DSCALE,DUP,DDN,DSEA,DSTR,DCHM,DBOT,DTOP,DGL)
-C ::::::: EHLQ SET 2 SUPERFAST DOUBLE PRECISION :::::::
-C
-* #include "pdf/expdp.inc"
- DOUBLE PRECISION
- + DX,DSCALE,DUP,DDN,DSEA,DSTR,DCHM,DBOT,DTOP,DGL
- PARAMETER (ALAM=0.29, Q02=5.0)
-* #include "pdf/w50511.inc"
- PARAMETER (BMAS=4.25,TMAS=178.1)
- DIMENSION XQ(9),TX(6),TT(6),TS(6),NEHLQ(8),CEHLQ(6,6,2,8)
-C
-C...The following data lines are coefficients needed in the
-C...Eichten, Hinchliffe, Lane, Quigg proton structure function
-C...parametrizations, see below.
-C...Powers of 1-x in different cases.
- DATA NEHLQ/3,4,7,6,7,7,7,7/
-C...Expansion coefficients for up valence quark distribution.
- DATA (((CEHLQ(IX,IT,NX,1),IX=1,6),IT=1,6),NX=1,2)/
- 1 7.237E-01,-2.189E-01,-2.995E-01,-1.909E-02,-1.477E-02, 2.500E-04,
- 2-5.314E-01,-2.425E-01, 3.283E-01, 1.119E-01, 2.223E-02, 7.070E-03,
- 3 2.289E-01, 1.890E-01,-9.859E-02,-6.900E-02,-1.747E-02,-5.080E-03,
- 4-1.041E-01,-1.084E-01, 2.108E-02, 2.975E-02, 9.830E-03, 2.830E-03,
- 5 4.394E-02, 5.116E-02,-1.410E-03,-1.055E-02,-4.230E-03,-1.270E-03,
- 6-1.991E-02,-2.539E-02,-2.780E-03, 3.430E-03, 1.720E-03, 5.500E-04,
- 1 2.410E-01, 2.884E-01, 9.369E-02, 1.900E-02, 2.530E-03, 2.400E-04,
- 2 1.765E-02,-9.220E-03,-3.037E-02,-2.085E-02,-8.440E-03,-2.810E-03,
- 3-6.450E-03,-5.260E-03, 1.720E-03, 3.110E-03, 1.830E-03, 8.700E-04,
- 4 2.120E-03, 2.320E-03, 2.600E-04,-4.900E-04,-3.900E-04,-2.300E-04,
- 5-6.900E-04,-8.200E-04,-2.000E-04, 7.000E-05, 9.000E-05, 6.000E-05,
- 6 2.400E-04, 3.100E-04, 1.100E-04, 0.000E+00,-2.000E-05,-2.000E-05/
-C...Expansion coefficients for down valence quark distribution.
- DATA (((CEHLQ(IX,IT,NX,2),IX=1,6),IT=1,6),NX=1,2)/
- 1 3.578E-01,-8.622E-02,-1.480E-01,-1.840E-02,-7.820E-03,-4.500E-04,
- 2-2.925E-01,-1.304E-01, 1.696E-01, 6.243E-02, 1.353E-02, 3.750E-03,
- 3 1.318E-01, 1.041E-01,-5.486E-02,-3.872E-02,-1.038E-02,-2.850E-03,
- 4-6.162E-02,-6.143E-02, 1.303E-02, 1.740E-02, 5.940E-03, 1.670E-03,
- 5 2.643E-02, 2.957E-02,-1.490E-03,-6.450E-03,-2.630E-03,-7.700E-04,
- 6-1.218E-02,-1.497E-02,-1.260E-03, 2.240E-03, 1.120E-03, 3.500E-04,
- 1 1.263E-01, 1.334E-01, 3.732E-02, 7.070E-03, 1.260E-03, 3.400E-04,
- 2 3.660E-03,-1.357E-02,-1.795E-02,-1.031E-02,-3.880E-03,-1.280E-03,
- 3-2.100E-03,-3.600E-04, 2.050E-03, 1.920E-03, 9.800E-04, 4.400E-04,
- 4 7.700E-04, 5.400E-04,-2.400E-04,-3.900E-04,-2.400E-04,-1.300E-04,
- 5-2.600E-04,-2.300E-04, 2.000E-05, 9.000E-05, 6.000E-05, 4.000E-05,
- 6 9.000E-05, 1.000E-04, 2.000E-05,-2.000E-05,-2.000E-05,-1.000E-05/
-C...Expansion coefficients for up and down sea quark distributions.
- DATA (((CEHLQ(IX,IT,NX,3),IX=1,6),IT=1,6),NX=1,2)/
- 1 1.008E-01,-7.100E-02, 1.973E-02,-5.710E-03, 2.930E-03,-9.900E-04,
- 2-5.271E-02,-1.823E-02, 1.792E-02,-6.580E-03, 1.750E-03,-1.550E-03,
- 3 1.220E-02, 1.763E-02,-8.690E-03,-8.800E-04,-1.160E-03,-2.100E-04,
- 4-1.190E-03,-7.180E-03, 2.360E-03, 1.890E-03, 7.700E-04, 4.100E-04,
- 5-9.100E-04, 2.040E-03,-3.100E-04,-1.050E-03,-4.000E-04,-2.400E-04,
- 6 1.190E-03,-1.700E-04,-2.000E-04, 4.200E-04, 1.700E-04, 1.000E-04,
- 1 1.081E+00,-1.189E+00, 3.868E-01,-8.617E-02, 1.115E-02,-1.180E-03,
- 2 9.917E-01,-1.396E+00, 4.998E-01,-1.159E-01, 1.674E-02,-1.720E-03,
- 3 5.099E-02,-1.338E-01, 9.173E-02,-2.885E-02, 5.890E-03,-6.500E-04,
- 4-3.178E-02, 5.703E-02,-2.070E-02, 2.440E-03, 1.100E-04,-9.000E-05,
- 5 8.970E-03,-1.392E-02, 2.050E-03, 6.500E-04,-2.300E-04, 2.000E-05,
- 6-2.340E-03, 3.010E-03, 5.000E-04,-3.900E-04, 6.000E-05, 1.000E-05/
-C...Expansion coefficients for gluon distribution.
- DATA (((CEHLQ(IX,IT,NX,4),IX=1,6),IT=1,6),NX=1,2)/
- 1 2.367E+00, 4.453E-01, 3.660E-01, 9.467E-02, 1.341E-01, 1.661E-02,
- 2-3.170E+00,-1.795E+00, 3.313E-02,-2.874E-01,-9.827E-02,-7.119E-02,
- 3 1.823E+00, 1.457E+00,-2.465E-01, 3.739E-02, 6.090E-03, 1.814E-02,
- 4-1.033E+00,-9.827E-01, 2.136E-01, 1.169E-01, 5.001E-02, 1.684E-02,
- 5 5.133E-01, 5.259E-01,-1.173E-01,-1.139E-01,-4.988E-02,-2.021E-02,
- 6-2.881E-01,-3.145E-01, 5.667E-02, 9.161E-02, 4.568E-02, 1.951E-02,
- 1 3.036E+01,-4.062E+01, 1.578E+01,-3.699E+00, 6.020E-01,-7.031E-02,
- 2 2.700E+01,-4.167E+01, 1.770E+01,-4.804E+00, 7.862E-01,-1.060E-01,
- 3-1.909E+00, 1.357E+00, 1.127E+00,-7.181E-01, 2.232E-01,-2.481E-02,
- 4-2.488E-01, 9.781E-01,-8.127E-01, 2.094E-01,-2.997E-02,-4.710E-03,
- 5 2.506E-01,-5.427E-01, 2.672E-01,-3.103E-02,-1.800E-03, 2.870E-03,
- 6-1.128E-01, 2.087E-01,-6.972E-02,-2.480E-03, 2.630E-03,-8.400E-04/
-C...Expansion coefficients for strange sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,5),IX=1,6),IT=1,6),NX=1,2)/
- 1 6.478E-02,-4.537E-02, 1.643E-02,-3.490E-03, 2.710E-03,-6.700E-04,
- 2-2.223E-02,-2.126E-02, 1.247E-02,-6.290E-03, 1.120E-03,-1.440E-03,
- 3-1.340E-03, 1.362E-02,-6.130E-03,-7.900E-04,-9.000E-04,-2.000E-04,
- 4 5.080E-03,-3.610E-03, 1.700E-03, 1.830E-03, 6.800E-04, 4.000E-04,
- 5-3.580E-03, 6.000E-05,-2.600E-04,-1.050E-03,-3.800E-04,-2.300E-04,
- 6 2.420E-03, 9.300E-04,-1.000E-04, 4.500E-04, 1.700E-04, 1.100E-04,
- 1 9.868E-01,-1.171E+00, 3.940E-01,-8.459E-02, 1.124E-02,-1.250E-03,
- 2 1.001E+00,-1.383E+00, 5.044E-01,-1.152E-01, 1.658E-02,-1.830E-03,
- 3 4.928E-02,-1.368E-01, 9.021E-02,-2.935E-02, 5.800E-03,-6.600E-04,
- 4-3.133E-02, 5.785E-02,-2.023E-02, 2.630E-03, 1.600E-04,-8.000E-05,
- 5 8.840E-03,-1.416E-02, 1.900E-03, 5.800E-04,-2.500E-04, 1.000E-05,
- 6-2.300E-03, 3.080E-03, 5.500E-04,-3.700E-04, 7.000E-05, 1.000E-05/
-C...Expansion coefficients for charm sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,6),IX=1,6),IT=1,6),NX=1,2)/
- 1 9.980E-03,-1.945E-02, 1.055E-02,-6.870E-03, 1.860E-03,-1.560E-03,
- 2 5.700E-03,-1.203E-02, 6.250E-03,-4.860E-03, 1.310E-03,-1.370E-03,
- 3-4.490E-03, 7.990E-03,-4.170E-03, 2.050E-03,-4.400E-04, 3.300E-04,
- 4 1.470E-03,-2.480E-03, 1.460E-03,-5.700E-04, 1.200E-04,-1.000E-05,
- 5-9.000E-05, 1.500E-04,-3.200E-04, 1.200E-04,-6.000E-05,-4.000E-05,
- 6-4.200E-04, 7.600E-04,-1.400E-04, 4.000E-05, 7.000E-05, 5.000E-05,
- 1 8.698E-01,-1.131E+00, 3.836E-01,-8.111E-02, 1.048E-02,-1.300E-03,
- 2 9.626E-01,-1.321E+00, 4.854E-01,-1.091E-01, 1.583E-02,-1.700E-03,
- 3 3.057E-02,-1.088E-01, 8.022E-02,-2.676E-02, 5.590E-03,-5.600E-04,
- 4-2.845E-02, 5.164E-02,-1.918E-02, 2.210E-03,-4.000E-05,-1.500E-04,
- 5 1.311E-02,-1.751E-02, 3.310E-03, 5.100E-04,-1.200E-04, 5.000E-05,
- 6-8.590E-03, 8.380E-03,-9.200E-04,-2.600E-04, 1.000E-05,-1.000E-05/
-C...Expansion coefficients for bottom sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,7),IX=1,6),IT=1,6),NX=1,2)/
- 1 8.980E-03,-1.459E-02, 7.510E-03,-4.410E-03, 1.310E-03,-1.070E-03,
- 2 5.970E-03,-9.440E-03, 4.800E-03,-3.020E-03, 9.100E-04,-8.500E-04,
- 3-3.050E-03, 4.440E-03,-2.100E-03, 8.500E-04,-2.400E-04, 1.400E-04,
- 4 5.300E-04,-1.300E-03, 5.600E-04,-2.700E-04, 3.000E-05,-2.000E-05,
- 5 2.000E-04, 1.400E-04,-1.100E-04, 1.000E-04, 0.000E+00, 0.000E+00,
- 6-2.600E-04, 3.200E-04, 0.000E+00,-3.000E-05, 1.000E-05,-1.000E-05,
- 1 8.672E-01,-1.174E+00, 4.265E-01,-9.252E-02, 1.244E-02,-1.460E-03,
- 2 8.500E-01,-1.194E+00, 4.630E-01,-1.083E-01, 1.614E-02,-1.830E-03,
- 3-2.241E-02,-5.630E-03, 2.815E-02,-1.425E-02, 3.520E-03,-4.300E-04,
- 4-7.300E-04, 8.030E-03,-5.780E-03, 1.380E-03,-1.300E-04,-4.000E-05,
- 5 3.460E-03,-5.380E-03, 1.960E-03,-2.100E-04, 1.000E-05, 1.000E-05,
- 6-4.850E-03, 5.950E-03,-1.890E-03, 2.600E-04,-3.000E-05, 0.000E+00/
-C...Expansion coefficients for top sea quark distribution.
- DATA (((CEHLQ(IX,IT,NX,8),IX=1,6),IT=1,6),NX=1,2)/
- 1 4.260E-03,-7.530E-03, 3.830E-03,-2.680E-03, 7.600E-04,-7.300E-04,
- 2 3.640E-03,-6.050E-03, 3.030E-03,-2.090E-03, 5.900E-04,-6.000E-04,
- 3-9.200E-04, 1.710E-03,-8.200E-04, 5.000E-04,-1.200E-04, 1.000E-04,
- 4-5.000E-05,-1.600E-04, 1.300E-04,-9.000E-05, 3.000E-05, 0.000E+00,
- 5 1.300E-04,-2.100E-04,-1.000E-05,-2.000E-05,-2.000E-05,-1.000E-05,
- 6-8.000E-05, 1.800E-04,-5.000E-05, 2.000E-05, 0.000E+00, 0.000E+00,
- 1 7.146E-01,-1.007E+00, 3.932E-01,-9.246E-02, 1.366E-02,-1.540E-03,
- 2 6.856E-01,-9.828E-01, 3.977E-01,-9.795E-02, 1.540E-02,-1.790E-03,
- 3-3.053E-02, 2.758E-02, 2.150E-03,-4.880E-03, 1.640E-03,-2.500E-04,
- 4 9.200E-04, 4.200E-04,-1.340E-03, 4.600E-04,-8.000E-05,-1.000E-05,
- 5 4.230E-03,-5.660E-03, 2.140E-03,-4.300E-04, 6.000E-05, 0.000E+00,
- 6-3.890E-03, 5.000E-03,-1.740E-03, 3.300E-04,-4.000E-05, 0.000E+00/
-C
-C
-C...Proton structure functions from Eichten, Hinchliffe, Lane, Quigg.
-C...Allowed variable range: 5 GeV^2 < Q^2 < 1E8 GeV^2; 1E-4 < x < 1
-C
-*...Set number of flavors - modify CERNLIB routine
- NFL=4
- IF (DSCALE.GT.BMAS) NFL=5
- IF (DSCALE.GT.TMAS) NFL=6
-C...Determine Lamdba and x and t expansion variables.
- X = DX
- SCALE = DSCALE
- KF = ABS(NFL)
- Q2 =SCALE*SCALE
- TMIN=LOG(Q02/ALAM**2)
- TMAX=LOG(1E8/ALAM**2)
- T=LOG(MAX(1.,Q2/ALAM**2))
- VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN)))
- NX=1
- IF(X.LE.0.1) NX=2
- IF(NX.EQ.1) VX=(2.*X-1.1)/0.9
- IF(NX.EQ.2) VX=MAX(-1.,(2.*LOG(X)+11.51293)/6.90776)
- CXS=1.
-
-C...Chebyshev polynomials for x and t expansion.
- TX(1)=1.
- TX(2)=VX
- TX(3)=2.*VX**2-1.
- TX(4)=4.*VX**3-3.*VX
- TX(5)=8.*VX**4-8.*VX**2+1.
- TX(6)=16.*VX**5-20.*VX**3+5.*VX
- TT(1)=1.
- TT(2)=VT
- TT(3)=2.*VT**2-1.
- TT(4)=4.*VT**3-3.*VT
- TT(5)=8.*VT**4-8.*VT**2+1.
- TT(6)=16.*VT**5-20.*VT**3+5.*VT
-
-C...Calculate structure functions.
- DO 120 KFL=1,6
- XQSUM=0.
- DO 110 IT=1,6
- DO 110 IX=1,6
- 110 XQSUM=XQSUM+CEHLQ(IX,IT,NX,KFL)*TX(IX)*TT(IT)
- 120 XQ(KFL)=XQSUM*(1.-X)**NEHLQ(KFL)*CXS
-C
-C...Put into output variables
- DUP = XQ(1)
- DDN = XQ(2)
- DSEA = XQ(3)
- DGL = XQ(4)
- DSTR = XQ(5)
- DCHM = XQ(6)
- DBOT = 0.D0
- DTOP = 0.D0
-C
- IF (KF.GE.5) THEN
-C
-C...Special expansion for bottom (threshold effects).
- BOT = 0.
- TMIN=7.4474
- IF(T.LE.TMIN) GOTO 140
- VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN)))
- TT(1)=1.
- TT(2)=VT
- TT(3)=2.*VT**2-1.
- TT(4)=4.*VT**3-3.*VT
- TT(5)=8.*VT**4-8.*VT**2+1.
- TT(6)=16.*VT**5-20.*VT**3+5.*VT
- XQSUM=0.
- DO 130 IT=1,6
- DO 130 IX=1,6
- 130 XQSUM=XQSUM+CEHLQ(IX,IT,NX,7)*TX(IX)*TT(IT)
- BOT=XQSUM*(1.-X)**NEHLQ(7)*CXS
- 140 CONTINUE
- ENDIF
- DBOT = BOT
-C--
- IF (KF.GE.6) THEN
-C
-C...Special expansion for top (threshold effects).
- TOP = 0.
- TMIN=10.8097
- TMIN=TMIN+2.*LOG(TMAS/30.)
- TMAX=TMAX+2.*LOG(TMAS/30.)
- IF(T.LE.TMIN) GOTO 160
- VT=MAX(-1.,MIN(1.,(2.*T-TMAX-TMIN)/(TMAX-TMIN)))
- TT(1)=1.
- TT(2)=VT
- TT(3)=2.*VT**2-1.
- TT(4)=4.*VT**3-3.*VT
- TT(5)=8.*VT**4-8.*VT**2+1.
- TT(6)=16.*VT**5-20.*VT**3+5.*VT
- XQSUM=0.
- DO 150 IT=1,6
- DO 150 IX=1,6
- 150 XQSUM=XQSUM+CEHLQ(IX,IT,NX,8)*TX(IX)*TT(IT)
- TOP=XQSUM*(1.-X)**NEHLQ(8)*CXS
- 160 CONTINUE
- ENDIF
- DTOP = TOP
-C
- RETURN
- END
diff --git a/gonsalves/mrs.f b/gonsalves/mrs.f
deleted file mode 100644
index 2bfa9e3..0000000
--- a/gonsalves/mrs.f
+++ /dev/null
@@ -1,459 +0,0 @@
-C-----------------------------------------------------------------------
-C Martin-Roberts-Stirling Parton Distributions
-C
- SUBROUTINE mrs (nset,x,qsq,pden)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION pden(-6:6)
- scale = dsqrt(qsq)
- mode = nset
- IF (mode.LE.3) THEN
- CALL mrs123 (x,scale,mode,upv,dnv,sea,str,chm,bot,glu)
- ELSEIF (mode.EQ.4.OR.mode.EQ.5) THEN
- CALL mrseb (x,scale,mode,upv,dnv,sea,str,chm,bot,glu)
- ELSE
- PRINT *,'Bad MRS mode ',nset
- STOP
- ENDIF
- pden(0) = glu/x
- pden(1) = (upv+sea)/x
- pden(2) = (dnv+sea)/x
- pden(-1) = sea/x
- pden(-2) = pden(-1)
- pden(3) = str/x
- pden(-3) = pden(3)
- pden(4) = chm/x
- pden(-4) = pden(4)
- pden(5) = bot/x
- pden(-5) = pden(5)
- pden(6) = 0d0
- pden(-6) = pden(6)
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE MRS123(X,SCALE,MODE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C***************************************************************C
-C C
-C MODE 1 CORRESPONDS TO C
-C MARTIN, ROBERTS, STIRLING (SOFT GLUE) WITH LAMBDA=.107 GEV C
-C C
-C MODE 2 CORRESPONDS TO C
-C MARTIN, ROBERTS, STIRLING (HARD GLUE) WITH LAMBDA=.250 GEV C
-C C
-C MODE 3 CORRESPONDS TO C
-C MARTIN, ROBERTS, STIRLING (1/RTX GLUE) WITH LAMBDA=.178 GEV C
-C C
-C -*- C
-C C
-C (NOTE THAT X TIMES THE PARTON DISTRIBUTION FUNCTION C
-C IS RETURNED I.E. G(X) = GLU/X ETC, AND THAT "SEA" C
-C IS THE LIGHT QUARK SEA I.E. UBAR(X)=DBAR(X)= ... C
-C = SEA/X FOR A PROTON. IF IN DOUBT, CHECK THE C
-C MOMENTUM SUM RULE! NOTE ALSO THAT SCALE=Q IN GEV) C
-C C
-C -*- C
-C C
-C (THE RANGE OF APPLICABILITY IS CURRENTLY: C
-C 10**-4 < X < 1 AND 5 < Q**2 < 1.31 * 10**6 C
-C HIGHER Q**2 VALUES CAN BE SUPPLIED ON REQUEST C
-C - PROBLEMS, COMMENTS ETC TO SRG$T3@GEN C
-C C
-C C
-C***************************************************************C
- IMPLICIT REAL*8(A-H,O-Z)
- IF(MODE.EQ.1) CALL STRUC1(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
- IF(MODE.EQ.2) CALL STRUC2(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
- IF(MODE.EQ.3) CALL STRUC3(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
- RETURN
- END
-C
- SUBROUTINE STRUC1(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C :::::::::::: MARTIN ROBERTS STIRLING :::::SOFT GLUE::::107 MEV::::
- IMPLICIT REAL*8(A-H,O-Z)
- DIMENSION F(6,42,19),G(6),XX(42),N0(6)
- SAVE
- DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
- . 1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
- . 1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
- . .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
- . .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
- . .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
- . .8D0,.9D0,1.D0/
- DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
- DATA N0/2,5,4,5,0,0/
- DATA INIT/0/
- IF(INIT.NE.0) GOTO 10
- OPEN(UNIT=31,FILE='mrs1.dat'
- & ,STATUS='OLD') !,READONLY)
- READ(31,*)
- INIT=1
- DO 20 N=1,41
- DO 20 M=1,19
- READ(31,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
- DO 25 I=1,6
- 25 F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
- 20 CONTINUE
- DO 30 J=1,15
- XX(J)=DLOG10(XX(J))+1.1D0
- DO 30 I=1,5
- DO 30 K=1,19
- 30 F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
- 50 FORMAT(6F10.5)
- DO 40 I=1,6
- DO 40 M=1,19
- 40 F(I,42,M)=0.D0
- CLOSE(31)
- 10 CONTINUE
- IF(X.LT.XMIN) X=XMIN
- IF(X.GT.XMAX) X=XMAX
- QSQ=SCALE**2
- IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
- IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
- XXX=X
- IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
- N=0
- 70 N=N+1
- IF(XXX.GT.XX(N+1)) GOTO 70
- A=(XXX-XX(N))/(XX(N+1)-XX(N))
- RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
- B=RM-DINT(RM)
- M=1+IDINT(RM)
- DO 60 I=1,6
- G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
- . + A*(1.D0-B)*F(I,N+1,M) + A*B*F(I,N+1,M+1)
- IF(N.GE.16) GOTO 65
- IF(I.EQ.6) GOTO 65
- FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
-
- G(I)=FAC**(G(I)/FAC)
- 65 CONTINUE
- G(I)=G(I)*(1.D0-X)**N0(I)
- 60 CONTINUE
- UPV=G(1)
- DNV=G(2)
- SEA=G(4)
- STR=G(4)
- CHM=G(5)
- GLU=G(3)
- BOT=G(6)
- RETURN
- END
- SUBROUTINE STRUC2(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C :::::::::::: MARTIN ROBERTS STIRLING :::::HARD GLUE::::250 MEV::::
- IMPLICIT REAL*8(A-H,O-Z)
- DIMENSION F(6,42,19),G(6),XX(42),N0(6)
- SAVE
- DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
- . 1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
- . 1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
- . .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
- . .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
- . .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
- . .8D0,.9D0,1.D0/
- DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
- DATA N0/2,5,4,5,0,0/
- DATA INIT/0/
- IF(INIT.NE.0) GOTO 10
- OPEN(UNIT=32,FILE='mrs2.dat'
- & ,STATUS='OLD') !,READONLY)
- READ(32,*)
- INIT=1
- DO 20 N=1,41
- DO 20 M=1,19
- READ(32,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
- DO 25 I=1,6
- 25 F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
- 20 CONTINUE
- DO 30 J=1,15
- XX(J)=DLOG10(XX(J))+1.1D0
- DO 30 I=1,5
- DO 30 K=1,19
- 30 F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
- 50 FORMAT(6F10.5)
- DO 40 I=1,6
- DO 40 M=1,19
- 40 F(I,42,M)=0.D0
- CLOSE(32)
- 10 CONTINUE
- IF(X.LT.XMIN) X=XMIN
- IF(X.GT.XMAX) X=XMAX
- QSQ=SCALE**2
- IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
- IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
- XXX=X
- IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
- N=0
- 70 N=N+1
- IF(XXX.GT.XX(N+1)) GOTO 70
- A=(XXX-XX(N))/(XX(N+1)-XX(N))
- RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
- B=RM-DINT(RM)
- M=1+IDINT(RM)
- DO 60 I=1,6
- G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
- . + A*(1.D0-B)*F(I,N+1,M) + A*B*F(I,N+1,M+1)
- IF(N.GE.16) GOTO 65
- IF(I.EQ.6) GOTO 65
- FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
-
- G(I)=FAC**(G(I)/FAC)
- 65 CONTINUE
- G(I)=G(I)*(1.D0-X)**N0(I)
- 60 CONTINUE
- UPV=G(1)
- DNV=G(2)
- SEA=G(4)
- STR=G(4)
- CHM=G(5)
- GLU=G(3)
- BOT=G(6)
- RETURN
- END
- SUBROUTINE STRUC3(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C :::::::::::: MARTIN ROBERTS STIRLING :::::1/RX GLUE::::178 MEV::::
- IMPLICIT REAL*8(A-H,O-Z)
- DIMENSION F(6,42,19),G(6),XX(42),N0(6)
- SAVE
- DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
- . 1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
- . 1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
- . .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
- . .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
- . .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
- . .8D0,.9D0,1.D0/
- DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
- DATA N0/3,4,5,5,0,0/
- DATA INIT/0/
- IF(INIT.NE.0) GOTO 10
- OPEN(UNIT=33,FILE='mrs3.dat'
- & ,STATUS='OLD') !,READONLY)
- READ(33,*)
- INIT=1
- DO 20 N=1,41
- DO 20 M=1,19
- READ(33,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
- DO 25 I=1,6
- 25 F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
- 20 CONTINUE
- DO 30 J=1,15
- XX(J)=DLOG10(XX(J))+1.1D0
- DO 30 I=1,5
- DO 30 K=1,19
- 30 F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
- 50 FORMAT(6F10.5)
- DO 40 I=1,6
- DO 40 M=1,19
- 40 F(I,42,M)=0.D0
- CLOSE(33)
- 10 CONTINUE
- IF(X.LT.XMIN) X=XMIN
- IF(X.GT.XMAX) X=XMAX
- QSQ=SCALE**2
- IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
- IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
- XXX=X
- IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
- N=0
- 70 N=N+1
- IF(XXX.GT.XX(N+1)) GOTO 70
- A=(XXX-XX(N))/(XX(N+1)-XX(N))
- RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
- B=RM-DINT(RM)
- M=1+IDINT(RM)
- DO 60 I=1,6
- G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
- . + A*(1.D0-B)*F(I,N+1,M) + A*B*F(I,N+1,M+1)
- IF(N.GE.16) GOTO 65
- IF(I.EQ.6) GOTO 65
- FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
-
- G(I)=FAC**(G(I)/FAC)
- 65 CONTINUE
- G(I)=G(I)*(1.D0-X)**N0(I)
- 60 CONTINUE
- UPV=G(1)
- DNV=G(2)
- SEA=G(4)
- STR=G(4)
- CHM=G(5)
- GLU=G(3)
- BOT=G(6)
- RETURN
- END
-C Received from J. Stirling February 1989
- SUBROUTINE MRSEB(X,SCALE,MODE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C***************************************************************C
-C C
-C C
-C NEW VERSIONS !!!! JANUARY 1989 (AS DESCRIBED IN C
-C "IMPROVED PARTON DISTRIBUTIONS ... " A.D. MARTIN, C
-C R.G. ROBERTS AND W.J. STIRLING PREPRINT RAL-88-113 ) C
-C C
-C MODE 1 CORRESPONDS TO C
-C MARTIN, ROBERTS, STIRLING (EMC FIT) WITH LAMBDA= 100 MEV C
-C C
-C MODE 2 CORRESPONDS TO C
-C MARTIN, ROBERTS, STIRLING (BCDMS FIT) WITH LAMBDA= 200 MEV C
-C C
-C ( SOFT GLUE : X G(X,Q0) = A (1-X)**4.4 ) C
-C C
-C -*- C
-C C
-C (NOTE THAT X TIMES THE PARTON DISTRIBUTION FUNCTION C
-C IS RETURNED I.E. G(X) = GLU/X ETC, AND THAT "SEA" C
-C IS THE LIGHT QUARK SEA I.E. UBAR(X)=DBAR(X)= ... C
-C = SEA/X FOR A PROTON. IF IN DOUBT, CHECK THE C
-C MOMENTUM SUM RULE! NOTE ALSO THAT SCALE=Q IN GEV) C
-C C
-C -*- C
-C C
-C (THE RANGE OF APPLICABILITY IS CURRENTLY: C
-C 10**-4 < X < 1 AND 5 < Q**2 < 1.31 * 10**6 C
-C HIGHER Q**2 VALUES CAN BE SUPPLIED ON REQUEST C
-C - PROBLEMS, COMMENTS ETC TO SRG$T3@GEN C
-C C
-C C
-C***************************************************************C
- IMPLICIT REAL*8(A-H,O-Z)
- IF(MODE.EQ.4) CALL STRUCE(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
- IF(MODE.EQ.5) CALL STRUCB(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
- RETURN
- END
-C
- SUBROUTINE STRUCE(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C :::::::::::: MARTIN ROBERTS STIRLING :::::SOFT GLUE:::: 91 MEV:::
- IMPLICIT REAL*8(A-H,O-Z)
- DIMENSION F(6,42,19),G(6),XX(42),N0(6)
- DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
- . 1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
- . 1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
- . .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
- . .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
- . .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
- . .8D0,.9D0,1.D0/
- DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
- DATA N0/2,5,4,5,0,0/
- DATA INIT/0/
- IF(INIT.NE.0) GOTO 10
- OPEN(UNIT=34,FILE='mrse.dat'
- & ,STATUS='OLD') !,READONLY)
- READ(34,*)
- INIT=1
- DO 20 N=1,41
- DO 20 M=1,19
- READ(34,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
- DO 25 I=1,6
- 25 F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
- 20 CONTINUE
- DO 30 J=1,15
- XX(J)=DLOG10(XX(J))+1.1D0
- DO 30 I=1,5
- DO 30 K=1,19
- 30 F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
- 50 FORMAT(6F10.5)
- DO 40 I=1,6
- DO 40 M=1,19
- 40 F(I,42,M)=0.D0
- CLOSE (34)
- 10 CONTINUE
- IF(X.LT.XMIN) X=XMIN
- IF(X.GT.XMAX) X=XMAX
- QSQ=SCALE**2
- IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
- IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
- XXX=X
- IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
- N=0
- 70 N=N+1
- IF(XXX.GT.XX(N+1)) GOTO 70
- A=(XXX-XX(N))/(XX(N+1)-XX(N))
- RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
- B=RM-DINT(RM)
- M=1+IDINT(RM)
- DO 60 I=1,6
- G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
- . + A*(1.D0-B)*F(I,N+1,M) + A*B*F(I,N+1,M+1)
- IF(N.GE.16) GOTO 65
- IF(I.EQ.6) GOTO 65
- FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
- G(I)=FAC**(G(I)/FAC)
- 65 CONTINUE
- G(I)=G(I)*(1.D0-X)**N0(I)
- 60 CONTINUE
- UPV=G(1)
- DNV=G(2)
- SEA=G(4)
- STR=G(4)
- CHM=G(5)
- GLU=G(3)
- BOT=G(6)
- RETURN
- END
- SUBROUTINE STRUCB(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
-C :::::::::::: MARTIN ROBERTS STIRLING :::::SOFT GLUE::::228 MEV:::
- IMPLICIT REAL*8(A-H,O-Z)
- DIMENSION F(6,42,19),G(6),XX(42),N0(6)
- DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
- . 1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
- . 1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
- . .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
- . .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
- . .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
- . .8D0,.9D0,1.D0/
- DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
- DATA N0/2,5,4,5,0,0/
- DATA INIT/0/
- IF(INIT.NE.0) GOTO 10
- OPEN(UNIT=35,FILE='mrsb.dat'
- & ,STATUS='OLD') !,READONLY)
- READ(35,*)
- INIT=1
- DO 20 N=1,41
- DO 20 M=1,19
- READ(35,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
- DO 25 I=1,6
- 25 F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
- 20 CONTINUE
- DO 30 J=1,15
- XX(J)=DLOG10(XX(J))+1.1D0
- DO 30 I=1,5
- DO 30 K=1,19
- 30 F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
- 50 FORMAT(6F10.5)
- DO 40 I=1,6
- DO 40 M=1,19
- 40 F(I,42,M)=0.D0
- CLOSE (35)
- 10 CONTINUE
- IF(X.LT.XMIN) X=XMIN
- IF(X.GT.XMAX) X=XMAX
- QSQ=SCALE**2
- IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
- IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
- XXX=X
- IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
- N=0
- 70 N=N+1
- IF(XXX.GT.XX(N+1)) GOTO 70
- A=(XXX-XX(N))/(XX(N+1)-XX(N))
- RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
- B=RM-DINT(RM)
- M=1+IDINT(RM)
- DO 60 I=1,6
- G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
- . + A*(1.D0-B)*F(I,N+1,M) + A*B*F(I,N+1,M+1)
- IF(N.GE.16) GOTO 65
- IF(I.EQ.6) GOTO 65
- FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
- G(I)=FAC**(G(I)/FAC)
- 65 CONTINUE
- G(I)=G(I)*(1.D0-X)**N0(I)
- 60 CONTINUE
- UPV=G(1)
- DNV=G(2)
- SEA=G(4)
- STR=G(4)
- CHM=G(5)
- GLU=G(3)
- BOT=G(6)
- RETURN
- END
diff --git a/gonsalves/mrst.f b/gonsalves/mrst.f
deleted file mode 100644
index 364c408..0000000
--- a/gonsalves/mrst.f
+++ /dev/null
@@ -1,34 +0,0 @@
-* Martin Roberts Stirling Thorne Parton Densities
- SUBROUTINE mrst(nset,x,qsq,pden)
- IMPLICIT NONE
- INTEGER nset,mode
- DOUBLE PRECISION x,qsq,pden(-6:6)
- DOUBLE PRECISION q,upv,dnv,usea,dsea,str,chm,bot,glu
-
- q=dsqrt(qsq)
- IF (nset.EQ.0) THEN
- mode = 1
- CALL mrstlo(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
- ELSEIF (nset.EQ.1.OR.nset.EQ.2) THEN
- mode=nset
- CALL mrst2002(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
- ELSE
- PRINT *,' bad MRST mode number',nset
- STOP
- ENDIF
- pden(0) = glu/x
- pden(1) = (upv+usea)/x
- pden(-1) = usea/x
- pden(2) = (dnv+dsea)/x
- pden(-2) = dsea/x
- pden(3) = str/x
- pden(-3) = pden(3)
- pden(4) = chm/x
- pden(-4) = pden(4)
- pden(5) = bot/x
- pden(-5) = pden(5)
- pden(6) = 0d0
- pden(-6) = pden(6)
-
- RETURN
- END
diff --git a/gonsalves/mrst2001lo.f b/gonsalves/mrst2001lo.f
deleted file mode 100644
index 7a162c5..0000000
--- a/gonsalves/mrst2001lo.f
+++ /dev/null
@@ -1,307 +0,0 @@
- subroutine mrstlo(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
-C***************************************************************C
-C C
-C This is a package for the new MRST 2001 LO parton C
-C distributions. C
-C Reference: A.D. Martin, R.G. Roberts, W.J. Stirling and C
-C R.S. Thorne, hep-ph/0201xxx C
-C C
-C There is 1 pdf set corresponding to mode = 1 C
-C C
-C Mode=1 gives the default set with Lambda(MSbar,nf=4) = 0.220 C
-C corresponding to alpha_s(M_Z) of 0.130 C
-C This set reads a grid whose first number is 0.02868 C
-C C
-C This subroutine uses an improved interpolation procedure C
-C for extracting values of the pdf's from the grid C
-C C
-C Comments to : W.J.Stirling@durham.ac.uk C
-C C
-C***************************************************************C
- implicit real*8(a-h,o-z)
- data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
- q2=q*q
- if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
- if(x.lt.xmin.or.x.gt.xmax) print 98,x
- if(mode.eq.1) then
- call mrst1lo(x,q2,upv,dnv,usea,dsea,str,chm,bot,glu)
- endif
- 99 format(' WARNING: Q^2 VALUE IS OUT OF RANGE ','q2= ',e10.5)
- 98 format(' WARNING: X VALUE IS OUT OF RANGE ','x= ',e10.5)
- return
- end
-
- subroutine mrst1lo(x,qsq,upv,dnv,usea,dsea,str,chm,bot,glu)
- implicit real*8(a-h,o-z)
- parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26)
- real*8 f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq),
- .f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb)
- real*8 qq(nq),xx(nx),cc1(nx,nq,4,4),cc2(nx,nq,4,4),
- .cc3(nx,nq,4,4),cc4(nx,nq,4,4),cc6(nx,nq,4,4),cc8(nx,nq,4,4),
- .ccc(nx,nqc,4,4),ccb(nx,nqb,4,4)
- real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
- data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
- . 1d-4,2d-4,4d-4,6d-4,8d-4,
- . 1d-3,2d-3,4d-3,6d-3,8d-3,
- . 1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
- . .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
- . .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
- . .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
- . .8d0,.9d0,1d0/
- data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
- . 1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
- . 1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
- . 1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
- . 1.8d6,3.2d6,5.6d6,1d7/
- data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
- data init/0/
- save
- xsave=x
- q2save=qsq
- if(init.ne.0) goto 10
- open(unit=33,file='lo2002.dat',status='old')
- do 20 n=1,nx-1
- do 20 m=1,nq
- read(33,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
- . f5(n,m),f7(n,m),f6(n,m),f8(n,m)
-c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
- 20 continue
- do 40 m=1,nq
- f1(nx,m)=0.d0
- f2(nx,m)=0.d0
- f3(nx,m)=0.d0
- f4(nx,m)=0.d0
- f5(nx,m)=0.d0
- f6(nx,m)=0.d0
- f7(nx,m)=0.d0
- f8(nx,m)=0.d0
- 40 continue
- do n=1,nx
- xxl(n)=dlog(xx(n))
- enddo
- do m=1,nq
- qql(m)=dlog(qq(m))
- enddo
-
- call jeppe1lo(nx,nq,xxl,qql,f1,cc1)
- call jeppe1lo(nx,nq,xxl,qql,f2,cc2)
- call jeppe1lo(nx,nq,xxl,qql,f3,cc3)
- call jeppe1lo(nx,nq,xxl,qql,f4,cc4)
- call jeppe1lo(nx,nq,xxl,qql,f6,cc6)
- call jeppe1lo(nx,nq,xxl,qql,f8,cc8)
-
- emc2=2.045
- emb2=18.5
-
- do 44 m=1,nqc
- qqlc(m)=qql(m+nqc0)
- do 44 n=1,nx
- fc(n,m)=f5(n,m+nqc0)
- 44 continue
- qqlc(1)=dlog(emc2)
- call jeppe1lo(nx,nqc,xxl,qqlc,fc,ccc)
-
- do 45 m=1,nqb
- qqlb(m)=qql(m+nqb0)
- do 45 n=1,nx
- fb(n,m)=f7(n,m+nqb0)
- 45 continue
- qqlb(1)=dlog(emb2)
- call jeppe1lo(nx,nqb,xxl,qqlb,fb,ccb)
-
-
- init=1
- 10 continue
-
- xlog=dlog(x)
- qsqlog=dlog(qsq)
-
- call jeppe2lo(xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
- call jeppe2lo(xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
- call jeppe2lo(xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
- call jeppe2lo(xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
- call jeppe2lo(xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
- call jeppe2lo(xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
-
- chm=0.d0
- if(qsq.gt.emc2) then
- call jeppe2lo(xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
- endif
-
- bot=0.d0
- if(qsq.gt.emb2) then
- call jeppe2lo(xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
- endif
-
- x=xsave
- qsq=q2save
- return
- 50 format(8f10.5)
- end
-
-c subroutine jeppe1lo(nx,my,xx,yy,ff,cc)
-c implicit real*8(a-h,o-z)
-c dimension xx(nx),yy(my),ff(nx,my),ff1(nx,my),ff2(nx,my),
-c xff12(nx,my),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
-c xcl(16),cc(nx,my,4,4),iwt(16,16)
-
- subroutine jeppe1lo(nx,my,xx,yy,ff,cc)
- implicit real*8(a-h,o-z)
- PARAMETER(NNX=49,MMY=37)
- dimension xx(nx),yy(my),ff(nx,my),ff1(NNX,MMY),ff2(NNX,MMY),
- xff12(NNX,MMY),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
- xcl(16),cc(nx,my,4,4),iwt(16,16)
-
- data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
- x -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
- x 2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
- x 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
- x 0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
- x 0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
- x -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
- x 9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
- x -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
- x 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
- x -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
- x 4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
-
-
- do 42 m=1,my
- dx=xx(2)-xx(1)
- ff1(1,m)=(ff(2,m)-ff(1,m))/dx
- dx=xx(nx)-xx(nx-1)
- ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
- do 41 n=2,nx-1
- ff1(n,m)=polderivlo(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
- xff(n+1,m))
- 41 continue
- 42 continue
-
- do 44 n=1,nx
- dy=yy(2)-yy(1)
- ff2(n,1)=(ff(n,2)-ff(n,1))/dy
- dy=yy(my)-yy(my-1)
- ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
- do 43 m=2,my-1
- ff2(n,m)=polderivlo(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
- xff(n,m+1))
- 43 continue
- 44 continue
-
- do 46 m=1,my
- dx=xx(2)-xx(1)
- ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
- dx=xx(nx)-xx(nx-1)
- ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
- do 45 n=2,nx-1
- ff12(n,m)=polderivlo(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
- xff2(n+1,m))
- 45 continue
- 46 continue
-
- do 53 n=1,nx-1
- do 52 m=1,my-1
- d1=xx(n+1)-xx(n)
- d2=yy(m+1)-yy(m)
- d1d2=d1*d2
-
- yy0(1)=ff(n,m)
- yy0(2)=ff(n+1,m)
- yy0(3)=ff(n+1,m+1)
- yy0(4)=ff(n,m+1)
-
- yy1(1)=ff1(n,m)
- yy1(2)=ff1(n+1,m)
- yy1(3)=ff1(n+1,m+1)
- yy1(4)=ff1(n,m+1)
-
- yy2(1)=ff2(n,m)
- yy2(2)=ff2(n+1,m)
- yy2(3)=ff2(n+1,m+1)
- yy2(4)=ff2(n,m+1)
-
- yy12(1)=ff12(n,m)
- yy12(2)=ff12(n+1,m)
- yy12(3)=ff12(n+1,m+1)
- yy12(4)=ff12(n,m+1)
-
- do 47 k=1,4
- z(k)=yy0(k)
- z(k+4)=yy1(k)*d1
- z(k+8)=yy2(k)*d2
- z(k+12)=yy12(k)*d1d2
- 47 continue
-
- do 49 l=1,16
- xxd=0.
- do 48 k=1,16
- xxd=xxd+iwt(k,l)*z(k)
- 48 continue
- cl(l)=xxd
- 49 continue
- l=0
- do 51 k=1,4
- do 50 j=1,4
- l=l+1
- cc(n,m,k,j)=cl(l)
- 50 continue
- 51 continue
- 52 continue
- 53 continue
- return
- end
-
- subroutine jeppe2lo(x,y,nx,my,xx,yy,cc,z)
- implicit real*8(a-h,o-z)
- dimension xx(nx),yy(my),cc(nx,my,4,4)
-
- n=locxlo(xx,nx,x)
- m=locxlo(yy,my,y)
-
- t=(x-xx(n))/(xx(n+1)-xx(n))
- u=(y-yy(m))/(yy(m+1)-yy(m))
-
- z=0.
- do 1 l=4,1,-1
- z=t*z+((cc(n,m,l,4)*u+cc(n,m,l,3))*u
- . +cc(n,m,l,2))*u+cc(n,m,l,1)
- 1 continue
- return
- end
-
- integer function locxlo(xx,nx,x)
- implicit real*8(a-h,o-z)
- dimension xx(nx)
- if(x.le.xx(1)) then
- locxlo=1
- return
- endif
- if(x.ge.xx(nx)) then
- locxlo=nx-1
- return
- endif
- ju=nx+1
- jl=0
- 1 if((ju-jl).le.1) go to 2
- jm=(ju+jl)/2
- if(x.ge.xx(jm)) then
- jl=jm
- else
- ju=jm
- endif
- go to 1
- 2 locxlo=jl
- return
- end
-
-
- real*8 function polderivlo(x1,x2,x3,y1,y2,y3)
- implicit real*8(a-h,o-z)
- polderivlo=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
- .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
- return
- end
diff --git a/gonsalves/mrst2002.f b/gonsalves/mrst2002.f
deleted file mode 100644
index 8ccef09..0000000
--- a/gonsalves/mrst2002.f
+++ /dev/null
@@ -1,409 +0,0 @@
- subroutine mrst2002(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
-C***************************************************************C
-C C
-C This is a package for the new MRST 2002 updated NLO and C
-C NNLO parton distributions. C
-C Reference: A.D. Martin, R.G. Roberts, W.J. Stirling and C
-C R.S. Thorne, hep-ph/0211080 C
-C C
-C There are 2 pdf sets corresponding to mode = 1, 2 C
-C C
-C Mode=1 gives the NLO set with alpha_s(M_Z,NLO) = 0.1197 C
-C This set reads a grid whose first number is 0.00949 C
-C C
-C Mode=2 gives the NNLO set with alpha_s(M_Z,NNLO) = 0.1154 C
-C This set reads a grid whose first number is 0.00685 C
-C C
-C Comments to : W.J.Stirling@durham.ac.uk C
-C C
-C***************************************************************C
- implicit real*8(a-h,o-z)
- data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
- q2=q*q
- if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
- if(x.lt.xmin.or.x.gt.xmax) print 98,x
- if(mode.eq.1) then
- call mrst1(x,q2,upv,dnv,usea,dsea,str,chm,bot,glu)
- elseif(mode.eq.2) then
- call mrst2(x,q2,upv,dnv,usea,dsea,str,chm,bot,glu)
- endif
- 99 format(' WARNING: Q^2 VALUE IS OUT OF RANGE ','q2= ',e10.5)
- 98 format(' WARNING: X VALUE IS OUT OF RANGE ','x= ',e10.5)
- return
- end
-
- subroutine mrst1(x,qsq,upv,dnv,usea,dsea,str,chm,bot,glu)
- implicit real*8(a-h,o-z)
- parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26)
- real*8 f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq),
- .f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb)
- real*8 qq(nq),xx(nx),cc1(nx,nq,4,4),cc2(nx,nq,4,4),
- .cc3(nx,nq,4,4),cc4(nx,nq,4,4),cc6(nx,nq,4,4),cc8(nx,nq,4,4),
- .ccc(nx,nqc,4,4),ccb(nx,nqb,4,4)
- real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
- data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
- . 1d-4,2d-4,4d-4,6d-4,8d-4,
- . 1d-3,2d-3,4d-3,6d-3,8d-3,
- . 1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
- . .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
- . .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
- . .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
- . .8d0,.9d0,1d0/
- data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
- . 1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
- . 1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
- . 1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
- . 1.8d6,3.2d6,5.6d6,1d7/
- data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
- data init/0/
- save
- xsave=x
- q2save=qsq
- if(init.ne.0) goto 10
- open(unit=33,file='mrst2002nlo.dat',status='old')
- do 20 n=1,nx-1
- do 20 m=1,nq
- read(33,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
- . f5(n,m),f7(n,m),f6(n,m),f8(n,m)
-c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
- 20 continue
- do 40 m=1,nq
- f1(nx,m)=0.d0
- f2(nx,m)=0.d0
- f3(nx,m)=0.d0
- f4(nx,m)=0.d0
- f5(nx,m)=0.d0
- f6(nx,m)=0.d0
- f7(nx,m)=0.d0
- f8(nx,m)=0.d0
- 40 continue
- do n=1,nx
- xxl(n)=dlog(xx(n))
- enddo
- do m=1,nq
- qql(m)=dlog(qq(m))
- enddo
-
- call jeppe1(nx,nq,xxl,qql,f1,cc1)
- call jeppe1(nx,nq,xxl,qql,f2,cc2)
- call jeppe1(nx,nq,xxl,qql,f3,cc3)
- call jeppe1(nx,nq,xxl,qql,f4,cc4)
- call jeppe1(nx,nq,xxl,qql,f6,cc6)
- call jeppe1(nx,nq,xxl,qql,f8,cc8)
-
- emc2=2.045
- emb2=18.5
-
- do 44 m=1,nqc
- qqlc(m)=qql(m+nqc0)
- do 44 n=1,nx
- fc(n,m)=f5(n,m+nqc0)
- 44 continue
- qqlc(1)=dlog(emc2)
- call jeppe1(nx,nqc,xxl,qqlc,fc,ccc)
-
- do 45 m=1,nqb
- qqlb(m)=qql(m+nqb0)
- do 45 n=1,nx
- fb(n,m)=f7(n,m+nqb0)
- 45 continue
- qqlb(1)=dlog(emb2)
- call jeppe1(nx,nqb,xxl,qqlb,fb,ccb)
-
-
- init=1
- 10 continue
-
- xlog=dlog(x)
- qsqlog=dlog(qsq)
-
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
-
- chm=0.d0
- if(qsq.gt.emc2) then
- call jeppe2(xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
- endif
-
- bot=0.d0
- if(qsq.gt.emb2) then
- call jeppe2(xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
- endif
-
- x=xsave
- qsq=q2save
- return
- 50 format(8f10.5)
- end
-
- subroutine mrst2(x,qsq,upv,dnv,usea,dsea,str,chm,bot,glu)
- implicit real*8(a-h,o-z)
- parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26)
- real*8 f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq),
- .f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb)
- real*8 qq(nq),xx(nx),cc1(nx,nq,4,4),cc2(nx,nq,4,4),
- .cc3(nx,nq,4,4),cc4(nx,nq,4,4),cc6(nx,nq,4,4),cc8(nx,nq,4,4),
- .ccc(nx,nqc,4,4),ccb(nx,nqb,4,4)
- real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
- data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
- . 1d-4,2d-4,4d-4,6d-4,8d-4,
- . 1d-3,2d-3,4d-3,6d-3,8d-3,
- . 1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
- . .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
- . .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
- . .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
- . .8d0,.9d0,1d0/
- data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
- . 1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
- . 1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
- . 1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
- . 1.8d6,3.2d6,5.6d6,1d7/
- data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
- data init/0/
- save
- xsave=x
- q2save=qsq
- if(init.ne.0) goto 10
- open(unit=33,file='mrst2002nnlo.dat',status='old')
- do 20 n=1,nx-1
- do 20 m=1,nq
- read(33,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
- . f5(n,m),f7(n,m),f6(n,m),f8(n,m)
-c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
- 20 continue
- do 40 m=1,nq
- f1(nx,m)=0.d0
- f2(nx,m)=0.d0
- f3(nx,m)=0.d0
- f4(nx,m)=0.d0
- f5(nx,m)=0.d0
- f6(nx,m)=0.d0
- f7(nx,m)=0.d0
- f8(nx,m)=0.d0
- 40 continue
- do n=1,nx
- xxl(n)=dlog(xx(n))
- enddo
- do m=1,nq
- qql(m)=dlog(qq(m))
- enddo
-
- call jeppe1(nx,nq,xxl,qql,f1,cc1)
- call jeppe1(nx,nq,xxl,qql,f2,cc2)
- call jeppe1(nx,nq,xxl,qql,f3,cc3)
- call jeppe1(nx,nq,xxl,qql,f4,cc4)
- call jeppe1(nx,nq,xxl,qql,f6,cc6)
- call jeppe1(nx,nq,xxl,qql,f8,cc8)
-
- emc2=2.045
- emb2=18.5
-
- do 44 m=1,nqc
- qqlc(m)=qql(m+nqc0)
- do 44 n=1,nx
- fc(n,m)=f5(n,m+nqc0)
- 44 continue
- qqlc(1)=dlog(emc2)
- call jeppe1(nx,nqc,xxl,qqlc,fc,ccc)
-
- do 45 m=1,nqb
- qqlb(m)=qql(m+nqb0)
- do 45 n=1,nx
- fb(n,m)=f7(n,m+nqb0)
- 45 continue
- qqlb(1)=dlog(emb2)
- call jeppe1(nx,nqb,xxl,qqlb,fb,ccb)
-
-
- init=1
- 10 continue
-
- xlog=dlog(x)
- qsqlog=dlog(qsq)
-
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
- call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
-
- chm=0.d0
- if(qsq.gt.emc2) then
- call jeppe2(xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
- endif
-
- bot=0.d0
- if(qsq.gt.emb2) then
- call jeppe2(xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
- endif
-
- x=xsave
- qsq=q2save
- return
- 50 format(8f10.5)
- end
- subroutine jeppe1(nx,my,xx,yy,ff,cc)
- implicit real*8(a-h,o-z)
- parameter(nnx=49,mmy=37)
- dimension xx(nx),yy(my),ff(nnx,mmy),ff1(nnx,mmy),ff2(nnx,mmy),
- xff12(nnx,mmy),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
- xcl(16),cc(nx,my,4,4),iwt(16,16)
-
- data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
- x -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
- x 2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
- x 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
- x 0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
- x 0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
- x -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
- x 9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
- x -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
- x 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
- x 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
- x -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
- x 4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
-
-
- do 42 m=1,my
- dx=xx(2)-xx(1)
- ff1(1,m)=(ff(2,m)-ff(1,m))/dx
- dx=xx(nx)-xx(nx-1)
- ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
- do 41 n=2,nx-1
- ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
- xff(n+1,m))
- 41 continue
- 42 continue
-
- do 44 n=1,nx
- dy=yy(2)-yy(1)
- ff2(n,1)=(ff(n,2)-ff(n,1))/dy
- dy=yy(my)-yy(my-1)
- ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
- do 43 m=2,my-1
- ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
- xff(n,m+1))
- 43 continue
- 44 continue
-
- do 46 m=1,my
- dx=xx(2)-xx(1)
- ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
- dx=xx(nx)-xx(nx-1)
- ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
- do 45 n=2,nx-1
- ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
- xff2(n+1,m))
- 45 continue
- 46 continue
-
- do 53 n=1,nx-1
- do 52 m=1,my-1
- d1=xx(n+1)-xx(n)
- d2=yy(m+1)-yy(m)
- d1d2=d1*d2
-
- yy0(1)=ff(n,m)
- yy0(2)=ff(n+1,m)
- yy0(3)=ff(n+1,m+1)
- yy0(4)=ff(n,m+1)
-
- yy1(1)=ff1(n,m)
- yy1(2)=ff1(n+1,m)
- yy1(3)=ff1(n+1,m+1)
- yy1(4)=ff1(n,m+1)
-
- yy2(1)=ff2(n,m)
- yy2(2)=ff2(n+1,m)
- yy2(3)=ff2(n+1,m+1)
- yy2(4)=ff2(n,m+1)
-
- yy12(1)=ff12(n,m)
- yy12(2)=ff12(n+1,m)
- yy12(3)=ff12(n+1,m+1)
- yy12(4)=ff12(n,m+1)
-
- do 47 k=1,4
- z(k)=yy0(k)
- z(k+4)=yy1(k)*d1
- z(k+8)=yy2(k)*d2
- z(k+12)=yy12(k)*d1d2
- 47 continue
-
- do 49 l=1,16
- xxd=0.
- do 48 k=1,16
- xxd=xxd+iwt(k,l)*z(k)
- 48 continue
- cl(l)=xxd
- 49 continue
- l=0
- do 51 k=1,4
- do 50 j=1,4
- l=l+1
- cc(n,m,k,j)=cl(l)
- 50 continue
- 51 continue
- 52 continue
- 53 continue
- return
- end
-
- subroutine jeppe2(x,y,nx,my,xx,yy,cc,z)
- implicit real*8(a-h,o-z)
- dimension xx(nx),yy(my),cc(nx,my,4,4)
-
- n=locx(xx,nx,x)
- m=locx(yy,my,y)
-
- t=(x-xx(n))/(xx(n+1)-xx(n))
- u=(y-yy(m))/(yy(m+1)-yy(m))
-
- z=0.
- do 1 l=4,1,-1
- z=t*z+((cc(n,m,l,4)*u+cc(n,m,l,3))*u
- . +cc(n,m,l,2))*u+cc(n,m,l,1)
- 1 continue
- return
- end
-
- integer function locx(xx,nx,x)
- implicit real*8(a-h,o-z)
- dimension xx(nx)
- if(x.le.xx(1)) then
- locx=1
- return
- endif
- if(x.ge.xx(nx)) then
- locx=nx-1
- return
- endif
- ju=nx+1
- jl=0
- 1 if((ju-jl).le.1) go to 2
- jm=(ju+jl)/2
- if(x.ge.xx(jm)) then
- jl=jm
- else
- ju=jm
- endif
- go to 1
- 2 locx=jl
- return
- end
-
-
- real*8 function polderiv(x1,x2,x3,y1,y2,y3)
- implicit real*8(a-h,o-z)
- polderiv=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
- .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
- return
- end
diff --git a/gonsalves/qt.f b/gonsalves/qt.f
deleted file mode 100644
index 306dcd1..0000000
--- a/gonsalves/qt.f
+++ /dev/null
@@ -1,1739 +0,0 @@
-************************************************************************
-* Numerical Program for Inclusive gamma*, W & Z production *
-* at large transverse momentum in hadron-hadron collisions *
-* Includes next-to-leading order QCD corrections given in *
-* R.J. Gonsalves, J. Pawlowski, C.-F. Wai, Phys. Rev. D40, *
-* 2245 (1989). *
-*----------------------------------------------------------------------*
-* UBVM FILE : QT FORTRAN *
-* Created : 1989 *
-* Revised : June 30, 1990 *
-* Updated : July-August 2004: *
-* Renamed qt.f *
-* Updated for latest MRST, CTEQ parton distributions *
-* Changed alpha 1/137 -> alpha(M_Z) *
-* Updated top quark mass *
-* September 2011: *
-* Corrected xlu,xld in setcon *
-* Corrected triangle loop code in plumin *
-*----------------------------------------------------------------------*
-* Example input file "qt.inp" *
-* computes dsigma/dQ_T**2 in p-pbar collisions at 1.8 TeV *
-* for W+ + W- production using MRST2002 NLO parton densities *
-* factorization/renormalization scales are set to Q_T *
-* Q_T values: 10,30,50,70,90,110,130,150,170,190 GeV *
-*---- Begin: Cut here and remove leading and trailing * ---------------*
-* Input parameters for qt.f: Inclusive gamma*, W, Z production *
-* Hadrons: lppb lns e pdset *
-* t f 1800 mrs02nlo *
-* Bosons: npwz lwpm q qt y *
-* 2 t 80 50 0 *
-* Partons: lallpr 1 2 3 4 5 6 7 8 9 10 11 12 *
-* t f t t f f f f f f f f f *
-* QCD: nscale nalpha lord1 lord2 *
-* 1 1 t t *
-* Scales: lscaleq lscalm lscalmu scalm scalmu *
-* t f f 0.0 0.0 *
-* Thresholds: linctri *
-* f *
-* Integrate: ndim ntot begin step npoints *
-* 3 0 10.0 20.0 10 *
-* Vegas: ncall itmx0 itmx ndev nprn alph0 alph *
-* 40000 3 7 2 -1 1.5 1.5 *
-* Ran Seed: iseed *
-* 13579 *
-* Output: lterm lauto ofname comment *
-* t f qt.out # *
-*---- End: Example input file -----------------------------------------*
-* Input file parameters: *
-* ------ line 3 ------ *
-* lppb: t = proton proton, f = proton antiproton *
-* lns: t = non-singlet, f = non-singlet and singlet *
-* e = hadron-hadron center of mass energy in GeV *
-* pdset: parton density set (8 character string) *
-* See BLOCK DATA pddata for sets implemented *
-* ----- line 5 ------ *
-* npwz: vector boson V produced *
-* 1 = gamma*, 2 = W-, 3 = W+, 4 = Z *
-* lwpm: t = sum W+ and W- f = single W- (npwz=2) or W+ (npwz=3) *
-* q = virtual photon mass in GeV [program sets q=M_(W/Z) for W/Z] *
-* qt = vector boson transverse momentum in GeV *
-* y = vector boson rapidity *
-* ------ line 7 ------ *
-* lallpr: t = include all subprocesses *
-* f = include subprocesses marked t *
-* nproc: 1 2 3 4 5 6 7 8 9 10 11 12 in function xcfin *
-* ------ line 9 ------ *
-* linctri: t = include virtual triangle diagram, f = exclude *
-* ------ line 11 ------ *
-* nscale: energy scale used for renormalization/factorization *
-* 1 = q_transverse, 2 = Q=M_V, 3 = sqrt(q_t**2 + Q**2) *
-* nalpha: use of NLO and LO QCD coupling alpha_s *
-* 1 = NLO + NLO**2, 2 = NLO + LO*NLO, 3 = NLO + LO**2 *
-* lord1: t = include, f = exclude leading order *
-* lord2: t = include, f = exclude next-leading order *
-* ------ line 13 ------ *
-* lscaleq: t = set renormalization factorization scales equal *
-* lscalm: t = cross section as function of factorization scale *
-* lscalmu: t = cross section as function of renormalization scale *
-* scalm: factorization scale is multiplied by sqrt(10**scalm) *
-* scalmu: renormalization scale is multiplied by sqrt(10**scalmu) *
-* ------ line 15 ------ *
-* ndim: dimension of phase space integration *
-* 2 = dsigma/dq_t**2/dy *
-* 3 = dsigma/dq_t**2 *
-* 4 = sigma_total with (q_t > q_t_min) *
-* ntot: normalize by total cross section *
-* 1 = compute and divide by total cross section *
-* begin: starting value for chosen variable *
-* step: step in chosen variable *
-* npoints: number of points in chosen variable *
-* ------ line 17 ------ *
-* ncall = number of integrand evaluations *
-* itmx0 = maximum number of warmup (discarded) iterations *
-* itmx = maximum number of iterations *
-* ndev = Fortran device for output from Vegas *
-* nprn : Vegas prints following on ndev for each iteration *
-* > 0 : integral, std dev, chi^2, grid information *
-* = 0 : integral, std dev, chi^2 *
-* < 0 : nothing *
-* alph0 = rate at which grid is modified during warmup *
-* alph = rate at which grid is modified during data taking *
-* ------ line 19 ------ *
-* iseed = random number seed for ran2 *
-* ------ line 21 ------ *
-* lterm: t = write numerical values on term, f = don't write *
-* lauto: t = derive output file name from input file name *
-* by replacing .inp with .out *
-* f = use ofname for output file *
-* ofname: name of output file *
-* comment: comment character for output file *
-************************************************************************
-C-----------------------------------------------------------------------
- BLOCK DATA pddata
- IMPLICIT NONE
- INTEGER maxset,i
- PARAMETER (maxset=13)
- CHARACTER pdset(maxset)*8,pdname(maxset)*45
- COMMON /pdsets/ pdset,pdname
- DOUBLE PRECISION lambda(maxset)
- COMMON /lambdas/ lambda
- DATA (pdset(i),i=1,maxset)/
- & 'do84set1',
- & 'do84set2',
- & 'ehlqset1',
- & 'ehlqset2',
- & 'mrs123m1',
- & 'mrs123m2',
- & 'mrs123m3',
- & 'mrs89em1',
- & 'mrs89bm2',
- & 'mrs01lo1',
- & 'mrs02nlo',
- & 'mrs02nnl',
- & 'ctq61nlo'
- & /
- DATA (pdname(i),i=1,maxset)/
- & 'Duke-Owens 1984 Set 1 Lambda=0.2 GeV ',
- & 'Duke-Owens 1984 Set 2 Lambda=0.4 GeV ',
- & 'EHLQ Set 1 Lambda=0.2 GeV ',
- & 'EHLQ Set 2 Lambda=0.29 GeV ',
- & 'MRS123 Mode 1 (soft glue) Lambda=0.107 GeV ',
- & 'MRS123 Mode 2 (hard glue) Lambda=0.250 GeV ',
- & 'MRS123 Mode 3 (1/RTX glue) Lambda=0.178 GeV ',
- & 'MRSEB Mode 1 (EMC FIT) Lambda=0.091 GeV ',
- & 'MRSEB Mode 2 (BCDMS FIT) Lambda=0.228 GeV ',
- & 'MRST2001 Mode 1 LO Lambda_4=0.220 GeV ',
- & 'MRST2002 Mode 1 NLO Lambda_4=0.323 GeV ',
- & 'MRST2002 Mode 2 NNLO Lambda_4=0.234 GeV ',
- & 'Cteq61 Mode 1 NLO Lambda_4=0.326 GeV '
- & /
- DATA (lambda(i),i=1,maxset)/
- & 0.200d0, 0.400d0, 0.200d0, 0.290d0,
- & 0.107d0, 0.250d0, 0.178d0, 0.091d0, 0.228d0,
- & 0.220d0, 0.323d0, 0.234d0, 0.326d0
- & /
- END
-C-----------------------------------------------------------------------
- PROGRAM qtmain
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- CHARACTER today*30,proces*50,numb(12)*3,string*50,ifname*50
- & ,ofname*50
- LOGICAL yes,no,lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- & ,lscaleq,lscalm,lscalmu
- PARAMETER (yes=.true.,no=.false.)
- DIMENSION result(5,1000)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- LOGICAL lasnnlo
- COMMON /nnlo/ asnnlo,lasnnlo
- COMMON /bveg1/ ncall,itmx,nprn,ndev,xl(10),xu(10),acc
- COMMON /bveg3/ alph,ndmx,mds
- COMMON /counter/ iread
- COMMON /iquad/ alph0,itmx0
- DATA (numb(i),i=1,12) /' 1 ',' 2 ',' 3 ',' 4 ',' 5 ',' 6 ',
- & ' 7 ',' 8 ',' 9 ','10 ','11 ','12 '/
- INTEGER maxset,i
- PARAMETER (maxset=13)
- CHARACTER pdset(maxset)*8,pdname(maxset)*45
- COMMON /pdsets/ pdset,pdname
- DOUBLE PRECISION lambda(maxset)
- COMMON /lambdas/ lambda
- CHARACTER pdsetname*8,comment*2
- COMMON /pdsetn/ pdsetname
- LOGICAL lterm,lauto
- LOGICAL linctri
- COMMON /triangle/ linctri
-
- CALL ctime(time(),today)
-* use iargc and getarg to get input file name from command line
- IF (iargc().GT.0) THEN
- CALL getarg(1,ifname)
- ELSE
- ifname='qt.inp'
- ENDIF
-
- ninput=1
- OPEN (unit=ninput,file=ifname,status='old')
- PRINT *,'Reading input from file ',ifname
- READ (ninput,*)
- READ (ninput,*)
- READ (ninput,*) lppb,lns,e,string
- READ (ninput,*)
- READ (ninput,*) npwz,lwpm,q,qt,y
- READ (ninput,*)
- READ (ninput,*) lallpr,(lproc(i),i=1,12)
- READ (ninput,*)
- READ (ninput,*) nscale,nalpha,lord1,lord2
- READ (ninput,*)
- READ (ninput,*) lscaleq,lscalm,lscalmu,scalm,scalmu
- READ (ninput,*)
- READ (ninput,*) linctri
- READ (ninput,*)
- READ (ninput,*) ndim,ntot,begin,step,npoints
- READ (ninput,*)
- READ (ninput,*) ncall,itmx0,itmx,ndev,nprn,alph0,alph
- READ (ninput,*)
- READ (ninput,*) iseed
- READ (ninput,*)
- READ (ninput,*) lterm,lauto,ofname,comment
-
- nset=0
- DO i=1,maxset
- IF (string(1:8).EQ.pdset(i)) THEN
- nset=i
- xlambd=lambda(i)
- pdsetname=pdset(i)
- ENDIF
- ENDDO
- IF (pdsetname.EQ.'mrs02nnl') THEN
- lasnnlo=.TRUE.
- ELSE
- lasnnlo=.FALSE.
- ENDIF
- IF (nset.EQ.0) THEN
- PRINT *,'Bad parton density set ',string
- STOP
- ENDIF
- IF (lallpr) THEN
- DO i=1,12
- lproc(i)=yes
- ENDDO
- ENDIF
- scalm=10d0**scalm
- scalmu=10d0**scalmu
- idum=iseed
- CALL setseed(idum)
- IF (lauto) THEN
- i=index(ifname,'.inp')
- IF (i.EQ.0) i=index(ifname,' ')
- ofname=ifname(1:i-1)//'.out'
- ENDIF
- OPEN (unit=2,file=ofname,status='unknown')
- IF (lterm) PRINT *,'Output will go to file ',ofname
-
- IF (npwz.EQ.1) lproc(12)=no
- IF (npwz.EQ.2.OR.npwz.EQ.3) THEN
- lproc(10)=no
- lproc(11)=no
- lproc(12)=no
- ENDIF
- IF (lwpm.AND.npwz.EQ.3) npwz=2
- IF (lns) THEN
- lproc(2)=no
- lproc(3)=no
- lproc(5)=no
- ENDIF
- IF (lord1.AND..NOT.lord2) THEN
- DO i=3,12
- lproc(i)=no
- ENDDO
- ENDIF
- proces=' '
- DO i=1,12
- IF (lproc(13-i)) proces=(numb(13-i)//proces)
- ENDDO
- IF (lord2) proces=('nlo '//proces)
- IF (lord1) proces=('lo '//proces)
-
- CALL setcon
- IF (ntot.EQ.1) THEN
- CALL settot (totxc,totsd,totchi)
- ENDIF
- nres=0
- IF (lscalm.OR.lscalmu) THEN
- IF (lscalm.AND.lscalmu.AND..NOT.lscaleq) THEN
- DO i=1,npoints
- scalm=10d0**(begin+step*(i-1))
- DO j=1,npoints
- scalmu=10d0**(begin+step*(j-1))
- nres=nres+1
- result(1,nres)=scalm
- result(2,nres)=scalmu
- CALL setup(result(3,nres),result(4,nres)
- & ,result(5,nres))
- IF (ntot.EQ.1) THEN
- result(3,nres)=result(3,nres)/totxc*2.
- result(4,nres)=result(4,nres)/totxc*2.
- ENDIF
- IF (lterm) WRITE (6,105) (result(k,nres),k=1,5)
- ENDDO
- ENDDO
- ELSE
- DO i=1,npoints
- scale=10d0**(begin+step*(i-1))
- nres=nres+1
- result(1,nres)=scale
- IF (lscalm.OR.lscaleq) scalm=scale
- IF (lscalmu.OR.lscaleq) scalmu=scale
- CALL setup (result(2,nres),result(3,nres)
- & ,result(4,nres))
- IF (ntot.EQ.1) THEN
- result(2,nres)=result(2,nres)/totxc*2.
- result(3,nres)=result(3,nres)/totxc*2.
- ENDIF
- IF (lterm) WRITE (6,100) (result(j,nres),j=1,4)
- ENDDO
- ENDIF
- ELSE
- DO qtp=begin,begin+step*(npoints-1),step
- qt=qtp
- nres=nres+1
- result(1,nres)=qt
- CALL setup (result(2,nres),result(3,nres),result(4,nres))
- IF (ntot.EQ.1) THEN
- result(2,nres)=result(2,nres)/totxc*2.
- result(3,nres)=result(3,nres)/totxc*2.
- ENDIF
- IF (lterm) WRITE (6,100) (result(j,nres),j=1,4)
- ENDDO
- ENDIF
-
- WRITE (ndev,110) comment,today,ifname
- WRITE (ndev,120) comment
- IF (lppb) string='p pbar collision'
- IF (.NOT. lppb) string='p p collision'
- IF (lns) string='Non-singlet: ppb - pp'
- WRITE (ndev,130) comment,string,e
- WRITE (ndev,140) comment,pdname(nset)
- IF (npwz.EQ.1) string='Virtual photon'
- IF (npwz.EQ.2) THEN
- string='W^- Production'
- IF (lwpm) string=' (W^-) + (W^+)'
- ENDIF
- IF (npwz.EQ.3) string='W^+ Production'
- IF (npwz.EQ.4) string='Z_0 Production'
- WRITE (ndev,150) comment,string,q
- WRITE (ndev,160) comment,proces
- WRITE (ndev,170) comment,nscale,nalpha
- IF (lscalm.OR.lscalmu) THEN
- IF (lscalm) WRITE (ndev,180) comment,qt,scalmu
- IF (lscalmu) WRITE (ndev,190) comment,qt,scalm
- ELSE
- WRITE (ndev,200) comment,scalm,scalmu
- ENDIF
- IF (npwz.EQ.4) THEN
- IF (linctri) THEN
- WRITE (ndev,205) comment,' included'
- ELSE
- WRITE (ndev,205) comment,' not included'
- ENDIF
- ENDIF
- IF (ndim.EQ.2) string='d sigma/(dqt**2 dy)'
- IF (ndim.EQ.3) string='d sigma/dqt**2'
- IF (ndim.EQ.4) string='sigma(qt > qtmin)'
- WRITE (ndev,210) comment,string
- IF (ntot.EQ.1) WRITE (ndev,220) comment,totxc,totsd,totchi
- WRITE (ndev,230) comment,ncall,itmx0,itmx,alph0,alph
- WRITE (ndev,240) comment,iseed
- WRITE (ndev,250) comment
- IF (lscalm.AND.lscalmu.AND..NOT.lscaleq) THEN
- DO i=1,nres
- IF (i.GT.1.AND.MOD(i-1,npoints).EQ.0) WRITE(ndev,106)
- WRITE (ndev,105) (result(j,i),j=1,5)
- ENDDO
- ELSE
- DO i=1,nres
- WRITE (ndev,100) (result(j,i),j=1,4)
- ENDDO
- ENDIF
-
- CLOSE (UNIT=ninput)
- CLOSE (UNIT=2)
-
- 100 FORMAT (1x,4(g9.3,3x))
- 105 FORMAT (1x,5(g9.3,3x))
- 106 FORMAT (1x)
- 110 FORMAT (a2,a30,'Input file: ',a30)
- 120 FORMAT (a2,'Electroweak Boson Production at Large Q_T')
- 130 FORMAT (a2,a30,'sqrt(S) = ',f7.1,' GeV')
- 140 FORMAT (a2,'Parton Densities: ',a45)
- 150 FORMAT (a2,'Vector Boson: ',a20,'Q = ',f7.2,' GeV')
- 160 FORMAT (a2,'Parton Processes: ',a50)
- 170 FORMAT (a2,'QCD: nscale = ',i2,2x,'nalpha = ',i2)
- 180 FORMAT (a2,'Scalm dependence at Q_T = ',f7.2,' scalmu = ',f7.2)
- 190 FORMAT (a2,'Scalmu dependence at Q_T = ',f7.2,' scalm = ',f7.2)
- 200 FORMAT (a2,'scalm = ',f7.2,' scalmu = ',f7.2)
- 205 FORMAT (a2,'Z_0 virtual triangle diagrams',a14)
- 210 FORMAT (a2,'Cross section computed: ',a30)
- 220 FORMAT (a2,'Divided by totxc/2 = ',3(g9.3,3x))
- 230 FORMAT (a2,'Vegas: ncall = ',i8,' itmx0 = ',i4,' itmx = ',i4,
- & ' alph0 = ',f4.2,' alph = ',f4.2)
- 240 FORMAT (a2,'Seed for random number generator = ',i12)
- 250 FORMAT (a2,'Variable Integral Std. Dev. ChiSqd/dof')
- END
-C-----------------------------------------------------------------------
- SUBROUTINE setcon
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
-
- pi=3.141592654
- tonbs=.389386d6
-* alpha=1./137.0360
- alpha=1./127.918
- sin2tw=0.2312
- sw=dsqrt(sin2tw)
- cw=dsqrt(1.d0-sin2tw)
- cf=4./3.
- ca=3.
- xnc=3.
- xmz=91.1876
-* xmw=xmz*cw
- xmw=80.425
- xmc=1.25
- xmb=4.25
-* xmt=178.1
- xmt=1e6
- qu=2./3.
- qd=-1./3.
- qw=1./dsqrt(2.d0)/sw
- xlu=1./2./sw/cw-qu*sw/cw
- xld=-1./2./sw/cw-qd*sw/cw
- xru=-qu*sw/cw
- xrd=-qd*sw/cw
- v2(1,1)=(0.9738)**2
- v2(1,2)=(0.2200)**2
- v2(1,3)=1.-v2(1,1)-v2(1,2)
- v2(2,1)=(0.224)**2
- v2(2,2)=(0.996)**2
- v2(2,3)=1.-v2(2,1)-v2(2,2)
- v2(3,1)=1.-v2(1,1)-v2(2,1)
- v2(3,2)=1.-v2(1,2)-v2(2,2)
- v2(3,3)=1.-v2(1,3)-v2(2,3)
- RETURN
- END
-C----------------------------------------------------------------------
- SUBROUTINE setup (ans,sd,chi2a)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- LOGICAL lasnnlo
- COMMON /nnlo/ asnnlo,lasnnlo
- COMMON /bveg1/ ncall,itmx,nprn,ndev,xl(10),xu(10),acc
- COMMON /bveg3/ alph,ndmx,mds
- COMMON /iquad/ alph0,itmx0
- EXTERNAL f2veg,f2dim,f3veg,f3dim,s2lims,s2llim,s2ulim
- EXTERNAL f4dim,f4veg
-
- sh=e**2
- IF (npwz.EQ.2.OR.npwz.EQ.3) q=xmw
- IF (npwz.EQ.4) q=xmz
- qq=q**2
- qt2=qt**2
- eqt=dsqrt(qt2+qq)
- IF (nscale.EQ.1) escal2=qt2
- IF (nscale.EQ.2) escal2=qq
- IF (nscale.EQ.3) escal2=qt2+qq
- fm=dlog(scalm*escal2/qq)
- fmu=dlog(scalmu*escal2/qq)
- CALL flavor
- IF (nset.LE.9) THEN
- aslo=alphas_old(1)/2./pi
- ELSE
- scale=dsqrt(escal2*scalmu)
- aslo=alphas(scale,xlambd,0)/2d0/pi
- ENDIF
- asnlo=aslo
- asnnlo=aslo
- IF (lord2) THEN
- IF (nset.LE.9) THEN
- asnlo=alphas_old(2)/2./pi
- ELSE
- scale=dsqrt(escal2*scalmu)
- asnlo=alphas(scale,xlambd,1)/2d0/pi
- asnnlo=alphas(scale,xlambd,2)/2d0/pi
- ENDIF
- ENDIF
- surd=dsqrt((sh-qq)**2-4.*sh*qt2)
- ymax=1./2.*dlog((sh+qq+surd)/(sh+qq-surd))
- qtmax=(sh-qq)/2d0/e
- qtmin=qt
- xl(1)=0d0
- xu(1)=sh-qq
- xl(2)=dlog(qq/sh)
- xu(2)=0d0
- xl(3)=-ymax
- xu(3)=ymax
- xl(4)=dlog(qtmin**2/sh)
- xu(4)=dlog(qtmax**2/sh)
- npr=nprn
- itm=itmx
- alp=alph
- nprn=-1
- itmx=itmx0
- alph=alph0
- IF (ndim.EQ.2) CALL vegas (2,f2veg,ans,sd,chi2a)
- IF (ndim.EQ.3) CALL vegas (3,f3veg,ans,sd,chi2a)
- IF (ndim.EQ.4) CALL vegas (4,f4veg,ans,sd,chi2a)
- nprn=npr
- itmx=itm
- alph=alp
- IF (ndim.EQ.2) CALL vegas1 (2,f2veg,ans,sd,chi2a)
- IF (ndim.EQ.3) CALL vegas (3,f3veg,ans,sd,chi2a)
- IF (ndim.EQ.4) CALL vegas1 (4,f4veg,ans,sd,chi2a)
- omult=4.*pi*alpha*cf*xnc*tonbs
- ans=ans*omult
- sd=sd*omult
- IF (npwz.EQ.1) THEN
- emult=alpha/3./pi/qq
- ans=ans*emult
- sd=sd*emult
- ENDIF
- ans=ans*pi
- sd=sd*pi
- RETURN
- END
-C----------------------------------------------------------------------
- SUBROUTINE settot (ans,sd,chi2a)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- COMMON /bveg1/ ncall,itmx,nprn,ndev,xl(10),xu(10),acc
- COMMON /bveg3/ alph,ndmx,mds
- COMMON /iquad/ alph0,itmx0
- EXTERNAL totveg
-
- sh=e**2
- IF (npwz.EQ.2.OR.npwz.EQ.3) q=xmw
- IF (npwz.EQ.4) q=xmz
- qq=q**2
- escal2=qq
- CALL flavor
- saveit=scalmu
- scalmu=1.
- IF (nset.LE.9) THEN
- aslo=alphas_old(1)/2./pi
- ELSE
- aslo=alphas(q,xlambd,1)/2d0/pi
- ENDIF
- scalmu=saveit
- ymax=1./2.*dlog(sh/qq)
- xl(1)=0d0
- xu(1)=1.
- xl(2)=-ymax
- xu(2)=ymax
- npr=nprn
- itm=itmx
- alp=alph
- nprn=-1
- itmx=itmx0
- alph=alph0
- CALL vegas (2,totveg,ans,sd,chi2a)
- nprn=npr
- itmx=itm
- alph=alp
- CALL vegas1 (2,totveg,ans,sd,chi2a)
- omult=4.*pi**2*alpha*xnc*tonbs
- ans=ans*omult
- sd=sd*omult
- IF (npwz.EQ.1) THEN
- emult=alpha/3./pi/qq
- ans=ans*emult
- sd=sd*emult
- ENDIF
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE flavor
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- CHARACTER pdsetname*8
- COMMON /pdsetn/ pdsetname
-
- qsq=escal2
- xnf=3.
- sqf2=qu**2+2.*qd**2
- sxrml=xrd-xld
- DO 1 i=1,3
- sv2u(i)=v2(i,1)+v2(i,2)
- sv2d(i)=v2(1,i)
-1 CONTINUE
- sv2=sv2u(1)
- szf2=(xlu**2+xru**2)/2.+xld**2+xrd**2
- IF (qsq.GT.4*xmc**2) THEN
- xnf=4.
- sqf2=sqf2+qu**2
- DO 2 i=1,3
- sv2d(i)=sv2d(i)+v2(2,i)
-2 CONTINUE
- sv2=sv2u(1)+sv2u(2)
- szf2=szf2+(xlu**2+xru**2)/2.
- sxrml=0d0
- ENDIF
- IF (pdsetname(1:4).EQ.'do84') RETURN
- IF (qsq.GT.4*xmb**2) THEN
- xnf=5.
- sqf2=sqf2+qd**2
- DO 3 i=1,3
- sv2u(i)=sv2u(i)+v2(i,3)
-3 CONTINUE
- sv2=sv2u(1)+sv2u(2)
- szf2=szf2+(xld**2+xrd**2)/2.
- sxrml=xrd-xld
- ENDIF
- IF (pdsetname(1:4).NE.'ehlq') RETURN
- IF (qsq.GT.4*xmt**2) THEN
- xnf=6.
- sqf2=sqf2+qu**2
- DO 4 i=1,3
- sv2d(i)=sv2d(i)+v2(3,i)
-4 CONTINUE
- sv2=sv2u(1)+sv2u(2)+sv2u(3)
- szf2=szf2+(xlu**2+xru**2)/2.
- sxrml=0d0
- ENDIF
- RETURN
- END
-C-----------------------------------------------------------------------
-C QCD Effective coupling - Bardeen et al. definition
-C
- DOUBLE PRECISION FUNCTION alphas_old (norder)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lproc(12)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- qsq=escal2*scalmu
- t=dlog(qsq/xlambd**2)
- b1=12.*pi/(33.-2.*xnf)
- alphas=b1/t
- IF (norder.EQ.2) THEN
- b2=24.*pi**2/(153.-19.*xnf)
- alphas=alphas*(1.-b1**2*dlog(t)/b2/t)
- ENDIF
- alphas_old=alphas
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f2dim (xlogx1,s2p)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- x1=dexp(xlogx1)
- s2=s2p
- f2dim=preigd ()
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f3dim (yp,xlogx1,s2p)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- y=yp
- x1=dexp(xlogx1)
- s2=s2p
- f3dim=preigd ()
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f4dim (xlqt2,yp,xlogx1,s2p)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- qt=e*dexp(xlqt2/2.)
- qt2=qt*qt
- eqt=dsqrt(qt2+qq)
- y=yp
- x1=dexp(xlogx1)
- s2=s2p
- f4dim=qt2*preigd ()
- RETURN
- END
-C----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f2tot (y,z)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- f2tot=pretot(y,z)
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f2veg (x,wgt)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION x(2)
- f2veg=f2dim(x(2),x(1))
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f3veg (x,wgt)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION x(3)
- f3veg=f3dim(x(3),x(2),x(1))
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION f4veg (x,wgt)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION x(4)
- f4veg=f4dim(x(4),x(3),x(2),x(1))
- RETURN
- END
-C----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION totveg (x,wgt)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION x(2)
- totveg=f2tot(x(2),x(1))
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION s2lims (xlogx1)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- ENTRY s2ulim (xlogx1)
- s2ulim = uh+dexp(xlogx1)*(sh+th-qq)
- RETURN
- ENTRY s2llim (xlogx1)
- s2llim = 0d0
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION preigd ()
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- DIMENSION pdd(13,2),pdd0(13,2)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- LOGICAL lasnnlo
- COMMON /nnlo/ asnnlo,lasnnlo
- LOGICAL linctri
- COMMON /triangle/ linctri
-
- th=dexp(y)
- uh=qq-e*eqt*th
- th=qq-e*eqt/th
-
-C Phase space boundaries for Vegas
- IF (ndim.EQ.4) THEN
- surd=dsqrt((sh-qq)**2-4.*sh*qt2)
- ymax=1./2.*dlog((sh+qq+surd)/(sh+qq-surd))
- IF(dabs(y).GT.ymax) GOTO 100
- ENDIF
- b=-uh/(sh+th-qq)
- IF (x1.LT.b) GOTO 100
- a=uh+x1*(sh+th-qq)
- IF (s2.GT.a) GOTO 100
-
- xjac=x1*sh+uh-qq
- x20=(-x1*th-qq*(1.-x1))/xjac
- x2=s2/xjac+x20
- fs2=dlog(s2/qq)
- fa=dlog(a/qq)
- fscale=escal2*scalm
- CALL compdd (x1,x2,x20,fscale,pdd,pdd0)
- alow=0d0
- ads2=0d0
- as2a=0d0
- afin=0d0
- CALL comvar (0,x20)
- check=(x1+x20)/2.*e-eqt*dcosh(y)
- IF (check.LT.0.D0) STOP 'PREIGD error: Insufficient Energy'
- DO 2 i=1,2
- IF (lord1) THEN
- IF (lproc(1).AND.i.EQ.1) alow=alow+xclow(1)*pdd0(1,i)
- IF (lproc(2)) alow=alow+xclow(2)*pdd0(2,i)
- ENDIF
- IF (lord2) THEN
- IF (lproc(1).AND.i.EQ.1) THEN
- ads2=ads2+xcds2(1)*pdd0(1,i)
- IF (linctri) THEN
- ads2=ads2+triang(1)*pdd0(12,i)
- ENDIF
- as2a=as2a-xcs2a(1)*pdd0(1,i)
- ENDIF
- IF (lproc(2)) THEN
- ads2=ads2+xcds2(2)*pdd0(2,i)
- IF (linctri) THEN
- ads2=ads2+triang(2)*pdd0(13,i)
- ENDIF
- as2a=as2a-xcs2a(2)*pdd0(2,i)
- ENDIF
- ENDIF
- CALL exchtu
-2 CONTINUE
- alow=alow/s
- ads2=ads2/s
- as2a=as2a/s
- IF (.NOT.lord2) GOTO 5
- CALL comvar (1,x2)
- DO 3 i=1,2
- IF (lproc(1).AND.i.EQ.1) as2a=as2a+xcs2a(1)*pdd(1,i)/s
- IF (lproc(2)) as2a=as2a+xcs2a(2)*pdd(2,i)/s
- DO 4 j=1,12
- IF (lproc(j)) afin=afin+xcfin(j)*pdd(j,i)
-4 CONTINUE
- CALL exchtu
-3 CONTINUE
- afin=afin/s+as2a
-5 CONTINUE
- IF (lasnnlo) THEN
- as23=asnnlo
- ELSE
- as23=asnlo
- ENDIF
- IF (nalpha.EQ.1) preigd=as23*alow/a+as23**2*(ads2/a+afin)
- IF (nalpha.EQ.2) preigd=as23*((alow+aslo*ads2)/a+aslo*afin)
- IF (nalpha.EQ.3) preigd=as23*alow/a+aslo**2*(ads2/a+afin)
- preigd=preigd*x1/xjac
- RETURN
-100 preigd=0d0
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION pretot (yq,z)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- DIMENSION pdd(13,2),pdd0(13,2),pdd1(13,2)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- tau=qq/sh
-
-C Phase space boundary for Vegas
- zmin=tau*dexp(2.*dabs(yq))
- IF (z.LT.zmin) GOTO 100
-
- x11=dsqrt(tau)
- x21=x11
- dexpyq=dexp(yq)
- x11=x11*dexpyq
- x21=x21/dexpyq
- x2=dsqrt(z)
- x1=x11/x2
- x2=x21/x2
- CALL compdd (x11,x21,x21,qq,pdd1,pdd0)
- CALL compdd (x1,x2,x2,qq,pdd,pdd0)
-
- pretot=0d0
- IF (lproc(1)) THEN
- qqb=1./(1.-zmin)*(1.+aslo*cf*(2.*pi**2/3.-8.))*pdd1(1,1)
- & +aslo*cf*(4.*dlog(1.-z)/(1.-z)*((1.+z**2)/z
- & *pdd(1,1)-2.*pdd1(1,1))
- & -2.*(1.+z**2)/z*dlog(z)*pdd(1,1))
- pretot=pretot+qqb
- ENDIF
- IF (lproc(2)) THEN
- qg=aslo*cf*((z**2+(1.-z)**2)*(dlog((1.-z)**2/z)-3./2.)
- & +2.-z**2/2.)*(pdd(2,1)+pdd(2,2))
- pretot=pretot+qg
- ENDIF
- pretot=pretot/sh
- RETURN
-100 pretot=0d0
- END
-C-----------------------------------------------------------------------
- SUBROUTINE comvar (ns2,x2)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- s=x1*x2*sh
- t=x1*th+(1.-x1)*qq
- u=x2*uh+(1.-x2)*qq
- s2p=s2*ns2
- d1=1./(s2p-t)
- d2=1./(s2p-u)
- d3=1./(s+t-s2p)
- d4=1./(s+u-s2p)
- d5=1./(s+qq-s2p)
- d6=1./(u*t-s2p*qq)
- d10=1./((u+t)**2-4.*s2p*qq)
- d9=dsqrt(d10)
- fs=dlog(s/qq)
- ft=dlog(-t/qq)
- fu=dlog(-u/qq)
- fst=fs-2.*dlog(-1./d1/t)
- fsu=fs-2.*dlog(-1./d2/u)
- ftu=dlog(d1*d2/d6)
- fstu=dlog(s*qq*d1*d2)
- fla=dlog((d9+d5)/(d9-d5))
- flt=dlog(s*qq*(1./(s2p*(2.*qq-u)-qq*t)/d1)**2)
- flu=dlog(s*qq*(1./(s2p*(2.*qq-t)-qq*u)/d2)**2)
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE exchtu
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- CALL exch (t,u)
- CALL exch (d1,d2)
- CALL exch (d3,d4)
- CALL exch (ft,fu)
- CALL exch (fst,fsu)
- CALL exch (flt,flu)
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE exch (a,b)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- temp=a
- a=b
- b=temp
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE compdd (x1,x2,x20,qsq,pdd,pdd0)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns,lwpm,lord1,lord2,lallpr,lproc(12)
- DIMENSION pdd(13,2), pdd0(13,2), a(-6:6), b(-6:6), b0(-6:6)
- DIMENSION qdd(13,2), qdd0(13,2)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /switch/ npwz,nscale,nalpha,nset,ndim,ntot
- & ,lppb,lns,lwpm,lord1,lord2,lproc
- COMMON /counter/ iread
- CHARACTER pdsetname*8
- COMMON /pdsetn/ pdsetname
-
- IF (pdsetname(1:4).EQ.'do84') THEN
- IF (pdsetname(8:8).EQ.'1') mode=1
- IF (pdsetname(8:8).EQ.'2') mode=2
- CALL dopden (mode,x1,qsq,a)
- CALL dopden (mode,x2,qsq,b)
- CALL dopden (mode,x20,qsq,b0)
- ELSEIF (pdsetname(1:4).EQ.'ehlq') THEN
- IF (pdsetname(8:8).EQ.'1') mode=1
- IF (pdsetname(8:8).EQ.'2') mode=2
- CALL ehlq (mode,x1,qsq,a)
- CALL ehlq (mode,x2,qsq,b)
- CALL ehlq (mode,x20,qsq,b0)
- ELSEIF (pdsetname(1:6).EQ.'mrs123') THEN
- IF (pdsetname(8:8).EQ.'1') mode=1
- IF (pdsetname(8:8).EQ.'2') mode=2
- IF (pdsetname(8:8).EQ.'3') mode=3
- CALL mrs (mode,x1,qsq,a)
- CALL mrs (mode,x2,qsq,b)
- CALL mrs (mode,x20,qsq,b0)
- ELSEIF (pdsetname(1:5).EQ.'mrs89') THEN
- IF (pdsetname(6:6).EQ.'e') mode=4
- IF (pdsetname(6:6).EQ.'b') mode=5
- CALL mrs (mode,x1,qsq,a)
- CALL mrs (mode,x2,qsq,b)
- CALL mrs (mode,x20,qsq,b0)
- ELSEIF (pdsetname(1:5).EQ.'mrs01') THEN
- IF (pdsetname(6:8).EQ.'lo1') mode=0
- CALL mrst (mode,x1,qsq,a)
- CALL mrst (mode,x2,qsq,b)
- CALL mrst (mode,x20,qsq,b0)
- ELSEIF (pdsetname(1:5).EQ.'mrs02') THEN
- IF (pdsetname(6:8).EQ.'nlo') mode=1
- IF (pdsetname(6:8).EQ.'nnl') mode=2
- CALL mrst (mode,x1,qsq,a)
- CALL mrst (mode,x2,qsq,b)
- CALL mrst (mode,x20,qsq,b0)
- ELSEIF (pdsetname(1:5).EQ.'ctq61') THEN
- IF (pdsetname(6:8).EQ.'nlo') mode=1
- CALL cteq (mode,x1,qsq,a)
- CALL cteq (mode,x2,qsq,b)
- CALL cteq (mode,x20,qsq,b0)
- ENDIF
-C Initial color averages
- DO 1 i=-6,6
- a(i) = a(i)/xnc
- b(i) = b(i)/xnc
- b0(i) = b0(i)/xnc
-1 CONTINUE
- glue = xnc/(xnc**2-1d0)
- a(0) = a(0)*glue
- b(0) = b(0)*glue
- b0(0) = b0(0)*glue
-C
- CALL plumin (lppb,lns,npwz,a,b,pdd)
- CALL plumin (lppb,lns,npwz,a,b0,pdd0)
- IF (lwpm.AND.npwz.EQ.2.AND..NOT.lns) THEN
- CALL plumin (lppb,lns,3,a,b,qdd)
- CALL plumin (lppb,lns,3,a,b0,qdd0)
- DO 2 i=1,13
- DO 2 j=1,2
- pdd(i,j)=pdd(i,j)+qdd(i,j)
-2 pdd0(i,j)=pdd0(i,j)+qdd0(i,j)
- ENDIF
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE plumin (lppb,lns,npwz,a,b,pdd)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL lppb,lns
- DIMENSION a(-6:6),b(-6:6),pdd(13,2),h1u(3),h1ub(3),h2u(3),h2ub(3)
- & ,h1d(3),h1db(3),h2d(3),h2db(3)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
- EXTERNAL sum
-
- qu2=qu**2
- qd2=qd**2
- w2=qw**2/2.
- zu2=(xlu**2+xru**2)/2.
- zd2=(xld**2+xrd**2)/2.
-
- CALL vector (a(1),a(4),a(6),h1u)
- CALL vector (a(-1),a(-4),a(-6),h1ub)
- CALL vector (a(2),a(3),a(5),h1d)
- CALL vector (a(-2),a(-3),a(-5),h1db)
- IF (lppb) THEN
- CALL vector (b(1),b(4),b(6),h2ub)
- CALL vector (b(-1),b(-4),b(-6),h2u)
- CALL vector (b(2),b(3),b(5),h2db)
- CALL vector (b(-2),b(-3),b(-5),h2d)
- ELSE
- CALL vector (b(1),b(4),b(6),h2u)
- CALL vector (b(-1),b(-4),b(-6),h2ub)
- CALL vector (b(2),b(3),b(5),h2d)
- CALL vector (b(-2),b(-3),b(-5),h2db)
- ENDIF
- sh1u=sum(h1u)+sum(h1ub)
- sh2u=sum(h2u)+sum(h2ub)
- sh1d=sum(h1d)+sum(h1db)
- sh2d=sum(h2d)+sum(h2db)
- sh1=sh1u+sh1d
- sh2=sh2u+sh2d
-
- IF (lns) THEN
- pdd(2,1)=0d0
- pdd(3,1)=0d0
- pdd(5,1)=0d0
- pdd(13,1)=0d0
- uv1=a(1)-a(-1)
- dv1=a(2)-a(-2)
- uv2=b(1)-b(-1)
- dv2=b(2)-b(-2)
- IF (npwz.EQ.1) THEN
- pdd(1,1)=qu2*uv1*uv2+qd2*dv1*dv2
- pdd(4,1)=sqf2*(uv1*uv2+dv1*dv2)
- pdd(7,1)=pdd(1,1)
- pdd(10,1)=2.*(qu*uv1+qd*dv1)*(qu*uv2+qd*dv2)
- pdd(12,1)=0d0
- ELSEIF (npwz.EQ.2.OR.npwz.EQ.3) THEN
- pdd(1,1)=w2*v2(1,1)*(uv1*dv2+dv1*uv2)
- pdd(4,1)=w2*sv2*(uv1*uv2+dv1*dv2)
- pdd(7,1)=w2*(sv2u(1)*dv1*dv2+sv2d(1)*uv1*uv2)
- pdd(10,1)=0d0
- pdd(12,1)=0d0
- ELSEIF (npwz.EQ.4) THEN
- pdd(1,1)=zu2*uv1*uv2+zd2*dv1*dv2
- pdd(4,1)=szf2*(uv1*uv2+dv1*dv2)
- pdd(7,1)=pdd(1,1)
- pdd(10,1)=((xlu+xru)*uv1+(xld+xrd)*dv1)/2.
- pdd(10,1)=pdd(10,1)*((xlu+xru)*uv2+(xld+xrd)*dv2)
- pdd(12,1)=(uv1*uv2+dv1*dv2)*(xru-xlu)*sxrml/2.
- ENDIF
- pdd(6,1)=pdd(1,1)
- pdd(8,1)=-pdd(1,1)/2.
- pdd(9,1)=-pdd(7,1)/2.
- pdd(11,1)=pdd(10,1)
- DO 55 i=1,13
-55 pdd(i,2)=pdd(i,1)
- RETURN
- ENDIF
-
-C 1. ann + qqb12
- IF (npwz.EQ.1) THEN
- pdd(1,1)=qu2*(dot(h1u,h2ub)+dot(h1ub,h2u))
- & +qd2*(dot(h1d,h2db)+dot(h1db,h2d))
- ELSEIF (npwz.EQ.2) THEN
- pdd(1,1)=w2*(prod(h2ub,v2,h1d)+prod(h1ub,v2,h2d))
- ELSEIF (npwz.EQ.3) THEN
- pdd(1,1)=w2*(prod(h1u,v2,h2db)+prod(h2u,v2,h1db))
- ELSEIF (npwz.EQ.4) THEN
- pdd(1,1)=zu2*(dot(h1u,h2ub)+dot(h1ub,h2u))
- & +zd2*(dot(h1d,h2db)+dot(h1db,h2d))
- ENDIF
- pdd(1,2)=pdd(1,1)
-
-C 2. com
- IF (npwz.EQ.1) THEN
- pdd(2,1)=qu2*sh1u+qd2*sh1d
- pdd(2,2)=qu2*sh2u+qd2*sh2d
- ELSEIF (npwz.EQ.2) THEN
- pdd(2,1)=w2*(dot(sv2u,h1d)+dot(h1ub,sv2d))
- pdd(2,2)=w2*(dot(sv2u,h2d)+dot(h2ub,sv2d))
- ELSEIF (npwz.EQ.3) THEN
- pdd(2,1)=w2*(dot(h1u,sv2d)+dot(sv2u,h1db))
- pdd(2,2)=w2*(dot(h2u,sv2d)+dot(sv2u,h2db))
- ELSEIF (npwz.EQ.4) THEN
- pdd(2,1)=zu2*sh1u+zd2*sh1d
- pdd(2,2)=zu2*sh2u+zd2*sh2d
- ENDIF
- pdd(2,1)=pdd(2,1)*b(0)
- pdd(2,2)=a(0)*pdd(2,2)
-
-C 3. fus
- IF (npwz.EQ.1) pdd(3,1)=sqf2
- IF (npwz.EQ.2.OR.npwz.EQ.3) pdd(3,1)=w2*sv2
- IF (npwz.EQ.4) pdd(3,1)=szf2
- pdd(3,1)=pdd(3,1)*a(0)*b(0)
- pdd(3,2)=pdd(3,1)
-
-C 4. qqb34 squared
- dd=dot(h1u,h2ub)+dot(h1d,h2db)+dot(h1ub,h2u)+dot(h1db,h2d)
- IF (npwz.EQ.1) pdd(4,1)=sqf2*dd
- IF (npwz.EQ.2.OR.npwz.EQ.3) pdd(4,1)=w2*sv2*dd
- IF (npwz.EQ.4) pdd(4,1)=szf2*dd
- pdd(4,2)=pdd(4,1)
-
-C 5. qqb56, qqb78, qq12, qq34, qq56, qq78, squared
- IF (npwz.EQ.1) THEN
- pdd(5,1)=qu2*sh1u+qd2*sh1d
- pdd(5,2)=qu2*sh2u+qd2*sh2d
- ELSEIF (npwz.EQ.2) THEN
- pdd(5,1)=w2*(dot(sv2u,h1d)+dot(h1ub,sv2d))
- pdd(5,2)=w2*(dot(sv2u,h2d)+dot(h2ub,sv2d))
- ELSEIF (npwz.EQ.3) THEN
- pdd(5,1)=w2*(dot(h1u,sv2d)+dot(sv2u,h1db))
- pdd(5,2)=w2*(dot(h2u,sv2d)+dot(sv2u,h2db))
- ELSEIF (npwz.EQ.4) THEN
- pdd(5,1)=zu2*sh1u+zd2*sh1d
- pdd(5,2)=zu2*sh2u+zd2*sh2d
- ENDIF
- pdd(5,1)=pdd(5,1)*sh2
- pdd(5,2)=sh1*pdd(5,2)
-
-C 6. qqb1256 and qqb1278
- pdd(6,1)=pdd(1,1)
- pdd(6,2)=pdd(1,2)
-
-C 7. qqb3456 and qqb3478
- IF (npwz.EQ.1.OR.npwz.EQ.4) THEN
- pdd(7,1)=pdd(1,1)
- pdd(7,2)=pdd(1,2)
- ELSEIF (npwz.EQ.2) THEN
- pdd(7,1)=w2*(dot3(sv2u,h1d,h2db)+dot3(h1ub,h2u,sv2d))
- pdd(7,2)=w2*(dot3(sv2u,h1db,h2d)+dot3(h1u,h2ub,sv2d))
- ELSEIF (npwz.EQ.3) THEN
- pdd(7,1)=w2*(dot3(sv2u,h1db,h2d)+dot3(h1u,h2ub,sv2d))
- pdd(7,2)=w2*(dot3(sv2u,h1d,h2db)+dot3(h1ub,h2u,sv2d))
- ENDIF
-
-C 8. qq1256 and qq3478
- IF (npwz.EQ.1) THEN
- pdd(8,1)=(qu2*(dot(h1u,h2u)+dot(h1ub,h2ub))
- & +qd2*(dot(h1d,h2d)+dot(h1db,h2db)))/2.
- pdd(8,2)=pdd(8,1)
- ELSEIF (npwz.EQ.2) THEN
- pdd(8,1)=w2*(prod(h2u,v2,h1d)+prod(h1ub,v2,h2db))/2.
- pdd(8,2)=w2*(prod(h1u,v2,h2d)+prod(h2ub,v2,h1db))/2.
- ELSEIF (npwz.EQ.3) THEN
- pdd(8,1)=w2*(prod(h1u,v2,h2d)+prod(h2ub,v2,h1db))/2.
- pdd(8,2)=w2*(prod(h2u,v2,h1d)+prod(h1ub,v2,h2db))/2.
- ELSEIF (npwz.EQ.4) THEN
- pdd(8,1)=(zu2*(dot(h1u,h2u)+dot(h1ub,h2ub))
- & +zd2*(dot(h1d,h2d)+dot(h1db,h2db)))/2.
- pdd(8,2)=pdd(8,1)
- ENDIF
-
-C 9. qq1278 and qq3456
- IF (npwz.EQ.1.OR.npwz.EQ.4) THEN
- pdd(9,1)=pdd(8,1)
- ELSEIF (npwz.EQ.2) THEN
- pdd(9,1)=w2*(dot3(sv2u,h1d,h2d)+dot3(h1ub,h2ub,sv2d))/2.
- ELSEIF (npwz.EQ.3) THEN
- pdd(9,1)=w2*(dot3(h1u,h2u,sv2d)+dot3(sv2u,h1db,h2db))/2.
- ENDIF
- pdd(9,2)=pdd(9,1)
-
-C 10. qqb5678LL and qq1234LR
- IF (npwz.EQ.1) THEN
- tmp1=qu*(sum(h1u)-sum(h1ub))+qd*(sum(h1d)-sum(h1db))
- tmp2=qu*(sum(h2ub)-sum(h2u))+qd*(sum(h2db)-sum(h2d))
- pdd(10,1)=tmp1*tmp2
- ELSEIF (npwz.EQ.2.OR.npwz.EQ.3) THEN
- pdd(10,1)=0d0
- ELSEIF (npwz.EQ.4) THEN
- tmp1=xlu*sum(h1u)+xld*sum(h1d)-xru*sum(h1ub)-xrd*sum(h1db)
- tmp2=xlu*sum(h2ub)+xld*sum(h2db)-xru*sum(h2u)-xrd*sum(h2d)
- pdd(10,1)=tmp1*tmp2/2.
- tmp1=xlu*sum(h1ub)+xld*sum(h1db)-xru*sum(h1u)-xrd*sum(h1d)
- tmp2=xlu*sum(h2u)+xld*sum(h2d)-xru*sum(h2ub)-xrd*sum(h2db)
- pdd(10,1)=pdd(10,1)+tmp1*tmp2/2.
- ENDIF
- pdd(10,2)=pdd(10,1)
-
-C 11. qqb5678LR and qq1234LL
- IF (npwz.LT.4) THEN
- pdd(11,1)=pdd(10,1)
- ELSEIF (npwz.EQ.4) THEN
- tmp1=xlu*sum(h1u)+xld*sum(h1d)-xru*sum(h1ub)-xrd*sum(h1db)
- tmp2=xru*sum(h2ub)+xrd*sum(h2db)-xlu*sum(h2u)-xld*sum(h2d)
- pdd(11,1)=tmp1*tmp2/2.
- tmp1=xlu*sum(h1ub)+xld*sum(h1db)-xru*sum(h1u)-xrd*sum(h1d)
- tmp2=xru*sum(h2u)+xrd*sum(h2d)-xlu*sum(h2ub)-xld*sum(h2db)
- pdd(11,1)=pdd(11,1)+tmp1*tmp2/2.
- ENDIF
- pdd(11,2)=pdd(11,1)
-
-C 12. qqb1234
- IF (npwz.LT.4) THEN
- pdd(12,1)=0d0
- ELSEIF (npwz.EQ.4) THEN
- pdd(12,1)=dot(h1u,h2ub)+dot(h1ub,h2u)
- pdd(12,1)=pdd(12,1)-dot(h1d,h2db)-dot(h1db,h2d)
- pdd(12,1)=pdd(12,1)*(xru-xlu)*sxrml/2.
- ENDIF
- pdd(12,2)=pdd(12,1)
-
-C 13. COM -- triangle loop contributions -- corrected 2011
- IF (npwz.LT.4) THEN
- pdd(13,1)=0d0
- pdd(13,2)=0d0
- ELSEIF (npwz.EQ.4) THEN
- pdd(13,1)=sum(h1u)+sum(h1ub)-sum(h1d)-sum(h1db)
- pdd(13,1)=pdd(13,1)*(xru-xlu)*sxrml*b(0)/2.
- pdd(13,2)=sum(h2u)+sum(h2ub)-sum(h2d)-sum(h2db)
- pdd(13,2)=a(0)*pdd(13,2)*(xru-xlu)*sxrml/2.
- ENDIF
-
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE vector (a,b,c,vec)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION vec(3)
- vec(1)=a
- vec(2)=b
- vec(3)=c
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION sum (a)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION a(3)
- sum=a(1)+a(2)+a(3)
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION dot (a,b)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION a(3),b(3)
- dot=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION dot3 (a,b,c)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION a(3),b(3),c(3)
- dot3=a(1)*b(1)*c(1)+a(2)*b(2)*c(2)+a(3)*b(3)*c(3)
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION prod (a,b,c)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION a(3),b(3,3),c(3)
- prod=0d0
- DO 1 i=1,3
- DO 1 j=1,3
-1 prod=prod+a(i)*b(i,j)*c(j)
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION xclow (nproc)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- GOTO (1,2) nproc
-C q qbar -> v G
-1 xclow=u/t+t/u+2.*s*qq/t/u
- RETURN
-C q G -> V q
-2 xclow=-s/t-t/s-2.*u*qq/s/t
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION xcds2 (nproc)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- sti=d3
- sui=d4
- tui=-d9
- GOTO (1,2) nproc
-1 CONTINUE
-C q qbar -> v G
- t0ut = u/t+t/u+2.*qq*s/u/t
- realp =
- & t0ut*(-(11./6.*ca-xnf/3.)*fa+67./18.*ca-5./9.*xnf
- & +ca*fa**2+2.*cf*(ft+fu-2.*fa)*fm
- & +(cf-ca/2.)*(pi**2./3.+(2.*fa+fs-fu-ft)**2)-3.*cf*fm)
- virtp =
- & t0ut*(-8.*cf-cf*fs**2+(2.*cf/3.-ca/6.)*pi**2
- & +ca/2.*(fs**2-(ft+fu)**2)+(11./6.*ca-xnf/3.)*fmu)
- & +cf*(s*(sti+sui)+(s+t)/u+(s+u)/t)
- & +ft*(cf*(4.*s**2+2.*s*t+4.*s*u+u*t)*sui**2+ca*t*sui)
- & +fu*(cf*(4.*s**2+2.*s*u+4.*s*t+t*u)*sti**2+ca*u*sti)
- & +(2.*cf-ca)*(2.*fs*(s**2*tui**2+2.*s*tui)
- & -qq*(u**2+t**2)/u/t*tui)
- & -(2.*cf-ca)*((s**2+(s+u)**2.)/u/t*fr1(qq,s,t)
- & +(s**2+(s+t)**2)/t/u*fr1(qq,s,u))
- & +ca*t0ut*fr2(qq,t,u) + 0.0
- xcds2 = realp + virtp
- RETURN
-2 CONTINUE
-C q G -> v G
- t0st = s/t+t/s+2.*qq*u/s/t
- realp =
- & -t0st*(cf*(7./2.+2*fm*(fu-fa)+fa**2-3./2.*(fm+fa))
- & +ca*(pi**2/6.+(fs-ft-fu)**2/2.+2.*ft*(fm-fa)
- & +2.*fa*(fs-fu+fa-fm))-fm*(11./6.*ca-xnf/3.))
- virtp =
- & t0st*(-8.*cf-cf*fu**2-1./3.*(cf-ca)*pi**2
- & +ca/2.*(fu**2-fs**2-ft**2)+(11./6.*ca-xnf/3.)*fmu)
- & +cf*(u*(tui+sui)+(u+t)/s+(s+u)/t)
- & +ft*(cf*(4.*u**2+2.*u*t+4.*s*u+s*t)*sui**2+ca*t*sui)
- & +fs*(cf*(4.*u**2+2.*s*u+4.*u*t+t*s)*tui**2+ca*s*tui)
- & +(2.*cf-ca)*(2.*fu*(u**2*sti**2+2.*u*sti)
- & -qq*(s**2+t**2)/s/t*sti)
- & -(2.*cf-ca)*((u**2+(u+s)**2)/s/t*(ft*fu-pi**2
- & +pi**2/2.-fr2(qq,t,u))+(u**2+(t+u)**2)/s/t*fr1(qq,s,u))
- & -ca*t0st*fr1(qq,s,t)
- virtp = - virtp
- xcds2 = realp + virtp
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION triang (nproc)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- GOTO (1,2) nproc
-1 CONTINUE
-C q qbar -> v G
- triang = -(s+qq)/(s-qq)*(1.-qq/(s-qq)*fs)
- RETURN
-
-2 CONTINUE
-C q G -> v G
- triang = (u+qq)/(u-qq)*(1.-qq/(u-qq)*fu)
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION fr1 (qq,s,t)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- fr1 = dlog(s/qq)*dlog(t/(qq-s))+0.5*dlog(qq/s)**2
- & -0.5*dlog((qq-t)/qq)**2+dilog(qq/s)-dilog(qq/(qq-t))
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION fr2 (qq,t,u)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- fr2 = 0.5*dlog((qq-t)/qq)**2+0.5*dlog((qq-u)/qq)**2
- & +dilog(qq/(qq-t))+dilog(qq/(qq-u))
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION xcs2a (nproc)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- GOTO (1,2) nproc
-1 CONTINUE
- t0ut = u/t+t/u+2.*qq*s/u/t
- anns2a =
- & t0ut*(xnf/3.-11./6.*ca+2.*cf*(ftu-2.*fm)
- & +(cf-ca/2.)*(4.*fstu-2.*ftu)
- & +fs2*(8.*cf-2.*ca))
- xcs2a = anns2a/s2
- RETURN
-2 CONTINUE
- t0st = s/t+t/s+2.*qq*u/s/t
- coms2a =
- & -t0st*(cf*(-3./2.+fsu+2.*ftu-2.*fm+(t+u)*fla*d9)
- & +ca*(2.*fstu+(fst-fsu)/2.-ftu-2.*fm))
- & -fs2*t0st*(2.*cf+4.*ca)
- & +(cf-ca/2)*((fla*(t+u)*d9+fsu)*(s+2.*u)/t
- & +2.*ftu*(t+2.*u)/s)
- xcs2a = coms2a/s2
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION xcfin (nproc)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- COMMON /const/ pi,tonbs,alpha,sin2tw,cf,ca,xnc,qw
- & ,xmw,xmz,xmc,xmb,xmt,qu,qd,xlu,xld,xru,xrd,v2(3,3)
- COMMON /variab/ e,q,qt,y,scalm,scalmu,qq,qt2,eqt,xlambd,escal2
- & ,sh,th,uh,aslo,asnlo,xnf,sqf2,sv2,sxrml,szf2,sv2d(3),sv2u(3)
- & ,x1,s2,s,t,u,d1,d2,d3,d4,d5,d6,d9,d10
- & ,fa,fm,fmu,fs,ft,fu,fs2,fst,fsu,ftu,fstu,fla,flt,flu
-
- GOTO (1,2,3,4,5,6,7,8,9,10,11,12) nproc
-
-1 CONTINUE
-C 1. ann and qqb12
- xcfin =
- & cf*(s*d1**2+2.*d1-s/u/t+(2.*qq*u+t*s2-2.*qq*s2)*d6/t
- & -2.*qq*(t-u)*d1/t/u)
- & +ca*(-11./6.*s/u/t+s**2*(3.*s2-4.*t)*d1**2/2./t/u
- & +2.*s*d1/u+qq/3./t**2)
- & +xnf/3.*(s/u/t-qq/t**2)
- & +2.*(cf-ca/2.)*(fstu+fs2-ftu)
- & *(qq-u)**2*d1/u*(d2-1/t)
- & +2.*(cf-ca/2)*(fstu+fs2)*(s-qq)/t/u
- & +cf*ftu*(4.*qq*(qq-t)**2*d6/u/t+(2.*(qq+s)-s2)*d6
- & +(s2-2.*s)/u/t)
- & -ca*ftu*qq/u/t
- & +cf*(fm-fs2)*(d6/u/t*(4.*u*t*(u-qq)-4.*qq*(u-qq)**2
- & -u*t*s2)+(qq-u)*d1**2-(2.*qq-u)*d1/t-2./t+qq/t**2
- & +4.*qq/u/t) + 0.0
- RETURN
-
-2 CONTINUE
-C 2. com
- xcfin =
- & fla*d9*d10*(d10*3./4.*cf*s*(s+qq-s2)*(u**2-t**2)*(1-u/t)
- & +cf*((u**2-t**2)*((qq-s2)/4./s+11./4.+s2/t-u/2./t)
- & -2.*s*(t+u**2/t)-s*(s-s2)*(u/t-3.)/2.+4.*s2*(t-u)+5.*s*u)
- & +ca*(1.-u/t)*((t+u)*(s+qq-s2)/4.+s*(s-s2)))
- & +fla*d9*(cf*((t-u-2.*s2)/4./s+(11.*s-2.*u-4.*s2)/4./t
- & +(2.*(1.+u/t)*(5.*qq+s2)-16.*qq*s2/t)/s)
- & -ca/2.*((1.+u/t)/2.+5.*(u/s-(u+s2)/t)-7.*(s+u*s2/s)/t
- & +3.*(t-3.*s2)/s+2.*(u**2+3.*s2**2)/s/t))
- & +d10*(cf*(3./2.*d10*s*(1.+u/t)*(t-u)**2+s*(u-3.*t)/2./t
- & +(1.-u/t)*(s2-7.*u/4.)+(t-u)*(11./4.+(3./2.*(t+u)-2.*s2)/s))
- & +ca*(u/t-1.)*(s+s2))
- & +2.*d1**3*s*t*(4.*cf-3.*ca-(fs2-fm)*(cf-ca))
- & +d1**2*(cf*(3.*s-8.*t-4.*u)-ca*(9.*s-4.*(t-u))
- & +2.*(fs2-fm)*(cf*(t+u)+ca*(3.*s+2.*u)))
- xcfin=xcfin
- & +d1*((cf*(fs2-fm)-ca*(fs2-fm+1./2.*(fstu+fs2
- & -ftu+(fst+flt)/2.)))*(4.*(1.+u/s)*(1.-u/t)-2.*(s/t+t/s))
- & -cf*(fs2-fm)*(s+2.*u)/t
- & +ca*((2.*flt+2.*fst)*(u+s)/t+2.*(fs2-fm)*((5.*s+4.*u)/t-2.))
- & -cf*(4.+(s-2.*u)/t+(2.*u-3.*t)/s-u**2/s/t)
- & +ca*(7.-(3.*s+2.*u)/t))
- & +d2*(2.*cf*(fm-fs2-s/t)+ca*(s/t-1.))
- & +d3*cf*(fsu-2.*ftu-flt)*(1.-s2/t)*(s2**2-2.*u*(s2-u))/s**2
- & +d6*cf*(2.*(fm-fs2-ftu)*(2.*s/t*(qq+2.*u)+t-s2+4.*u/t*(2.*u-s2)
- & +2.*(u-s2)*(t-2.*u+2.*u**2/t)/s)
- & +s2-t+4.*u*(qq/t+(1.-u/t)*(t-qq)/s))
- xcfin=xcfin
- & +cf*((fsu-2.*ftu-flt)*((2.*u*(s2-u)-s2**2)/s-s2)/s/t
- & -2.*ftu*(1./t-2./s+4.*u/s/t)+flt*(u-3.*s2+5.*qq)/s/t
- & +fsu*(5./s+(1.+2.*t/s-3.*s2/s)/u-s2*qq/s/u**2)
- & +(fs2-fm)*(8./s+3./t+1./u-4.*u/s/t+(2.*t-3.*s2)/s/u
- & -qq*(1./t**2+s2/s/u**2))-2./s+3./4./t-1./u+3.*u/s/t
- & -(t+s2)/s/u+qq*(s2/s/u**2-1./2./t**2))
- & +ca*(2.*(fs2-fm)*(6./s-2./t-2.*qq/s/t-3.*s2/s/t)
- & -2.*fst*(-15./4./s+9.*s2/4./s/t+qq/t**2)
- & +2.*(fs2-fm+fst)*((s+s2**2/s)/t**2
- & +2.*s2*qq/t**3-4.*qq/t*(1.-qq/t)/s
- & +2.*s2/s/t*qq/t*(1.-s2/t)*(1.-qq/t))
- & +(fs2+fstu-(fst+fsu)/2.)*(4.-2.*u/t-s2/t)/s+fsu*(1.-u/t)/s
- & +ftu*(2.*u/t-1.)/s-2.*flt*(3./t+1./s+2.*(u-s2)/s/t)
- & -1./s-1./t+2.*s2/s/t+(3.*qq+4.*s2*(1.-3.*qq/t))/t**2
- & -2.*(s+s2**2./s)/t**2-12.*s2/s/t*qq/t*(1.-s2/t)*(1.-qq/t)
- & +2./s*qq/t*(1.-qq/t)) + 0.0
- RETURN
-
-3 CONTINUE
-C 3. fus
- xcfin =
- & ca*d5*d9*fla*(d10*t**2*(3./2.*d10*(2.*s-t-u)*(u**2-t**2)
- & +(2.*s-3.*t+u+(t**2-u**2)/s))-3.*s/2.-t*(9./2.+5.*t/s+3.*u/s))
- &+ca*d9*fla*(d10*t**2*(3./2.*d10*(u**2-t**2)+1.)*(qq-s2)/s
- & -11.*(qq-s2)/s/4.-4.*(1.+t/s))
- &+cf*d9*fla*(2.*d5*(t+u-2.*s)+6.)
- &+ca*d10*t*(3.*d10*(t-u)**2*(1.+(t+u)/2./s-2.*s*d5)
- & +d5*(2.*s-3.*t+3.*u)-1.-u/s+s2*(t-u)/s**2
- & +(d5-1./2./s)*(t-u)**2/s)
- &+cf*4.*d6*((fs2-fm+ftu)*(s/2.+2.*t*(1.+t/s))-t*(1.+(t+u)/s))
- &+d5*2*(d1*(ca-2.*cf)*t**2*(fst+flt)/s
- & +d3*cf*(fsu-flt-2.*ftu)*(s+2.*u*(1.+u/s)))
- &+d5*(2.*cf*(fst-(3.+4.*(t+u)/s)*flt-(4.+8.*t/s)*ftu)
- & +ca/2.*((fst+flt)*(2.+(u+3.*t)/s)+1.+t*(t-u)/s**2))
- xcfin = xcfin
- &+2.*cf*d3*(2.*ftu+flt-fsu)*(2.+(t+u)/s)
- &+8.*t*d1**2*(cf*(1.-fs2+fm)+ca)
- &+4.*(cf-ca/2.)*d1*d2*(fs2-ftu+fstu)*t*u/s
- &+d1*((ca-2.*cf)*((2.+(t+u)/s)*(fst+flt)+2.*(u-t)/s
- & *(fs2+fstu-ftu))-8.*cf*(fs2-fm-1.)+4.*ca*(2.+(u-t)/s))
- &+cf*(2.*(fs2-fm)*(2./s-2./t-(u+s2*qq/t)/s/t)
- & +4.*(fs2+fstu)/s-2.*fst/t*(2.+(u+s2*qq/t)/s)
- & +(4.*qq*(1.+s2/t)-2.*u)/s/t)
- &-ca/s*(2.*(fs2+fstu)+fst/2.+5./2.*flt+15./4.) + 0.0
- RETURN
-
-4 CONTINUE
-C 4. qqb34
- xcfin =
- & d5/2.*(u/s-1.)-5./4./s
- & +d10*(u*d5*(-2.*s+3./2.*u/s*(t-u)+4.*u-2.*t)
- & +u/2./s*(2.*s+2.*s2+t-u))
- & +d10**2*3.*u**2*(u-t)*(d5*(2.*s-u-t)-(s+s2)/s)
- & +d9*fla*(d5*(2.*u**2/s+5./2.*u+3./2.*s)+(3./4.*s+u-s2/2.)/s)
- & +d9*d10*fla*(u**2*d5*(3.*u-t-u**2/s+t**2/s-2.*s)
- & +u/s*(2.*s2*t-u*t-2.*s2**2+4.*u*s2-3.*u**2+2.*s*s2-u*s))
- & +d9*d10**2*fla*(3.*d5*(u**2-t**2)*u**2*(s-qq)
- & +3.*u**2*qq/s*(u-t)*(u+t-2.*s2)) + 0.0
- xcfin=xcfin/2.
- RETURN
-
-5 CONTINUE
-C 5. qqb56
- xcfin =
- & d9*d10*s*fla*(u**2-t**2)/t
- & +d9/t*fla*(3.*s+2.*u)
- & +d10*(t+u-2.*s2)*(1.-u/t)
- & +(fs2-fm)*d1*(s*d1+4.*s/t+2.*u/t)
- & +(fst-fm+fs2)*((s+s2-qq)**2/s/t**2
- & +(u/t-2.*s2*qq/t**2)**2/s+4.*(qq/t-1.)/t)
- & +fst*2*(1.-qq/t)/t
- & -d1*(s*d1-1.-u/t)+2.*(u/t+s2*(1./t+1./s-s2/t/s))/t
- & +(qq/t-1.)/t+12.*s2*qq*(u-s2*qq/t)/s/t**3 + 0.0
- xcfin=xcfin/2.
- RETURN
-
-6 CONTINUE
-C 6. qqb1256
- xcfin =
- & d9*fla*(d10**2*3.*s**2*qq*(t-u)**2*(1./t+1./u)
- & +d10*s*qq*(5.*s/t-7.*s/u+4.*u/t-4.)
- & +(2.*s**2/u+(qq*(2.*s+u)-u*s)/t)/s2+(2.*s2-u)/t-2.)
- & +d10**2*3.*s*qq*(2.*s2-u-t)*(t-u)**2/u/t
- & +d10*(s*(s-s2)*(7./u-5./t)+(t/u-u/t)*(7.*s+t+3.*u-2.*s2)/2.
- & +3.*s*(t/u-u/t)/2.+2.*s-2.*s2*(1.-u/t))
- & +d1*s**2*(d1+3./t)/u+fst*((2.*s*(1.+s/u)+u)/s2-2.*qq/t)/t
- & -1./2./u+3./2./t-3.*qq/t**2
- xcfin=xcfin*(cf-ca/2.)
- RETURN
-
-7 CONTINUE
-C 7. qqb3456
- xcfin =
- & d9*fla*(6.*s*s2*qq*d10**2*(t-u)**2/t
- & +d10*(2.*s*u+(s2*(1.+t/s)-t*(t+u)/2./s)
- & *(t+u)*(1.-u/t)+3.*s/2.*(t-u)*(1.-u/t)
- & -4.*s2*qq*(1.+(s-u)/t))
- & +d5*(6.*s+2.*u+(t**2-u**2)/2./s)
- & -1.+(s/2.+2.*u-s2)/t+(3.*t+2.*u-6.*s2)/s
- & +(u**2-2*u*s2+4.*s2**2)/s/t)
- & +d10**2*3.*s2*(t-u)**2*(1.+u/t-2.*qq/t)
- & +d10*(-t/2.+3.*u+s2+(t**2-u**2)/s-u*(u+6.*s2)/2/t)
- & +d5*(fst+flt)*(2.*(s+u)/t+(t+u**2/t)/2./s) +2.*d1
- & +fst/2.*(2./t+3./s+(u-2.*s2)/s/t)
- & +flt/2.*(-2./t+1./s-(u+2.*s2)/s/t)
- & +1./2./t-1./s+(2.*u-4.*s2*qq/t)/s/t + 0.0
- xcfin=xcfin*(cf-ca/2.)
- RETURN
-
-8 CONTINUE
-C 8. qq1256
- xcfin =
- & (fst+flt)*(4.*d1*d5*s**2/t-2.*d1*(1.-u/t)-4./t)
- & +4.*qq*fst/t**2 + 0.0
- xcfin=xcfin*(cf-ca/2.)
- RETURN
-
-9 CONTINUE
-C 9. qq1278
- xcfin =
- & (2.*fs2-fst-fsu+2.*fstu)*(-2.-u/t-t/u+2.*s2*(1./u+1./t
- & -s2/u/t))/s
- & -4.*(1.+(u/t+t/u)/2.-s2*qq*(1./t**2+1./u**2))/s
- xcfin=xcfin*(cf-ca/2.)
- RETURN
-
-10 CONTINUE
-C 10. qqb5678LL
- xcfin =
- & d9*fla*(-d10*s/t*(t-u)**2+2.*d5*(s-t)+3.*s/t+2.*s2/t-1.)
- & +d10*(u/t-1.)*(u-2.*s2)
- & +d3*d5*(fsu-flt-2.*ftu)*(2.*s*(s+u)+u**2)/t
- & +d5*((fst-ftu)*(1.+u/t+2.*t/u+4.*s*(1./t+1./u+s/t/u))
- & -(flt+ftu)*(1+u/t))
- & +d3*(2.*ftu-fsu+flt)*(1.+2.*s/t+u/t) + 2.*d1*s/u
- & +ftu*(s+s2)/t/u-(fst+flt+1.-s/u)/t
- xcfin = xcfin/2.
- RETURN
-
-11 CONTINUE
-C 11. qqb5678LR
- xcfin =
- & 2.*d5*(d3*(fsu-flt-2.*ftu)+2.*(fst-ftu)/u)*s**2/t
- & +d3*(1.+2.*s/t-u/t)*(2.*ftu-fsu+flt)
- & +2.*(ftu-fsu)*(s+s2-u)/t/u - 2.*(ftu+flt)/t
- xcfin = xcfin/2.
- RETURN
-
-12 CONTINUE
-C 12. qqb1234
- xcfin =
- & d9*fla*(d10*(d10*3.*s*qq*(2.*s+t+u)*(t-u)**2/t
- & +qq*(2.*s+t+u)*(u-t-s)/t) - qq/t)
- & +d10*(-d10*3./2./t*(s+qq-s2)*(2.*s+t+u)*(t-u)**2
- & +(2.*s+t+u)*(3./2.+(s-s2)/t-u/2./t))
- xcfin = xcfin/2.
- RETURN
-
- END
-C-----------------------------------------------------------------------
-C Dilogarithm Li2(x) for -1 <= x <= 1
-C
- DOUBLE PRECISION FUNCTION dilog (x)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- EXTERNAL dilitg
- pi=3.141592654
- IF (DABS(x).GT.1.d0) PAUSE
- & 'Dilog: Argument outside range [-1,1]'
- IF (x.EQ.1.d0) THEN
- dilog = pi**2/6.
- RETURN
- ENDIF
- IF (x.EQ.0.d0) THEN
- dilog=0.
- RETURN
- ENDIF
- xx=x
- IF (x.GT..5d0) xx=1.d0-x
- IF (x.LT.0.d0) xx=-x/(1.-x)
- CALL dilgau (dilitg,0.d0,xx,dilog)
- IF (x.GT..5d0) THEN
- dilog=-dilog+pi**2/6.-dlog(x)*dlog(xx)
- ELSEIF (x.LT.0.d0) THEN
- dilog=-dilog-.5d0*(dlog(1.-x))**2
- ENDIF
- RETURN
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION dilitg (x)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- dilitg = -dlog(1.d0-x)/x
- RETURN
- END
-C-----------------------------------------------------------------------
-C Gaussian Integrator - from Numerical Recipes - page 122
-C
- SUBROUTINE dilgau (func,a,b,ss)
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- DIMENSION x(5), w(5)
- DATA x /.1488743389d0,.4333953941d0,.6794095682d0,
- & .8650633666d0,.9739065285d0/
- DATA w /.2955242247d0,.2692667193d0,.2190863625d0,
- & .1494513491d0,.0666713443d0/
- xm=0.5d0*(b+a)
- xr=0.5d0*(b-a)
- ss=0.d0
- DO 11 j=1,5
- dx=xr*x(j)
- ss=ss+w(j)*(func(xm+dx)+func(xm-dx))
-11 CONTINUE
- ss=xr*ss
- RETURN
- END
diff --git a/gonsalves/qt.inp b/gonsalves/qt.inp
deleted file mode 100644
index 0cad6d3..0000000
--- a/gonsalves/qt.inp
+++ /dev/null
@@ -1,21 +0,0 @@
- Input parameters for qt.f: Inclusive gamma*, W, Z production
- Hadrons: lppb lns e pdset
- t f 1800 mrs02nlo
- Bosons: npwz lwpm q qt y
- 2 t 80 50 0
- Partons: lallpr 1 2 3 4 5 6 7 8 9 10 11 12
- t f t t f f f f f f f f f
- QCD: nscale nalpha lord1 lord2
- 1 1 t t
- Scales: lscaleq lscalm lscalmu scalm scalmu
- t f f 0.0 0.0
- Thresholds: linctri
- f
- Integrate: ndim ntot begin step npoints
- 3 0 10.0 20.0 10
- Vegas: ncall itmx0 itmx ndev nprn alph0 alph
- 40000 3 7 2 -1 1.5 1.5
- Ran Seed: iseed
- 13579
- Output: lterm lauto ofname comment
- t f qt.out #
diff --git a/gonsalves/vegas.f b/gonsalves/vegas.f
deleted file mode 100644
index 094db7e..0000000
--- a/gonsalves/vegas.f
+++ /dev/null
@@ -1,328 +0,0 @@
-C-----------------------------------------------------------------------
-C Vegas
-C
- BLOCK DATA vegdat
-C
-C makes default parameter assignments for vegas
- IMPLICIT DOUBLE PRECISION(a-h,o-z)
- COMMON/bveg1/ncall,itmx,nprn,ndev,xl(10),xu(10),acc
- COMMON/bveg2/it,ndo,si,swgt,schi,xi(50,10)
- COMMON/bveg3/alph,ndmx,mds
- COMMON/rnsd/iseed
- DATA ncall/5000/,itmx/5/,nprn/5/,acc/-1./,
- 1 xl/0.,0.,0.,0.,0.,0.,0.,0.,0.,0./,
- 2 xu/1.,1.,1.,1.,1.,1.,1.,1.,1.,1./,
- 3 alph/1.5/,ndmx/50/,mds/1/,ndev/6/,
- 4 ndo/1/,xi/500*1./,it/0/,si,swgt,schi/3*0./
-
- END
-C-----------------------------------------------------------------------
- SUBROUTINE VEGAS(ndim,fxn,avgi,sd,chi2a)
-C
-C Subroutine performs ndim-dimensional Monte Carlo integ'n
-C - by G.P. Lepage Sept 1976/(rev)Aug 1979
-C - algorithm described in J Comp Phys 27,192(1978)
-C
- IMPLICIT DOUBLE PRECISION(a-h,o-z)
- COMMON/bveg1/ncall,itmx,nprn,ndev,xl(10),xu(10),acc
- COMMON/bveg2/it,ndo,si,swgt,schi,xi(50,10)
- COMMON/bveg3/alph,ndmx,mds
- COMMON/bveg4/calls,ti,tsi
- DIMENSION d(50,10),di(50,10),xin(50),r(50),dx(10),ia(10),
- 1 kg(10),dt(10),x(10)
-C DIMENSION rand(10)
- double precision rand(10)
- DATA one/1d0/
- sqrt(a)=dsqrt(a)
- alog(a)=dlog(a)
- abs(a)=dabs(a)
-C
- ndo=1
- DO 1 j=1,ndim
-1 xi(1,j)=one
-C
- ENTRY vegas1(ndim,fxn,avgi,sd,chi2a)
-C - initializes cumulative variables, but not grid
- it=0
- si=0d0
- swgt=si
- schi=si
-C
- ENTRY vegas2(ndim,fxn,avgi,sd,chi2a)
-C - no initialization
- nd=ndmx
- ng=1
- IF(mds.EQ.0) GO TO 2
- ng=(ncall/2d0)**(1d0/ndim)
- mds=1
- IF((2*ng-ndmx).LT.0) GO TO 2
- mds=-1
- npg=ng/ndmx+1
- nd=ng/npg
- ng=npg*nd
-2 k=ng**ndim
- npg=ncall/k
- IF(npg.LT.2) npg=2
- calls=npg*k
- dxg=one/ng
- dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-one)
- xnd=nd
- ndm=nd-1
- dxg=dxg*xnd
- xjac=one/calls
- DO 3 j=1,ndim
- dx(j)=xu(j)-xl(j)
-3 xjac=xjac*dx(j)
-C
-C rebin, preserving bin density
- IF(nd.EQ.ndo) GO TO 8
- rc=ndo/xnd
- DO 7 j=1, ndim
- k=0
- xn=0.
- dr=xn
- i=k
-4 k=k+1
- dr=dr+one
- xo=xn
- xn=xi(k,j)
-5 IF(rc.GT.dr) GO TO 4
- i=i+1
- dr=dr-rC
- xin(i)=xn-(xn-xo)*dr
- IF(i.LT.ndm) GO TO 5
- DO 6 i=1,ndm
-6 xi(i,j)=xin(i)
-7 xi(nd,j)=one
- ndo=nd
-C
-8 IF(nprn.ge.0) WRITE(ndev,200) ndim,calls,it,itmx,acc,nprn,
- 1 alph,mds,nd,(xl(j),xu(j),j=1,ndim)
- IF(nprn.ge.0) WRITE(ndev,222)
-C
- ENTRY vegas3(ndim,fxn,avgi,sd,chi2a)
-C - main integration loop
-9 it=it+1
- ti=0d0
- tsi=ti
- DO 10 j=1,ndim
- kg(j)=1
- DO 10 i=1,nd
- d(i,j)=ti
-10 di(i,j)=ti
-C
-11 fb=0d0
- f2b=fb
- k=0
-12 k=k+1
- CALL randa(ndim,rand)
- wgt=xjaC
- DO 15 j=1,ndim
- xn=(kg(j)-rand(j))*dxg+one
- ia(j)=xn
- IF(ia(j).GT.1) GO TO 13
- xo=xi(ia(j),j)
- rc=(xn-ia(j))*xo
- GO TO 14
-13 xo=xi(ia(j),j)-xi(ia(j)-1,j)
- rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
-14 x(j)=xl(j)+rc*dx(j)
-15 wgt=wgt*xo*xnd
-C
- f=wgt
- f=f*fxn(x,wgt)
- f2=f*f
- fb=fb+f
- f2b=f2b+f2
- DO 16 j=1,ndim
- di(ia(j),j)=di(ia(j),j)+f
-16 IF(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2
- IF(k.LT.npg) GO TO 12
-C
- f2b=sqrt(f2b*npg)
- f2b=(f2b-fb)*(f2b+fb)
- ti=ti+fb
- tsi=tsi+f2b
- IF(mds.ge.0) GO TO 18
- DO 17 j=1,ndim
-17 d(ia(j),j)=d(ia(j),j)+f2b
-18 k=ndim
-19 kg(k)=mod(kg(k),ng)+1
- IF(kg(k).ne.1) GO TO 11
- k=k-1
- IF(k.GT.0) GO TO 19
-C
-C compute final results for this iteration
- tsi=tsi*dv2g
- ti2=ti*ti
- wgt=one/tsi
- si=si+ti*wgt
- swgt=swgt+wgt
- schi=schi+ti2*wgt
- avgi=si/swgt
- chi2a=(schi-si*avgi)/(it-.9999d0)
- sd=sqrt(one/swgt)
-C
- IF(nprn.LT.0) GO TO 21
- tsi=sqrt(tsi)
- WRITE(ndev,201) it,ti,tsi,avgi,sd,chi2a
- IF(nprn.EQ.0) GO TO 21
- DO 20 j=1,ndim
-20 WRITE(ndev,202) j,(xi(i,j),di(i,j),i=1+nprn/2,nd,nprn)
-C
-c refine grid
-21 DO 23 j=1,ndim
- xo=d(1,j)
- xn=d(2,j)
- d(1,j)=(xo+xn)/2.
- dt(j)=d(1,j)
- DO 22 i=2,ndm
- d(i,j)=xo+xn
- xo=xn
- xn=d(i+1,j)
- d(i,j)=(d(i,j)+xn)/3.
-22 dt(j)=dt(j)+d(i,j)
- d(nd,j)=(xo+xn)/2.
-23 dt(j)=dt(j)+d(nd,j)
-C
- DO 28 j=1,ndim
- rc=0.
- DO 24 i=1,nd
- r(i)=0d0
- IF(d(i,j).le.0.) GO TO 24
- xo=dt(j)/d(i,j)
- r(i)=((xo-one)/xo/alog(xo))**alph
-24 rc=rc+r(i)
- rc=rc/xnd
- k=0
- xn=0d0
- dr=xn
- i=k
-25 k=k+1
- dr=dr+r(k)
- xo=xn
- xn=xi(k,j)
-26 IF(rc.GT.dr) GO TO 25
- i=i+1
- dr=dr-rC
- xin(i)=xn-(xn-xo)*dr/r(k)
- IF(i.LT.ndm) GO TO 26
- DO 27 i=1,ndm
-27 xi(i,j)=xin(i)
-28 xi(nd,j)=one
-C
- IF(it.LT.itmx.AND.acc*abs(avgi).LT.sd) GO TO 9
-200 FORMAT(/35h Input parameters for Vegas: ndim=,i3,8h ncall=,f8.0
- 1 /28x,5h it+,i5,7h itmx=,i5/28x,6h acc=,g9.3
- 2 /28x,7h nprn=,i3,7h alph=,f5.2/28x,6h mds=,i3,6h nd=,i4
- 3 /28x,10h (xl,xu)=,(t40,2h( ,g12.6,3h , ,g12.6,2h )))
-222 FORMAT(
- 4 //1x,'Iter.#. Integral Std.-Dev. Accum.-Int. Std.-Dev. ',
- 5 'Chi2a',/)
-201 FORMAT(1x,i3,5x,5(g9.3,3x))
-*201 FORMAT(///21h integration by vegas//14h iteration no.,i3,
-* 1 14h: integral =,g14.8/24x,10hstd dev =,g10.4/
-* 2 34h accumulated results: integral =,g14.8/
-* 3 24x,10hstd dev =,g10.4/24x,17hchi**2 per it'n =,g10.4)
-202 FORMAT(/15h DATA for axis ,i2/25h x delta i ,
- 1 24h x delta i ,18h x delta i,
- 2 /(1h ,f7.6,1x,g11.4,5x,f7.6,1x,g11.4,5x,f7.6,1x,g11.4))
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE randa(n,rand)
-C subroutine generates uniformly distributed random no's x(i),i=1,n
- IMPLICIT NONE
- INTEGER n
- DOUBLE PRECISION rand(n),ran2,ranf
- EXTERNAL ran2
- INTEGER i,idum
- COMMON /ranseed/ idum
- DO i=1,n
-C rand(i)=ran2(idum)
- rand(i)=ranf(idum)
- ENDDO
- RETURN
- END
-C-----------------------------------------------------------------------
- BLOCK DATA seeddata
- IMPLICIT NONE
- INTEGER idum
- COMMON /ranseed/ idum
- DATA idum /-123456789/
- END
-C-----------------------------------------------------------------------
- SUBROUTINE getseed(seed)
- IMPLICIT NONE
- INTEGER seed,idum
- COMMON /ranseed/ idum
- seed=idum
- RETURN
- END
-C-----------------------------------------------------------------------
- SUBROUTINE setseed(seed)
- IMPLICIT NONE
- INTEGER seed,idum
- COMMON /ranseed/ idum
- IF (seed.GT.0) THEN
- idum=-seed
- ELSE
- idum=seed
- ENDIF
- RETURN
- END
-C-----------------------------------------------------------------------
- double precision FUNCTION ran2(idum)
-* from Numerical recipes
-* converted to double precision
- implicit none
- INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
- double precision AM,EPS,RNMX
- PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
- *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
- *NTAB=32,NDIV=1+IMM1/NTAB,EPS=3d-16,RNMX=1d0-EPS)
- INTEGER idum2,j,k,iv(NTAB),iy
- SAVE iv,iy,idum2
- DATA idum2/123456789/, iv/NTAB*0/, iy/0/
- if (idum.le.0) then
- idum=max(-idum,1)
- idum2=idum
- do 11 j=NTAB+8,1,-1
- k=idum/IQ1
- idum=IA1*(idum-k*IQ1)-k*IR1
- if (idum.lt.0) idum=idum+IM1
- if (j.le.NTAB) iv(j)=idum
-11 continue
- iy=iv(1)
- endif
- k=idum/IQ1
- idum=IA1*(idum-k*IQ1)-k*IR1
- if (idum.lt.0) idum=idum+IM1
- k=idum2/IQ2
- idum2=IA2*(idum2-k*IQ2)-k*IR2
- if (idum2.lt.0) idum2=idum2+IM2
- j=1+iy/NDIV
- iy=iv(j)-idum2
- iv(j)=idum
- if(iy.lt.1)iy=iy+IMM1
- ran2=min(AM*iy,RNMX)
- return
- END
-C-----------------------------------------------------------------------
- DOUBLE PRECISION FUNCTION RANF(ISEED)
-C Park-Miller generator using Schrage's algorithm
- IMPLICIT NONE
- INTEGER ISEED
- INTEGER IA,IC,IQ,IR
- PARAMETER (IA=16807,IC=2147483647,IQ=127773,IR=2836)
- INTEGER IH,IL,IT
- IH=ISEED/IQ
- IL=MOD(ISEED,IQ)
- IT=IA*IL-IR*IH
- IF(IT.GT.0) THEN
- ISEED=IT
- ELSE
- ISEED=IC+IT
- ENDIF
- RANF=ISEED/DBLE(IC)
- END