LHAPDF 5.2.2 source code.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / inputPDF.f
1       
2       subroutine weightPDF(x)
3       implicit none
4       integer nset,imem,nfmax
5       real*8 x,Q
6       nset=1
7       call weightPDFM(nset,x)
8       return
9
10       entry GetNF(nfmax)
11       nset=1
12       call GetNFM(nset,nfmax)
13       return
14
15       entry GetThreshold(imem,Q)
16       nset=1
17       call GetThresholdM(nset,imem,Q)
18       return
19       
20       end
21       
22       subroutine parmPDF(nset,x,pdf)
23       implicit none
24       include 'parmsetup.inc'
25       integer nset
26       character*16 name(nmxset)
27       integer nmem(nmxset),ndef(nmxset),mem
28       common/NAME/name,nmem,ndef,mem
29       real*4 xp(40),x4,fgdis4
30       character*16 s1,s2
31       integer i,j,imem,nop,parmN,Nfunc(nmxset),Fw(nmxset),nfmax,M
32       real*8 x,N,b0,Poly,pdf(-6:6),Fparm(nopmax),F(nofmax)
33       real*8 Ccoef(nmxset,-6:6,nofmax),Fpow(nmxset,nofmax)
34       real*8 Q,Treshold(nmxset,-6:6)
35 c      data Treshold/39*0d0/
36       integer Fmap(nmxset,nofmax,npfmax)
37       integer Ftype(nmxset,nofmax),Fn(nmxset,nofmax),Ctype(nmxset,-6:6)
38       integer lhasilent
39       common/lhasilent/lhasilent
40       logical first
41       data first/.true./
42       save Nfunc,Fn,Fw,Fpow,Fmap,Ccoef,Fparm,Ftype,Ctype,first,Treshold
43             
44       do i=1,Nfunc(nset)
45          if (Ftype(nset,i).eq.1) then
46             Poly=1.0
47
48             do j=4,Fn(nset,i)
49          Poly=Poly+Fparm(Fmap(nset,i,j))*x**(float(j-3)/Fpow(nset,i))
50             enddo
51             Poly=Fparm(Fmap(nset,i,1))*Poly
52        F(i)=x**Fparm(Fmap(nset,i,2))*(1.0-x)**Fparm(Fmap(nset,i,3))*Poly
53          endif
54          if (Ftype(nset,i).eq.2) then
55             if (x.lt.0.9999999) then
56       Poly=Fparm(Fmap(nset,i,2))*log(x)+Fparm(Fmap(nset,i,3))*log(1.0-x)
57      .      +Fparm(Fmap(nset,i,4))*x
58      .      +Fparm(Fmap(nset,i,6))*log(1.0+x*exp(Fparm(Fmap(nset,i,5))))
59          F(i)=Fparm(Fmap(nset,i,1))*exp(Poly)
60             else
61                F(i)=0d0
62             endif
63          endif
64          if (Ftype(nset,i).eq.101) then
65        Poly=exp(Fparm(Fmap(nset,i,1)))*x**(Fparm(Fmap(nset,i,2))-1)
66      .       *(1d0-x)**Fparm(Fmap(nset,i,3))
67        Poly=
68      . Poly+(1d0+Fparm(Fmap(nset,i,4))*x)*(1d0-x)**Fparm(Fmap(nset,i,5))
69             b0=10d0
70             if (Poly.gt.b0) then
71                F(i)=Poly
72             elseif (Poly.lt.-b0) then
73                F(i)=0d0
74             else
75                F(i)=Poly+log(1d0+exp(-b0*Poly)-exp(-b0))/b0
76             endif
77          endif
78 c - to add the mrst2004 gluon convolution
79          if (Ftype(nset,i).eq.201) then
80            xp(2) = Fparm(Fmap(nset,i,1))
81            xp(3) = Fparm(Fmap(nset,i,2))
82            xp(23)= Fparm(Fmap(nset,i,3))
83            xp(16)= Fparm(Fmap(nset,i,4))
84            xp(5) = Fparm(Fmap(nset,i,5))
85            xp(40)= Fparm(Fmap(nset,i,6))
86            xp(24)= Fparm(Fmap(nset,i,7)) 
87            xp(20)= Fparm(Fmap(nset,i,8)) 
88            xp(4) = Fparm(Fmap(nset,i,9)) 
89            xp(1) = Fparm(Fmap(nset,i,10))
90            x4 = x
91            call gconv(x4,xp,fgdis4)
92           F(i) = FGDIS4
93          endif
94 c--
95       enddo
96       do i=-6,6
97          pdf(i)=0.0
98          if (Ctype(nset,i).gt.0) then
99             if (Ctype(nset,i).eq.1) then
100                do j=1,Nfunc(nset)
101                   pdf(i)=pdf(i)+Ccoef(nset,i,j)*F(j)
102 c                 print *,i,j,Ccoef(i,j),F(j)
103                enddo
104             endif
105             if (Ctype(nset,i).eq.101) then
106                if (i.eq.-2) then
107        pdf(i)=F(int(Ccoef(nset,i,1)))/(F(int(Ccoef(nset,i,2)))+1d0)
108                endif
109                if (i.eq.-1) then
110       pdf(i)=
111      .F(int(Ccoef(nset,i,1)))/(1d0/F(int(Ccoef(nset,i,2)))+1d0)
112                endif
113                if (i.eq.1) then
114                   pdf(i)=F(1)+pdf(-1)
115                endif
116                if (i.eq.2) then
117                   pdf(i)=F(2)+pdf(-2)
118                endif
119             endif
120          endif
121       enddo
122 c      print *,pdf
123       return
124 *
125       entry weightPDFM(nset,x)
126       if (Fw(nset).ge.0) then
127          x=Fparm(Fw(nset))
128       else
129          call numberPDF(nop)
130          x=1.0/float(nop)
131       endif
132       return
133 *
134 c      entry GetParmPDFM(nset,imem,x)
135       entry GetParmPDF(nset,imem,x)
136       x=Fparm(imem)
137       return
138 *
139       entry GetNfM(nset,nfmax)
140       nfmax=0
141       do i=1,6
142          if (Treshold(nset,-i).ge.0d0) nfmax=nfmax+1
143          if (Treshold(nset,i).ge.0d0) nfmax=nfmax+1
144       enddo
145       nfmax=nfmax/2
146       return
147 *
148       entry GetThresholdM(nset,imem,Q)
149       Q=Treshold(nset,imem)
150       return
151 *
152 c      entry InitEvolvePDFM(nset,imem)
153       entry InitEvolvePDF(nset,imem)
154 c      print*, 'calling listPDF',nset,imem
155       call listPDF(nset,imem,Fparm)
156 c      print *,Fparm
157       return
158 *
159 c      entry initInputPDFM(nset)
160       entry initInputPDF(nset)
161
162       if(first) then
163         do i=1,nmxset
164           do j=-6,6
165              Treshold(i,j)=0.0d0
166           enddo
167         enddo
168         first=.false.
169       endif
170       read(1,*) s1,Fw(nset),Nfunc(nset)
171 c      print *,s1,Fw,Nfunc
172       if(lhasilent.eq.0) then
173         write(*,*) 'Parametrization: ',s1
174         write(*,*)
175       endif
176       do i=1,Nfunc(nset)
177          Ftype(nset,i)=-1
178          read(1,*) s1,s2
179          if (index(s2,'x-taylor').eq.1) then
180             Ftype(nset,i)=1
181             read(1,*) FPow(nset,i),Fn(nset,i)
182             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
183          endif
184          if (index(s2,'log-pade').eq.1) then
185             Ftype(nset,i)=2
186             FPow(nset,i)=0d0
187             read(1,*) Fn(nset,i)
188             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
189          endif
190          if (index(s2,'cteq6-ratio').eq.1) then
191             Ftype(nset,i)=101
192             Fpow(nset,i)=0d0
193             Fn(nset,i)=5
194             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
195          endif
196          if (index(s2,'convol').eq.1) then
197             Ftype(nset,i)=201
198             read(1,*) Fn(nset,i)
199             read(1,*) (Fmap(nset,i,j),j=1,Fn(nset,i))
200          endif
201          if (Ftype(nset,i).lt.0) then
202             write(*,*) 'File description error:'
203             write(*,*) 'Unknown functional ',s2
204             stop
205          endif
206       enddo
207       read(1,*) s1
208       do i=-6,6
209          Ctype(nset,i)=-1
210          read(1,*) s1,s2
211 c        print *,s1,s2
212          if (index(s2,'none').eq.1) then
213             Ctype(nset,i)=0
214             Treshold(nset,i)=-1d0
215          endif
216          if (index(s2,'treshold').eq.1) then
217             Ctype(nset,i)=0
218             read(1,*) Treshold(nset,i)
219          endif
220          if (index(s2,'composite').eq.1) then
221             Ctype(nset,i)=1
222             Treshold(nset,i)=0d0
223             read(1,*) (Ccoef(nset,i,j),j=1,Nfunc(nset))
224          endif
225          if (index(s2,'cteq6-ratio').eq.1) then
226             Ctype(nset,i)=101
227             Treshold(nset,i)=0d0
228             read(1,*) (Ccoef(nset,i,j),j=1,3)
229          endif
230          if (Ctype(nset,i).lt.0) then
231             write(*,*) 'File description error:'
232             write(*,*) 'Unknown composit type ',s2
233             stop
234          endif
235       enddo
236       if (Fw(nset).ge.0) then
237          write(*,*) '***********************************************'
238          write(*,*) '* Note that this is a weigthed PDF set.       *'
239          write(*,*) '* See manual for proper use.                  *'
240          write(*,*) '***********************************************'
241       endif
242       return
243 *
244       end
245
246       FUNCTION BETA_LHA(X1,X2)
247       CALL GAMMA_LHA(X1,G1,IER)
248       IF(IER.NE.0) write(16,*) 'GAMMA_LHA ERROR: IER= ',IER,X1,X2
249       CALL GAMMA_LHA(X2,G2,IER)
250       IF(IER.NE.0) write(16,*) 'GAMMA_LHA ERROR: IER= ',IER,X1,X2
251       X3=X1+X2
252       CALL GAMMA_LHA(X3,G3,IER)
253       IF(IER.NE.0) write(16,*) 'GAMMA_LHA ERROR: IER= ',IER,X1,X2
254       BETA_LHA=G1*G2/G3
255       RETURN
256       END
257
258       SUBROUTINE GAMMA_LHA(XX,GX,IER)
259       IF(XX-34.5)6,6,4
260     4 IER=2
261       GX=1.E38
262       RETURN
263     6 X=XX
264       ERR=1.0E-6
265       IER=0
266       GX=1.0
267       IF(X-2.0)50,50,15
268    10 IF(X-2.0)110,110,15
269    15 X=X-1.0
270       GX=GX*X
271       GO TO 10
272    50 IF(X-1.0)60,120,110
273 C        SEE IF X IS NEAR NEGATIVE INTEGER OR ZERO
274    60 IF(X-ERR)62,62,80
275    62 K=X
276       Y=FLOAT(K)-X
277       IF(ABS(Y)-ERR)130,130,64
278    64 IF(1.0-Y-ERR)130,130,70
279 C        X NOT NEAR A NEGATIVE INTEGER OR ZERO
280    70 IF(X-1.0)80,80,110
281    80 GX=GX/X
282       X=X+1.0
283       GO TO 70
284   110 Y=X-1.0
285       GY=1.0+Y*(-0.5771017+Y*(+0.9858540+Y*(-0.8764218+Y*(+0.8328212+
286      1Y*(-0.5684729+Y*(+0.2548205+Y*(-0.05149930)))))))
287       GX=GX*GY
288   120 RETURN
289   130 IER=1
290       RETURN
291       END
292
293       FUNCTION  ALPHA(T,AL)
294       COMMON/AINPUT/IORD,QSCT,QSDT
295       COMMON/PARAM/PARA(40)
296       DATA PI/3.14159/
297       DATA TOL/.0005/
298       ITH=0
299       TT=T
300       qsctt=qsct/4.
301       qsdtt=qsdt/4.
302 c      AL=para(1)
303       AL2=AL*AL
304       FLAV=4.
305       QS=AL2*EXP(T)
306
307       if(qs.lt.0.5d0) then   !!  running stops below 0.5
308           qs=0.5d0
309           t=alog(qs/al2)
310           tt=t
311       endif
312
313       IF(QS.gt.QSCTT) GO        TO 12  
314       IF(QS.lt.QSDTT) GO        TO 312  
315    11 CONTINUE
316       B0=11-2.*FLAV/3. 
317       IF(IORD)1,1,2
318 c     IF(IORD)2,2,2     !TAKE CARE !!
319     1 CONTINUE
320       ALPHA=4.*PI/B0/T
321       RETURN
322     2 CONTINUE
323       X1=4.*PI/B0
324       B1=102.-38.*FLAV/3.
325       X2=B1/B0**2
326       AS=X1/T*(1.-X2*aLOG(T)/T)
327     5 CONTINUE
328       F=-T+X1/AS-X2*aLOG(X1/AS+X2)
329       FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
330       AS2=AS-F/FP
331       DEL=ABS(F/FP/AS)
332       IF(DEL-TOL)3,3,4
333     3 CONTINUE
334       ALPHA=AS2
335       IF(ITH.EQ.0) RETURN
336       GO TO (13,14,15) ITH
337     4 CONTINUE
338       AS=AS2
339       GO TO 5
340    12 ITH=1
341       T=aLOG(QSCTT/AL2)
342       GO TO 11
343    13 ALFQC4=ALPHA
344       FLAV=5.   
345       ITH=2
346       GO TO 11
347    14 ALFQC5=ALPHA
348       ITH=3
349       T=TT
350       GO TO 11
351    15 ALFQS5=ALPHA
352       ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
353       ALPHA=1./ALFINV
354       RETURN
355
356   311 CONTINUE
357       B0=11-2.*FLAV/3. 
358       IF(IORD)31,31,32
359 c     IF(IORD)32,32,32  !TAKE CARE !!
360    31 CONTINUE
361       ALPHA=4.*PI/B0/T
362       RETURN
363    32 CONTINUE
364       X1=4.*PI/B0
365       B1=102.-38.*FLAV/3.
366       X2=B1/B0**2
367       AS=X1/T*(1.-X2*aLOG(T)/T)
368    35 CONTINUE
369       F=-T+X1/AS-X2*aLOG(X1/AS+X2)
370       FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
371       AS2=AS-F/FP
372       DEL=ABS(F/FP/AS)
373       IF(DEL-TOL)33,33,34
374    33 CONTINUE
375       ALPHA=AS2
376       IF(ITH.EQ.0) RETURN
377       GO TO (313,314,315) ITH
378    34 CONTINUE
379       AS=AS2
380       GO TO 35
381   312 ITH=1
382       T=aLOG(QSDTT/AL2)
383       GO TO 311
384   313 ALFQC4=ALPHA
385       FLAV=3.   
386       ITH=2
387       GO TO 311
388   314 ALFQC3=ALPHA
389       ITH=3
390       T=TT
391       GO TO 311
392   315 ALFQS3=ALPHA
393       ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
394       ALPHA=1./ALFINV
395       RETURN
396       END
397
398       SUBROUTINE WATE96
399 C*******************************************************************
400 C*****                                                         *****
401 C***** THE X(I) AND W(I) ARE THE DIRECT OUTPUT FROM A PROGRAM  *****
402 C***** USING NAG ROUTINE D01BCF TO CALCULATE THE               *****
403 C***** GAUSS-LEGENDRE WEIGHTS FOR 96 POINT INTEGRATION.        *****
404 C***** THEY AGREE TO TYPICALLY 14 DECIMAL PLACES WITH THE      *****
405 C***** TABLE IN ABRAMOWITZ & STEGUN, PAGE 919.                 *****
406 C*****                                                         *****
407 C***** ---->   PETER HARRIMAN, APRIL 3RD 1990.                 *****
408 C*****                                                         *****
409 C*******************************************************************
410       DIMENSION X(48),W(48)
411       COMMON/GAUS96/XI(96),WI(96),nterms,XX(97)
412       NTERMS=96
413
414       X( 1)=   0.01627674484960183561
415       X( 2)=   0.04881298513604856015
416       X( 3)=   0.08129749546442434360
417       X( 4)=   0.11369585011066471632
418       X( 5)=   0.14597371465489567682
419       X( 6)=   0.17809688236761733390
420       X( 7)=   0.21003131046056591064
421       X( 8)=   0.24174315616383866556
422       X( 9)=   0.27319881259104774468
423       X(10)=   0.30436494435449495954
424       X(11)=   0.33520852289262397655
425       X(12)=   0.36569686147231213885
426       X(13)=   0.39579764982890709712
427       X(14)=   0.42547898840729897474
428       X(15)=   0.45470942216774136446
429       X(16)=   0.48345797392059470382
430       X(17)=   0.51169417715466604391
431       X(18)=   0.53938810832435567233
432       X(19)=   0.56651041856139533470
433       X(20)=   0.59303236477757022282
434       X(21)=   0.61892584012546672523
435       X(22)=   0.64416340378496526886
436       X(23)=   0.66871831004391424358
437       X(24)=   0.69256453664216964528
438       X(25)=   0.71567681234896561582
439       X(26)=   0.73803064374439816819
440       X(27)=   0.75960234117664555964
441       X(28)=   0.78036904386743123629
442       X(29)=   0.80030874413913884180
443       X(30)=   0.81940031073792957139
444       X(31)=   0.83762351122818502758
445       X(32)=   0.85495903343459936363
446       X(33)=   0.87138850590929436968
447       X(34)=   0.88689451740241818933
448       X(35)=   0.90146063531585023110
449       X(36)=   0.91507142312089592706
450       X(37)=   0.92771245672230655266
451       X(38)=   0.93937033975275308073
452       X(39)=   0.95003271778443564022
453       X(40)=   0.95968829144874048809
454       X(41)=   0.96832682846326217918
455       X(42)=   0.97593917458513455843
456       X(43)=   0.98251726356301274934
457       X(44)=   0.98805412632962202890
458       X(45)=   0.99254390032376081654
459       X(46)=   0.99598184298720747465
460       X(47)=   0.99836437586317963722
461       X(48)=   0.99968950388322870559
462       W( 1)=   0.03255061449236316962
463       W( 2)=   0.03251611871386883307
464       W( 3)=   0.03244716371406427668
465       W( 4)=   0.03234382256857594104
466       W( 5)=   0.03220620479403026124
467       W( 6)=   0.03203445623199267876
468       W( 7)=   0.03182875889441101874
469       W( 8)=   0.03158933077072719007
470       W( 9)=   0.03131642559686137819
471       W(10)=   0.03101033258631386231
472       W(11)=   0.03067137612366917839
473       W(12)=   0.03029991542082762553
474       W(13)=   0.02989634413632842385
475       W(14)=   0.02946108995816795100
476       W(15)=   0.02899461415055528410
477       W(16)=   0.02849741106508543861
478       W(17)=   0.02797000761684838950
479       W(18)=   0.02741296272602931385
480       W(19)=   0.02682686672559184485
481       W(20)=   0.02621234073567250055
482       W(21)=   0.02557003600534944960
483       W(22)=   0.02490063322248370695
484       W(23)=   0.02420484179236479915
485       W(24)=   0.02348339908592633665
486       W(25)=   0.02273706965832950717
487       W(26)=   0.02196664443874448477
488       W(27)=   0.02117293989219144572
489       W(28)=   0.02035679715433347898
490       W(29)=   0.01951908114014518992
491       W(30)=   0.01866067962741165898
492       W(31)=   0.01778250231604547316
493       W(32)=   0.01688547986424539715
494       W(33)=   0.01597056290256253144
495       W(34)=   0.01503872102699521608
496       W(35)=   0.01409094177231515264
497       W(36)=   0.01312822956696188190
498       W(37)=   0.01215160467108866759
499       W(38)=   0.01116210209983888144
500       W(39)=   0.01016077053500880978
501       W(40)=   0.00914867123078384552
502       W(41)=   0.00812687692569928101
503       W(42)=   0.00709647079115442616
504       W(43)=   0.00605854550423662775
505       W(44)=   0.00501420274292825661
506       W(45)=   0.00396455433844564804
507       W(46)=   0.00291073181793626202
508       W(47)=   0.00185396078894924657
509       W(48)=   0.00079679206555731759
510       DO 1 I=1,48
511       XI(I)=-X(49-I)
512       WI(I)=W(49-I)
513       XI(I+48)=X(I)
514       WI(I+48)=W(I)
515     1 CONTINUE
516       DO 2 I=1,96
517     2 XX(I)=0.5*(XI(I)+1.)
518       XX(97)=1.0
519       EXPON=1.0
520       DO 3 I=1,96
521       YI=2.*(0.5*(1.+XI(I)))**EXPON-1.
522       WI(I)=WI(I)/(1.+XI(I))*(1.+YI)*EXPON
523       XI(I)=YI
524       XX(I)=0.5*(1.+YI)
525     3 CONTINUE
526       RETURN
527       END
528       
529       subroutine gconv(x,xp,fgdis)
530       COMMON/AINPUT/IORD,QSCT,QSDT
531       common/GAUS96/XI(96),WI(96),NTERMS,XX(97)
532       dimension xp(40)
533       logical first
534       data first/.true./
535       if(first) then
536          call wate96
537          first=.false.
538       endif
539       PI = 3.14159
540       PI2 = PI*PI
541       iord = 1
542       qsdt=8.18    !!  This is the value of 4m_c^2
543       qsct=74.0    !!  This is the value of 4m_b^2
544       cf = 4./3.
545       eta4 = xp(40)
546       T=alog(1/xp(1)**2)
547 c      AL = 0.550/(4.* 3.14159)
548       AL=ALPHA(T,xp(1))/(4.* pi)
549       rx=sqrt(x)
550       FF1=BETA_LHA(XP(2),XP(3)+1.)+XP(16)*BETA_LHA(XP(2)+1.,XP(3)+1.)+
551      XXP(23)*BETA_LHA(XP(2)+0.5,XP(3)+1.)
552       FF2=BETA_LHA(XP(2)+1.,XP(3)+1.)+
553      XXP(16)*BETA_LHA(XP(2)+2.,XP(3)+1.)+
554      XXP(23)*BETA_LHA(XP(2)+1.5,XP(3)+1.)
555       FF3=BETA_LHA(XP(5),ETA4+1.)+XP(20)*BETA_LHA(XP(5)+1.,ETA4+1.)+
556      XXP(24)*BETA_LHA(XP(5)+0.5,ETA4+1.)
557       FF4=BETA_LHA(XP(5)+1.,ETA4+1.)+XP(20)*BETA_LHA(XP(5)+2.,ETA4+1.)+
558      XXP(24)*BETA_LHA(XP(5)+1.5,ETA4+1.)
559       COEFU=2.*XP(4)/FF1
560       COEFD=   XP(4)/FF3
561 c          print *,'coefu ',coefu
562 c          print *,'coefd ',coefd
563       UV=coefu*X**XP(2)*(1.-X)**XP(3)*(1.+XP(16)*X+XP(23)*SQRT(X))
564       DV=coefd*X**XP(5)*(1.-X)**ETA4*(1.+XP(20)*X+XP(24)*SQRT(X))
565       FGDIS=al*CF*(-9.-2.*PI2/3.+alog(1.-x)*(-3.+
566      .2.*alog(1.-x)))*(UV+DV)
567      
568       DO 23 M=1,NTERMS
569       Y=0.5*(1.-X)*XI(M)+0.5*(1.+X)
570       XY=X/Y
571       UVXY=coefu*XY**XP(2)*(1.-XY)**XP(3)*(1.+XP(16)*XY
572      .+XP(23)*SQRT(XY))
573       DVXY=coefd*XY**XP(5)*(1.-XY)**ETA4*(1.+XP(20)*XY
574      .+XP(24)*SQRT(XY))
575       AL1=ALOG(1.-Y)
576       C22=CF*(6.+4.*Y-2.*(1.+Y*Y)/(1.-Y)*ALOG(Y)-2.*(1.+Y)*ALOG(1.-Y))
577       C23=CF*(-3.+4.*ALOG(1.-Y))/(1.-Y)
578
579       FGDIS=FGDIS+.5*(1.-X)*WI(M)*al*
580      .(C22*(uvxy+dvxy)+C23*(uvxy+dvxy-uv-dv))
581
582    23 CONTINUE
583       
584       return
585       end