COVERITY reports many FORWARD_NULL, corrected; AliEMCALv2: initialization of fHits...
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapQCDNUM3.f
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
24 c      
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)
59 c      read(1,*) s1
60 c      print *,s1
61 c         gridname='large.grid'
62 c         nx=400
63 c         xmin=1d-6
64 c         xmax=1d0
65 c         nq=112
66 c         qmin=1d0
67 c         qmax=1d10
68 c         read(1,*) dummy
69 *      else
70          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
71 c*      endif
72 c      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)
82 c      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)
128 c      print *,'entering QCDNUMpdf',nset
129       call GetAlfas(nset,alfas0,scale0)
130 c      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)
136 c        print *,'calling parmPDF',ix
137          call parmPDF(nset,xx,f)
138 c        if(ix.lt.6) print *,nset,xx,f
139          singlet=0.d0
140          do i=1,nf0
141             singlet=singlet+f(i)+f(-i)
142 c           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)
161 c      
162 c      CALL EVPLUS('DP',iq0,1,nq)
163 c      CALL EVPLUS('UP',iq0,1,nq)
164 c      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)
168 c
169       call getnset(iset)
170       call save_pdfqcd(iset)
171       return
172 *
173       end
174