Reading QA thresholds from external file II
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapQCDNUM3.f
CommitLineData
4e9e3152 1 subroutine QCDNUM3evolve(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 call getnset(iset)
20 if (iset.ne.isetlast) then
21 call get_pdfqcd(iset)
22 isetlast = iset
23 endif
24c
25 Q2=Q*Q
26 nf=3
27 pdf(0)=QPDFXQ('GLUON' ,x,Q2,iflag)
28 singlet=QPDFXQ('SINGLET',x,Q2,iflag)
29 dm= QPDFXQ('DM',x,Q2,IFLAG)
30 um= QPDFXQ('UM',x,Q2,IFLAG)
31 dp= QPDFXQ('DP',x,Q2,IFLAG)
32 up= QPDFXQ('UP',x,Q2,IFLAG)
33 sp= QPDFXQ('SP',x,Q2,IFLAG)
34 ub=0.5d0*(up-um+singlet/dble(nf))
35 db=0.5d0*(dp-dm+singlet/dble(nf))
36 sb=0.5d0*(sp+singlet/dble(nf))
37
38 pdf(1)=dm+db
39 pdf(2)=um+ub
40 pdf(3)=sb
41 pdf(4)=0d0
42 pdf(5)=0d0
43 pdf(6)=0d0
44 pdf(-1)=db
45 pdf(-2)=ub
46 pdf(-3)=sb
47 pdf(-4)=0d0
48 pdf(-5)=0d0
49 pdf(-6)=0d0
50 return
51*
52 entry QCDNUM3alfa(alfas,Q)
53 Q2=Q*Q
54 nf=3
55 alfas=QALFAS(Q2,Qlam,nf,iflag)
56 return
57*
58 entry QCDNUM3read(nset)
59c read(1,*) s1
60c print *,s1
61c gridname='large.grid'
62c nx=400
63c xmin=1d-6
64c xmax=1d0
65c nq=112
66c qmin=1d0
67c qmax=1d10
68c read(1,*) dummy
69* else
70 read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
71c* endif
72c iset = nset
73 return
74*
75 entry QCDNUM3init(nset,Eorder,Q2fit)
76 call qninit
77 call QNISET('ORDER',Eorder+1)
78 call grxdef(nx,xmin)
79 call grqdef(nq,qmin,qmax)
80 call grqinp(Q2fit,1)
81 iq0=IQFROMQ(Q2fit)
82c nf0= NFLGET(iq0)
83 nf0=3
84 CALL QNBOOK(2,'dm')
85 CALL QNBOOK(3,'um')
86 CALL QNBOOK(4,'dp')
87 CALL QNBOOK(5,'up')
88 CALL QNBOOK(6,'sp')
89 CALL QNLSET('W1ANA',.TRUE.)
90 CALL QNLSET('W2NUM',.TRUE.)
91 CALL QNLSET('W2STF',.FALSE.)
92 call qthres(1d10,2d10)
93 if (index(gridname,'none').eq.1) then
94 call qnfilw(0,0)
95 else
96 qnerr=-1
97 open(unit=2,status='old',file=gridname,
98 . form='unformatted',err=1)
99 call QNREAD(2,1,qnerr)
100 1 close(2)
101 if (qnerr.ne.0) then
102 write(*,*) 'Grid file problem: ',gridname
103 if (qnerr.lt.0) then
104 write(*,*) 'Grid file does not exist'
105 write(*,*) 'Calculating and creating grid file'
106 call qnfilw(0,0)
107 open(unit=2,status='unknown',file=gridname,
108 . form='unformatted')
109 call QNDUMP(2)
110 close(2)
111 else
112 write(*,*) 'Existing grid file is inconsistent'
113 if (qnerr.eq.1)
114 . write(*,*) 'Defined grid different'
115 if (qnerr.eq.2)
116 . write(*,*) 'Heavy quark weight table different'
117 if (qnerr.eq.3)
118 . write(*,*) 'Charm mass different'
119 if (qnerr.eq.4)
120 . write(*,*) 'Bottom mass different'
121 stop
122 endif
123 endif
124 endif
125 return
126*
127 entry QCDNUM3pdf(nset)
128c print *,'entering QCDNUMpdf',nset
129 call GetAlfas(nset,alfas0,scale0)
130c print *,alfas0,scale0
131 Q2=scale0*scale0
132 CALL QNRSET('ALFAS',alfas0)
133 CALL QNRSET('ALFQ0',Q2)
134 DO ix = 1,nx
135 xx = XFROMIX(ix)
136c print *,'calling parmPDF',ix
137 call parmPDF(nset,xx,f)
138c if(ix.lt.6) print *,nset,xx,f
139 singlet=0.d0
140 do i=1,nf0
141 singlet=singlet+f(i)+f(-i)
142c print *,i,singlet
143 end do
144 gl=f(0)
145 dm=f(1)-f(-1)
146 um=f(2)-f(-2)
147 dp=f(1)+f(-1)-singlet/dble(nf0)
148 up=f(2)+f(-2)-singlet/dble(nf0)
149 sp=f(3)+f(-3)-singlet/dble(nf0)
150 CALL QNPSET('SINGLET',ix,iq0,singlet)
151 CALL QNPSET('GLUON',ix,iq0,gl)
152 CALL QNPSET('DM',ix,iq0,DM)
153 CALL QNPSET('UM',ix,iq0,UM)
154 CALL QNPSET('DP',ix,iq0,DP)
155 CALL QNPSET('UP',ix,iq0,UP)
156 CALL QNPSET('SP',ix,iq0,SP)
157 ENDDO
158 CALL EVOLSG(iq0,1,nq)
159 CALL EVOLNM('DM',iq0,1,nq)
160 CALL EVOLNM('UM',iq0,1,nq)
161c
162c CALL EVPLUS('DP',iq0,1,nq)
163c CALL EVPLUS('UP',iq0,1,nq)
164c CALL EVPLUS('SP',iq0,1,nq)
165 CALL EVOLNP('DP',iq0,1,nq)
166 CALL EVOLNP('UP',iq0,1,nq)
167 CALL EVOLNP('SP',iq0,1,nq)
168c
169 call getnset(iset)
170 call save_pdfqcd(iset)
171 return
172*
173 end
174