leptox.F

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

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

      SUBROUTINE leptox

      IMPLICIT NONE

C...Select process and choose kinematical variables (x,y; x,Q2; x,W2)
C...according to the differential cross section.

      COMMON /lintrl/ psave(3,4,5),ksave(4),xmin,xmax,ymin,ymax,
     &q2min,q2max,w2min,w2max,ilep,inu,ig,iz
      REAL psave,xmin,xmax,ymin,ymax,q2min,q2max,w2min,w2max
      INTEGER ksave,ilep,inu,ig,iz
      SAVE /lintrl/

*
* 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 /loptim/ optx(4),opty(4),optq2(4),optw2(4),comfac
      REAL optx,opty,optq2,optw2,comfac
      SAVE /loptim/

* this Common block appears only in S LEPTOX 
      COMMON /flinfo/ rflq,rflg,rflm,rflt
      REAL rflq,rflg,rflm,rflt
      SAVE /flinfo/ 

      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/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
      INTEGER  kchg
      REAL  pmas,parf,vckm
      SAVE /ludat2/


      INTEGER ih,i,ncut,ii
      INTEGER lqcd,ltm,lht,lkinem
      REAL w2low,w2upp,ylow,yupp,q2low,q2upp,pnt,pqh,s,pm2,
     +xpq,gfq2,aemcor,ulalem,zlep,a,b,flg,flq,
     +flm,flt,f2em,pqh17,weight,hx,xfact,q2fact,hy,yfact,hq2,
     +w2fact,sigl,sigr,sigma,viol,yq,yqb,fyq,hw2,which
      REAL rlu
*...Added array XDPQ to store delta parton distributions
      dimension pqh(17,2),pnt(2,2),xpq(-6:6)
*---
CGI-001122...d- & u-quark weights for p & n target & helicity state
      REAL pntdu(2,2,2)
CMM-010306...anti d- & anti u-quark weights for p & n target & helicity state
      REAL pntdus(2,2,2)     
      DOUBLE PRECISION dari27,dari28
      DATA dari27,dari28/2*0.d0/
      DATA w2low,w2upp,ylow,yupp,q2low,q2upp/6*0./

*PEPSI>>
      IF (lst(40).NE.0) THEN
        CALL polleptox
        RETURN
      ENDIF 
*PEPSI<<
      DO 10 ih=1,2
      DO 5 i=1,2
CGI-001122
      pntdu(i,1,ih)=0.
      pntdu(i,2,ih)=0.
CMM-010306
      pntdus(i,1,ih)=0.
      pntdus(i,2,ih)=0.      
    5 pnt(i,ih)=0.
      DO 6 i=1,8
      ewqc(1,ih,i)=0.
    6 ewqc(2,ih,i)=0.
      DO 10 i=1,17
   10 pqh(i,ih)=0.
      DO 20 i=1,17
   20 pq(i)=0.

      lst(21)=0
      ncut=0
      s=parl(21)
      pm2=psave(3,2,5)**2
      IF(lst(2).NE.1) THEN
       q2low=max(q2min,x*ymin*s,(w2min-pm2)*x/max(1.-x,1.e-22))
       q2upp=min(q2max,x*ymax*s,(w2max-pm2)*x/max(1.-x,1.e-22))
       ylow=max(ymin,q2min/max(s*x,1.e-22),
     & (w2min-pm2)/max(s*(1.-x),1.e-22))
       yupp=min(ymax,q2max/max(s*x,1.e-22),
     & (w2max-pm2)/max(s*(1.-x),1.e-22))
       w2low=max(w2min,(1.-x)*ymin*s+pm2,q2min*(1.-x)/max(x,1.e-22)+pm2)
       w2upp=min(w2max,(1.-x)*ymax*s+pm2,q2max*(1.-x)/max(x,1.e-22)+pm2)
       goto 110
      ENDIF

      IF(pari(28).LT.0.5) THEN
C...For first call, reset double precision counters.
        dari27=0.d0
        dari28=0.d0
      ENDIF
  100 dari28=dari28+1.d0
      pari(28)=dari28
  101 CONTINUE
C...Choose x according to the distribution
C...hx(x) =  a + b/x + c/x**2 + d/x**3. In detail
C...hq=OPTX(1)/(XMAX-XMIN) + 1/ln(XMAX/XMIN)*OPTX(2)/X
C...   +XMIN*XMAX/(XMAX-XMIN)*OPTX(3)/X**2
C...   +2*(XMIN*XMAX)**2/(XMAX**2-XMIN**2)*OPTX(4)/X**3
      which=(optx(1)+optx(2)+optx(3)+optx(4))*rlu(0)
      IF(which.LE.optx(1)) THEN
        x=xmin+rlu(0)*(xmax-xmin)
      ELSEIF(which.LE.(optx(1)+optx(2))) THEN
        x=xmin*(xmax/xmin)**rlu(0)
      ELSEIF(which.LE.(optx(1)+optx(2)+optx(3))) THEN
        x=xmin*xmax/(xmax+rlu(0)*(xmin-xmax))
      ELSE
        x=sqrt((xmin*xmax)**2/(xmax**2+rlu(0)*(xmin**2-xmax**2)))
      ENDIF
      IF(lst(31).EQ.1) THEN
C...Choose Q**2 according to the distribution
C...hq(Q2) =  a + b/(Q2) + c/(Q2)**2 + d/(Q2)**3. In detail
C...hq=OPTQ2(1)/(Q2MAX-Q2MIN) + 1/ln(Q2MAX/Q2MIN)*OPTQ2(2)/Q2
C...   +Q2MIN*Q2MAX/(Q2MAX-Q2MIN)*OPTQ2(3)/Q2**2
C...   +2*(Q2MIN*Q2MAX)**2/(Q2MAX**2-Q2MIN**2)*OPTQ2(4)/Q2**3
        q2low=max(q2min,x*ymin*s,(w2min-pm2)*x/(1.-x))
        q2upp=min(q2max,x*ymax*s,(w2max-pm2)*x/(1.-x))
        IF(q2upp.LT.q2low) goto 101
        which=(optq2(1)+optq2(2)+optq2(3)+optq2(4))*rlu(0)
        IF(which.LE.optq2(1)) THEN
          q2=q2low+rlu(0)*(q2upp-q2low)
        ELSEIF(which.LE.(optq2(1)+optq2(2))) THEN
          q2=q2low*(q2upp/q2low)**rlu(0)
        ELSEIF(which.LE.(optq2(1)+optq2(2)+optq2(3))) THEN
          q2=q2low*q2upp/(q2upp+rlu(0)*(q2low-q2upp))
        ELSE
         q2=sqrt((q2low*q2upp)**2/(q2upp**2+rlu(0)*(q2low**2-q2upp**2)))
        ENDIF
        y=q2/(parl(21)*x)
        IF(y.LT.ymin.OR.y.GT.ymax) goto 100
      ELSEIF(lst(31).EQ.2) THEN
C...Choose y according to the distribution
C...hy(y) =  a + b/y + c/y**2 + d/y**3. In detail
C...hy=OPTY(1)/(YMAX-YMIN) + 1/ln(YMAX/YMIN)*OPTY(2)/Y
C...   +YMIN*YMAX/(YMAX-YMIN)*OPTY(3)/Y**2
C...   +2*(YMIN*YMAX)**2/(YMAX**2-YMIN**2)*OPTY(4)/Y**3
        ylow=max(ymin,q2min/(s*x),(w2min-pm2)/(s*(1.-x)))
        yupp=min(ymax,q2max/(s*x),(w2max-pm2)/(s*(1.-x)))
        IF(yupp.LT.ylow) goto 101
        which=(opty(1)+opty(2)+opty(3)+opty(4))*rlu(0)
        IF(which.LE.opty(1)) THEN
          y=ylow+rlu(0)*(yupp-ylow)
        ELSEIF(which.LE.(opty(1)+opty(2))) THEN
          y=ylow*(yupp/ylow)**rlu(0)
        ELSEIF(which.LE.(opty(1)+opty(2)+opty(3))) THEN
          y=ylow*yupp/(yupp+rlu(0)*(yupp-ylow))
        ELSE
          y=sqrt((ylow*yupp)**2/(yupp**2+rlu(0)*(ylow**2-yupp**2)))
        ENDIF
        q2=x*y*parl(21)
        IF(q2.LT.q2min.OR.q2.GT.q2max) goto 100
      ELSEIF(lst(31).EQ.3) THEN
C...Choose W**2 according to the distribution
C...hw(W2) =  a + b/(W2) + c/(W2)**2 + d/(W2)**3. In detail
C...hw=OPTW2(1)/(W2MAX-W2MIN) + 1/ln(W2MAX/W2MIN)*OPTW2(2)/W2
C...   +W2MIN*W2MAX/(W2MAX-W2MIN)*OPTW2(3)/W2**2
C...   +2*(W2MIN*W2MAX)**2/(W2MAX**2-W2MIN**2)*OPTW2(4)/W2**3
        w2low=max(w2min,(1.-x)*ymin*s+pm2,q2min*(1.-x)/x+pm2)
        w2upp=min(w2max,(1.-x)*ymax*s+pm2,q2max*(1.-x)/x+pm2)
        IF(w2upp.LT.w2low) goto 101
        which=(optw2(1)+optw2(2)+optw2(3)+optw2(4))*rlu(0)
        IF(which.LE.optw2(1)) THEN
          w2=w2low+rlu(0)*(w2upp-w2low)
        ELSEIF(which.LE.(optw2(1)+optw2(2))) THEN
          w2=w2low*(w2upp/w2low)**rlu(0)
        ELSEIF(which.LE.(optw2(1)+optw2(2)+optw2(3))) THEN
          w2=w2low*w2upp/(w2upp+rlu(0)*(w2low-w2upp))
        ELSE
         w2=sqrt((w2low*w2upp)**2/(w2upp**2+rlu(0)*(w2low**2-w2upp**2)))
        ENDIF
        y=(w2-pm2)/((1.-x)*parl(21))
        q2=x*y*parl(21)
        IF(y.LT.ymin.OR.y.GT.ymax) goto 100
        IF(q2.LT.q2min.OR.q2.GT.q2max) goto 100
      ENDIF

  110 IF(lkinem(lst(2)).NE.0) THEN
        ncut=ncut+1
        IF(lst(2).EQ.1) THEN
          IF(ncut.LE.9999) goto 100
          IF(lst(3).GE.1) WRITE(6,1200)
        ENDIF
        lst(21)=2
        RETURN
      ENDIF
      pari(24)=(1.+(1.-y)**2)/2.
      pari(25)=1.-y
      pari(26)=(1.-(1.-y)**2)/2.
      CALL lnstrf(x,q2,xpq)
C...Lepton helicity state, only one contributes in some cases.
      ih=1
      IF(parl(6).GT.+0.99) ih=2
  200 lst(30)=sign(1.,ih-1.5)
      pqh(17,ih)=0.
      pnt(1,ih)=0.
      pnt(2,ih)=0.
CGI-001122
      pntdu(1,1,ih)=0.
      pntdu(1,2,ih)=0.
      pntdu(2,1,ih)=0.
      pntdu(2,2,ih)=0.
CMM-010306
      pntdus(1,1,ih)=0.
      pntdus(1,2,ih)=0.
      pntdus(2,1,ih)=0.
      pntdus(2,2,ih)=0.
      IF(lst(23).EQ.2) THEN
C...Charged current: zero cross-section for one helicity state.
        IF(ksave(1).LT.0.AND.ih.EQ.1
     &  .OR.ksave(1).GT.0.AND.ih.EQ.2) goto 240
        yq=pari(24)-lst(30)*pari(26)
        yqb=pari(24)+lst(30)*pari(26)
        IF(pari(11).GT.1.e-06) THEN
          IF(k(3,2).LT.0) THEN
            pnt(1,ih)=(1.-pari(11))*pari(13)*yq
            pnt(2,ih)=pari(11)*pari(12)*yq
CGI-001122...Save u-quark weights for p & n target
            pntdu(1,2,ih)=pari(13)*yq
            pntdu(2,2,ih)=pari(12)*yq
CMM-010306...Save anti u-quark weights for p & n target
            pntdus(1,2,ih)=pari(43)*yq
            pntdus(2,2,ih)=pari(42)*yq
          ELSE
            pnt(1,ih)=(1.-pari(11))*pari(12)*yq
            pnt(2,ih)=pari(11)*pari(13)*yq
CGI-001122...Save d-quark weights for p & n target
            pntdu(1,1,ih)=pari(12)*yq
            pntdu(2,1,ih)=pari(13)*yq
CMM-010306...Save anti d-quark weights for p & n target
            pntdus(1,1,ih)=pari(42)*yq
            pntdus(2,1,ih)=pari(43)*yq
           ENDIF
        ENDIF
        DO 220 i=1,lst(12)
        IF(k(3,2)*qc(i).LT.0) THEN
          pqh(i,ih)=xpq(i)*yq
        ELSE
          pqh(i+lst(12),ih)=xpq(-i)*yqb
        ENDIF
  220   CONTINUE
      ELSE
C...Neutral current: electromagnetic or weak or both with interference.
        gfq2=q2/(pmas(23,1)**2+q2)*sqrt(2.)*parl(17)*pmas(23,1)**2/
     &  (3.1415927*parl(16))
C...Correction to obtain Q**2 dependent alpha-em, if desired.
        aemcor=1.
        IF(lst(18).GE.2) aemcor=ulalem(q2)/parl(16)
        ii=3-ih
        zlep=zl(ih,ilep+2*inu)
        DO 230 i=1,max(lst(12),lst(13))
        a=(-ig*qc(i)*aemcor+iz*gfq2*zlep*zq(ih,i))**2
        b=(-ig*qc(i)*aemcor+iz*gfq2*zlep*zq(ii,i))**2
C...Save helicity-dependent electroweak quark couplings for later use.
        ewqc(1,ih,i)=a
        ewqc(2,ih,i)=b
        IF(i.GT.lst(12)) goto 230
        fyq=(a+b)*pari(24)+(a-b)*pari(26)
        pqh(i,ih)=xpq(i)*fyq
        IF(i.LE.2.AND.pari(11).GT.1.e-06) THEN
          pnt(1,ih)=pnt(1,ih)+(1.-pari(11))*pari(11+i)*fyq
          pnt(2,ih)=pnt(2,ih)+pari(11)*pari(14-i)*fyq
CGI-001122...Save d- & u-quark weights for p & n target
          pntdu(1,i,ih)=pari(11+i)*fyq
          pntdu(2,i,ih)=pari(14-i)*fyq
CMM-010306...Save anti d- & u-quark weights for p & n target
          pntdus(1,i,ih)=pari(41+i)*fyq
          pntdus(2,i,ih)=pari(44-i)*fyq
        ENDIF
        pqh(i+lst(12),ih)=xpq(-i)*((a+b)*pari(24)-(a-b)*pari(26))
  230   CONTINUE
      ENDIF
  240 CONTINUE
      DO 300 i=1,lst(12)
  300 pqh(17,ih)=pqh(17,ih)+pqh(i,ih)+pqh(i+lst(12),ih)

      IF(abs(parl(6)).LT.0.99.AND.ih.EQ.1) THEN
        ih=2
        goto 200
      ENDIF

      flq=0.
      flg=0.
      flm=0.
      flt=0.
      IF(lst(11).NE.0.AND.(lst(23).EQ.1.OR.lst(23).EQ.4)
     &.AND.lst(2).NE.-3) THEN
C...Include F_L for photon exchange (unless QCD grid being set up)
        lqcd=mod(lst(11),10)
        ltm=mod(lst(11)/10,10)
        lht=lst(11)/100
C...Include QCD, target mass and/or higher twist contr. to long. str fcn
C...FL from interpolation.
        IF(lqcd.EQ.1.OR.ltm.EQ.1) CALL flipol(flq,flg,flm)
C...Event simulation: if requested, get FL by event-by-event integration
        IF(lst(2).GT.0.AND.
     &  (lqcd.EQ.2.OR.ltm.EQ.2)) CALL flintg(flq,flg,flm)
        IF(ltm.GE.1.OR.lht.GE.1) THEN
          f2em=0.
          DO 301 i=1,lst(12)
  301     f2em=f2em+qc(i)**2*(xpq(i)+xpq(-i))
          IF(ltm.GE.1) flm=flm-2.*x**2*psave(3,2,5)**2/q2*f2em
          IF(lht.GE.1) flt=8.*parl(19)/q2*f2em
        ENDIF
        DO 305 ih=1,2
        pqh17=pqh(17,ih)
C...Note factor 2 at the end, since PQH(IH,17) contains overall factor 2
        pqh(17,ih)=pqh(17,ih)-y**2*(flq+flg+flm+flt)
        DO 305 i=1,16
  305   pqh(i,ih)=pqh(i,ih)*pqh(17,ih)/pqh17
      ENDIF

      DO 310 i=1,17
  310 pq(i)=(1.-parl(6))/2.*pqh(i,1)+(1.+parl(6))/2.*pqh(i,2)

C...Relative contribution from longitudinal str. fcn. and higher twist.
       rflq=-y**2*flq/max(pq(17),1.e-33)
       rflg=-y**2*flg/max(pq(17),1.e-33)
       rflm=-y**2*flm/max(pq(17),1.e-33)
       rflt=-y**2*flt/max(pq(17),1.e-33)

C...Common factor for matrix elements.
      IF(lst(31).EQ.1) THEN
        IF(lst(23).EQ.2) THEN
          comfac=1./x/(1.+q2/pmas(24,1)**2)**2
        ELSE
          comfac=1./x/q2**2
        ENDIF
      ELSEIF(lst(31).EQ.2) THEN
        IF(lst(23).EQ.2) THEN
          comfac=1./(1.+q2/pmas(24,1)**2)**2*parl(21)
        ELSE
          comfac=1./q2**2*parl(21)
        ENDIF
      ELSEIF(lst(31).EQ.3) THEN
        IF(lst(23).EQ.2) THEN
          comfac=1./x/(1.+q2/pmas(24,1)**2)**2  * x/(1.-x)
        ELSE
          comfac=1./x/q2**2 * x/(1.-x)
        ENDIF
      ENDIF
C-check: Move change of COMFAC to below??
C...Prepare for Q2 weighting.
C     WEIGHT=1/Q2**2
      weight=1.d0
      comfac=comfac/weight
      IF(lst(2).LE.-2) RETURN
      hx=optx(1)/(xmax-xmin) + 1./alog(xmax/xmin)*optx(2)/x
     &+xmin*xmax/(xmax-xmin)*optx(3)/x**2
     &+2*(xmin*xmax)**2/(xmax**2-xmin**2)*optx(4)/x**3
      xfact=optx(1)+optx(2)+optx(3)+optx(4)
      IF(lst(31).EQ.1) THEN
        hq2=optq2(1)/(q2upp-q2low)
     &  +1./alog(q2upp/q2low)*optq2(2)/q2
     &  +q2low*q2upp/(q2upp-q2low)*optq2(3)/q2**2
     &  +2*(q2low*q2upp)**2/(q2upp**2-q2low**2)*optq2(4)/q2**3
        q2fact=optq2(1)+optq2(2)+optq2(3)+optq2(4)
        comfac=comfac*xfact*q2fact/hx/hq2
      ELSEIF(lst(31).EQ.2) THEN
        hy=opty(1)/(yupp-ylow)+1./alog(yupp/ylow)*opty(2)/y
     &  +ylow*yupp/(yupp-ylow)*opty(3)/y**2
     &  +2*(ylow*yupp)**2/(yupp**2-ylow**2)*opty(4)/y**3
        yfact=opty(1)+opty(2)+opty(3)+opty(4)
        comfac=comfac*xfact*yfact/hx/hy
      ELSEIF(lst(31).EQ.3) THEN
        hw2=optw2(1)/(w2upp-w2low)
     &  +1./alog(w2upp/w2low)*optw2(2)/w2
     &  +w2low*w2upp/(w2upp-w2low)*optw2(3)/w2**2
     &  +2*(w2low*w2upp)**2/(w2upp**2-w2low**2)*optw2(4)/w2**3
        w2fact=optw2(1)+optw2(2)+optw2(3)+optw2(4)
        comfac=comfac*xfact*w2fact/hx/hw2
      ENDIF
      IF(lst(2).LE.0) RETURN

C-check: Move change of COMFAC to here?
      sigl=(1.-parl(6))/2.*pqh(17,1)
      sigr=(1.+parl(6))/2.*pqh(17,2)
      sigma=sigl+sigr
      IF(lst(2).EQ.1) THEN
C...When chosing (x,y), reject according to maximum of "cross-section",
C...update cross-section estimate.
        dari27=dari27+dble(sigma)*dble(comfac)*weight
        pari(27)=dari27
        viol=sigma*comfac/pari(lst(23))
        IF(viol.GT.pari(32)) THEN
          pari(32)=viol
          IF(pari(32).GT.1.) THEN
            pari(lst(23))=pari(lst(23))*pari(32)
            IF(lst(3).GE.1) WRITE(6,1300) pari(32),int(pari(30)+1),
     &      pari(lst(23)),x,y,q2,w2
            pari(32)=1.
          ENDIF
        ENDIF
        IF(viol.LT.rlu(0)) goto 100
        parl(24)=pari(31)*dari27/dari28
      ENDIF

      IF(abs(parl(6)).LT.0.99) THEN
C...Choose helicity of incoming lepton.
        ih=1
        IF(rlu(0)*sigma.GT.sigl) ih=2
      ENDIF
      lst(30)=sign(1.,ih-1.5)
CGI-001122...Save parton density weights for chosen helicity
      DO 510 i=1,17
  510 pq(i)=pqh(i,ih)

C...Choose target nucleon, proton or neutron.
      lst(22)=1
      k(2,2)=2212
      IF(pari(11).GT.1.e-06) THEN
        IF(rlu(0).LT.(pari(11)*(pqh(17,ih)-pnt(1,ih)-pnt(2,ih))+
     &  pnt(2,ih))/pqh(17,ih)) THEN
          lst(22)=2
          k(2,2)=2112
        ENDIF
CGI-001122...Save parton density weights for chosen target nucleon
CMM-010306...and take also care of the anti u- and d-quarks.
        pq(17)=pq(17)-pq(1)-pq(2)-pq(1+lst(12))-pq(2+lst(12))
        pq(1)=pntdu(lst(22),1,ih)
        pq(2)=pntdu(lst(22),2,ih)
        pq(1+lst(12))=pntdus(lst(22),1,ih)
        pq(2+lst(12))=pntdus(lst(22),2,ih)      
        pq(17)=pq(17)+pq(1)+pq(2)+pq(1+lst(12))+pq(2+lst(12))
      ENDIF

      RETURN
 1200 FORMAT(' Warning: LEPTOX is looping, cannot find allowed ',
     &'phase space point due to cuts,',/,
     &10x,'check, in particular, CUT(11) to CUT(14)')
 1300 FORMAT(' Warning: maximum violated by a factor ',f7.3,
     &' in event ',i7,/,' maximum increased by this factor to ',e12.3,
     &/,' Point of violation: x, y, Q**2, W**2 = ',4g10.3)
      END