[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / 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
      call CtLhbldat2
 
      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)		!error corrected (jcp)
c        imem = mem
        call getnset(iset)
        call setnmem(iset,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
      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)

c note, this routine doesn't have the features of not recalculating 
c the x or q point values if they have not changed since the last call.
c that makes it slower than cteq6, but since cteq5 is already obsolete,
c speed is not so important for it. 
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
c crude treatment if outside the defined interval.
      If (Jx .LT. 0) then
         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
4e9e3152	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
	52	call CtLhbldat2
	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) !error corrected (jcp)
98	c imem = mem
99	call getnset(iset)
100	call setnmem(iset,mem)
101	return
102	*
103	end
104
105	FUNCTION CtLhPartonX5 (IPRTN, X, Q)
106	C
107	C Given the parton distribution function in the array UPD in
108	C COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
109	C x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
110	C
111	IMPLICIT DOUBLE PRECISION (A-H, O-Z)
112	C
113	PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
114	PARAMETER (MXPQX = (MXF 2 +2) MXQ * MXX)
115	PARAMETER (M= 2, M1 = M + 1)
116	C
117	Common
118	> / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
119	> / CtqPar2 / Nx, Nt, NfMx
120	> / XQrange / Qini, Qmax, Xmin
121	C
122	Dimension Fq(M1), Df(M1)
123
124	c note, this routine doesn't have the features of not recalculating
125	c the x or q point values if they have not changed since the last call.
126	c that makes it slower than cteq6, but since cteq5 is already obsolete,
127	c speed is not so important for it.
128	C Work with Log (Q)
129	QG = LOG (Q/AL)
130
131	C Find lower end of interval containing X
132	JL = -1
133	JU = Nx+1
134	11 If (JU-JL .GT. 1) Then
135	JM = (JU+JL) / 2
136	If (X .GT. XV(JM)) Then
137	JL = JM
138	Else
139	JU = JM
140	Endif
141	Goto 11
142	Endif
143
144	Jx = JL - (M-1)/2
145	c crude treatment if outside the defined interval.
146	If (Jx .LT. 0) then
147	Jx = 0
148	Elseif (Jx .GT. Nx-M) Then
149	Jx = Nx - M
150	Endif
151	C Find the interval where Q lies
152	JL = -1
153	JU = NT+1
154	12 If (JU-JL .GT. 1) Then
155	JM = (JU+JL) / 2
156	If (QG .GT. QL(JM)) Then
157	JL = JM
158	Else
159	JU = JM
160	Endif
161	Goto 12
162	Endif
163
164	Jq = JL - (M-1)/2
165	If (Jq .LT. 0) Then
166	Jq = 0
167	c If (Q .lt. Qini) Print '(A, 2(1pE12.4))',
168	c > ' WARNING: Q << Qini, extrapolation used; Q,Qini=',Q,Qini
169	Elseif (Jq .GT. Nt-M) Then
170	Jq = Nt - M
171	c If (Q .gt. Qmax) Print '(A, 2(1pE12.4))',
172	c > ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
173	Endif
174
175	If (Iprtn .GE. 3) Then
176	Ip = - Iprtn
177	Else
178	Ip = Iprtn
179	EndIf
180	C Find the off-set in the linear array Upd
181	JFL = Ip + NfMx
182	J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx
183	C
184	C Now interpolate in x for M1 Q's
185	Do 21 Iq = 1, M1
186	J1 = J0 + (Nx+1)*(Iq-1) + 1
187	Call CtLhPolint3 (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
188	21 Continue
189	C Finish off by interpolating in Q
190	Call CtLhPolint3 (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)
191
192	CtLhPartonX5 = Ftmp
193	C
194	RETURN
195
196	C ****************************
197	END
198
199	Function CtLhCtq5Pdf (Iparton, X, Q)
200	Implicit Double Precision (A-H,O-Z)
201	Logical Warn
202	Common
203	> / CtqPar2 / Nx, Nt, NfMx
204	> / QCDtable / Alambda, Nfl, Iorder
205
206	Data Warn /.true./
207	save Warn
208
209	If (X .lt. 0D0 .or. X .gt. 1D0) Then
210	Print *, 'X out of range in CtLhCtq5Pdf: ', X
211	Stop
212	Endif
213	If (Q .lt. Alambda) Then
214	Print *, 'Q out of range in CtLhCtq5Pdf: ', Q
215	Stop
216	Endif
217	If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
218	If (Warn) Then
219	C put a warning for calling extra flavor.
220	Warn = .false.
221	Print *, 'Warning: Iparton out of range in CtLhCtq5Pdf: '
222	> , Iparton
223	Endif
224	CtLhCtq5Pdf = 0D0
225	Return
226	Endif
227
228	CtLhCtq5Pdf = CtLhPartonX5 (Iparton, X, Q)
229	if(CtLhCtq5Pdf.lt.0.D0) CtLhCtq5Pdf = 0.D0
230
231	Return
232
233	C ********************
234	End