Removing obsolete constants.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / wrapQCDNUM4.f
1       subroutine QCDNUM4evolve(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
18 *
19 c      print *,'QCDNUM4evolve'
20       call getnset(iset)
21       if (iset.ne.isetlast) then
22         call get_pdfqcd(iset)
23         isetlast = iset
24       endif
25 c      
26       Q2=Q*Q
27       nf=4
28       if (Q2.lt.tc2) nf=3
29       pdf(0)=QPDFXQ('GLUON'  ,x,Q2,iflag)
30       singlet=QPDFXQ('SINGLET',x,Q2,iflag)
31       dm= QPDFXQ('DM',x,Q2,IFLAG)
32       um= QPDFXQ('UM',x,Q2,IFLAG)
33       dp= QPDFXQ('DP',x,Q2,IFLAG)
34       up= QPDFXQ('UP',x,Q2,IFLAG)
35       sp= QPDFXQ('SP',x,Q2,IFLAG)
36       ub=0.5d0*(up-um+singlet/dble(nf))
37       db=0.5d0*(dp-dm+singlet/dble(nf))
38       sb=0.5d0*(sp+singlet/dble(nf))
39       cb=0.d0
40       if (nf.ge.4) then
41          cp= QPDFXQ('CP',X,Q2,IFLAG)
42          cb=0.5d0*(cp+singlet/dble(nf))
43       end if
44
45       pdf(1)=dm+db
46       pdf(2)=um+ub
47       pdf(3)=sb
48       pdf(4)=cb
49       pdf(5)=0d0
50       pdf(6)=0d0
51       pdf(-1)=db
52       pdf(-2)=ub
53       pdf(-3)=sb
54       pdf(-4)=cb
55       pdf(-5)=0d0
56       pdf(-6)=0d0
57       return
58 *
59       entry QCDNUM4alfa(alfas,Q)
60       Q2=Q*Q
61       nf=4
62       if (Q2.lt.mc2) nf=3
63       alfas=QALFAS(Q2,Qlam,nf,iflag)
64       return
65 *
66       entry QCDNUM4read(nset)
67 c      read(1,*) s1
68 c      print *,s1
69 c         gridname='large.grid'
70 c         nx=400
71 c         xmin=1d-6
72 c         xmax=1d0
73 c         nq=112
74 c         qmin=1d0
75 c         qmax=1d10
76 c         read(1,*) dummy
77 *      else
78          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
79 c*      endif
80 c      iset = nset
81       return
82 *
83       entry QCDNUM4init(nset,Eorder,Q2fit)
84       call qninit
85       call QNISET('ORDER',Eorder+1)
86       call grxdef(nx,xmin)
87       call grqdef(nq,qmin,qmax)
88       call grqinp(Q2fit,1)
89       call getQMassM(nset,4,mc)
90       mc2=mc*mc
91 c      print *,mc
92       CALL QNRSET('CMASS',mc)
93       CALL QNRSET('MCALF',mc)
94       call getThresholdM(nset,4,tc)
95       tc2=tc*tc
96       call GRQINP(tc2,1)
97       CALL GRQINP(tc2-1.0d-4,1)
98 c      print *,'setting iq0',Q2fit
99       call qthres(tc2,2d10)
100       iq0=IQFROMQ(Q2fit)
101       iqc=IQFROMQ(tc2)
102       nf0= NFLGET(iq0)
103 c      print *,iq0,iqc,nf0,q2fit,tc2
104 c      nf0=3
105 c      print *,nf0
106       CALL QNBOOK(2,'dm')
107       CALL QNBOOK(3,'um')
108       CALL QNBOOK(4,'dp')
109       CALL QNBOOK(5,'up')
110       CALL QNBOOK(6,'sp')
111       CALL QNBOOK(7,'cp')
112       CALL QNLSET('W1ANA',.TRUE.)
113       CALL QNLSET('W2NUM',.TRUE.)
114       CALL QNLSET('W2STF',.FALSE.)
115       if (index(gridname,'none').eq.1) then
116          call qnfilw(0,0)
117       else
118          qnerr=-1
119          open(unit=2,status='old',file=gridname,
120      .        form='unformatted',err=1)
121          call QNREAD(2,1,qnerr)
122  1       close(2)
123          if (qnerr.ne.0) then
124             write(*,*) 'Grid file problem: ',gridname
125             if (qnerr.lt.0) then 
126                write(*,*) 'Grid file does not exist'
127                write(*,*) 'Calculating and creating grid file'
128                call qnfilw(0,0)
129                open(unit=2,status='unknown',file=gridname,
130      .              form='unformatted')
131                call QNDUMP(2)
132                close(2)
133             else
134                write(*,*) 'Existing grid file is inconsistent'
135                if (qnerr.eq.1)
136      .              write(*,*) 'Defined grid different'
137                if (qnerr.eq.2)
138      .              write(*,*) 'Heavy quark weight table different'
139                if (qnerr.eq.3)
140      .              write(*,*) 'Charm mass different'
141                if (qnerr.eq.4)
142      .              write(*,*) 'Bottom mass different'
143                stop
144             endif
145          endif
146       endif
147       return
148 *
149       entry QCDNUM4pdf(nset)
150 c      print *,'entering QCDNUMpdf',nset
151       call GetAlfas(nset,alfas0,scale0)
152 c      print *,alfas0,scale0
153       Q2=scale0*scale0
154       CALL QNRSET('ALFAS',alfas0)      
155       CALL QNRSET('ALFQ0',Q2)
156       DO ix = 1,nx
157          xx = XFROMIX(ix)
158 c        print *,'calling parmPDF',ix
159          call parmPDF(nset,xx,f)
160 c        if(ix.lt.6) print *,nset,xx,f
161          singlet=0.d0
162          do i=1,nf0
163             singlet=singlet+f(i)+f(-i)
164 c           print  *,i,singlet
165          end do
166          gl=f(0)
167          dm=f(1)-f(-1)
168          um=f(2)-f(-2)
169          dp=f(1)+f(-1)-singlet/dble(nf0)
170          up=f(2)+f(-2)-singlet/dble(nf0)
171          sp=f(3)+f(-3)-singlet/dble(nf0)
172 c         print *,ix,iq0
173          CALL QNPSET('SINGLET',ix,iq0,singlet)
174          CALL QNPSET('GLUON',ix,iq0,gl)
175          CALL QNPSET('DM',ix,iq0,DM)
176          CALL QNPSET('UM',ix,iq0,UM)
177          CALL QNPSET('DP',ix,iq0,DP)
178          CALL QNPSET('UP',ix,iq0,UP)
179          CALL QNPSET('SP',ix,iq0,SP)
180       ENDDO
181       CALL EVOLSG(iq0,1,nq)
182       CALL EVOLNM('DM',iq0,1,nq)
183       CALL EVOLNM('UM',iq0,1,nq)
184       CALL EVLSEA4('dp',iq0,iqc,nq)
185       CALL EVLSEA4('up',iq0,iqc,nq)
186       CALL EVLSEA4('sp',iq0,iqc,nq)
187
188 c      print *,'calling evols - heavy...'
189
190 C--   Heavy quark evolution
191
192       CALL EVOLCP4('cp',iqc,nq)
193 c
194       call getnset(iset)
195       call save_pdfqcd(iset)
196       return
197 *
198       end
199       subroutine EVLSEA4(name,IQ0,IQC,NQGRI)
200       implicit none
201       
202       CHARACTER*(*) name
203       integer iq0,iqc,nqgri
204       real*8 f34,f45,f43,f54,factor
205       parameter(f34=1.d0/12.d0,f43=-1.d0/12.d0)
206
207 c      print *,iq0,iqc,nqgri
208       If(IQ0.le.IQC)then
209          CAll EVPLUS(name,IQ0,1,IQC)
210          CALL QADDSI(name,IQC,f34)
211          CAll EVPLUS(name,IQC,IQC,NQGRI)
212       else if(IQ0.gt.IQC)then
213          CAll EVPLUS(name,IQ0,IQC,NQGRI)
214          CALL QADDSI(name,IQC,f43)
215          CAll EVPLUS(name,IQC,1,IQC)
216       end if         
217
218       end
219
220
221       subroutine EVOLCP4(name,IQC,NQGRI)
222       implicit none
223       
224       CHARACTER*(*) name
225       integer iq0,iqc,nqgri
226       real*8 f4
227       parameter(f4=-1.d0/4.d0)
228
229 c
230 c     First set to zero to avoid adding -1/4Singl at each iteration
231 c
232
233       CAll QNPNUL(name)
234       CALL QADDSI(name,IQC,f4)      
235       CAll EVPLUS(name,IQC,IQC,NQGRI)
236       
237       end