lxp.F

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

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

      SUBROUTINE lxp(XP,IFAIL)

      IMPLICIT NONE

C...Choose value of XP according to QCD matrix elements weighted by
C...structure functions.

*
* 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/

      COMMON /lgrid/ nxx,nww,xx(31),ww(21),pqg(31,21,3),pqqb(31,21,2), 
     &qgmax(31,21,3),qqbmax(31,21,2),ycut(31,21),xtot(31,21),np
      REAL    xx,ww,pqg,pqqb,qgmax,qqbmax,ycut,xtot
      INTEGER nxx,nww,np
      SAVE /lgrid/


      common/debug1/iw,ix,ip
      integer iw,ix,ip
      
      INTEGER ifail,loop
      REAL xp,xpmin,xpmax,ap,bp,qxpmax,xpweit
      REAL rlu,dsigma
      DOUBLE PRECISION dxpmax

      ifail=1
CAE.Scheme for ME cutoff: W2, Q2, mixed, z-shat
      IF(lst(20).LE.1) THEN
        xpmin=dble(x)/(1.d0-2.d0*(1.d0-dble(x))*dble(parl(27)))
        dxpmax=dble(x)/(dble(x)+(1.d0-dble(x))*dble(parl(27)))
        xpmax=sngl(dxpmax)
      ELSEIF(lst(20).EQ.2) THEN
        xpmin=dble(x)/(1.d0-2.d0*dble(x)*dble(parl(27)))
        dxpmax=1.d0/(1.d0+dble(parl(27)))
        xpmax=sngl(dxpmax)
      ELSEIF(lst(20).EQ.3.OR.lst(20).EQ.4) THEN
        xpmin=x
        dxpmax=1.d0/(1.d0+dble(parl(9)))
        xpmax=dxpmax
      ELSEIF(lst(20).EQ.5) THEN
        xpmin=x
        dxpmax=dble(q2)/(dble(q2)+dble(parl(9)))
        xpmax=dxpmax
      ELSEIF(lst(20).EQ.6) THEN
        xpmin=x
        dxpmax=dble(q2)/(dble(q2)+dble(parl(27)))
        xpmax=dxpmax
      ELSE
        WRITE(6,*) 'LXP: No such jet scheme!'
      ENDIF
CAE
      IF(xpmin.LT.x.OR.xpmin.GT.1.) RETURN
      IF(xpmin.GE.xpmax) RETURN
      ap=1.-xpmin
      bp=(1.d0-dxpmax)/ap
      IF(lst(24).EQ.2) THEN
        qxpmax=pari(15)
        IF(lst(17).NE.0.AND.lst(19).GE.0) qxpmax=
     &  pari(24)*pari(15)+pari(25)*pari(16)+pari(26)*pari(17)
      ELSE
        qxpmax=pari(18)
        IF(lst(17).NE.0.AND.lst(19).GE.0) qxpmax=
     &  pari(24)*pari(18)+pari(25)*pari(19)
      ENDIF

C...Safety factor on max value.
      qxpmax=qxpmax*1.05
      loop=0
  100 loop=loop+1
      IF(loop.GT.1000) RETURN
      xp=1.-ap*bp**rlu(0)
      xpweit=dsigma(xp)/qxpmax
      IF(xpweit.LT.rlu(0)) goto 100
      ifail=0
      RETURN
      END