attrad.f

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

C
C
C
C****************************************************************
C conduct soft radiation according to dipole approxiamtion
C****************************************************************
  SUBROUTINE attrad(IERROR)
C
  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
  SAVE  /hiparnt/
        common/hijdat/hidat0(10,10),hidat(10)
  SAVE  /hijdat/
  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
  SAVE  /lujets/
  ierror=0

C.....S INVARIANT MASS-SQUARED BETWEEN PARTONS I AND I+1......
C.....SM IS THE LARGEST MASS-SQUARED....

40  sm=0.
        jl=1
  DO 30 i=1,n-1
     s=2.*(p(i,4)*p(i+1,4)-p(i,1)*p(i+1,1)-p(i,2)*p(i+1,2)
     &    -p(i,3)*p(i+1,3))+p(i,5)**2+p(i+1,5)**2
     IF(s.LT.0.) s=0.
     wp=sqrt(s)-1.5*(p(i,5)+p(i+1,5))
     IF(wp.GT.sm) THEN
        pbt1=p(i,1)+p(i+1,1)
        pbt2=p(i,2)+p(i+1,2)
        pbt3=p(i,3)+p(i+1,3)
        pbt4=p(i,4)+p(i+1,4)
        btt=(pbt1**2+pbt2**2+pbt3**2)/pbt4**2
        IF(btt.GE.1.0-1.0e-10) go to 30
        IF((i.NE.1.OR.i.NE.n-1).AND.
     &             (k(i,2).NE.21.AND.k(i+1,2).NE.21)) go to 30
        jl=i
        sm=wp
     ENDIF
30  CONTINUE
  s=(sm+1.5*(p(jl,5)+p(jl+1,5)))**2
        IF(sm.LT.hipr1(5)) goto 2
     
C.....MAKE PLACE FOR ONE GLUON.....
        IF(jl+1.EQ.n) goto 190
        DO 160 j=n,jl+2,-1
          k(j+1,1)=k(j,1)
    k(j+1,2)=k(j,2)
          DO 150 m=1,5
C+++BAC
       v(j+1,m) = v(j,m)
C---BAC
150         p(j+1,m)=p(j,m)
160       CONTINUE
190     n=n+1
     
C.....BOOST TO REST SYSTEM FOR PARTICLES JL AND JL+1.....
        p1=p(jl,1)+p(jl+1,1)
        p2=p(jl,2)+p(jl+1,2)
        p3=p(jl,3)+p(jl+1,3)
        p4=p(jl,4)+p(jl+1,4)
        bex=-p1/p4
        bey=-p2/p4
        bez=-p3/p4
  imin=jl
  imax=jl+1
        CALL atrobo(0.,0.,bex,bey,bez,imin,imax,ierror)
  IF(ierror.NE.0) RETURN
C.....ROTATE TO Z-AXIS....
        cth=p(jl,3)/sqrt(p(jl,4)**2-p(jl,5)**2)
        IF(abs(cth).GT.1.0)  cth=max(-1.,min(1.,cth))
        theta=acos(cth)
        phi=ulangl(p(jl,1),p(jl,2))
        CALL atrobo(0.,-phi,0.,0.,0.,imin,imax,ierror)
        CALL atrobo(-theta,0.,0.,0.,0.,imin,imax,ierror)
     
C.....CREATE ONE GLUON AND ORIENTATE.....
     
1 CALL ar3jet(s,x1,x3,jl)
        CALL arorie(s,x1,x3,jl)   
  IF(hidat(2).GT.0.0) THEN
           ptg1=sqrt(p(jl,1)**2+p(jl,2)**2)
           ptg2=sqrt(p(jl+1,1)**2+p(jl+1,2)**2)
           ptg3=sqrt(p(jl+2,1)**2+p(jl+2,2)**2)
     ptg=max(ptg1,ptg2,ptg3)
     IF(ptg.GT.hidat(2)) THEN
        fmfact=exp(-(ptg**2-hidat(2)**2)/hipr1(2)**2)
        IF(atl_ran(nseed).GT.fmfact) go to 1
     ENDIF
  ENDIF
C.....ROTATE AND BOOST BACK.....
  imin=jl
  imax=jl+2
        CALL atrobo(theta,phi,-bex,-bey,-bez,imin,imax,ierror)
  IF(ierror.NE.0) RETURN
C.....ENUMERATE THE GLUONS.....
        k(jl+2,1)=k(jl+1,1)
  k(jl+2,2)=k(jl+1,2)
  k(jl+2,3)=k(jl+1,3)
  k(jl+2,4)=k(jl+1,4)
  k(jl+2,5)=k(jl+1,5)
        p(jl+2,5)=p(jl+1,5)
        k(jl+1,1)=2
  k(jl+1,2)=21
  k(jl+1,3)=0
  k(jl+1,4)=0
  k(jl+1,5)=0
        p(jl+1,5)=0.

  v(jl+1,1) = 0.
  v(jl+1,2) = 0.
  v(jl+1,3) = 0.
  v(jl+1,4) = 0.
  v(jl+1,5) = 0.

C----THETA FUNCTION DAMPING OF THE EMITTED GLUONS. FOR HADRON-HADRON.
C----R0=VFR(2)
C       IF(VFR(2).GT.0.) THEN
C       PTG=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
C       PTGMAX=WSTRI/2.
C       DOPT=SQRT((4.*PAR(71)*VFR(2))/WSTRI)
C       PTOPT=(DOPT*WSTRI)/(2.*VFR(2))
C       IF(PTG.GT.PTOPT) IORDER=IORDER-1
C       IF(PTG.GT.PTOPT) GOTO 1
C       ENDIF
C-----
      IF(sm.GE.hipr1(5)) goto 40

2       k(1,1)=2
  k(1,3)=0
  k(1,4)=0
  k(1,5)=0
        k(n,1)=1
  k(n,3)=0
  k(n,4)=0
  k(n,5)=0

        RETURN
        END