bin/BCVEGPY/evntinit.F

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

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

0004       subroutine evntinit
0005       implicit double precision(a-h, o-z)
0006         implicit integer(i-n)
0007
0008 #include "invegas.h"
0009 #include "bcvegpy_set_par.inc"
0010
0011 c...three pythia functions return integers, so need declaring.

0012       integer pycomp
0013 c...external statement.

0014       external pydata,totfun
0015
0016 c****************************************************

0017 c*************** define most of the common block used

0018 c****************************************************

0019
0020 c...pythia common block.

0021       common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
0022       parameter (maxnup=500)
0023       common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0024      &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0025      &vtimup(maxnup),spinup(maxnup)
0026       save /hepeup/
0027
0028       parameter (maxpup=100)
0029       integer pdfgup,pdfsup,lprup
0030       common/heprup/idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
0031      &idwtup,nprup,xsecup(maxpup),xerrup(maxpup),xmaxup(maxpup),
0032      &lprup(maxpup)
0033       save /heprup/
0034
0035 c...user process event common block.

0036       common/pypars/mstp(200),parp(200),msti(200),pari(200)
0037       common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
0038       common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
0039         double complex colmat,bundamp
0040       common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0041      &  colmat(10,64),bundamp(4),pmomzero(5,8)
0042         common/counter/ibcstate,nev
0043         common/grade/xi(NVEGBIN,10)
0044         common/vegcross/vegsec,vegerr,iveggrade
0045         common/colflow/amp2cf(10),smatval
0046
0047 c...transform of running information.

0048         common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0049      &  smin,smax,ymin,ymax,psetamin,psetamax
0050         common/confine/ptcut,etacut
0051         common/histcol/inx
0052
0053 c...generate---switch for full events. vegasopen---switch for using

0054 c...vegas.

0055         logical generate,vegasopen,usegrade
0056         common/genefull/generate
0057
0058 c...parameters transformtion.

0059       common/funtrans/nq2,npdfu
0060         common/usertran/ishower,idpp
0061
0062 c...transform of the vegas information

0063       common/vegasinf/number,nitmx
0064
0065 c...to get the subprocess cross-section.

0066       common/subopen/subfactor,subenergy,isubonly
0067
0068 c...to get the distribution of an extra factor z=(2(k1+k2).p_bc)/shat.

0069       common/extraz/zfactor,zmin,zmax
0070         common/outpdf/ioutpdf,ipdfnum
0071
0072 c...used in grv98l

0073         common/intinip/iinip
0074         common/intinif/iinif
0075
0076 c...transform some variables

0077       common/loggrade/ievntdis,igenerate,ivegasopen,igrade
0078
0079 c...for transform the subprocess information, i.e.,  whether using

0080 c...the subprocess q\bar{q}->bc+b+\bar{c} to generate events.

0081       common/qqbar/iqqbar,iqcode
0082
0083         dimension nfin(20),alamin(20)
0084       common/coloct/ioctet
0085         common/octmatrix/coeoct
0086 c...transform the first derivative of the wavefunction at the zero or the

0087 c...wavefunction at the zero.

0088         common/wavezero/fbc
0089
0090 c...data:lambda and n_f values for parton distributions. this is obtained

0091 c...from pythia and used for vegas running.

0092       data alamin/0.177d0,0.239d0,0.247d0,0.2322d0,0.248d0,0.248d0,
0093      &0.192d0,0.326d0,2*0.2d0,0.2d0,0.2d0,0.29d0,0.2d0,0.4d0,5*0.2d0/,
0094      &nfin/20*4/
0095
0096 c...xbcsec(8) records the total differential cross-sections for different

0097 c...states: 1---singlet 1s0; 2---singlet 3s1; 7---octet 1s0; 8---octet 3s1;

0098 c...3---singlet 1p1; 4---singlet 3p0; 5---singlet 3p1; 6---singlet 3p2.

0099         common/mixevnt/xbcsec(8),imix,imixtype
0100         common/mixevnt2/xbcsum,ibclimit
0101         character*1 pfile
0102
0103 c************************************************

0104 c... open recording files.--record mainly the

0105 c... vegas running information

0106 c************************************************

0107         call upopenfile(imix,imixtype,ibcstate,ioctet)
0108
0109 c************************************************

0110
0111         generate =.false.
0112         if(igenerate.eq.1) generate =.true.
0113
0114         vegasopen=.false.
0115         if(ivegasopen.eq.1) vegasopen=.true.
0116
0117         usegrade=.false.
0118         if(igrade.eq.1) usegrade=.true.
0119
0120         if(isubonly.eq.1) subfactor=subenergy/ecm
0121
0122 c--------------------------------------------

0123 c...reference light-like momentum, which can be choosing arbitrarily.

0124 c...this can be used to check tht rightness of the program.

0125 c...--used only for s-wave gluon fusion mechanism.

0126         pmomup(1,8)=3.0d0
0127       pmomup(2,8)=4.0d0
0128         pmomup(3,8)=5.0d0
0129         pmomup(4,8)=dsqrt(pmomup(1,8)**2+pmomup(2,8)**2+pmomup(3,8)**2)
0130         pmomup(5,8)=0.0d0
0131
0132 c--------------------------------------------

0133 c...output program logo.

0134       call uplogo
0135
0136 c--------------------------------------------

0137 c...setting the alphas.

0138 c   mstu(111)=0---fixed; d=1---fir order; =2---sec order

0139         mstp(2)  =mstu(111)
0140
0141 c---------------------------------------------------------

0142 c...choose lambda value to be used in case of using vegas.

0143         if(vegasopen) then
0144            if(ioutpdf.eq.0) then
0145 c...this part is from the inner part of pythia.

0146               if(mstp(51).ge.1 .and. mstp(51).le.20) then
0147               paru(112)=alamin(mstp(51))
0148               mstu(112)=nfin(mstp(51))
0149             end if
0150            else
0151               if(ipdfnum.eq.100) then
0152                          paru(112)=0.1750d0    !grv98l
0153                      mstu(112)=4
0154               end if
0155                   if(ipdfnum.eq.200) then
0156                          paru(112)=0.220d0     !msrt2001l
0157                  mstu(112)=4
0158               end if
0159                   if(ipdfnum.eq.300) then
0160                          paru(112)=0.215d0     !cteq6l1
0161                  mstu(112)=4
0162               end if
0163            end if
0164         end if
0165
0166 c---------------------------------

0167 c...logo for vegas.

0168         call vegaslogo(vegasopen)
0169
0170 c******************************************************************

0171 c...using vegas to get the importance function and then the

0172 c...sampling stratage: xi(), which will be stored in grade files.

0173 c******************************************************************

0174
0175       if(vegasopen) then
0176           if(isubonly.eq.0) then
0177           if(imix.eq.0) then
0178 c--------------------------------------

0179 c... using existed grade to generate more precise grade.

0180 c... this way is only used for imix=0 (preferable);

0181 c... for imix=1, the way is the same,

0182 c... but is tedious to write all of them here.

0183 c--------------------------------------

0184                  if(iveggrade.eq.1) then
0185                write(*,'(a)')'existed grade to gene. more precise grade.'
0186                    write(3,'(a)')'existed grade to gene. more precise grade.'
0187                    do i=1,NVEGBIN
0188                   read(11,*) (xi(i,j),j=1,7)
0189                end do
0190 c              open(unit=29,file='data/newgrade.grid',status='unknown')

0191                open(unit=29,file='newgrade.grid',status='unknown')
0192              end if
0193                  call vegas(totfun,7,number,nitmx,2)
0194                  do i=1,NVEGBIN
0195                if(iveggrade.eq.1) then
0196                           write(29,*) (xi(i,j),j=1,7)
0197                else
0198                   write(11,*) (xi(i,j),j=1,7)
0199                    end if
0200                  end do
0201                 end if
0202
0203 c******************************************************

0204
0205           if(imix.eq.1) then
0206               if(imixtype.eq.1.or.imixtype.eq.2) then
0207                    ibcstate=1                         !1s0 and singlet
0208              write(*,'(a)') '...generating cross-sec and grade (1s0)...'
0209              write(3,'(a)') '...generating cross-sec and grade (1s0)...'
0210                    call paraswave(ibcstate)
0211                    call vegas(totfun,7,number,nitmx,2)
0212              xbcsec(1)=vegsec
0213                    do i=1,NVEGBIN
0214                    write(36,*) (xi(i,j),j=1,7)
0215                    end do
0216                write(36,*) '+'
0217              write(36,*) vegsec,crossmax
0218                    rewind(36)
0219               end if
0220
0221 c********************************

0222               if(imixtype.eq.1.or.imixtype.eq.2) then
0223                    ibcstate=2                         !3s1 and singlet
0224              write(*,'(a)') '...generating cross-sec and grade (3s1)...'
0225              write(3,'(a)') '...generating cross-sec and grade (3s1)...'
0226              call paraswave(ibcstate)
0227              call vegas(totfun,7,number,nitmx,2)
0228              xbcsec(2)=vegsec
0229              do i=1,NVEGBIN
0230                    write(37,*) (xi(i,j),j=1,7)
0231                    end do
0232                write(37,*) '+'
0233              write(37,*) vegsec,crossmax
0234                rewind(37)
0235             end if
0236
0237 c********************************

0238               if(imixtype.eq.1.or.imixtype.eq.3) then
0239                    ibcstate=3
0240              write(*,'(a)') '...generating cross-sec and grade (1p1)...'
0241              write(3,'(a)') '...generating cross-sec and grade (1p1)...'
0242              call parapwave
0243              call vegas(totfun,7,number,nitmx,2)
0244              xbcsec(3)=vegsec
0245              do i=1,NVEGBIN
0246                    write(38,*) (xi(i,j),j=1,7)
0247                    end do
0248                write(38,*) '+'
0249              write(38,*) vegsec,crossmax
0250                rewind(38)
0251               end if
0252
0253 c********************************

0254               if(imixtype.eq.1.or.imixtype.eq.3) then
0255                    ibcstate=4
0256              write(*,'(a)') '...generating cross-sec and grade (3p0)...'
0257              write(3,'(a)') '...generating cross-sec and grade (3p0)...'
0258              call parapwave
0259              call vegas(totfun,7,number,nitmx,2)
0260              xbcsec(4)=vegsec
0261              do i=1,NVEGBIN
0262                    write(39,*) (xi(i,j),j=1,7)
0263                    end do
0264                write(39,*) '+'
0265              write(39,*) vegsec,crossmax
0266                rewind(39)
0267              end if
0268
0269 c********************************

0270              if(imixtype.eq.1.or.imixtype.eq.3) then
0271                    ibcstate=5
0272              write(*,'(a)') '...generating cross-sec and grade (3p1)...'
0273              write(3,'(a)') '...generating cross-sec and grade (3p1)...'
0274              call parapwave
0275              call vegas(totfun,7,number,nitmx,2)
0276              xbcsec(5)=vegsec
0277              do i=1,NVEGBIN
0278                    write(46,*) (xi(i,j),j=1,7)
0279                    end do
0280                write(46,*) '+'
0281              write(46,*) vegsec,crossmax
0282                rewind(46)
0283               end if
0284
0285 c********************************

0286               if(imixtype.eq.1.or.imixtype.eq.3) then
0287                    ibcstate=6
0288              write(*,'(a)') '...generating cross-sec and grade (3p2)...'
0289              write(3,'(a)') '...generating cross-sec and grade (3p2)...'
0290              call parapwave
0291              call vegas(totfun,7,number,nitmx,2)
0292              xbcsec(6)=vegsec
0293              do i=1,NVEGBIN
0294                    write(47,*) (xi(i,j),j=1,7)
0295                    end do
0296                write(47,*) '+'
0297              write(47,*) vegsec,crossmax
0298                rewind(47)
0299               end if
0300
0301 c********************************

0302               if(imixtype.eq.1.or.imixtype.eq.3) then
0303                    ibcstate=7                         !1s0 and octet
0304              write(*,'(a)') '.generating cross-sec and grade (1s0)_8...'
0305              write(3,'(a)') '.generating cross-sec and grade (1s0)_8...'
0306              call paraswave(ibcstate)
0307              ibcstate=1                   !return to original definition
0308                    call vegas(totfun,7,number,nitmx,2)
0309              xbcsec(7)=vegsec
0310              do i=1,NVEGBIN
0311                    write(48,*) (xi(i,j),j=1,7)
0312                    end do
0313                write(48,*) '+'
0314              write(48,*) vegsec,crossmax
0315                rewind(48)
0316               end if
0317
0318 c********************************

0319               if(imixtype.eq.1.or.imixtype.eq.3) then
0320                    ibcstate=8                         !3s1 and octet
0321              write(*,'(a)') '.generating cross-sec and grade (3s1)_8...'
0322              write(3,'(a)') '.generating cross-sec and grade (3s1)_8...'
0323              call paraswave(ibcstate)
0324              ibcstate=2                   !return to original definition
0325              call vegas(totfun,7,number,nitmx,2)
0326              do i=1,NVEGBIN
0327                    write(49,*) (xi(i,j),j=1,7)
0328                    end do
0329              xbcsec(8)=vegsec
0330                write(49,*) '+'
0331              write(49,*) vegsec,crossmax
0332                rewind(49)
0333               end if
0334                 end if
0335           else
0336            write(*,'(a,a,1x,g12.5)')'getting subprocess info....',
0337      &           'at c.m. energy(gev)',subenergy
0338            write(3,'(a,a,1x,g12.5)')'getting subprocess info....',
0339      &           'at c.m. energy(gev)',subenergy
0340            write(*,'(a)') 'the info. of hadron collider is no use!!!'
0341            write(3,'(a)') 'the info. of hadron collider is no use!!!'
0342           if(iveggrade.eq.1) then
0343            write(*,'(a)')'using existed grade to gene. precise grade.'
0344            write(3,'(a)')'using existed grade to gene. precise grade.'
0345                    do i=1,NVEGBIN
0346                   read(11,*) (xi(i,j),j=1,7)
0347                end do
0348 c              open(unit=29,file='data/newgrade.grid',status='unknown')

0349                open(unit=29,file='newgrade.grid',status='unknown')
0350             end if
0351                 ndim=5
0352                 call vegas(totfun,ndim,number,nitmx,2)
0353 c...record the (new) generated grade.

0354             do i=1,NVEGBIN
0355                if(iveggrade.eq.1) then
0356                             write(29,*) (xi(i,j),j=1,ndim)
0357                else
0358                     write(11,*) (xi(i,j),j=1,ndim)
0359                    end if
0360                 end do
0361           end if
0362           call vegasend(vegasopen,ievntdis,usegrade)
0363       else
0364           call vegasend(vegasopen,ievntdis,usegrade)
0365           if(isubonly.eq.0) then
0366              ndim=7
0367           else
0368              ndim=5
0369              write(*,'(a)')'getting the info. of the subporcess....'
0370                  write(3,'(a)')'getting the info. of the subporcess....'
0371           end if
0372
0373 c...initialize the grade.

0374           rc=1.0d0/NVEGBIN
0375         do 77 j=1,ndim
0376            xi(NVEGBIN,j)=1.0d0
0377            dr=0.0d0
0378         do 77 i=1,NVEGBIN-1
0379            dr=dr+rc
0380            xi(i,j)=dr
0381 77      continue
0382
0383 c*****************************************************************

0384 c...using the existed grade. one thing should be care here is that

0385 c...the existed grade should be formed under the same parameters.

0386         if(imix.eq.0) then
0387             if(usegrade) then
0388              write(*,'(a)') 'using the existed vegas grade.'
0389                  write(3,'(a)') 'using the existed vegas grade.'
0390                  do i=1,NVEGBIN
0391                 read(11,*) (xi(i,j),j=1,7)
0392              end do
0393             end if
0394           end if
0395       end if
0396
0397 c********************************************

0398 c***  find the total-corssection for each state

0399 c***  from the previous generated grade, used

0400 c***  when vegas not used.

0401 c********************************************

0402
0403       if(imix.eq.1 .and. (.not.vegasopen)) then
0404
0405          if(imixtype.eq.1) then
0406           do i=1,4
0407            do
0408              read(36+i-1,*) pfile
0409              if(pfile.eq.'+') exit
0410            enddo
0411            read(36+i-1,*) xbcsec(i)
0412            rewind(36+i-1)
0413           enddo
0414           do i=1,4
0415            do
0416              read(46+i-1,*) pfile
0417              if(pfile.eq.'+') exit
0418            enddo
0419            read(46+i-1,*) xbcsec(4+i)
0420            rewind(46+i-1)
0421           enddo
0422          end if
0423 c---------------------------------------

0424          if(imixtype.eq.2) then
0425           do i=1,2
0426            do
0427              read(36+i-1,*) pfile
0428              if(pfile.eq.'+') exit
0429            enddo
0430            read(36+i-1,*) xbcsec(i)
0431            rewind(36+i-1)
0432           enddo
0433          end if
0434 c----------------------------------------

0435          if(imixtype.eq.3) then
0436           do i=3,4
0437            do
0438              read(36+i-1,*) pfile
0439              if(pfile.eq.'+') exit
0440            enddo
0441            read(36+i-1,*) xbcsec(i)
0442            rewind(36+i-1)
0443           enddo
0444           do i=1,4
0445            do
0446              read(46+i-1,*) pfile
0447              if(pfile.eq.'+') exit
0448            enddo
0449            read(46+i-1,*) xbcsec(4+i)
0450            rewind(46+i-1)
0451           enddo
0452          end if
0453         end if
0454
0455 c*******************************************

0456
0457         if(imix.eq.1) then
0458          if(imixtype.eq.1) then
0459           xbcsum=0.0d0
0460           do i=1,8
0461            xbcsum=xbcsum+xbcsec(i)
0462           end do
0463           ibclimit=8
0464          end if
0465
0466          if(imixtype.eq.2) then
0467           xbcsum=0.0d0
0468           do i=1,2
0469            xbcsum=xbcsum+xbcsec(i)
0470           end do
0471           ibclimit=2
0472          end if
0473
0474          if(imixtype.eq.3) then
0475           xbcsum=0.0d0
0476           do i=3,8
0477            xbcsum=xbcsum+xbcsec(i)
0478           end do
0479           ibclimit=8
0480          end if
0481         end if
0482
0483         crossmax=0.0d0
0484
0485 c***************************************************

0486 c***    some switches for event generation     *****

0487 c***************************************************

0488 c...switch off aspects: initial and final state

0489 c...showers, multiple interactions, hadronization.

0490       if(ishower.eq.0) then
0491           mstp(61) =0
0492         mstp(71) =0
0493         mstp(81) =0
0494         mstp(111)=0
0495         end if
0496
0497 c...expanded event listing (required for histogramming).

0498       mstp(125)=2
0499
0500 c...pyinit--->to initialize the generation procedure.

0501       call pyinit('user',' ',' ',0d0)
0502
0503 c...set bc+ stable.

0504       mdcy(pycomp(541),1)=0
0505 c***************************************************

0506
0507       return
0508         end