|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
|
c**********************************************************************
c**********************************************************************
c*** vegas program
c*** vegas()---- integration over (0,1) integratal variables
c*** generand--- generate random number according to grade function
c*** randa------ generate a group of random number with pyr(0)
c**********************************************************************
c**********************************************************************
c***********************************************************************
c...changes are made which make the program only to get the importance
c...function.
subroutine vegas(fxn,ndim,ncall,itmx,nprn)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c arguments: c
c fxn = function to be integrated/mapped c
c ndim = # dimensions c
c ncall = maximum total # of calls to the function c
c itmx = maximum # of iterations allowed c
c nprn = printout level: c
c >=2 additionally inf. about accumulated values c
c per iteration. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit double precision(a-h,o-z)
implicit integer (i-n)
c...pythia common block.
parameter (maxnup=500)
common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
&istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
&vtimup(maxnup),spinup(maxnup)
save /hepeup/
#include "invegas.h"
#include "bcvegpy_set_par.inc"
c...user process event common block.
common/grade/xi(NVEGBIN,10)
common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
& smin,smax,ymin,ymax,psetamin,psetamax
common/vegcross/vegsec,vegerr,iveggrade
dimension x(10),xin(NVEGBIN),r(NVEGBIN),ia(10),d(NVEGBIN,10)
c data alph/1.5d0/
alph=1.50d0
c...number of bins.
nd=NVEGBIN
c===============================================================
c a)) initializing some variables
c===============================================================
si =0.0d0
si2 =0.0d0
swgt =0.0d0
schi =0.0d0
scalls=0.0d0
calls =ncall
xnd =nd
ndm =nd-1
c...............................................................
c defining the initial intervals distribution
c...............................................................
if(iveggrade.eq.0) then
rc=1.0d0/xnd
do 7 j=1,ndim
xi(nd,j)=1.0d0
dr=0.0d0
do 7 i=1,ndm
dr=dr+rc
xi(i,j)=dr
7 continue
end if
c===============================================================
c b)) iterating loop
c===============================================================
do 1000 it=1,itmx
c if(it.eq.1) then
c call time(begin_time)
c write(*,'(a)') begin_time
c write(3,'(a)') begin_time
c end if
c...............................................................
c initializing iteration variables
c...............................................................
ti =0.0d0
sfun2=0.0d0
do 10 j=1,ndim
do 10 i=1,nd
d(i,j)=0.0d0
10 continue
do 11 jj=1,ncall
call generand(ndim,xnd,x,ia,wgt)
call phpoint(x,wt)
if(wt.lt.1.0d-15) then
tempfxn=0.0d0
else
tempfxn=fxn(x,wt)*wgt
end if
c...record the maximum differential cross-section.
if(tempfxn.gt.crossmax) then
crossmax=tempfxn
end if
fun=tempfxn/calls
fun2=fun*fun
weight=wgt/calls
ti=ti+fun
sfun2=sfun2+fun2
do j=1,ndim
iaj=ia(j)
d(iaj,j)=d(iaj,j)+fun2
end do
11 continue
c...............................................................
c computing the integral and error values
c...............................................................
ti2=ti*ti
tsi=dsqrt((sfun2*calls-ti2)/(calls-1.0d0))
wgta=ti2/tsi**2
si=si+ti*wgta
si2=si2+ti2
swgt=swgt+wgta
schi=schi+ti2*wgta
scalls=scalls+calls
avgi=si/swgt
sd=swgt*it/si2
chi2a=0.0d0
if(it.gt.1)chi2a=sd*(schi/swgt-avgi*avgi)/(it-1)
sd=1.0d0/dsqrt(sd)
err=sd*100.0d0/avgi
c...record the cross-section and the corresponding err obtained.
c...may be used for pythia running(initialization).
vegsec=avgi !pb
vegerr=err
c...............................................................
c printing
c...............................................................
c...cross unit pb.
if(nprn.ge.2) then
print 201,it,avgi,sd,err
write(3,201) it,avgi,sd,err
c call time(end_time)
c write(*,'(a)') end_time
c write(3,'(a)') end_time
end if
c778 format('integral value =',g10.4,'+/-',
c & g10.4,3x,' chi**2=',g10.4, /' ')
201 format('iter.no',i3,' acc.result(pb)==>',g14.5,'+/-',g10.4,
. '... % err=',g10.2)
c===============================================================
c c)) redefining the grid
c===============================================================
c...............................................................
c smoothing the f**2 valued stored for each interval
c...............................................................
do 23 j=1,ndim
xo=d(1,j)
xn=d(2,j)
d(1,j)=(xo+xn)/2.0d0
x(j)=d(1,j)
do 22 i=2,ndm
d(i,j)=xo+xn
xo=xn
xn=d(i+1,j)
d(i,j)=(d(i,j)+xn)/3.0d0
x(j)=x(j)+d(i,j)
22 continue
d(nd,j)=(xn+xo)/2.0d0
x(j)=x(j)+d(nd,j)
23 continue
c...............................................................
c computing the 'importance function' of each interval
c...............................................................
do 28 j=1,ndim
rc=0.0d0
do i=1,nd
r(i)=0.0d0
if(d(i,j).le.0.) go to 224
xo=x(j)/d(i,j)
r(i)=((xo-1.0d0)/xo/dlog(xo))**alph
224 rc=rc+r(i)
end do
c...............................................................
c redefining the size of each interval
c...............................................................
rc=rc/xnd
kk=0
xn=0.0d0
dr=0.0d0
i=0
25 kk=kk+1
dr=dr+r(kk)
xo=xn
xn=xi(kk,j)
26 if(rc.gt.dr) go to 25
i=i+1
dr=dr-rc
xin(i)=xn-(xn-xo)*dr/r(kk)
if(i.lt.ndm) go to 26
do i=1,ndm
xi(i,j)=xin(i)
end do
xi(nd,j)=1.0d0
28 continue
1000 continue
return
end
c**************************************************************
subroutine generand(ndim,xnd,x,ia,weight)
implicit double precision(a-h,o-z)
implicit integer (i-n)
#include "invegas.h"
#include "bcvegpy_set_par.inc"
common/grade/xi(NVEGBIN,10)
dimension ia(10),x(10),randx(10)
c...to get the subprocess cross-section.
common/subopen/subfactor,subenergy,isubonly
weight=1.0d0
call randa(ndim,randx)
c...............................................................
c computing the point position
c...............................................................
do 15 j=1,ndim
xn=randx(j)*xnd+1.0d0
ia(j)=xn
xim1=0.0d0
if(ia(j).gt.1) xim1=xi(ia(j)-1,j)
xo=xi(ia(j),j)-xim1
x(j)=xim1+(xn-ia(j))*xo
weight=weight*xo*xnd
15 continue
if(isubonly.eq.1) then
x(6)=subfactor
x(7)=subfactor
end if
return
end
c*****************************************************************
subroutine randa(n,randx)
implicit double precision(a-h,o-z)
implicit integer (i-n)
dimension randx(n)
do i=1,n
randx(i)=pyr(0)
end do
return
end
|