Change needed for G4
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapQCDNUM4.f
CommitLineData
4e9e3152 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*
19c print *,'QCDNUM4evolve'
20 call getnset(iset)
21 if (iset.ne.isetlast) then
22 call get_pdfqcd(iset)
23 isetlast = iset
24 endif
25c
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)
67c read(1,*) s1
68c print *,s1
69c gridname='large.grid'
70c nx=400
71c xmin=1d-6
72c xmax=1d0
73c nq=112
74c qmin=1d0
75c qmax=1d10
76c read(1,*) dummy
77* else
78 read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
79c* endif
80c 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
91c 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)
98c print *,'setting iq0',Q2fit
99 call qthres(tc2,2d10)
100 iq0=IQFROMQ(Q2fit)
101 iqc=IQFROMQ(tc2)
102 nf0= NFLGET(iq0)
103c print *,iq0,iqc,nf0,q2fit,tc2
104c nf0=3
105c 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)
150c print *,'entering QCDNUMpdf',nset
151 call GetAlfas(nset,alfas0,scale0)
152c 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)
158c print *,'calling parmPDF',ix
159 call parmPDF(nset,xx,f)
160c if(ix.lt.6) print *,nset,xx,f
161 singlet=0.d0
162 do i=1,nf0
163 singlet=singlet+f(i)+f(-i)
164c 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)
172c 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
188c print *,'calling evols - heavy...'
189
190C-- Heavy quark evolution
191
192 CALL EVOLCP4('cp',iqc,nq)
193c
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
207c 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
229c
230c First set to zero to avoid adding -1/4Singl at each iteration
231c
232
233 CAll QNPNUL(name)
234 CALL QADDSI(name,IQC,f4)
235 CAll EVPLUS(name,IQC,IQC,NQGRI)
236
237 end