lzp.F

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

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

      SUBROUTINE lzp(XP,ZP,IFAIL)

      IMPLICIT NONE

C...Choose value of ZP according to QCD matrix elements.
*
* 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 ifail,loop,ipart
      REAL xp,ap,bp
      REAL rlu,dzpmax
      REAL d,e,dp,dqcd,zpweit
      REAL dz,dx,da,db,dc,dd,de
      INTEGER ih,i
      REAL zp,c1,c2,szp,fqg,fqq,zpmin,zpmax,a,b,
     +csign,c,cp,aa,bb,z1

      DATA c1,c2/0.2122066,0.0795775/,dzpmax,szp,cp/3*0./
      fqg(dz,dx,da,db,dc)=da*(dz**2+dx**2)/(1.-dx)+2.*da*dx*dz*(1.-dz)
     &+2.*da*(1.-dz)+4.*db*dx*dz*(1.-dz)+dc*(dz**2+dx**2)/(1.-dx)+
     &2.*dc*(dx+dz)*(1.-dz)
      fqq(dz,dx,da,db,dc,dd,de)=da*dd*(dz**2+(1.-dz)**2)+db*de*dz*
     &(1.-dz)+dc*dd*(2.*dz-1.)

*PEPSI>>
      IF(lst(40).NE.0) THEN
       CALL pollzp(xp,zp,ifail)
       RETURN
      ENDIF
*PEPSI<<

      ifail=1
      ih=1
      IF(lst(30).EQ.1) ih=2
CAE.Scheme for ME cutoff: W2, Q2, mixed, z-shat
      IF(lst(20).LE.1) THEN
        zpmin=(1.-x)*xp/(xp-x)*parl(27)
      ELSEIF(lst(20).EQ.2) THEN
        zpmin=x*xp/(xp-x)*parl(27)
      ELSEIF(lst(20).GE.3.AND.lst(20).LE.5) THEN
        zpmin=parl(27)
      ELSEIF(lst(20).EQ.6) THEN
        zpmin=parl(8)
      ELSE
        WRITE(6,*) 'LZP: No such jet scheme!'
      ENDIF
CAE
      IF(zpmin.LE.0..OR.zpmin.GE.0.5) RETURN
      zpmax=1.-zpmin
      i=iabs(lst(25))
      ap=1.-zpmin
      bp=zpmin/ap
      IF(lst(23).EQ.2) THEN
        a=pari(24)
        b=pari(25)
        csign=-lst(30)*isign(1,lst(25))*pari(26)
      ELSE
        a=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(24)
        b=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(25)
        c=(ewqc(1,ih,i)-ewqc(2,ih,i))*pari(26)
        csign=-c*lst(30)*isign(1,lst(25))
      ENDIF
      IF(lst(24).EQ.2) THEN
        dzpmax=max(fqg(zpmin,xp,a,b,csign),fqg(zpmax,xp,a,b,csign))
        aa=2.*(a+csign)/(1.-xp)-4.*a*xp-8.*b*xp-4.*csign
        IF(abs(aa).GT.1.e-20) THEN
          bb=2.*a*(xp-1.)+4.*b*xp+2.*csign*(1.-xp)
          z1=-bb/aa
          IF(z1.GT.zpmin.AND.z1.LT.zpmax)
     &    dzpmax=max(dzpmax,fqg(z1,xp,a,b,csign))
        ENDIF
        dzpmax=dzpmax*c1*1.05
      ELSEIF(lst(24).EQ.3) THEN
        cp=1./bp**2
        d=xp**2+(1.-xp)**2
        e=8.*xp*(1-xp)
        dzpmax=max(fqq(zpmin,xp,a,b,csign,d,e),
     &  fqq(zpmax,xp,a,b,csign,d,e))
        aa=4.*a*d-2.*b*e
        IF(abs(aa).GT.1.e-20) THEN
          bb=b*e-2.*a*d+2.*csign*d
          z1=-bb/aa
          IF(z1.GT.zpmin.AND.z1.LT.zpmax)
     &    dzpmax=max(dzpmax,fqq(z1,xp,a,b,csign,d,e))
        ENDIF
        dzpmax=dzpmax*c2*1.05
      ENDIF
      ipart=lst(24)-1
      loop=0
  100 loop=loop+1
      IF(loop.GT.1000) RETURN
      IF(lst(24).EQ.2) THEN
        zp=1.-ap*bp**rlu(0)
        szp=1.-zp
      ELSEIF(lst(24).EQ.3) THEN
        dp=bp*cp**rlu(0)
        zp=dp/(1.+dp)
        szp=zp*(1.-zp)
      ENDIF
      zpweit=szp*(a*dqcd(0,ipart,1,xp,zp,0.)+b*dqcd(0,ipart,2,xp,zp,0.)
     &+csign*dqcd(0,ipart,3,xp,zp,0.))/dzpmax
      IF(zpweit.LT.rlu(0)) goto 100
      ifail=0
      RETURN
      END