bin/BCVEGPY/vegas.F

0001
0002 c**********************************************************************

0003 c**********************************************************************

0004 c*** vegas program

0005 c*** vegas()---- integration over (0,1) integratal variables

0006 c*** generand--- generate random number according to grade function

0007 c*** randa------ generate a group of random number with pyr(0)

0008 c**********************************************************************

0009 c**********************************************************************

0010
0011 c***********************************************************************

0012 c...changes are made which make the program only to get the importance

0013 c...function.

0014       subroutine vegas(fxn,ndim,ncall,itmx,nprn)
0015 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0016 c    arguments:                                                        c

0017 c       fxn    = function to be integrated/mapped                      c

0018 c       ndim   = # dimensions                                          c

0019 c       ncall  = maximum total # of calls to the function              c

0020 c       itmx   = maximum # of iterations allowed                       c

0021 c       nprn   = printout level:                                       c

0022 c              >=2  additionally inf. about accumulated values         c

0023 c                     per iteration.                                   c

0024 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0025       implicit double precision(a-h,o-z)
0026         implicit integer (i-n)
0027
0028 c...pythia common block.

0029       parameter (maxnup=500)
0030       common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0031      &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0032      &vtimup(maxnup),spinup(maxnup)
0033       save /hepeup/
0034
0035 #include "invegas.h"
0036 #include "bcvegpy_set_par.inc"
0037
0038 c...user process event common block.

0039         common/grade/xi(NVEGBIN,10)
0040         common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0041      &  smin,smax,ymin,ymax,psetamin,psetamax
0042         common/vegcross/vegsec,vegerr,iveggrade
0043
0044
0045       dimension x(10),xin(NVEGBIN),r(NVEGBIN),ia(10),d(NVEGBIN,10)
0046 c      data alph/1.5d0/

0047
0048       alph=1.50d0
0049 c...number of bins.

0050       nd=NVEGBIN
0051 c===============================================================

0052 c a))   initializing some variables

0053 c===============================================================

0054       si    =0.0d0
0055       si2   =0.0d0
0056       swgt  =0.0d0
0057       schi  =0.0d0
0058       scalls=0.0d0
0059         calls =ncall
0060       xnd   =nd
0061       ndm   =nd-1
0062 c...............................................................

0063 c  defining the initial intervals distribution

0064 c...............................................................

0065       if(iveggrade.eq.0) then
0066           rc=1.0d0/xnd
0067         do 7 j=1,ndim
0068           xi(nd,j)=1.0d0
0069           dr=0.0d0
0070         do 7 i=1,ndm
0071           dr=dr+rc
0072           xi(i,j)=dr
0073 7       continue
0074       end if
0075
0076 c===============================================================

0077 c b))  iterating loop

0078 c===============================================================

0079       do 1000 it=1,itmx
0080 c        if(it.eq.1) then

0081 c           call time(begin_time)

0082 c           write(*,'(a)') begin_time

0083 c           write(3,'(a)') begin_time

0084 c         end if

0085 c...............................................................

0086 c  initializing iteration variables

0087 c...............................................................

0088         ti   =0.0d0
0089         sfun2=0.0d0
0090         do 10 j=1,ndim
0091         do 10 i=1,nd
0092          d(i,j)=0.0d0
0093 10      continue
0094           do 11 jj=1,ncall
0095           call generand(ndim,xnd,x,ia,wgt)
0096             call phpoint(x,wt)
0097                 if(wt.lt.1.0d-15) then
0098               tempfxn=0.0d0
0099                 else
0100                   tempfxn=fxn(x,wt)*wgt
0101             end if
0102 c...record the maximum differential cross-section.

0103                 if(tempfxn.gt.crossmax) then
0104                crossmax=tempfxn
0105                 end if
0106           fun=tempfxn/calls
0107           fun2=fun*fun
0108           weight=wgt/calls
0109           ti=ti+fun
0110           sfun2=sfun2+fun2
0111           do j=1,ndim
0112             iaj=ia(j)
0113             d(iaj,j)=d(iaj,j)+fun2
0114           end do
0115 11      continue
0116 c...............................................................

0117 c  computing the integral and error values

0118 c...............................................................

0119         ti2=ti*ti
0120         tsi=dsqrt((sfun2*calls-ti2)/(calls-1.0d0))
0121         wgta=ti2/tsi**2
0122         si=si+ti*wgta
0123         si2=si2+ti2
0124         swgt=swgt+wgta
0125         schi=schi+ti2*wgta
0126         scalls=scalls+calls
0127         avgi=si/swgt
0128         sd=swgt*it/si2
0129         chi2a=0.0d0
0130         if(it.gt.1)chi2a=sd*(schi/swgt-avgi*avgi)/(it-1)
0131         sd=1.0d0/dsqrt(sd)
0132         err=sd*100.0d0/avgi
0133
0134 c...record the cross-section and the corresponding err obtained.

0135 c...may be used for pythia running(initialization).

0136           vegsec=avgi  !pb
0137           vegerr=err
0138 c...............................................................

0139 c  printing

0140 c...............................................................

0141 c...cross unit pb.

0142           if(nprn.ge.2) then
0143             print 201,it,avgi,sd,err
0144                 write(3,201) it,avgi,sd,err
0145 c           call time(end_time)

0146 c               write(*,'(a)') end_time

0147 c               write(3,'(a)') end_time

0148           end if
0149
0150 c778     format('integral value =',g10.4,'+/-',

0151 c     &   g10.4,3x,' chi**2=',g10.4, /' ')

0152 201     format('iter.no',i3,' acc.result(pb)==>',g14.5,'+/-',g10.4,
0153      .   '... % err=',g10.2)
0154 c===============================================================

0155 c  c))   redefining the grid

0156 c===============================================================

0157 c...............................................................

0158 c  smoothing the f**2 valued stored for each interval

0159 c...............................................................

0160         do 23 j=1,ndim
0161           xo=d(1,j)
0162           xn=d(2,j)
0163           d(1,j)=(xo+xn)/2.0d0
0164           x(j)=d(1,j)
0165           do 22 i=2,ndm
0166             d(i,j)=xo+xn
0167             xo=xn
0168             xn=d(i+1,j)
0169             d(i,j)=(d(i,j)+xn)/3.0d0
0170             x(j)=x(j)+d(i,j)
0171 22        continue
0172             d(nd,j)=(xn+xo)/2.0d0
0173             x(j)=x(j)+d(nd,j)
0174 23      continue
0175 c...............................................................

0176 c   computing the 'importance function' of each interval

0177 c...............................................................

0178         do 28 j=1,ndim
0179          rc=0.0d0
0180          do i=1,nd
0181           r(i)=0.0d0
0182           if(d(i,j).le.0.) go to 224
0183           xo=x(j)/d(i,j)
0184           r(i)=((xo-1.0d0)/xo/dlog(xo))**alph
0185 224       rc=rc+r(i)
0186          end do
0187 c...............................................................

0188 c  redefining the size of each interval

0189 c...............................................................

0190         rc=rc/xnd
0191         kk=0
0192         xn=0.0d0
0193         dr=0.0d0
0194         i=0
0195 25      kk=kk+1
0196         dr=dr+r(kk)
0197         xo=xn
0198         xn=xi(kk,j)
0199 26      if(rc.gt.dr) go to 25
0200         i=i+1
0201         dr=dr-rc
0202         xin(i)=xn-(xn-xo)*dr/r(kk)
0203         if(i.lt.ndm) go to 26
0204         do i=1,ndm
0205           xi(i,j)=xin(i)
0206         end do
0207         xi(nd,j)=1.0d0
0208 28      continue
0209 1000  continue
0210
0211         return
0212       end
0213
0214 c**************************************************************

0215       subroutine generand(ndim,xnd,x,ia,weight)
0216         implicit double precision(a-h,o-z)
0217       implicit integer (i-n)
0218
0219 #include "invegas.h"
0220 #include "bcvegpy_set_par.inc"
0221
0222         common/grade/xi(NVEGBIN,10)
0223         dimension ia(10),x(10),randx(10)
0224
0225 c...to get the subprocess cross-section.

0226       common/subopen/subfactor,subenergy,isubonly
0227
0228       weight=1.0d0
0229
0230       call randa(ndim,randx)
0231 c...............................................................

0232 c  computing the point position

0233 c...............................................................

0234         do 15 j=1,ndim
0235          xn=randx(j)*xnd+1.0d0
0236          ia(j)=xn
0237          xim1=0.0d0
0238          if(ia(j).gt.1) xim1=xi(ia(j)-1,j)
0239          xo=xi(ia(j),j)-xim1
0240          x(j)=xim1+(xn-ia(j))*xo
0241          weight=weight*xo*xnd
0242 15    continue
0243
0244       if(isubonly.eq.1) then
0245          x(6)=subfactor
0246            x(7)=subfactor
0247       end if
0248
0249       return
0250         end
0251
0252 c*****************************************************************

0253
0254       subroutine randa(n,randx)
0255       implicit double precision(a-h,o-z)
0256         implicit integer (i-n)
0257       dimension randx(n)
0258
0259       do i=1,n
0260          randx(i)=pyr(0)
0261       end do
0262
0263       return
0264       end