bin/BCVEGPY/genevnt.F

0001 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0002 c upinit----initial parameters for pythia.                           c

0003 c upevnt----call event process and possible error messages.          c

0004 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0005 c...       bc in color-singlet and color-octet states.               c

0006 c copyright (c) z.x zhang, chafik driouich, paula eerola and x.g. wu c

0007 c reference: comput.phys.commun. 159,192(2004);                      c

0008 c reference: comput.phys.commun. 174,241(2006);                      c

0009 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

0010
0011       subroutine upinit
0012
0013 c...double precision and integer declarations.

0014       implicit double precision(a-h, o-z)
0015       implicit integer(i-n)
0016
0017 c...selection of hard scattering subprocesses.

0018       common/pypars/mstp(200),parp(200),msti(200),pari(200)
0019       save /pypars/
0020
0021       INTEGER MAXPUP
0022       PARAMETER (MAXPUP=100)
0023       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0024       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0025       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0026      &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0027      &LPRUP(MAXPUP)
0028       SAVE /HEPRUP/
0029
0030 c...user process event common block.

0031       parameter (maxnup=500)
0032       common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0033      &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0034      &vtimup(maxnup),spinup(maxnup)
0035       save /hepeup/
0036
0037 c...extra commonblock to transfer run info.

0038         common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0039      &  smin,smax,ymin,ymax,psetamin,psetamax
0040         common/histcol/inx
0041
0042 c...the user own transfer of information.

0043       double complex colmat,bundamp
0044       common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0045      &  colmat(10,64),bundamp(4),pmomzero(5,8)
0046         common/counter/ibcstate,nev
0047         logical generate
0048         common/genefull/generate
0049         common/vegcross/vegsec,vegerr,iveggrade
0050
0051 c...parameters transformtion.

0052       common/funtrans/nq2,npdfu
0053         common/usertran/ishower,idpp
0054 c...transform some variables

0055       common/loggrade/ievntdis,igenerate,ivegasopen,igrade
0056
0057 C...Lines to read in assumed never longer than 200 characters.

0058       PARAMETER (MAXLEN=200)
0059
0060
0061
0062 c...set up incoming beams. tevotran

0063       if(npdfu.eq.1) then
0064            idbmup(1) = 2212
0065          idbmup(2) = -2212
0066         end if
0067
0068 c...set up incoming beams. lhc

0069         if(npdfu.eq.2) then
0070            idbmup(1) = 2212
0071            idbmup(2) = 2212
0072         end if
0073
0074       ebmup(1)  = 0.5d0*ecm
0075       ebmup(2)  = 0.5d0*ecm
0076
0077 c...set up the external process.

0078       idwtup   = idpp
0079       nprup    = 1
0080       lprup(1) = 1001
0081         idprup   = lprup(1)
0082
0083 c...set up g+g --> b_c + b + c~ : maximum cross section in pb.

0084 c...using the default xmaxup(1)=0 to make pyevnt accept almost

0085 c...all the upevnt events. crossmax is the maximum differential

0086 c...cross-section.

0087         if(idwtup.eq.1) then
0088           if(generate) then
0089          if(ivegasopen.eq.1) then
0090             xmaxup(1)=crossmax
0091            else
0092             write(*,'(a)')
0093      &   'warning: here should input a maximum differential cross-sec'
0094             write(*,'(a)')
0095      &   '!stop here! input a proper value in (subroutine upinit)!!!!'
0096         stop 'or running vegas to get the correct value!program stop!'
0097 c...      xmaxup(1)=100! the value added here! find the value in old runs.

0098 c...note: this value should be the one obtained under the same condition.

0099            end if
0100           else
0101              xmaxup(1)=0.0d0
0102           end if
0103         end if
0104
0105 c...the value of xsecup(1) can be given arbitrarily. all of them are

0106 c...only used for pystat(1) to produce a sensible cross-section table.

0107 c...the actural total cross-sec of xsecup is (appcross).

0108         if (idwtup.eq.3) then
0109           xsecup(1)=vegsec*1.0d+3  ! vegas value (nb) for initialization.
0110           xmaxup(1)=crossmax       ! maximum differential cross-section
0111         end if
0112
0113       return
0114       end
0115
0116 c**************************************************************************

0117
0118       subroutine upevnt
0119
0120 c...all real arithmetic in double precision.

0121       implicit double precision(a-h, o-z)
0122       implicit integer(i-n)
0123
0124 C...PYTHIA commonblock: only used to provide read unit MSTP(162).

0125       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0126       SAVE /PYPARS/
0127
0128
0129 c...user process event common block.

0130
0131       INTEGER MAXNUP
0132       PARAMETER (MAXNUP=500)
0133       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0134       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0135       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0136      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0137      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0138       SAVE /HEPEUP/
0139
0140       parameter (maxpup=100)
0141       integer pdfgup,pdfsup,lprup
0142       common/heprup/idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
0143      &idwtup,nprup,xsecup(maxpup),xerrup(maxpup),xmaxup(maxpup),
0144      &lprup(maxpup)
0145       save /heprup/
0146
0147 C...Lines to read in assumed never longer than 200 characters.

0148       PARAMETER (MAXLEN=200)
0149
0150 #include "invegas.h"
0151 #include "bcvegpy_set_par.inc"
0152
0153        common/grade/xi(NVEGBIN,10)
0154
0155         common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0156      &  smin,smax,ymin,ymax,psetamin,psetamax
0157         common/counter/ibcstate,nev
0158         logical generate
0159         common/genefull/generate
0160
0161 c...get the approximate total corss-section.

0162       common/totcross/appcross
0163
0164 c...to get the subprocess cross-section.

0165       common/subopen/subfactor,subenergy,isubonly
0166
0167 c...for transform the subprocess information, i.e.,  whether using

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

0169       common/qqbar/iqqbar,iqcode
0170       common/mixevnt/xbcsec(8),imix,imixtype
0171
0172       dimension x(10),ia(10)
0173
0174 cc#include "invegas.h"

0175
0176 c...if to get mixed results, then one should redefine grade

0177 c...(the sampling importance function), the old grade files

0178 c...must be exist and be generated under the same condition.

0179       if(imix.eq.1 .and. isubonly.eq.0) call initmixgrade
0180
0181 c...call the respective routine to generate event.

0182       if(idprup.eq.1001) then
0183          xnd  =NVEGBIN*1.0d0
0184          if(isubonly.eq.0) then
0185             ndim =7
0186          else
0187             ndim =5
0188          end if
0189
0190 c...using the generated grade to generate the events points.

0191          call generand(ndim,xnd,x,ia,wgt)
0192          call phpoint(x,wt)
0193          if(wt.lt.1.0d-16) then
0194             xwgtup=0.0d0
0195          else
0196             xwgtup=totfun(x,wt)*wgt
0197          end if
0198
0199          if(idwtup.eq.1.and.generate) then
0200 * a dirty trick: limiting the diff xsection to current xwgtup

0201 * the current xwgtup is high enough, so the total xsection

0202 * does not deviate much. this trick is from yu gouz's program

0203             if(xwgtup.gt.xmaxup(1)) then
0204                xwgtup = xmaxup(1)*0.9999999d0
0205             end if
0206          end if
0207
0208 c...calculate approximate crossection.

0209          appcross=appcross+xwgtup/nev*1.0d-3 !nb
0210
0211 c...record the maximum differential cross-section.

0212          if(xwgtup.gt.crossmax) then
0213             crossmax=xwgtup
0214          end if
0215 c...gluon-gluon fusion. for all the s- and p- wave states.

0216          if(iqqbar.eq.0) then
0217             call bcpythia(21)
0218          end if
0219 c...  quark-antiquark annihilation, only for color-singlet s-wave.

0220            if(iqqbar.eq.1 .and. (ibcstate.eq.1.or.ibcstate.eq.2)) then
0221              call bcpythia(iqcode)
0222            end if
0223       else
0224          write(*,*) 'fatal error! unknown process',idprup
0225          stop
0226       end if
0227
0228       return
0229       end
0230
0231 c**************************************************************************

0232
0233 C...bcvegpy_write_lhe

0234 C...write the Bc event in the format needed for the

0235 C...Les Houches event record.

0236
0237       SUBROUTINE bcvegpy_write_lhe
0238
0239 C...Double precision and integer declarations.

0240       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0241       IMPLICIT INTEGER(I-N)
0242
0243 C...Commonblocks.

0244       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0245       COMMON/PYCTAG/NCT,MCT(4000,2)
0246       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0247       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0248       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0249       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0250       COMMON/PYINT1/MINT(400),VINT(400)
0251       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0252       COMMON/PYINT4/MWID(500),WIDS(500,5)
0253       SAVE /PYJETS/,/PYCTAG/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,
0254      &/PYINT1/,/PYINT2/,/PYINT4/
0255
0256 C...HEPEUP for output.

0257       INTEGER MAXNUP
0258       PARAMETER (MAXNUP=500)
0259       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0260       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0261       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0262      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0263      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0264       SAVE /HEPEUP/
0265
0266 c IDWTUP = 3 requires XWGTUP = 1

0267       XWGTUP_wlhe = 1
0268
0269 C...Optionally write out event to disk. Minimal size for time/spin fields.

0270       IF(MSTP(162).GT.0) THEN
0271         WRITE(MSTP(162),5200) NUP,IDPRUP,XWGTUP_wlhe,SCALUP,AQEDUP
0272      +   ,AQCDUP
0273         DO 190 I=1,NUP
0274           IF(VTIMUP(I).EQ.0D0) THEN
0275             WRITE(MSTP(162),5300) IDUP(I),ISTUP(I),MOTHUP(1,I),
0276      &      MOTHUP(2,I),ICOLUP(1,I),ICOLUP(2,I),(PUP(J,I),J=1,5),
0277      &      ' 0. 9.'
0278           ELSE
0279             WRITE(MSTP(162),5400) IDUP(I),ISTUP(I),MOTHUP(1,I),
0280      &      MOTHUP(2,I),ICOLUP(1,I),ICOLUP(2,I),(PUP(J,I),J=1,5),
0281      &      VTIMUP(I),' 9.'
0282           ENDIF
0283   190   CONTINUE
0284
0285 C...Optional extra line with parton-density information.

0286         IF(MSTP(165).GE.1) WRITE(MSTP(162),5500) MSTI(15),MSTI(16),
0287      &  PARI(33),PARI(34),PARI(23),PARI(29),PARI(30)
0288       ENDIF
0289
0290 C...Print formats.

0291
0292  5200 FORMAT(1P,2I6,4E14.6)
0293  5300 FORMAT(1P,I8,5I5,5E18.10,A6)
0294  5400 FORMAT(1P,I8,5I5,5E18.10,E12.4,A3)
0295  5500 FORMAT(1P,'#pdf ',2I5,5E18.10)
0296
0297       RETURN
0298       END
0299
0300 C*********************************************************************

0301
0302 C...bcvegpy_PYUPIN

0303 C...Fills the HEPRUP commonblock with info on incoming beams and allowed

0304 C...processes, and optionally stores that information on file.

0305
0306       SUBROUTINE bcvegpy_PYUPIN
0307
0308 C...Double precision and integer declarations.

0309       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0310       IMPLICIT INTEGER(I-N)
0311
0312 C...Commonblocks.

0313       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0314       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0315       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0316       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0317       SAVE /PYJETS/,/PYSUBS/,/PYPARS/,/PYINT5/
0318
0319 C...User process initialization commonblock.

0320       INTEGER MAXPUP
0321       PARAMETER (MAXPUP=100)
0322       INTEGER IDBMUP,PDFGUP,PDFSUP,IDWTUP,NPRUP,LPRUP
0323       DOUBLE PRECISION EBMUP,XSECUP,XERRUP,XMAXUP
0324       COMMON/HEPRUP/IDBMUP(2),EBMUP(2),PDFGUP(2),PDFSUP(2),
0325      &IDWTUP,NPRUP,XSECUP(MAXPUP),XERRUP(MAXPUP),XMAXUP(MAXPUP),
0326      &LPRUP(MAXPUP)
0327       SAVE /HEPRUP/
0328
0329 C...Store info on incoming beams.

0330 c      IDBMUP(1)=K(1,2)

0331 c      IDBMUP(2)=K(2,2)

0332 c      EBMUP(1)=P(1,4)

0333 c      EBMUP(2)=P(2,4)

0334       PDFGUP(1)=0
0335       PDFGUP(2)=0
0336       PDFSUP(1)=MSTP(51)
0337       PDFSUP(2)=MSTP(51)
0338
0339 C...Event weighting strategy.

0340       IDWTUP=3
0341
0342 C...Info on individual processes.

0343 c      NPRUP=0

0344 c      DO 100 ISUB=1,500

0345 c        IF(MSUB(ISUB).EQ.1) THEN

0346 c          NPRUP=NPRUP+1

0347 c          XSECUP(NPRUP)=1D9*XSEC(ISUB,3)

0348 c          XERRUP(NPRUP)=XSECUP(NPRUP)/SQRT(MAX(1D0,DBLE(NGEN(ISUB,3))))

0349 c          XMAXUP(NPRUP)=1D0

0350 c          LPRUP(NPRUP)=ISUB

0351 c        ENDIF

0352 c  100 CONTINUE

0353
0354 C...Write info to file.

0355       IF(MSTP(161).GT.0) THEN
0356         WRITE(MSTP(161),5100) IDBMUP(1),IDBMUP(2),EBMUP(1),EBMUP(2),
0357      &  PDFGUP(1),PDFGUP(2),PDFSUP(1),PDFSUP(2),IDWTUP,NPRUP
0358         DO 110 IPR=1,NPRUP
0359           WRITE(MSTP(161),5200) XSECUP(IPR),XERRUP(IPR),XMAXUP(IPR),
0360      &    LPRUP(IPR)
0361   110   CONTINUE
0362       ENDIF
0363
0364 C...Formats for printout.

0365  5100 FORMAT(1P,2I8,2E14.6,6I6)
0366  5200 FORMAT(1P,3E14.6,I6)
0367
0368       RETURN
0369       END