dsigma.F

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

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

      FUNCTION dsigma(XP)
      
      IMPLICIT NONE
      
C...Differential cross section for first order QCD processes.

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

      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/

      
      REAL dsigma,xp,dqcdi,ulmass,poldsigma
      INTEGER i,ih,il,iu,ip
      REAL pqh,amu,xi,zpmin,zpmax,xpq,wq,wqb,tq,tqb,t1,s13,sgn,sig
      dimension xpq(-6:6),pqh(17,2)

*PEPSI>>
      IF (lst(40).NE.0) THEN
         dsigma = poldsigma(xp)
         RETURN
      ENDIF
*PEPSI<<


      dsigma=0.
      DO 10 i=1,17
      pqh(i,1)=0.
      pqh(i,2)=0.
   10 pq(i)=0.

      mstj(93)=1
      amu=ulmass(2)
      IF(lst(32).EQ.1.AND.lst(19).GE.0.AND.lst(17).EQ.1) THEN
        il=lst(32)
        iu=lst(32)
      ELSE
        il=1
        iu=3
        IF(lst(23).EQ.1.OR.lst(24).EQ.3) iu=2
      ENDIF
      xi=x/xp
C...Scheme for ME cutoff: W2, Q2, mixed
      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).GE.6) THEN
        zpmin=parl(8)
      ENDIF

      IF(zpmin.LE.0..OR.zpmin.GE.0.5) RETURN
      zpmax=1.d0-dble(zpmin)
      CALL lnstrf(xi,q2,xpq)
      IF(lst(24).EQ.3) goto 3000

C...Gluon bremsstrahlung process, i.e. qg-event.
 2000 DO 2500 ip=il,iu
      sig=dqcdi(1,ip,xp,zpmin,zpmax)
      sgn=sign(1.,5.-2.*ip)
      DO 2300 ih=1,2
      IF(ih.EQ.1) THEN
        IF(parl(6).GT.0.99) goto 2300
        IF(lst(32).EQ.0.AND.lst(30).NE.-1) goto 2300
      ELSEIF(ih.EQ.2) THEN
        IF(parl(6).LT.-0.99) goto 2300
        IF(lst(32).EQ.0.AND.lst(30).NE.1) goto 2300
      ENDIF
      IF(lst(32).NE.0) lst(30)=sign(1.,ih-1.5)
      IF(lst(23).NE.2) THEN
        DO 2100 i=1,lst(12)
        wq=xpq(i)*sig*(ewqc(1,ih,i)+sgn*ewqc(2,ih,i))
        wqb=xpq(-i)*sig*sgn*(ewqc(1,ih,i)+sgn*ewqc(2,ih,i))
C...Include y-dependence.
        wq=wq*pari(23+ip)
        wqb=wqb*pari(23+ip)
        pqh(i,ih)=pqh(i,ih)+wq
        pqh(i+lst(12),ih)=pqh(i+lst(12),ih)+wqb
        pqh(17,ih)=pqh(17,ih)+wq+wqb
 2100   CONTINUE
      ELSEIF(lst(23).EQ.2) THEN
C...Zero CC 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 2300
        IF(ip.EQ.3) THEN
          tq=-lst(30)
          tqb=-tq
        ELSE
          tq=1.
          tqb=1.
        ENDIF
        DO 2200 i=1,lst(12)
        t1=-k(3,2)*qc(i)
        IF(t1.GT.0) THEN
          wq=xpq(i)*sig*tq
          wqb=0.
        ELSE
          wqb=xpq(-i)*sig*tqb
          wq=0.
        ENDIF
C...Include y-dependence.
        wq=wq*pari(23+ip)
        wqb=wqb*pari(23+ip)
        pqh(i,ih)=pqh(i,ih)+wq
        pqh(i+lst(12),ih)=pqh(i+lst(12),ih)+wqb
        pqh(17,ih)=pqh(17,ih)+wq+wqb
 2200   CONTINUE
        ENDIF
 2300 CONTINUE
 2500 CONTINUE
      DO 2600 i=1,17
 2600 pq(i)=(1.-parl(6))/2.*pqh(i,1)+(1.+parl(6))/2.*pqh(i,2)
      ih=1
      IF(lst(30).EQ.1) ih=2
      IF(lst(32).EQ.0) THEN
C...Simulation: cross section for chosen helicity state only.
        dsigma=pqh(17,ih)
      ELSEIF(lst(19).EQ.-1) THEN
C...Integration event-by-event: normalize and include alpha_s*1/(1-xp)
        dsigma=pqh(17,ih)/pari(20)*parl(25)/(1.-xp)
C...Max of dsigma/dxp for L- and R-handed lepton.
        IF(pqh(17,1).GT.pari(15)) pari(15)=pqh(17,1)
        IF(pqh(17,2).GT.pari(16)) pari(16)=pqh(17,2)
      ELSE
C...Integration for grid: normalize and include alpha_s*1/(1-xp)
        dsigma=pq(17)/pari(20)*parl(25)/(1.-xp)
        IF(lst(17).EQ.0) THEN
C...Fixed beam energy, max of dsigma/dxp for L- and R-handed lepton.
          IF(pqh(17,1).GT.pari(15)) pari(15)=pqh(17,1)
          IF(pqh(17,2).GT.pari(16)) pari(16)=pqh(17,2)
        ELSE
C...Variable beam energy, max of dsigma/dxp for S, T, I contributions.
          IF(pq(17)/pari(23+lst(32)).GT.pari(14+lst(32)))
     &    pari(14+lst(32))=pq(17)/pari(23+lst(32))
        ENDIF
      ENDIF
      RETURN

C...Boson-gluon fusion, i.e. qq-event.
 3000 s13=q2*(1.-xp)/xp
      IF(s13.LT.4.*amu**2) RETURN
      DO 3500 ip=il,iu
      sig=xpq(0)*dqcdi(2,ip,xp,zpmin,zpmax)
      DO 3300 ih=1,2
      IF(ih.EQ.1) THEN
        IF(parl(6).GT.0.99) goto 3300
        IF(lst(32).EQ.0.AND.lst(30).NE.-1) goto 3300
      ELSEIF(ih.EQ.2) THEN
        IF(parl(6).LT.-0.99) goto 3300
        IF(lst(32).EQ.0.AND.lst(30).NE.1) goto 3300
      ENDIF
      IF(lst(32).NE.0) lst(30)=sign(1.,ih-1.5)
      IF(lst(23).NE.2) THEN
        DO 3100 i=1,lst(13)
        mstj(93)=1
        IF(s13.LT.4.*ulmass(i)**2) goto 3100
        wq=sig/2.*(ewqc(1,ih,i)+ewqc(2,ih,i))
        wqb=wq
C...Include y-dependence.
        wq=wq*pari(23+ip)
        wqb=wqb*pari(23+ip)
        pqh(i,ih)=pqh(i,ih)+wq
        pqh(i+lst(13),ih)=pqh(i+lst(13),ih)+wqb
        pqh(17,ih)=pqh(17,ih)+wq+wqb
 3100   CONTINUE
      ELSEIF(lst(23).EQ.2) THEN
C...Zero CC 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 3300
        DO 3200 i=1,lst(13)
        mstj(93)=1
        IF(s13.LT.(amu+ulmass(i))**2) goto 3200
        IF(k(3,2)*qc(i).LT.0) THEN
          wq=sig
          wqb=0.
        ELSE
          wqb=sig
          wq=0.
        ENDIF
C...Include y-dependence.
        wq=wq*pari(23+ip)
        wqb=wqb*pari(23+ip)
        pqh(i,ih)=pqh(i,ih)+wq
        pqh(i+lst(13),ih)=pqh(i+lst(13),ih)+wqb
        pqh(17,ih)=pqh(17,ih)+wq+wqb
 3200   CONTINUE
      ENDIF
 3300 CONTINUE
 3500 CONTINUE
      DO 3600 i=1,17
 3600 pq(i)=(1.-parl(6))/2.*pqh(i,1)+(1.+parl(6))/2.*pqh(i,2)
      ih=1
      IF(lst(30).EQ.1) ih=2
      IF(lst(32).EQ.0) THEN
C...Simulation: cross section for chosen helicity state only.
        dsigma=pqh(17,ih)
      ELSEIF(lst(19).EQ.-1) THEN
C...Integration event-by-event: normalize and include alpha_s*1/(1-xp)
        dsigma=pqh(17,ih)/pari(20)*parl(25)/(1.-xp)
C...Max of dsigma/dxp for L- and R-handed lepton.
        IF(pqh(17,1).GT.pari(18)) pari(18)=pqh(17,1)
        IF(pqh(17,2).GT.pari(19)) pari(19)=pqh(17,2)
      ELSE
C...Integration for grid: normalize and include alpha_s*1/(1-xp)
        dsigma=pq(17)/pari(20)*parl(25)/(1.-xp)
        IF(lst(17).EQ.0) THEN
C...Fixed beam energy, max of dsigma/dxp for L- and R-handed lepton.
          IF(pqh(17,1).GT.pari(18)) pari(18)=pqh(17,1)
          IF(pqh(17,2).GT.pari(19)) pari(19)=pqh(17,2)
        ELSE
C...Variable beam energy, max of dsigma/dxp for S, T, I contributions.
          IF(pq(17)/pari(23+lst(32)).GT.pari(17+lst(32)))
     &    pari(17+lst(32))=pq(17)/pari(23+lst(32))
        ENDIF
      ENDIF
      RETURN
      END