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