Reading QA thresholds from external file II
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapQCDNUM.f
1       subroutine QCDNUMevolve(x,Q,pdf)
2       implicit none
3       character*64 gridname
4       character*16 s1,dummy
5       integer Eorder,index
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
12       integer nx,nq
13       integer nset,iset,isetlast
14       data isetlast/-1/
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
18 *
19       call getnset(iset)
20 c fix 14.7.2006 to speed up sue of QCDNUM
21       if(iset.ne.isetlast) then
22         call get_pdfqcd(iset)
23         isetlast = iset
24       endif
25 c      
26       Q2=Q*Q
27       nf=5
28       if (Q2.lt.tb2) nf=4
29       if (Q2.lt.tc2) nf=3
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))
40       cb=0.d0
41       if (nf.ge.4) then
42          cp= QPDFXQ('CP',X,Q2,IFLAG)
43          cb=0.5d0*(cp+singlet/dble(nf))
44       end if
45       bb=0.d0
46       if (nf.ge.5) then
47          bp= QPDFXQ('BP',X,Q2,IFLAG)            
48          bb=0.5d0*(bp+singlet/dble(nf))
49       end if
50
51       pdf(1)=dm+db
52       pdf(2)=um+ub
53       pdf(3)=sb
54       pdf(4)=cb
55       pdf(5)=bb
56       pdf(6)=0d0
57       pdf(-1)=db
58       pdf(-2)=ub
59       pdf(-3)=sb
60       pdf(-4)=cb
61       pdf(-5)=bb
62       pdf(-6)=0d0
63       return
64 *
65       entry QCDNUMalfa(alfas,Q)
66       Q2=Q*Q
67       nf=6
68       if (Q2.lt.mt2) nf=5
69       if (Q2.lt.mb2) nf=4
70       if (Q2.lt.mc2) nf=3
71       alfas=QALFAS(Q2,Qlam,nf,iflag)
72       return
73 *
74       entry QCDNUMread(nset)
75 c      read(1,*) s1
76 c      print *,s1
77 c         gridname='large.grid'
78 c         nx=400
79 c         xmin=1d-6
80 c         xmax=1d0
81 c         nq=112
82 c         qmin=1d0
83 c         qmax=1d10
84 c         read(1,*) dummy
85 *      else
86          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
87 c*      endif
88 c      iset = nset
89       return
90 *
91       entry QCDNUMinit(nset,Eorder,Q2fit)
92       call qninit
93       call QNISET('ORDER',Eorder+1)
94       call grxdef(nx,xmin)
95       call grqdef(nq,qmin,qmax)
96       call grqinp(Q2fit,1)
97       call getQMassM(nset,4,mc)
98       mc2=mc*mc
99       CALL QNRSET('CMASS',mc)
100       CALL QNRSET('MCALF',mc)
101       call getThresholdM(nset,4,tc)
102       tc2=tc*tc
103       call GRQINP(tc2,1)
104       CALL GRQINP(tc2-1.0d-4,1)
105       call getQMassM(nset,5,mb)
106       mb2=mb*mb
107       CALL QNRSET('BMASS',mb)
108       CALL QNRSET('MBALF',mb)
109       call getThresholdM(nset,5,tb)
110       tb2=tb*tb
111       call GRQINP(tb2,1)
112       CALL GRQINP(tb2-1.0d-4,1)
113       call getQMassM(nset,6,mt)
114       mt2=mt*mt
115       CALL QNRSET('TMASS',mt)
116       CALL QNRSET('MTALF',mt)
117       call qthres(tc2,tb2)
118       iq0=IQFROMQ(Q2fit)
119       iqc=IQFROMQ(tc2)
120       iqb=IQFROMQ(tb2)
121       nf0= NFLGET(iq0)
122 c      print *,nf0
123       CALL QNBOOK(2,'dm')
124       CALL QNBOOK(3,'um')
125       CALL QNBOOK(4,'dp')
126       CALL QNBOOK(5,'up')
127       CALL QNBOOK(6,'sp')
128       CALL QNBOOK(7,'cp')
129       CALL QNBOOK(8,'bp')
130       CALL QNLSET('W1ANA',.TRUE.)
131       CALL QNLSET('W2NUM',.TRUE.)
132       CALL QNLSET('W2STF',.FALSE.)
133       if (index(gridname,'none').eq.1) then
134          call qnfilw(0,0)
135       else
136          qnerr=-1
137          open(unit=2,status='old',file=gridname,
138      .        form='unformatted',err=1)
139          call QNREAD(2,1,qnerr)
140  1       close(2)
141          if (qnerr.ne.0) then
142             write(*,*) 'Grid file problem: ',gridname
143             if (qnerr.lt.0) then 
144                write(*,*) 'Grid file does not exist'
145                write(*,*) 'Calculating and creating grid file'
146                call qnfilw(0,0)
147                open(unit=2,status='unknown',file=gridname,
148      .              form='unformatted')
149                call QNDUMP(2)
150                close(2)
151             else
152                write(*,*) 'Existing grid file is inconsistent'
153                if (qnerr.eq.1)
154      .              write(*,*) 'Defined grid different'
155                if (qnerr.eq.2)
156      .              write(*,*) 'Heavy quark weight table different'
157                if (qnerr.eq.3)
158      .              write(*,*) 'Charm mass different'
159                if (qnerr.eq.4)
160      .              write(*,*) 'Bottom mass different'
161                stop
162             endif
163          endif
164       endif
165       return
166 *
167       entry QCDNUMpdf(nset)
168 c      print *,'entering QCDNUMpdf',nset
169       call GetAlfas(nset,alfas0,scale0)
170 c      print *,alfas0,scale0
171       Q2=scale0*scale0
172       CALL QNRSET('ALFAS',alfas0)      
173       CALL QNRSET('ALFQ0',Q2)
174       DO ix = 1,nx
175          xx = XFROMIX(ix)
176 c        print *,'calling parmPDF',ix
177          call parmPDF(nset,xx,f)
178 c        if(ix.lt.6) print *,nset,xx,f
179          singlet=0.d0
180          do i=1,nf0
181             singlet=singlet+f(i)+f(-i)
182 c           print  *,i,singlet
183          end do
184          gl=f(0)
185          dm=f(1)-f(-1)
186          um=f(2)-f(-2)
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)
204       ENDDO
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)
217
218 c      print *,'calling evols - heavy...'
219
220 C--   Heavy quark evolution
221
222       CALL EVOLCP('cp',iqc,iqb,nq)
223       CALL EVOLBP('bp',iqb,nq)
224 c
225       call getnset(iset)
226       call save_pdfqcd(iset)
227       return
228 *
229       end
230       subroutine EVLSEA(name,IQ0,IQC,IQB,NQGRI)
231       implicit none
232       
233       CHARACTER*(*) name
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,
237      >          f54=-1.d0/20.d0)
238
239
240       If(IQ0.le.IQC)then
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)
261       end if         
262
263       end
264
265
266       subroutine EVOLCP(name,IQC,IQB,NQGRI)
267       implicit none
268       
269       CHARACTER*(*) name
270       integer iq0,iqc,iqb,nqgri
271       real*8 f4,f45
272       parameter(f4=-1.d0/4.d0,f45=1.d0/20.d0)
273
274 c
275 c     First set to zero to avoid adding -1/4Singl at each iteration
276 c
277
278       CAll QNPNUL(name)
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)
283       
284       end
285
286       subroutine EVOLBP(name,IQB,NQGRI)
287       implicit none
288       
289       CHARACTER*(*) name
290       integer iqb,nqgri
291       real*8 f5
292       parameter(f5=-1.d0/5.d0)
293
294 c
295 c     First set to zero to avoid adding -1/5Singl at each iteration
296 c
297
298       CAll QNPNUL(name)
299       CALL QADDSI(name,IQB,f5)      
300       CAll EVPLUS(name,IQB,IQB,NQGRI)
301       
302       end
303 c
304 c routine to swap in/out the array PDFQCD(MXX,MQ2,0:10) in QCDNUM as required 
305 c by nset
306 c
307       subroutine get_pdfqcd(nset)
308 c
309       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
310 c      
311       include 'parmsetup.inc'     
312 c
313       PARAMETER ( MXX = 410 )
314       PARAMETER ( MQ2 =  120 )
315 C--   Do not set the following parameter to zero!
316       PARAMETER ( NDFMAX = 20)
317 c
318       COMMON/QCPASS/
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)
329 c
330       real*8  PDFSAV(MXX,MQ2,0:10,nmxset)
331       save PDFSAV
332 c      
333       do k=0,10
334         do j=1,mq2
335           do i=1,mxx
336              PDFQCD(i,j,k)=PDFSAV(i,j,k,nset) 
337           enddo     
338         enddo     
339       enddo 
340 c
341       return
342 c
343       entry save_pdfqcd(nset)
344 c
345       do k=0,10
346         do j=1,mq2
347           do i=1,mxx
348              PDFSAV(i,j,k,nset)=PDFQCD(i,j,k) 
349           enddo     
350         enddo     
351       enddo 
352 c
353       return
354       end
355     
356