[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapQCDNUM4.f

      subroutine QCDNUM4evolve(x,Q,pdf)
      implicit none
      character*64 gridname
      character*16 s1,dummy
      integer Eorder,index
      real*8 x,Q,Qlam,Q2,Q2fit,alfas0,scale0,alfas
      real*8 mc,mc2,mb,mb2,mt,mt2,tc,tc2,tb,tb2
      real*8 singlet,dm,um,dp,up,sp,ub,db,sb,cb,bb,bp,cp,gl,xx
      real*8 QALFAS,QPDFXQ,XFROMIX
      real*8 f(-6:6),pdf(-6:6)
      integer iq0,iqc,iqb,nf,nf0,qnerr,NFLGET,iflag,ix,i,IQFROMQ
      integer nx,nq
      integer nset,iset,isetlast
      data isetlast/-1/
      real*8 xmin,xmax,qmin,qmax,S
      save iq0,iqc,iqb,nf0,mc2,mb2,tc2,tb2,mt2
      save nx,xmin,xmax,nq,qmin,qmax,gridname
*
c      print *,'QCDNUM4evolve'
      call getnset(iset)
      if (iset.ne.isetlast) then
        call get_pdfqcd(iset)
        isetlast = iset
      endif
c      
      Q2=Q*Q
      nf=4
      if (Q2.lt.tc2) nf=3
      pdf(0)=QPDFXQ('GLUON'  ,x,Q2,iflag)
      singlet=QPDFXQ('SINGLET',x,Q2,iflag)
      dm= QPDFXQ('DM',x,Q2,IFLAG)
      um= QPDFXQ('UM',x,Q2,IFLAG)
      dp= QPDFXQ('DP',x,Q2,IFLAG)
      up= QPDFXQ('UP',x,Q2,IFLAG)
      sp= QPDFXQ('SP',x,Q2,IFLAG)
      ub=0.5d0*(up-um+singlet/dble(nf))
      db=0.5d0*(dp-dm+singlet/dble(nf))
      sb=0.5d0*(sp+singlet/dble(nf))
      cb=0.d0
      if (nf.ge.4) then
         cp= QPDFXQ('CP',X,Q2,IFLAG)
         cb=0.5d0*(cp+singlet/dble(nf))
      end if

      pdf(1)=dm+db
      pdf(2)=um+ub
      pdf(3)=sb
      pdf(4)=cb
      pdf(5)=0d0
      pdf(6)=0d0
      pdf(-1)=db
      pdf(-2)=ub
      pdf(-3)=sb
      pdf(-4)=cb
      pdf(-5)=0d0
      pdf(-6)=0d0
      return
*
      entry QCDNUM4alfa(alfas,Q)
      Q2=Q*Q
      nf=4
      if (Q2.lt.mc2) nf=3
      alfas=QALFAS(Q2,Qlam,nf,iflag)
      return
*
      entry QCDNUM4read(nset)
c      read(1,*) s1
c      print *,s1
c         gridname='large.grid'
c         nx=400
c         xmin=1d-6
c         xmax=1d0
c         nq=112
c         qmin=1d0
c         qmax=1d10
c         read(1,*) dummy
*      else
         read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
c*      endif
c      iset = nset
      return
*
      entry QCDNUM4init(nset,Eorder,Q2fit)
      call qninit
      call QNISET('ORDER',Eorder+1)
      call grxdef(nx,xmin)
      call grqdef(nq,qmin,qmax)
      call grqinp(Q2fit,1)
      call getQMassM(nset,4,mc)
      mc2=mc*mc
c      print *,mc
      CALL QNRSET('CMASS',mc)
      CALL QNRSET('MCALF',mc)
      call getThresholdM(nset,4,tc)
      tc2=tc*tc
      call GRQINP(tc2,1)
      CALL GRQINP(tc2-1.0d-4,1)
c      print *,'setting iq0',Q2fit
      call qthres(tc2,2d10)
      iq0=IQFROMQ(Q2fit)
      iqc=IQFROMQ(tc2)
      nf0= NFLGET(iq0)
c      print *,iq0,iqc,nf0,q2fit,tc2
c      nf0=3
c      print *,nf0
      CALL QNBOOK(2,'dm')
      CALL QNBOOK(3,'um')
      CALL QNBOOK(4,'dp')
      CALL QNBOOK(5,'up')
      CALL QNBOOK(6,'sp')
      CALL QNBOOK(7,'cp')
      CALL QNLSET('W1ANA',.TRUE.)
      CALL QNLSET('W2NUM',.TRUE.)
      CALL QNLSET('W2STF',.FALSE.)
      if (index(gridname,'none').eq.1) then
         call qnfilw(0,0)
      else
         qnerr=-1
         open(unit=2,status='old',file=gridname,
     .        form='unformatted',err=1)
         call QNREAD(2,1,qnerr)
 1       close(2)
         if (qnerr.ne.0) then
            write(*,*) 'Grid file problem: ',gridname
            if (qnerr.lt.0) then 
               write(*,*) 'Grid file does not exist'
               write(*,*) 'Calculating and creating grid file'
               call qnfilw(0,0)
               open(unit=2,status='unknown',file=gridname,
     .              form='unformatted')
               call QNDUMP(2)
               close(2)
            else
               write(*,*) 'Existing grid file is inconsistent'
               if (qnerr.eq.1)
     .              write(*,*) 'Defined grid different'
               if (qnerr.eq.2)
     .              write(*,*) 'Heavy quark weight table different'
               if (qnerr.eq.3)
     .              write(*,*) 'Charm mass different'
               if (qnerr.eq.4)
     .              write(*,*) 'Bottom mass different'
               stop
            endif
         endif
      endif
      return
*
      entry QCDNUM4pdf(nset)
c      print *,'entering QCDNUMpdf',nset
      call GetAlfas(nset,alfas0,scale0)
c      print *,alfas0,scale0
      Q2=scale0*scale0
      CALL QNRSET('ALFAS',alfas0)      
      CALL QNRSET('ALFQ0',Q2)
      DO ix = 1,nx
         xx = XFROMIX(ix)
c	 print *,'calling parmPDF',ix
         call parmPDF(nset,xx,f)
c	 if(ix.lt.6) print *,nset,xx,f
         singlet=0.d0
         do i=1,nf0
            singlet=singlet+f(i)+f(-i)
c	    print  *,i,singlet
         end do
         gl=f(0)
         dm=f(1)-f(-1)
         um=f(2)-f(-2)
         dp=f(1)+f(-1)-singlet/dble(nf0)
         up=f(2)+f(-2)-singlet/dble(nf0)
         sp=f(3)+f(-3)-singlet/dble(nf0)
c         print *,ix,iq0
         CALL QNPSET('SINGLET',ix,iq0,singlet)
         CALL QNPSET('GLUON',ix,iq0,gl)
         CALL QNPSET('DM',ix,iq0,DM)
         CALL QNPSET('UM',ix,iq0,UM)
         CALL QNPSET('DP',ix,iq0,DP)
         CALL QNPSET('UP',ix,iq0,UP)
         CALL QNPSET('SP',ix,iq0,SP)
      ENDDO
      CALL EVOLSG(iq0,1,nq)
      CALL EVOLNM('DM',iq0,1,nq)
      CALL EVOLNM('UM',iq0,1,nq)
      CALL EVLSEA4('dp',iq0,iqc,nq)
      CALL EVLSEA4('up',iq0,iqc,nq)
      CALL EVLSEA4('sp',iq0,iqc,nq)

c      print *,'calling evols - heavy...'

C--   Heavy quark evolution

      CALL EVOLCP4('cp',iqc,nq)
c
      call getnset(iset)
      call save_pdfqcd(iset)
      return
*
      end
      subroutine EVLSEA4(name,IQ0,IQC,NQGRI)
      implicit none
      
      CHARACTER*(*) name
      integer iq0,iqc,nqgri
      real*8 f34,f45,f43,f54,factor
      parameter(f34=1.d0/12.d0,f43=-1.d0/12.d0)

c      print *,iq0,iqc,nqgri
      If(IQ0.le.IQC)then
         CAll EVPLUS(name,IQ0,1,IQC)
         CALL QADDSI(name,IQC,f34)
         CAll EVPLUS(name,IQC,IQC,NQGRI)
      else if(IQ0.gt.IQC)then
         CAll EVPLUS(name,IQ0,IQC,NQGRI)
         CALL QADDSI(name,IQC,f43)
         CAll EVPLUS(name,IQC,1,IQC)
      end if         

      end


      subroutine EVOLCP4(name,IQC,NQGRI)
      implicit none
      
      CHARACTER*(*) name
      integer iq0,iqc,nqgri
      real*8 f4
      parameter(f4=-1.d0/4.d0)

c
c     First set to zero to avoid adding -1/4Singl at each iteration
c

      CAll QNPNUL(name)
      CALL QADDSI(name,IQC,f4)      
      CAll EVPLUS(name,IQC,IQC,NQGRI)
      
      end
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