[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / wrapcteq5.f

      subroutine CTEQ5evolve(x,Q,pdf)
      implicit real*8(a-h,o-z)
      include 'parmsetup.inc'
      character*16 name(nmxset)
      integer nmem(nmxset),ndef(nmxset),mmem
      common/NAME/name,nmem,ndef,mmem
      integer nset
      real*8 pdf(-6:6)
      Character Line*80
      PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
      PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
      PARAMETER (M= 2, M1 = M + 1)
      Common
     > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
     > / CtqPar2 / Nx, Nt, NfMx
     > / XQrange / Qini, Qmax, Xmin
     > / QCDtable /  Alambda, Nfl, Iorder
     > / Masstbl / Amass(6)
      data pi / 3.141592653589793d0 /
      save
c
      U =         X * CtLhCtq5Pdf(1,X,Q)
      D =         X * CtLhCtq5Pdf(2,X,Q)
      USEA =      X * CtLhCtq5Pdf(-1,X,Q)
      DSEA =      X * CtLhCtq5Pdf(-2,X,Q)
      STR =       X * CtLhCtq5Pdf(3,X,Q)
      CHM =       X * CtLhCtq5Pdf(4,X,Q)
      BOT =       X * CtLhCtq5Pdf(5,X,Q)
      GLU  =      X * CtLhCtq5Pdf(0,X,Q)
      UPV=U-USEA
      DNV=D-DSEA
c      
      pdf(0)  = glu
      pdf(1)  = dnv+dsea
      pdf(-1) = dsea
      pdf(2)  = upv+usea
      pdf(-2) = usea
      pdf(3)  = str
      pdf(-3) = str
      pdf(4)  = chm
      pdf(-4) = chm
      pdf(5)  = bot
      pdf(-5) = bot
      pdf(6)  = 0.0d0
      pdf(-6) = 0.0d0

      return
*
      entry CTEQ5read(nset)

      call CtLhbldat1           !this line was missing in previous releases (jcp)
      call CtLhbldat2           !this line was missing in previous releases (jcp)
 
      read(1,*)nmem(nset),ndef(nset)
      Read  (1, '(A)') Line     
      Read  (1, '(A)') Line
      Read  (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
      Iorder = Nint(Dr)
      Nfl = Nint(Fl)
      Alambda = Al

      Read  (1, '(A)') Line 
      Read  (1, *) NX,  NT, NfMx

      Read  (1, '(A)') Line
      Read  (1, *) QINI, QMAX, (QL(I), I =0, NT)

      Read  (1, '(A)') Line
      Read  (1, *) XMIN, (XV(I), I =0, NX)

      Do 11 Iq = 0, NT
         QL(Iq) = Log (QL(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  (1, '(A)') Line
      Read  (1, *, IOSTAT=IRET) (UPD(I), I=1,Npts)
      Return
*    

      entry CTEQ5alfa(alfas,Qalfa)

       alfas = pi*CtLhALPI(Qalfa)

      return
*
      entry CTEQ5init(Eorder,Q2fit)

      return
*
      entry CTEQ5pdf(mem)
        imem = mem
      return
* 
      end

      FUNCTION CtLhPartonX5 (IPRTN, X, Q)
C
C   Given the parton distribution function in the array UPD in
C   COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
C   x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
      PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
      PARAMETER (M= 2, M1 = M + 1)
C
      Logical First
      Common 
     > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
     > / CtqPar2 / Nx, Nt, NfMx
     > / XQrange / Qini, Qmax, Xmin
C
      Dimension Fq(M1), Df(M1)

      Data First /.true./
      save First
C                                                 Work with Log (Q)
      QG  = LOG (Q/AL)

C                           Find lower end of interval containing X
      JL = -1
      JU = Nx+1
 11   If (JU-JL .GT. 1) Then
         JM = (JU+JL) / 2
         If (X .GT. XV(JM)) Then
            JL = JM
         Else
            JU = JM
         Endif
         Goto 11
      Endif

      Jx = JL - (M-1)/2
      If (X .lt. Xmin .and. First ) Then
         First = .false.
c         Print '(A, 2(1pE12.4))', 
c     >     ' WARNING: X << Xmin, extrapolation used; X,Xmin=',X,Xmin
         If (Jx .LT. 0) Jx = 0
      Elseif (Jx .GT. Nx-M) Then
         Jx = Nx - M
      Endif
C                                    Find the interval where Q lies
      JL = -1
      JU = NT+1
 12   If (JU-JL .GT. 1) Then
         JM = (JU+JL) / 2
         If (QG .GT. QL(JM)) Then
            JL = JM
         Else
            JU = JM
         Endif
         Goto 12
      Endif

      Jq = JL - (M-1)/2
      If (Jq .LT. 0) Then
         Jq = 0
c         If (Q .lt. Qini)  Print '(A, 2(1pE12.4))', 
c     >     ' WARNING: Q << Qini, extrapolation used; Q,Qini=',Q,Qini
      Elseif (Jq .GT. Nt-M) Then
         Jq = Nt - M
c         If (Q .gt. Qmax)  Print '(A, 2(1pE12.4))', 
c     >     ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
      Endif

      If (Iprtn .GE. 3) Then
         Ip = - Iprtn
      Else
         Ip = Iprtn
      EndIf
C                             Find the off-set in the linear array Upd
      JFL = Ip + NfMx
      J0  = (JFL * (NT+1) + Jq) * (NX+1) + Jx
C
C                                           Now interpolate in x for M1 Q's
      Do 21 Iq = 1, M1
         J1 = J0 + (Nx+1)*(Iq-1) + 1
         Call CtLhPolint3 (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
 21   Continue
C                                          Finish off by interpolating in Q
      Call CtLhPolint3 (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)

      CtLhPartonX5 = Ftmp
C
      RETURN
C                        ****************************
      END

      Function CtLhCtq5Pdf (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 CtLhCtq5Pdf: ', X
	Stop
      Endif
      If (Q .lt. Alambda) Then
	Print *, 'Q out of range in CtLhCtq5Pdf: ', 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 CtLhCtq5Pdf: '
     >              , Iparton
         Endif
         CtLhCtq5Pdf = 0D0
         Return
      Endif

      CtLhCtq5Pdf = CtLhPartonX5 (Iparton, X, Q)
      if(CtLhCtq5Pdf.lt.0.D0)  CtLhCtq5Pdf = 0.D0

      Return

C                             ********************
      End
Commit	Line	Data
3c5d1739	1	subroutine CTEQ5evolve(x,Q,pdf)
	2	implicit real*8(a-h,o-z)
	3	include 'parmsetup.inc'
	4	character*16 name(nmxset)
	5	integer nmem(nmxset),ndef(nmxset),mmem
	6	common/NAME/name,nmem,ndef,mmem
	7	integer nset
	8	real*8 pdf(-6:6)
	9	Character Line*80
	10	PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
	11	PARAMETER (MXPQX = (MXF 2 +2) MXQ * MXX)
	12	PARAMETER (M= 2, M1 = M + 1)
	13	Common
	14	> / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
	15	> / CtqPar2 / Nx, Nt, NfMx
	16	> / XQrange / Qini, Qmax, Xmin
	17	> / QCDtable / Alambda, Nfl, Iorder
	18	> / Masstbl / Amass(6)
	19	data pi / 3.141592653589793d0 /
	20	save
	21	c
	22	U = X * CtLhCtq5Pdf(1,X,Q)
	23	D = X * CtLhCtq5Pdf(2,X,Q)
	24	USEA = X * CtLhCtq5Pdf(-1,X,Q)
	25	DSEA = X * CtLhCtq5Pdf(-2,X,Q)
	26	STR = X * CtLhCtq5Pdf(3,X,Q)
	27	CHM = X * CtLhCtq5Pdf(4,X,Q)
	28	BOT = X * CtLhCtq5Pdf(5,X,Q)
	29	GLU = X * CtLhCtq5Pdf(0,X,Q)
	30	UPV=U-USEA
	31	DNV=D-DSEA
	32	c
	33	pdf(0) = glu
	34	pdf(1) = dnv+dsea
	35	pdf(-1) = dsea
	36	pdf(2) = upv+usea
	37	pdf(-2) = usea
	38	pdf(3) = str
	39	pdf(-3) = str
	40	pdf(4) = chm
	41	pdf(-4) = chm
	42	pdf(5) = bot
	43	pdf(-5) = bot
	44	pdf(6) = 0.0d0
	45	pdf(-6) = 0.0d0
	46
	47	return
	48	*
	49	entry CTEQ5read(nset)
	50
	51	call CtLhbldat1 !this line was missing in previous releases (jcp)
	52	call CtLhbldat2 !this line was missing in previous releases (jcp)
	53
	54	read(1,*)nmem(nset),ndef(nset)
	55	Read (1, '(A)') Line
	56	Read (1, '(A)') Line
	57	Read (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
	58	Iorder = Nint(Dr)
	59	Nfl = Nint(Fl)
	60	Alambda = Al
	61
	62	Read (1, '(A)') Line
	63	Read (1, *) NX, NT, NfMx
	64
65	Read (1, '(A)') Line
66	Read (1, *) QINI, QMAX, (QL(I), I =0, NT)
67
68	Read (1, '(A)') Line
69	Read (1, *) XMIN, (XV(I), I =0, NX)
70
71	Do 11 Iq = 0, NT
72	QL(Iq) = Log (QL(Iq) /Al)
73	11 Continue
74	C
75	C Since quark = anti-quark for nfl>2 at this stage,
76	C we Read out only the non-redundent data points
77	C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
78
79	Nblk = (NX+1) * (NT+1)
80	Npts = Nblk * (NfMx+3)
81
82	Read (1, '(A)') Line
83	Read (1, *, IOSTAT=IRET) (UPD(I), I=1,Npts)
84	Return
85	*
86
87	entry CTEQ5alfa(alfas,Qalfa)
88
89	alfas = pi*CtLhALPI(Qalfa)
90
91	return
92	*
93	entry CTEQ5init(Eorder,Q2fit)
94
95	return
96	*
97	entry CTEQ5pdf(mem)
98	imem = mem
99	return
100	*
101	end
102
103	FUNCTION CtLhPartonX5 (IPRTN, X, Q)
104	C
105	C Given the parton distribution function in the array UPD in
106	C COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
107	C x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
108	C
109	IMPLICIT DOUBLE PRECISION (A-H, O-Z)
110	C
111	PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
112	PARAMETER (MXPQX = (MXF 2 +2) MXQ * MXX)
113	PARAMETER (M= 2, M1 = M + 1)
114	C
115	Logical First
116	Common
117	> / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
118	> / CtqPar2 / Nx, Nt, NfMx
119	> / XQrange / Qini, Qmax, Xmin
120	C
121	Dimension Fq(M1), Df(M1)
122
123	Data First /.true./
124	save First
125	C Work with Log (Q)
126	QG = LOG (Q/AL)
127
128	C Find lower end of interval containing X
129	JL = -1
130	JU = Nx+1
131	11 If (JU-JL .GT. 1) Then
132	JM = (JU+JL) / 2
133	If (X .GT. XV(JM)) Then
134	JL = JM
135	Else
136	JU = JM
137	Endif
138	Goto 11
139	Endif
140
141	Jx = JL - (M-1)/2
142	If (X .lt. Xmin .and. First ) Then
143	First = .false.
144	c Print '(A, 2(1pE12.4))',
145	c > ' WARNING: X << Xmin, extrapolation used; X,Xmin=',X,Xmin
146	If (Jx .LT. 0) Jx = 0
147	Elseif (Jx .GT. Nx-M) Then
148	Jx = Nx - M
149	Endif
150	C Find the interval where Q lies
151	JL = -1
152	JU = NT+1
153	12 If (JU-JL .GT. 1) Then
154	JM = (JU+JL) / 2
155	If (QG .GT. QL(JM)) Then
156	JL = JM
157	Else
158	JU = JM
159	Endif
160	Goto 12
161	Endif
162
163	Jq = JL - (M-1)/2
164	If (Jq .LT. 0) Then
165	Jq = 0
166	c If (Q .lt. Qini) Print '(A, 2(1pE12.4))',
167	c > ' WARNING: Q << Qini, extrapolation used; Q,Qini=',Q,Qini
168	Elseif (Jq .GT. Nt-M) Then
169	Jq = Nt - M
170	c If (Q .gt. Qmax) Print '(A, 2(1pE12.4))',
171	c > ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
172	Endif
173
174	If (Iprtn .GE. 3) Then
175	Ip = - Iprtn
176	Else
177	Ip = Iprtn
178	EndIf
179	C Find the off-set in the linear array Upd
180	JFL = Ip + NfMx
181	J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx
182	C
183	C Now interpolate in x for M1 Q's
184	Do 21 Iq = 1, M1
185	J1 = J0 + (Nx+1)*(Iq-1) + 1
186	Call CtLhPolint3 (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
187	21 Continue
188	C Finish off by interpolating in Q
189	Call CtLhPolint3 (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)
190
191	CtLhPartonX5 = Ftmp
192	C
193	RETURN
194	C ****************************
195	END
196
197	Function CtLhCtq5Pdf (Iparton, X, Q)
198	Implicit Double Precision (A-H,O-Z)
199	Logical Warn
200	Common
201	> / CtqPar2 / Nx, Nt, NfMx
202	> / QCDtable / Alambda, Nfl, Iorder
203
204	Data Warn /.true./
205	save Warn
206
207	If (X .lt. 0D0 .or. X .gt. 1D0) Then
208	Print *, 'X out of range in CtLhCtq5Pdf: ', X
209	Stop
210	Endif
211	If (Q .lt. Alambda) Then
212	Print *, 'Q out of range in CtLhCtq5Pdf: ', Q
213	Stop
214	Endif
215	If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
216	If (Warn) Then
217	C put a warning for calling extra flavor.
218	Warn = .false.
219	Print *, 'Warning: Iparton out of range in CtLhCtq5Pdf: '
220	> , Iparton
221	Endif
222	CtLhCtq5Pdf = 0D0
223	Return
224	Endif
225
226	CtLhCtq5Pdf = CtLhPartonX5 (Iparton, X, Q)
227	if(CtLhCtq5Pdf.lt.0.D0) CtLhCtq5Pdf = 0.D0
228
229	Return
230
231	C ********************
232	End