3342cc77851a55609be5bafd0453da3d1cb32172
[u/mrichter/AliRoot.git] / EVGEN / PartEloss / integralgauss.f
1 *******************************************************************
2       SUBROUTINE qsimp(sr,sd)
3       REAL*8           xr,length,noc,qq,nnorm,rri,fracc
4       COMMON /input/   xr,length,noc,qq,nnorm,rri,fracc
5       INTEGER JMAX
6       REAL*8 a1,b1,diff1,d1,sr,sd,EPS,fu1r,fu1d
7       EXTERNAL func1, func2
8       PARAMETER (EPS=1.e-9, JMAX=12)
9       INTEGER j,j1
10       REAL*8 osr,ostr,str
11       REAL*8 del,sumr,x
12       ostr=-1.e10
13       osr=-1.e10
14       a1=0.0
15       b1=(0.99999-xr)*fracc
16       diff1 = b1-a1
17       d1 = 0.5*diff1
18 *
19       res=dgauss(func2,a1,b1,1.d-6)
20       call func1(a1,fu1r,fu1d)
21       sd=fu1d
22       sr=res
23       END
24
25       function func2(yy)
26       implicit double precision (a-h,o-z)
27       call func1(yy,fu1r,fu1d)
28       func2=fu1r
29       return
30       end
31
32
33 **************************************************************
34 *
35       SUBROUTINE func1(yy,funr,fund)
36 *
37       REAL*8           funr,yy,fund
38       REAL*8           xr,length,noc,qq,nnorm,rri,fracc
39       COMMON /input/   xr,length,noc,qq,nnorm,rri,fracc
40       REAL*8           fragm
41       EXTERNAL         lookup
42       REAL*8           cont, disc, wwt, tepsi
43 *
44       tepsi = yy
45       wwt = tepsi
46       if(wwt.ge.1.3) then
47          call lookup(rri,1.d0,cont,disc)
48          funr=0.0
49          fund=disc
50       else
51          call lookup(rri,wwt,cont,disc)
52          funr = cont*fragm(xr/(1.0-tepsi/fracc),qq)
53      .         /(1.0-tepsi/fracc)
54          fund = disc
55       endif
56       END
57 *******************************************************************
58       SUBROUTINE lookup(rrrr,xxxx,continuous,discrete)
59 *
60       REAL*8           xx(400), da(30), ca(30,260), rrr(30)
61       COMMON /data/    xx, da, ca, rrr
62       REAL*8           rrrr,xxxx, continuous, discrete
63       REAL*8           rrin, xxin
64       INTEGER          nrlow, nrhigh, nxlow, nxhigh
65       REAL*8           rrhigh, rrlow, rfraclow, rfrachigh
66       REAL*8           xfraclow, xfrachigh
67       REAL*8           clow, chigh
68 *
69       rrin = rrrr
70       xxin = xxxx
71 *
72 *    determine the tabulated values xx(nxlow), xx(nxhigh)
73 *    rrlow, rrhigh such that
74 *    xx(nxlow) < xxin <  xx(nxhigh)
75 *    rrlow < rrin < rrhigh
76 *
77       nxlow = int(xxin/0.005) + 1
78       nxhigh = nxlow + 1
79       xfraclow = (xx(nxhigh)-xxin)/0.005
80       xfrachigh = (xxin - xx(nxlow))/0.005
81 *
82       do 666, nr=1,30
83          if (rrin.lt.rrr(nr)) then
84             rrhigh = rrr(nr)
85          else
86             rrhigh = rrr(nr-1)
87             rrlow = rrr(nr)
88             nrlow = nr
89             nrhigh = nr-1
90             goto 665
91          endif
92  666     enddo
93  665     continue
94 *
95       rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
96       rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
97 *
98       clow = xfraclow*ca(nrlow,nxlow)+xfrachigh*ca(nrlow,nxhigh)
99       chigh = xfraclow*ca(nrhigh,nxlow)+xfrachigh*ca(nrhigh,nxhigh)
100       continuous = rfraclow*clow + rfrachigh*chigh
101       discrete = rfraclow*da(nrlow) + rfrachigh*da(nrhigh)
102 *
103       END
104
105 ***************************************************************
106
107 c       BKK FF
108
109       FUNCTION fragmbkk(xxx,qqq)
110       REAL*8   alphav, betav, gammav, nv
111       REAL*8   alphas, betas, gammas, ns
112       REAL*8   sbar, xx, qq, xxx, qqq, lambda, fragv, frags
113 *
114       xx = xxx
115       qq = qqq
116       lambda = 1.0
117       sbar=log(log(qq*qq/(lambda*lambda))/log(4.0/(lambda*lambda)))
118 *
119       alphav = -1.0 - 0.0272*sbar
120       betav = 1.2 + 0.67*sbar
121       gammav = -0.393*sbar
122       nv = 0.551 - 0.053*sbar - 0.032*sbar*sbar
123 *
124 *      alphav = -1.0 - 0.059*sbar
125 *      betav = 1.2 + 0.6*sbar
126 *      gammav = -0.163*sbar
127 *      nv = 0.338 - 0.064*sbar - 0.0105*sbar*sbar
128 *
129       alphas = -1.0 + 0.447*sbar - 0.266*sbar*sbar
130       betas = 4.7 - 2.88*sbar + 2.05*sbar*sbar
131       gammas = -9.01*sbar + 4.36*sbar*sbar
132       ns = 1.23 + 2.85*sbar - 1.6*sbar*sbar
133 *
134 *      alphas = -1.0 + 0.757*sbar - 0.537*sbar*sbar
135 *      betas = 5.26 - 5.22*sbar + 3.62*sbar*sbar
136 *      gammas = -13.6*sbar + 8.17*sbar*sbar
137 *      ns = 1.19 + 4.20*sbar - 2.86*sbar*sbar
138 *
139       fragv = nv*(xx**alphav)*((1.0-xx)**betav)*((1.0+xx)**gammav)
140       frags = ns*(xx**alphas)*((1.0-xx)**betas)*((1.0+xx)**gammas)
141 c      fragmbkk = (fragv+frags)
142       fragmbkk = fragv
143       END
144
145 ***************************************************************
146
147       FUNCTION fragm(xxx,qqq)
148       REAL*8   alpha, beta, gamma, n
149       REAL*8   sbar, xx, qq, xxx, qqq, lambda, fragv, frags
150 *
151       xx = xxx
152       qq = qqq
153       lambda = 0.088
154       sbar=log(log(qq*qq/(lambda*lambda))/log(2.0/(lambda*lambda)))
155 *
156
157 c       u or d -> pi+ + pi-
158
159       n=0.54610-0.22946*sbar-0.22594*sbar**2+0.21119*sbar**3
160       alpha=-1.46616-0.45404*sbar-0.12684*sbar**2+0.27646*sbar**3
161       beta=1.01864+0.95367*sbar-1.09835*sbar**2+0.74657*sbar**3
162       gamma=-0.01877*sbar+0.02949*sbar**2
163
164 c       g -> pi+ + pi-
165
166       n=6.04510-6.61523*sbar-1.64978*sbar**2+2.68223*sbar**3
167       alpha=-.71378+0.14705*sbar-1.08423*sbar**2-.43182*sbar**3
168       beta=2.92133+1.48429*sbar+1.32887*sbar**2-1.78696*sbar**3
169       gamma=0.23086*sbar-0.29182*sbar**2
170
171       fragm = n*xx**alpha*(1.-xx)**beta*(1.+gamma/xx)/2.
172
173       END
174