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