pollzp.F

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

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

      SUBROUTINE pollzp(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)
*
*
* 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,xi
      REAL rlu,dzpmax
      REAL d,e,f,dp,poldqcd,dqcd,zpweit
      REAL dz,dx,da,db,dc,dd,de,df
      INTEGER ih,i,ifl
      REAL zp,c1,c2,szp,fqg,fqqp,zpmin,zpmax,a,b,
     +csign,c,xpq,xdpq,cp,aa,bb,z1
*HI      REAL FQQ
      dimension xpq(-6:6),xdpq(-6:6)
*---
      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)
*HI      FQQ(DZ,DX,DA,DB,DC,DD,DE)=DA*DD*(DZ**2+(1.-DZ)**2)+DB*DE*DZ*
*HI     &(1.-DZ)+DC*DD*(2.*DZ-1.)
*...Polarized case:
      fqqp(dz,dx,da,db,dc,dd,de,df)=da*dd*(dz**2+(1.-dz)**2)+db*de*dz*
     &(1.-dz)+dc*df*((1.-dz)**2+dz**2)
*---
      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.OR.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
         
         xi=x/xp
         CALL lnstrf(xi,q2,xpq)
         CALL dnstrf(xi,q2,xdpq)
         
         IF (lst(24).EQ.2) THEN
            ifl=lst(25)
            a=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(24)*xpq(ifl)
            b=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(25)*xpq(ifl)
            c=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(26)*lst(40)*xdpq(ifl)
         ELSE IF (lst(24).EQ.3) THEN
            a=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(24)*xpq(0)
            b=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(25)*xpq(0)
            c=(ewqc(1,ih,i)+ewqc(2,ih,i))*pari(26)*lst(40)*xdpq(0)
         ENDIF
         csign=c 

      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
         e=8.*xp*(1-xp)
         d=xp**2+(1.-xp)**2
         f=(2.*xp-1.)
         dzpmax=max(fqqp(zpmin,xp,a,b,csign,d,e,f),
     &            fqqp(zpmax,xp,a,b,csign,d,e,f))
         aa=4.*(a*d+csign*f)-2.*b*e
         IF(abs(aa).GT.1.e-20) THEN
            bb=b*e-2.*(a*d+csign*f)
            z1=-bb/aa
            IF(z1.GT.zpmin.AND.z1.LT.zpmax) THEN
               dzpmax=max(dzpmax,fqqp(z1,xp,a,b,csign,d,e,f))
            ENDIF
         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*poldqcd(0,ipart,3,xp,zp,0.))/dzpmax
      IF(zpweit.LT.rlu(0)) goto 100
      ifail=0
      RETURN
      END