]>
Commit | Line | Data |
---|---|---|
4e9e3152 | 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 |