Reading QA thresholds from external file II
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapgrv.f
1       subroutine GRVevolve(xin,qin,pdf)
2       implicit real*8 (a-h,o-z)
3       include 'parmsetup.inc'
4       PARAMETER(ngrid=2)
5       PARAMETER (NPART=6, NX=68, NQ=27, NARG=2)
6       character*16 name(nmxset)
7       integer nmem(nmxset),ndef(nmxset),mmem
8       common/NAME/name,nmem,ndef,mmem
9       DIMENSION XXUVF(0:ngrid,NX,NQ), XXDVF(0:ngrid,NX,NQ), 
10      +          XXDEF(0:ngrid,NX,NQ), XXUDF(0:ngrid,NX,NQ),
11      1          XXSF(0:ngrid,NX,NQ),   XXGF(0:ngrid,NX,NQ),
12      +          XUVF(NX,NQ), XDVF(NX,NQ), 
13      +          XDEF(NX,NQ), XUDF(NX,NQ),
14      1          XSF(NX,NQ),   XGF(NX,NQ),
15      +          PARTON (NPART,NQ,NX-1), 
16      2          QS(NQ), XB(NX), XT(NARG), NA(NARG), ARRF(NX+NQ) 
17       CHARACTER*80 LINE
18       dimension pdf(-6:6)
19       save 
20       x=xin
21       q2=qin*qin      
22 *...CHECK OF X AND Q2 VALUES : 
23       IF ( (X.LT.0.99D-9) .OR. (X.GT.1.D0) ) THEN
24          WRITE(6,91) 
25   91     FORMAT (2X,'PARTON INTERPOLATION: X OUT OF RANGE')
26          STOP
27       ENDIF
28       IF ( (Q2.LT.0.799) .OR. (Q2.GT.1.01E6) ) THEN
29          WRITE(6,92) 
30   92     FORMAT (2X,'PARTON INTERPOLATION: Q2 OUT OF RANGE')
31          STOP
32       ENDIF
33 *...INTERPOLATION :
34       DO IQ=1,NQ
35       DO IX=1,NX
36       xuvf(ix,iq)=xxuvf(imem,ix,iq)
37       xdvf(ix,iq)=xxdvf(imem,ix,iq)
38       xdef(ix,iq)=xxdef(imem,ix,iq)
39       xudf(ix,iq)=xxudf(imem,ix,iq)
40       xsf(ix,iq)=xxsf(imem,ix,iq)
41       xgf(ix,iq)=xxgf(imem,ix,iq)
42       enddo
43       enddo
44       XT(1) = DLOG(X)
45       XT(2) = DLOG(Q2)
46       X1 = 1.- X
47       XV = X**0.5
48       XS = X**(-0.2)
49       UV = FINT(NARG,XT,NA,ARRF,XUVF) * X1**3 * XV
50       DV = FINT(NARG,XT,NA,ARRF,XDVF) * X1**4 * XV
51       DE = FINT(NARG,XT,NA,ARRF,XDEF) * X1**7 * XV
52       UD = FINT(NARG,XT,NA,ARRF,XUDF) * X1**7 * XS
53       US = 0.5 * (UD - DE)
54       DS = 0.5 * (UD + DE)
55       SS = FINT(NARG,XT,NA,ARRF,XSF)  * X1**7 * XS
56       GL = FINT(NARG,XT,NA,ARRF,XGF)  * X1**5 * XS 
57 *
58
59       pdf(-6) = 0.0d0
60        pdf(6) = 0.0d0
61       pdf(-5) = 0.0d0
62        pdf(5) = 0.0d0
63       pdf(-4) = 0.0d0
64        pdf(4) = 0.0d0
65       pdf(-3) = ss
66        pdf(3) = ss
67       pdf(-2) = us
68        pdf(2) = uv+us
69       pdf(-1) = ds
70        pdf(1) = dv+ds
71        pdf(0) = gl
72       return
73 c      
74 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
75       entry GRVread(nset)
76 c      
77 c      print *,'calling grvread'
78       read(1,*)nmem(nset),ndef(nset)
79
80       read(1,93)xb
81 c      print *,xb
82       read(1,93)qs
83 c      print *,qs
84   93  format(8e8.2)
85 c
86       do ng=0,nmem(nset)
87 c      
88       READ(1,89) LINE
89   89  FORMAT(A80)
90 c      print *,line
91       DO 15 M = 1, NX-1 
92       DO 15 N = 1, NQ
93       READ(1,90) PARTON(1,N,M), PARTON(2,N,M), PARTON(3,N,M), 
94      1           PARTON(4,N,M), PARTON(5,N,M), PARTON(6,N,M) 
95   90  FORMAT (6(1PE10.3))
96   15  CONTINUE     
97 *
98 *....ARRAYS FOR THE INTERPOLATION SUBROUTINE :
99       DO 10 IQ = 1, NQ
100       DO 20 IX = 1, NX-1
101         XB0V = XB(IX)**0.5 
102         XB0S = XB(IX)**(-0.2) 
103         XB1 = 1.-XB(IX)
104         xXUVF(ng,IX,IQ) = PARTON(1,IQ,IX) / (XB1**3 * XB0V)
105         xXDVF(ng,IX,IQ) = PARTON(2,IQ,IX) / (XB1**4 * XB0V)
106         xXDEF(ng,IX,IQ) = PARTON(3,IQ,IX) / (XB1**7 * XB0V) 
107         xXUDF(ng,IX,IQ) = PARTON(4,IQ,IX) / (XB1**7 * XB0S)
108         xXSF(ng,IX,IQ)  = PARTON(5,IQ,IX) / (XB1**7 * XB0S)
109         xXGF(ng,IX,IQ)  = PARTON(6,IQ,IX) / (XB1**5 * XB0S)
110   20  CONTINUE
111         xXUVF(ng,NX,IQ) = 0.E0
112         xXDVF(ng,NX,IQ) = 0.E0
113         xXDEF(ng,NX,IQ) = 0.E0
114         xXUDF(ng,NX,IQ) = 0.E0
115         xXSF(ng,NX,IQ)  = 0.E0
116         xXGF(ng,NX,IQ)  = 0.E0
117   10  CONTINUE  
118       NA(1) = NX
119       NA(2) = NQ
120       DO 30 IX = 1, NX
121         ARRF(IX) = DLOG(XB(IX))
122   30  CONTINUE
123       DO 40 IQ = 1, NQ
124         ARRF(NX+IQ) = DLOG(QS(IQ))
125   40  CONTINUE
126       
127       enddo
128
129       return
130 c
131 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
132       entry GRValfa(alfas,qalfa)
133       call getnset(iset)
134       q2alfa = qalfa*qalfa
135       call GetOrderAsM(iset,iord)
136       nord=iord+1
137       alfas=grvals(q2alfa,nord)
138       alfas = 4.0d0*3.14159d0*alfas
139       return
140 c
141 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
142       entry GRVinit(Eorder,Q2fit)
143       return
144 c
145 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
146       entry GRVpdf(mem)
147       imem = mem
148       return
149 c
150  1000 format(5e13.5)
151       end
152 c
153       FUNCTION FINT(NARG,ARG,NENT,ENT,TABLE)
154 *********************************************************************
155 *                                                                   *
156 *   THE INTERPOLATION ROUTINE (CERN LIBRARY ROUTINE E104)           *
157 *                                                                   *
158 *********************************************************************
159       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
160       DIMENSION ARG(5),NENT(5),ENT(10),TABLE(10)
161       DIMENSION D(5),NCOMB(5),IENT(5)
162       KD=1
163       M=1
164       JA=1
165          DO 5 I=1,NARG
166       NCOMB(I)=1
167       JB=JA-1+NENT(I)
168          DO 2 J=JA,JB
169       IF (ARG(I).LE.ENT(J)) GO TO 3
170     2 CONTINUE
171       J=JB
172     3 IF (J.NE.JA) GO TO 4
173       J=J+1
174     4 JR=J-1
175       D(I)=(ENT(J)-ARG(I))/(ENT(J)-ENT(JR))
176       IENT(I)=J-JA
177       KD=KD+IENT(I)*M
178       M=M*NENT(I)
179     5 JA=JB+1
180       FINT=0.
181    10 FAC=1.
182       IADR=KD
183       IFADR=1
184          DO 15 I=1,NARG
185       IF (NCOMB(I).EQ.0) GO TO 12
186       FAC=FAC*(1.-D(I))
187       GO TO 15
188    12 FAC=FAC*D(I)
189       IADR=IADR-IFADR
190    15 IFADR=IFADR*NENT(I)
191       FINT=FINT+FAC*TABLE(IADR)
192       IL=NARG
193    40 IF (NCOMB(IL).EQ.0) GO TO 80
194       NCOMB(IL)=0
195       IF (IL.EQ.NARG) GO TO 10
196       IL=IL+1
197          DO 50  K=IL,NARG
198    50 NCOMB(K)=1
199       GO TO 10
200    80 IL=IL-1
201       IF(IL.NE.0) GO TO 40
202       RETURN
203       END
204            
205 c      FUNCTION ALPHAS (Q2, NAORD)
206       FUNCTION grvals (Q2, NAORD)
207 *********************************************************************
208 *                                                                   *
209 *   THE ALPHA_S ROUTINE.                                            *
210 *                                                                   *
211 *   INPUT :  Q2    =  scale in GeV**2  (not too low, of course);    *
212 *            NAORD =  1 (LO),  2 (NLO).                             *
213 *                                                                   *
214 *   OUTPUT:  alphas_s/(4 pi) for use with the GRV(98) partons.      *  
215 *                                                                   *
216 *******************************************************i*************
217 *
218       IMPLICIT DOUBLE PRECISION (A - Z)
219       INTEGER NF, K, I, NAORD
220       DIMENSION LAMBDAL (3:6),  LAMBDAN (3:6), Q2THR (3)
221 *
222 *...HEAVY QUARK THRESHOLDS AND LAMBDA VALUES :
223       DATA Q2THR   /  1.960,  20.25,  30625. /
224       DATA LAMBDAL / 0.2041, 0.1750, 0.1320, 0.0665 /
225       DATA LAMBDAN / 0.2994, 0.2460, 0.1677, 0.0678 /
226 *
227 *...DETERMINATION OF THE APPROPRIATE NUMBER OF FLAVOURS :
228       NF = 3
229       DO 10 K = 1, 3
230       IF (Q2 .GT. Q2THR (K)) THEN
231          NF = NF + 1
232       ELSE
233           GO TO 20
234        END IF
235   10   CONTINUE
236 *
237 *...LO ALPHA_S AND BETA FUNCTION FOR NLO CALCULATION :
238   20   B0 = 11.- 2./3.* NF
239        B1 = 102.- 38./3.* NF
240        B10 = B1 / (B0*B0)
241        IF (NAORD .EQ. 1) THEN
242          LAM2 = LAMBDAL (NF) * LAMBDAL (NF)
243          ALP  = 1./(B0 * DLOG (Q2/LAM2))
244          GO TO 1
245        ELSE IF (NAORD .EQ. 2) then
246          LAM2 = LAMBDAN (NF) * LAMBDAN (NF)
247          B1 = 102.- 38./3.* NF
248          B10 = B1 / (B0*B0)
249        ELSE
250          WRITE (6,91)
251   91     FORMAT ('INVALID CHOICE FOR ORDER IN ALPHA_S')
252          STOP
253        END IF
254 *
255 *...START VALUE FOR NLO ITERATION :
256        LQ2 = DLOG (Q2 / LAM2)
257        ALP = 1./(B0*LQ2) * (1.- B10*DLOG(LQ2)/LQ2)
258 *
259 *...EXACT NLO VALUE, FOUND VIA NEWTON PROCEDURE :
260        DO 2 I = 1, 6
261        XL  = DLOG (1./(B0*ALP) + B10)
262        XLP = DLOG (1./(B0*ALP*1.01) + B10)
263        XLM = DLOG (1./(B0*ALP*0.99) + B10)
264        Y  = LQ2 - 1./ (B0*ALP) + B10 * XL
265        Y1 = (- 1./ (B0*ALP*1.01) + B10 * XLP
266      1       + 1./ (B0*ALP*0.99) - B10 * XLP) / (0.02D0*ALP)
267        ALP = ALP - Y/Y1
268   2    CONTINUE
269 *
270 *...OUTPUT :
271 c  1    ALPHAS = ALP
272   1    grvals = ALP
273        RETURN
274        END
275
276
277