lqqbev.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lqqbev.F

C **********************************************************************

      SUBROUTINE lqqbev

      IMPLICIT NONE

C...Generate boson-gluon fusion event, choose xp and zp according to
C...QCD matrix elements and apply cuts for softness and collinearness.


*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LPTOU2)
*
      COMMON /leptou/ cut(14),lst(40),parl(30),
     &                x,y,w2,q2,u
      REAL cut,parl,x,y,w2,q2,u
      INTEGER lst
      SAVE /leptou/

      COMMON /linter/ pari(50),ewqc(2,2,8),qc(8),zl(2,4),zq(2,8),pq(17)
      REAL  pari,ewqc,qc,zl,zq,pq
      SAVE /linter/

      INTEGER  nlupdm,nplbuf
      parameter(nlupdm=4000,nplbuf=5)
      common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
      INTEGER  n,k
      REAL  p,v
      SAVE /lujets/

      common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
      INTEGER  mstu,mstj
      REAL  paru,parj
      SAVE /ludat1/

      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
      INTEGER  kchg
      REAL  pmas,parf,vckm
      SAVE /ludat2/

      
      INTEGER j1,j2,j3,j4,ifail,lqmcut
      REAL w,xp,zp,xt,pt,phi,pt2,epz,
     +amifl1,amifl3,amr1,amr2,tm2r1,tm2r2
      REAL ulmass,rlu,plu
      INTEGER ifl1,iflr2,iflr1,ifl3,ifl1a,ifl3a,nremh,ir1,ir2,
     +     i,j,ifl1s
      INTEGER lucomp,kfifl1,kfifl2

      lst(24)=3
      w=sqrt(w2)
      j1=mstu(1)
      j2=mstu(1)+1
      j3=mstu(1)+2
      j4=mstu(1)+3
       
      CALL lxp(xp,ifail)
      IF(ifail.NE.0) goto 999

C...Choose flavour of produced quark-antiquark pair.
  200 CALL lflav(ifl1,ifl3)

      IF(lst(21).NE.0) RETURN
      IF(ifl1.LT.0) THEN
C...Put quark in first position
        ifl1s=ifl1
        ifl1=ifl3
        ifl3=ifl1s
      ENDIF
      CALL lzp(xp,zp,ifail)
      IF(ifail.NE.0) goto 999
      ifl1a=iabs(ifl1)
      ifl3a=iabs(ifl3)
      mstj(93)=1
      amifl1=ulmass(ifl1)
      mstj(93)=1
      amifl3=ulmass(ifl3)


      IF(lst(14).EQ.0.OR.(lst(8).GE.2.AND.mod(lst(8),10).NE.9)) THEN
C...If baryon production from target remnant is neglected the
C...target remnant is approximated by a gluon.
        IF(w.LT.amifl1+amifl3+parj(32)) goto 999
        IF(lqmcut(xp,zp,amifl1,0.,amifl3).NE.0) goto 999
C       Pass the (consituent quark) masses we are using to LU3ENT
        p(j1,5)=amifl1
        p(j2,5)=0
        p(j3,5)=amifl3
        mstu(10)=1
        CALL lu3ent(j1,ifl1,21,ifl3,w,pari(21),pari(23))
        mstu(10)=2
        k(mstu(1)+1,3)=2
C...Align target remnant (gluon) along -z axis
        CALL lurobo(-acos(-p(j2,3)/sqrt(p(j2,3)**2+p(j2,1)**2)),
     &  0.,0.,0.,0.)
C...Phi-rotation to bring quark to phi=0.
        CALL lurobo(0.,-plu(j1,15),0.,0.,0.)
      ELSE

        IF(w.LT.amifl1+amifl3+0.9+2.*parj(32)) goto 999
        IF(lqmcut(xp,zp,amifl1,1.,amifl3).NE.0) goto 999
        p(j1,5)=amifl1
        p(j3,5)=amifl3
C...Choose target valence quark/diquark to form jet system with
C...produced antiquark/quark.
        iflr2=int(1.+lst(22)/3.+rlu(0))
        IF(iflr2.EQ.lst(22)) THEN
          iflr1=2101
          IF(rlu(0).GT.parl(4)) iflr1=2103
        ELSE
          iflr1=1000*iflr2+100*iflr2+3
        ENDIF
        iflr2=3-iflr2
        mstj(93)=1
        amr1=ulmass(iflr1)
CJR--
        kfifl1=lucomp(iflr1)
        IF (kfifl1.EQ.90) THEN
          amr1=amr1-2*parl(20)
        ELSEIF (1.LE.kfifl1 .AND. kfifl1.LE.6) THEN
          amr1=amr1-parl(20)
        ENDIF
        mstj(93)=1
        amr2=ulmass(iflr2)
        kfifl2=lucomp(iflr2)
        IF (kfifl2.EQ.90) THEN
          amr2=amr2-2.*parl(20)
        ELSEIF (1.LE.kfifl2 .AND. kfifl2.LE.6) THEN
          amr2=amr2-parl(20)
        ENDIF
CJR--

        nremh=0
  310   nremh=nremh+1
        IF(nremh.GT.100) goto 999
        CALL lprikt(parl(14),pt,phi)
        CALL lremh(0,pt,iflr1,iflr2,xt)
        pt2=pt**2
        tm2r1=amr1**2+pt2
        tm2r2=amr2**2+pt2
        p(j2,5)=sqrt(tm2r1/(1.-xt)+tm2r2/xt)
        IF(lqmcut(xp,zp,amifl1,p(j2,5),amifl3).NE.0) goto 310
        mstu(10)=1
        CALL lu3ent(j1,ifl1,21,ifl3,w,pari(21),pari(23))
        mstu(10)=2
C...Align target remnant (effective gluon) along -z axis
        CALL lurobo(-acos(-p(j2,3)/sqrt(p(j2,3)**2+p(j2,1)**2)),
     &0.,0.,0.,0.)
C...Phi-rotation to bring quark to phi=0.
        CALL lurobo(0.,-plu(j1,15),0.,0.,0.)
        epz=p(j2,4)-p(j2,3)
        IF(ifl1.GT.0) THEN
          ir1=j2
          ir2=j4
        ELSE
          ir1=j4
          ir2=j2
        ENDIF
        p(ir1,1)=pt*cos(phi)
        p(ir1,2)=pt*sin(phi)
        p(ir1,3)=-0.5*((1.-xt)*epz-tm2r1/(1.-xt)/epz)
        p(ir1,4)= 0.5*((1.-xt)*epz+tm2r1/(1.-xt)/epz)
        p(ir1,5)=amr1
        p(ir2,1)=-p(ir1,1)
        p(ir2,2)=-p(ir1,2)
        p(ir2,3)=-0.5*(xt*epz-tm2r2/xt/epz)
        p(ir2,4)= 0.5*(xt*epz+tm2r2/xt/epz)
        p(ir2,5)=amr2
        k(ir1,1)=1
        k(ir1,2)=iflr1
        k(ir2,1)=1
        k(ir2,2)=iflr2
        k(j3,1)=2
        DO 320 i=j1,j4
          DO 320 j=3,5
  320       k(i,j)=0
        n=j4
        k(ir1,3)=2
        k(ir2,3)=2
       IF((p(j1,4)+p(j2,4))**2-(p(j1,1)+p(j2,1))**2-(p(j1,3)+p(j2,3))**2
     &  -p(j2,2)**2.LT.(p(j1,5)+p(j2,5)+parj(32))**2) goto 310
       IF((p(j3,4)+p(j4,4))**2-(p(j3,1)+p(j4,1))**2-(p(j3,3)+p(j4,3))**2
     &  -p(j4,2)**2.LT.(p(j3,5)+p(j4,5)+parj(32))**2) goto 310

      ENDIF

      CALL lazimu(xp,zp)
      lst(21)=0
      RETURN

  999 lst(21)=5
      RETURN
      END