1 subroutine QCDNUMevolve(x,Q,pdf)
6 real*8 x,Q,Qlam,Q2,Q2fit,alfas0,scale0,alfas
7 real*8 mc,mc2,mb,mb2,mt,mt2,tc,tc2,tb,tb2
8 real*8 singlet,dm,um,dp,up,sp,ub,db,sb,cb,bb,bp,cp,gl,xx
9 real*8 QALFAS,QPDFXQ,XFROMIX
10 real*8 f(-6:6),pdf(-6:6)
11 integer iq0,iqc,iqb,nf,nf0,qnerr,NFLGET,iflag,ix,i,IQFROMQ
13 integer nset,iset,isetlast
15 real*8 xmin,xmax,qmin,qmax,S
16 save iq0,iqc,iqb,nf0,mc2,mb2,tc2,tb2,mt2
17 save nx,xmin,xmax,nq,qmin,qmax,gridname,isetlast
20 c fix 14.7.2006 to speed up sue of QCDNUM
21 if(iset.ne.isetlast) then
30 pdf(0)=QPDFXQ('GLUON' ,x,Q2,iflag)
31 singlet=QPDFXQ('SINGLET',x,Q2,iflag)
32 dm= QPDFXQ('DM',x,Q2,IFLAG)
33 um= QPDFXQ('UM',x,Q2,IFLAG)
34 dp= QPDFXQ('DP',x,Q2,IFLAG)
35 up= QPDFXQ('UP',x,Q2,IFLAG)
36 sp= QPDFXQ('SP',x,Q2,IFLAG)
37 ub=0.5d0*(up-um+singlet/dble(nf))
38 db=0.5d0*(dp-dm+singlet/dble(nf))
39 sb=0.5d0*(sp+singlet/dble(nf))
42 cp= QPDFXQ('CP',X,Q2,IFLAG)
43 cb=0.5d0*(cp+singlet/dble(nf))
47 bp= QPDFXQ('BP',X,Q2,IFLAG)
48 bb=0.5d0*(bp+singlet/dble(nf))
65 entry QCDNUMalfa(alfas,Q)
71 alfas=QALFAS(Q2,Qlam,nf,iflag)
74 entry QCDNUMread(nset)
77 c gridname='large.grid'
86 read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
91 entry QCDNUMinit(nset,Eorder,Q2fit)
93 call QNISET('ORDER',Eorder+1)
95 call grqdef(nq,qmin,qmax)
97 call getQMassM(nset,4,mc)
99 CALL QNRSET('CMASS',mc)
100 CALL QNRSET('MCALF',mc)
101 call getThresholdM(nset,4,tc)
104 CALL GRQINP(tc2-1.0d-4,1)
105 call getQMassM(nset,5,mb)
107 CALL QNRSET('BMASS',mb)
108 CALL QNRSET('MBALF',mb)
109 call getThresholdM(nset,5,tb)
112 CALL GRQINP(tb2-1.0d-4,1)
113 call getQMassM(nset,6,mt)
115 CALL QNRSET('TMASS',mt)
116 CALL QNRSET('MTALF',mt)
130 CALL QNLSET('W1ANA',.TRUE.)
131 CALL QNLSET('W2NUM',.TRUE.)
132 CALL QNLSET('W2STF',.FALSE.)
133 if (index(gridname,'none').eq.1) then
137 open(unit=2,status='old',file=gridname,
138 . form='unformatted',err=1)
139 call QNREAD(2,1,qnerr)
142 write(*,*) 'Grid file problem: ',gridname
144 write(*,*) 'Grid file does not exist'
145 write(*,*) 'Calculating and creating grid file'
147 open(unit=2,status='unknown',file=gridname,
148 . form='unformatted')
152 write(*,*) 'Existing grid file is inconsistent'
154 . write(*,*) 'Defined grid different'
156 . write(*,*) 'Heavy quark weight table different'
158 . write(*,*) 'Charm mass different'
160 . write(*,*) 'Bottom mass different'
167 entry QCDNUMpdf(nset)
168 c print *,'entering QCDNUMpdf',nset
169 call GetAlfas(nset,alfas0,scale0)
170 c print *,alfas0,scale0
172 CALL QNRSET('ALFAS',alfas0)
173 CALL QNRSET('ALFQ0',Q2)
176 c print *,'calling parmPDF',ix
177 call parmPDF(nset,xx,f)
178 c if(ix.lt.6) print *,nset,xx,f
181 singlet=singlet+f(i)+f(-i)
187 dp=f(1)+f(-1)-singlet/dble(nf0)
188 up=f(2)+f(-2)-singlet/dble(nf0)
189 sp=f(3)+f(-3)-singlet/dble(nf0)
190 c print *,'calling qpnset sg...',ix,iq0,singlet
191 CALL QNPSET('SINGLET',ix,iq0,singlet)
192 c print *,'calling qpnset gl...',ix,iq0,gl
193 CALL QNPSET('GLUON',ix,iq0,gl)
194 c print *,'calling qpnset dm...',ix,iq0,dm
195 CALL QNPSET('DM',ix,iq0,DM)
196 c print *,'calling qpnset um...',ix,iq0,um
197 CALL QNPSET('UM',ix,iq0,UM)
198 c print *,'calling qpnset dp...',ix,iq0,dp
199 CALL QNPSET('DP',ix,iq0,DP)
200 c print *,'calling qpnset up...',ix,iq0,up
201 CALL QNPSET('UP',ix,iq0,UP)
202 c print *,'calling qpnset sp...',ix,iq0,sp
203 CALL QNPSET('SP',ix,iq0,SP)
205 c print *,'calling evols...'
206 CALL EVOLSG(iq0,1,nq)
207 c print *,'calling evols dm...'
208 CALL EVOLNM('DM',iq0,1,nq)
209 c print *,'calling evols um...'
210 CALL EVOLNM('UM',iq0,1,nq)
211 c print *,'calling evols dp...',iq0,iqc,iqb,nq
212 CALL EVLSEA('dp',iq0,iqc,iqb,nq)
213 c print *,'calling evols up...'
214 CALL EVLSEA('up',iq0,iqc,iqb,nq)
215 c print *,'calling evols sp...'
216 CALL EVLSEA('sp',iq0,iqc,iqb,nq)
218 c print *,'calling evols - heavy...'
220 C-- Heavy quark evolution
222 CALL EVOLCP('cp',iqc,iqb,nq)
223 CALL EVOLBP('bp',iqb,nq)
226 call save_pdfqcd(iset)
230 subroutine EVLSEA(name,IQ0,IQC,IQB,NQGRI)
234 integer iq0,iqc,iqb,nqgri
235 real*8 f34,f45,f43,f54,factor
236 parameter(f34=1.d0/12.d0,f45=1.d0/20.d0,f43=-1.d0/12.d0,
241 CAll EVPLUS(name,IQ0,1,IQC)
242 CALL QADDSI(name,IQC,f34)
243 CAll EVPLUS(name,IQC,IQC,IQB)
244 CALL QADDSI(name,IQB,f45)
245 CAll EVPLUS(name,IQB,IQB,NQGRI)
246 else if(IQ0.gt.IQC.and.IQ0.le.IQB)then
247 CAll EVPLUS(name,IQ0,IQC,IQB)
248 CALL QADDSI(name,IQC,f43)
249 CAll EVPLUS(name,IQC,1,IQC)
250 CALL QADDSI(name,IQC,f34)
251 CALL QADDSI(name,IQB,f45)
252 CAll EVPLUS(name,IQB,IQB,NQGRI)
253 else if(IQ0.gt.IQB)then
254 CAll EVPLUS(name,IQ0,IQB,NQGRI)
255 CALL QADDSI(name,IQB,f54)
256 CAll EVPLUS(name,IQB,IQC,IQB)
257 CALL QADDSI(name,IQB,f45)
258 CALL QADDSI(name,IQC,f43)
259 CAll EVPLUS(name,IQC,1,IQC)
260 CALL QADDSI(name,IQC,f34)
266 subroutine EVOLCP(name,IQC,IQB,NQGRI)
270 integer iq0,iqc,iqb,nqgri
272 parameter(f4=-1.d0/4.d0,f45=1.d0/20.d0)
275 c First set to zero to avoid adding -1/4Singl at each iteration
279 CALL QADDSI(name,IQC,f4)
280 CAll EVPLUS(name,IQC,IQC,IQB)
281 CALL QADDSI(name,IQB,f45)
282 CAll EVPLUS(name,IQB,IQB,NQGRI)
286 subroutine EVOLBP(name,IQB,NQGRI)
292 parameter(f5=-1.d0/5.d0)
295 c First set to zero to avoid adding -1/5Singl at each iteration
299 CALL QADDSI(name,IQB,f5)
300 CAll EVPLUS(name,IQB,IQB,NQGRI)
304 c routine to swap in/out the array PDFQCD(MXX,MQ2,0:10) in QCDNUM as required
307 subroutine get_pdfqcd(nset)
309 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
311 include 'parmsetup.inc'
313 PARAMETER ( MXX = 410 )
314 PARAMETER ( MQ2 = 120 )
315 C-- Do not set the following parameter to zero!
316 PARAMETER ( NDFMAX = 20)
319 +ALPHA0, Q0ALFA, ASLAST, QALAST,
320 +ALFASQ(MQ2), ALFAPQ(MQ2), ALFA2Q(MQ2),
321 +DELUP(MQ2), DELDN(MQ2), PDFQCD(MXX,MQ2,0:10),
322 +FNSQCD(MXX,MQ2),DNSQCD(MXX,MQ2),
323 +FSIQCD(MXX,MQ2),DSIQCD(MXX,MQ2),
324 +FGLQCD(MXX,MQ2),DGGQCD(MXX,MQ2),
325 +FSTORE(MXX,MQ2,31:30+NDFMAX),IDFAST(7,30),NDFAST,
326 +MARKFF(MXX,MQ2),MARKFH(MXX,MQ2),MARKQQ(MQ2),
327 +ISTFID(31:30+NDFMAX),IPDFID(31:30+NDFMAX),IEALFA(MQ2),
328 +IQL_LAST(10),IQ0_LAST(10),IQH_LAST(10)
330 real*8 PDFSAV(MXX,MQ2,0:10,nmxset)
336 PDFQCD(i,j,k)=PDFSAV(i,j,k,nset)
343 entry save_pdfqcd(nset)
348 PDFSAV(i,j,k,nset)=PDFQCD(i,j,k)