lmeps.F

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

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

      SUBROUTINE lmeps

      IMPLICIT NONE


*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LEPTO2)
*

      COMMON /leptou/ cut(14),lst(40),parl(30),
     &                xlp,ylp,w2lp,q2lp,ulp
      REAL cut,parl,xlp,ylp,w2lp,q2lp,ulp
      INTEGER lst
      SAVE /leptou/

      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 /lboost/ dbeta(2,3),stheta(2),sphi(2),pb(5),phir
      DOUBLE PRECISION dbeta
      REAL stheta,sphi,pb,phir
      SAVE /lboost/

      COMMON /pypara/ ipy(80),pypar(80),pyvar(80)
      REAL pypar,pyvar
      INTEGER ipy
      SAVE /pypara/ 

*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LYPRO2)
*
      COMMON /lyproc/ isub,kfl(3,2),x(2),sh,th,uh,q2,xsec(0:40)
      REAL    x,sh,th,uh,q2,xsec
      INTEGER isub,kfl
      SAVE /lyproc/ 

      COMMON /lyint1/ xq(2,-6:6),dsig(-6:6,-6:6,5),fsig(10,10,3)
      REAL    xq,dsig,fsig
      SAVE /lyint1/ 
 


      INTEGER i,j,ks,ip2,ns,ifl,it,ipu1,ipu2
      REAL ps,xr,robo,xpq,qmax,t1,t2
      REAL ulangl
      DOUBLE PRECISION dpq2,drobo(5)
      DOUBLE PRECISION deltap(4),dplong,dbtot,dgamma,droot
      dimension ks(9,5),ps(9,5),robo(5),xpq(-6:6)
**HI      SAVE KS,PS

C     CALL GULIST(100,2)
C...Save event record in hadronic cms
      DO 10 i=1,7
      DO 10 j=1,5
      ks(i,j)=k(i,j)
   10 ps(i,j)=p(i,j)
C...Rearrange event record to PYSSPA standard
      ip2=6
      IF(lst(24).EQ.3) ip2=7
      DO 20 j=1,5
      k(3,j)=0.
      p(3,j)=0.
      k(4,j)=0
      p(4,j)=0.
      k(5,j)=ks(3,j)
      p(5,j)=ps(3,j)
      k(7,j)=ks(4,j)
      p(7,j)=ps(4,j)
      k(8,j)=ks(5,j)
      p(8,j)=ps(5,j)
      k(9,j)=ks(4,j)
      p(9,j)=ps(4,j)
      k(10,j)=ks(ip2,j)
   20 p(10,j)=ps(ip2,j)
      k(5,3)=3
      k(6,3)=4
      k(7,3)=5
      k(8,3)=6
      k(9,3)=5
      k(10,3)=6
      DO 30 i=5,10
   30 k(i,1)=21
      k(9,1)=0
C...Incoming parton = outgoing 2 parton - boson fourvectors
      DO 40 j=1,4
   40 p(6,j)=p(8,j)+p(10,j)-p(5,j)
      p(6,5)=0.
      k(6,2)=lst(25)
      IF(lst(24).EQ.3) k(6,2)=21
      n=10
C     CALL GULIST(101,2)

      xr=xlp
      dpq2=dble(q2lp)
C...Partons with colour information in hadronic cms frame.
      DO 120 i=11,27
      DO 120 j=1,5
      k(i,j)=0
      p(i,j)=0.
  120 v(i,j)=0.
      ns=20
      DO 130 j=1,5
      k(ns+1,j)=k(5,j)
      p(ns+1,j)=p(5,j)
      k(ns+3,j)=k(6,j)
      p(ns+3,j)=p(6,j)
      k(ns+5,j)=k(8,j)
      p(ns+5,j)=p(8,j)
      k(ns+6,j)=k(10,j)
  130 p(ns+6,j)=p(10,j)
C...Old standard continuation lines
      k(ns+2,1)=-1
      k(ns+2,3)=ns+1
      k(ns+4,1)=-1
      k(ns+4,3)=ns+3
      p(ns+4,3)=27
      p(ns+4,4)=27
C...Origin and colour info for incoming parton
      k(ns+3,1)=13
      k(ns+3,3)=2
      k(ns+3,4)=27
      k(ns+3,5)=27
C...Colour info for two outgoing partons
      k(ns+5,1)=3
      k(ns+6,1)=3
      IF(k(ns+6,2).EQ.21) THEN
C...qg-event
        IF(k(ns+5,2).GT.0) THEN
          k(ns+5,4)=(ns+6)*mstu(5)
          k(ns+5,5)=(ns+7)*mstu(5)
          k(ns+6,4)=(ns+7)*mstu(5)
          k(ns+6,5)=(ns+5)*mstu(5)
        ELSE
          k(ns+5,4)=(ns+7)*mstu(5)
          k(ns+5,5)=(ns+6)*mstu(5)
          k(ns+6,4)=(ns+5)*mstu(5)
          k(ns+6,5)=(ns+7)*mstu(5)
        ENDIF
      ELSE
C...qqbar-event
        k(ns+5,4)=(ns+7)*mstu(5)
        k(ns+5,5)=(ns+7)*mstu(5)
        k(ns+6,4)=(ns+7)*mstu(5)
        k(ns+6,5)=(ns+7)*mstu(5)
      ENDIF
C...Effective outgoing parton = sum of both outgoing partons
      k(ns+7,1)=14
      k(ns+7,3)=3
      IF(lst(24).EQ.2) THEN
        k(ns+7,2)=k(ns+5,2)
        IF(k(ns+7,2).EQ.21) WRITE(6,*) ' Warning: K(NS+7,2)=',k(ns+7,2)
        IF(k(ns+7,2).GT.0) THEN
          k(ns+7,4)=(ns+3)*mstu(5)+26
          k(ns+7,5)=(ns+3)*mstu(5)+25
        ELSE
          k(ns+7,4)=(ns+3)*mstu(5)+25
          k(ns+7,5)=(ns+3)*mstu(5)+26
        ENDIF
      ELSE
        k(ns+7,2)=21
        IF(k(ns+5,2).GT.0) THEN
          k(ns+7,4)=(ns+3)*mstu(5)+25
          k(ns+7,5)=(ns+3)*mstu(5)+26
        ELSE
          k(ns+7,4)=(ns+3)*mstu(5)+26
          k(ns+7,5)=(ns+3)*mstu(5)+25
        ENDIF
      ENDIF
      DO 140 j=1,4
  140 p(ns+7,j)=p(8,j)+p(10,j)

      it=ns+7
      IF(abs(p(it,1)).GT.0.1.OR.abs(p(it,2)).GT.0.1) THEN
C       WRITE(6,*) 'Warning: non-zero pt on final shower initiator'
C       WRITE(6,*) '1:',IT,K(IT,2),P(IT,1),P(IT,2),P(IT,3),P(IT,4),P(IT,5)
C       WRITE(6,*) '1:',8 ,K( 8,2),P( 8,1),P( 8,2),P( 8,3),P( 8,4),P( 8,5)
C       WRITE(6,*) '1:',10,K(10,2),P(10,1),P(10,2),P(10,3),P(10,4),P(10,5)
        lst(21)=12
        RETURN
      ENDIF
      p(it,1)=0.
      p(it,2)=0.

      p(ns+7,5)=sqrt(max(0.,p(ns+7,4)**2-p(ns+7,1)**2-p(ns+7,2)**2-
     &p(ns+7,3)**2))
      n=ns+7
C     CALL GULIST(103,2)

C...Scale for bremsstrahlung etc.
      q2=q2lp
      ipy(40)=10
      ipy(47)=n
C...Save quantities for later use.
      x(1)=1.
      x(2)=xr
      CALL lystfu(k(2,2),xr,q2lp,xpq)
      DO 160 ifl=-6,6
  160 xq(2,ifl)=xpq(ifl)
      IF(lst(23).EQ.1) THEN
        isub=39
        ipy(11)=1
      ELSEIF(lst(23).EQ.3) THEN
        isub=39
        ipy(11)=2
      ELSEIF(lst(23).EQ.4) THEN
        isub=39
        ipy(11)=3
      ELSEIF(lst(23).EQ.2) THEN
        isub=40
      ENDIF
      kfl(2,1)=k(5,2)
      kfl(2,2)=k(6,2)
      kfl(1,1)=kfl(2,1)
      kfl(1,2)=kfl(2,2)
      IF(isub.EQ.39) kfl(3,1)=k(1,2)
      IF(isub.EQ.40) kfl(3,1)=k(1,2)+isign(1,k(1,2))
      kfl(3,2)=k(27,2)
      pyvar(2)=parl(21)
      pyvar(1)=sqrt(pyvar(2))
      pyvar(3)=p(1,5)
      pyvar(4)=p(2,5)
      pyvar(5)=pyvar(1)/2.
      ipy(41)=k(1,2)
      ipy(42)=k(2,2)
      ipy(48)=0

C...Generate timelike parton shower (if required)
      IF(ipy(13).EQ.1) THEN
        CALL lscale(1,qmax)
        CALL lushow(25,26,qmax)
      ENDIF
      it=25
      IF(n.GE.27) it=27
      ns=n
C     CALL GULIST(104,2)

C...Generate spacelike parton shower (if required)
      ipu1=0
      ipu2=23
      IF(x(2)*(1.+(p(it,5)**2+pypar(22))/p(21,5)**2).GT.0.999) THEN
        lst(21)=13
        RETURN
      ENDIF
      IF(ipy(14).GE.1) THEN
        CALL lysspa(ipu1,ipu2)
        IF(lst(21).NE.0) RETURN
      ENDIF
      IF(n.EQ.ns) THEN 
        DO 220 i=ns+1,ns+4
        DO 220 j=1,5
        k(i,j)=0
        p(i,j)=0.
  220   v(i,j)=0.
        k(ns+1,1)=11
        k(ns+1,2)=kfl(2,1)
        k(ns+1,3)=21
        DO 230 j=1,5
  230   p(ns+1,j)=p(21,j)
        k(ns+2,1)=-1
        k(ns+2,3)=ns+1
        k(ns+3,1)=13
        k(ns+3,2)=kfl(2,2)
        k(ns+3,3)=23
        k(ns+3,4)=23
        k(ns+3,5)=23
        p(ns+3,3)=(p(it,5)**2+q2lp)*(p(21,4)-p(21,3))/(2.*q2lp)
        p(ns+3,4)=-p(ns+3,3)
        k(ns+4,1)=-1
        k(ns+4,3)=ns+3
        p(ns+4,3)=23
        p(ns+4,4)=23
        p(24,1)=ns+3
        p(24,2)=ns+3
        k(23,4)=k(23,4)+(ns+3)*mstu(5)
        k(23,5)=k(23,5)+(ns+3)*mstu(5)
        ipu1=0
        ipu2=ns+3
        n=n+4
      ENDIF
C     CALL GULIST(105,2)

C...Rotate and boost outgoing parton shower
      IF(n.GT.31) THEN
        k(n+1,1)=0
        DO 210 j=1,4
  210   p(n+1,j)=p(ns+1,j)+p(ns+3,j)
        IF(p(n+1,4).LE.1.01*p(it,5)) THEN
          lst(21)=14
          RETURN
        ENDIF
        robo(1)=ulangl(p(it,3),sqrt(p(it,1)**2+p(it,2)**2))
        robo(2)=ulangl(p(it,1),p(it,2))
      IF(abs(robo(1)).GT.0.001.OR.abs(robo(2)).GT.0.001) THEN
      WRITE(6,*) '1:',it,k(it,2),p(it,1),p(it,2),p(it,3),p(it,4),p(it,5)
      WRITE(6,*) '   ROBO(1-2)=',robo(1),robo(2)
      ENDIF
        CALL ludbrb(25,ns,0.,-robo(2),0.d0,0.d0,0.d0)
        CALL ludbrb(25,ns,-robo(1),0.,0.d0,0.d0,0.d0)
C...Replace old rotation method with x,y,z-boost to preserve QCD phi-dep
        deltap(1)=dble(p(n+1,1))
        deltap(2)=dble(p(n+1,2))
        deltap(3)=dble(p(n+1,3)) - dble(p(it,3))
        deltap(4)=sqrt(deltap(1)**2+deltap(2)**2+deltap(3)**2)
        IF(deltap(4).LT.1.d-11) goto 410
        dplong=-(dble(p(it,3))*deltap(3))/deltap(4)
        droot=max(0.d0,dble(p(n+1,4))**2-dble(p(it,4))**2+dplong**2)
        dbtot=-(dplong*dble(p(it,4))-dble(p(n+1,4))*sqrt(droot))/
C    &  SQRT(DBLE(P(N+1,4))**2-DBLE(P(IT,4))**2+DPLONG**2))/
     &  (dplong**2+dble(p(n+1,4))**2)
        dgamma=1.d0/sqrt(1.d0-dbtot**2)
        DO 400 i = 1,3
  400   drobo(i+2)=deltap(i)/(dgamma/(dgamma+1.d0)*
     &  (dble(p(n+1,4))-dgamma*dble(p(it,4)))+dgamma*dble(p(it,4)))
        CALL ludbrb(25,ns,0.,0.,drobo(3),drobo(4),drobo(5))
  410   CONTINUE
C...End phi-correction
      ENDIF
C     CALL GULIST(106,2)

      q2=q2lp
C...Hadron remnant and primordial kt
      ipy(47)=n
      CALL lyremn(ipu1,ipu2)
      IF(ipy(48).EQ.1) THEN
        lst(21)=15
        RETURN
      ENDIF
C     CALL GULIST(107,2)

C...Rearrange partons along strings
      mstu(24)=0
      CALL luprep(0)
      IF(mstu(24).NE.0) THEN
C       CALL GULIST(188,2)
        IF(lst(3).GE.1) WRITE(6,*) ' LUPREP error MSTU(24)= ',mstu(24)
        lst(21)=16
        RETURN
      ENDIF
C     CALL GULIST(109,2)

C...Clean up event record -> order:
C...1=inc. lepton; 2=inc. nucleon; 3=exch boson; 4=scat. lepton;
C...5=inc. parton before initial shower; 6=inc. parton at hard scattering
C...after shower; 7,8=first,second parton from hard scattering
C...before final shower
      lst(26)=7
      DO 510 j=1,5
      k(n+1,j)=k(4,j)
  510 p(n+1,j)=p(4,j)
      DO 520 j=1,5
      k(3,j)=k(5,j)
      p(3,j)=p(5,j)
      k(4,j)=k(9,j)
      p(4,j)=p(9,j)
      k(5,j)=k(n+1,j)
      p(5,j)=p(n+1,j)
      k(6,j)=k(ns+3,j)
      p(6,j)=p(ns+3,j)
C     K(7,J)=K(IT,J)
C     P(7,J)=P(IT,J)
      k(7,j)=k(25,j)
      p(7,j)=p(25,j)
      k(8,j)=k(26,j)
      p(8,j)=p(26,j)
  520 CONTINUE
      k(3,3)=1
      k(4,3)=1
      k(6,1)=21
      k(6,3)=5
      k(6,4)=0
      k(6,5)=0
      k(7,1)=21
      k(7,3)=6
      k(7,4)=0
      k(7,5)=0
      k(8,1)=21
      k(8,3)=6
      k(8,4)=0
      k(8,5)=0
C...Activate line with scattered lepton.
      k(4,1)=1
C...Deactivate obsolete lines 9, 10, 21, NS+1 (extra lines with boson)
      k(9,1)=0
      k(10,1)=0
      k(21,1)=0
      IF(k(ns+1,2).EQ.k(3,2)) k(ns+1,1)=0
C...Zero irrelevant lines with K(I,1)<0
      DO 540 i=1,n
      IF(k(i,1).LT.0) THEN
        DO 530 j=1,5
        k(i,j)=0
  530   p(i,j)=0.
      ENDIF
  540 CONTINUE
C     CALL GULIST(110,2)
C...Delete internal parton lines, i.e. with K(I,1)=13,14
      IF(mod(lst(4)/10,10).EQ.0) THEN
        CALL ltimex(t1)
        CALL luedit(14)
        CALL ltimex(t2)
C       CALL GULIST(111,2)
      ENDIF
C...Delete empty lines
      CALL ltimex(t1)
      CALL luedit(12)
      CALL ltimex(t2)
C     CALL GULIST(112,2)

      RETURN
      END