hijels.f

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

C
C
C*******************************************************************
CThis subroutine performs elastic scattering between two nucleons
C
C*******************************************************************
        SUBROUTINE hijels(PSC1,PSC2)
        IMPLICIT DOUBLE PRECISION(d)
        dimension psc1(5),psc2(5)
        common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
        SAVE  /hiparnt/
        common/ranseed/nseed
        SAVE  /ranseed/
C
        cc=1.0-hint1(12)/hint1(13)
        rr=(1.0-cc)*hint1(13)/hint1(12)/(1.0-hipr1(33))-1.0
        bb=0.5*(3.0+rr+sqrt(9.0+10.0*rr+rr**2))
        ep=sqrt((psc1(1)-psc2(1))**2+(psc1(2)-psc2(2))**2
     &          +(psc1(3)-psc2(3))**2)
        IF(ep.LE.0.1) RETURN
        els0=98.0/ep+52.0*(1.0+rr)**2
        pcm1=psc1(1)+psc2(1)
        pcm2=psc1(2)+psc2(2)
        pcm3=psc1(3)+psc2(3)
        ecm=psc1(4)+psc2(4)
        am1=psc1(5)**2
        am2=psc2(5)**2
        amm=ecm**2-pcm1**2-pcm2**2-pcm3**2
        IF(amm.LE.psc1(5)+psc2(5)) RETURN
C               ********elastic scattering only when approaching
C                               to each other
        pmax=(amm**2+am1**2+am2**2-2.0*amm*am1-2.0*amm*am2
     &                  -2.0*am1*am2)/4.0/amm
        pmax=abs(pmax)
20      tt=atl_ran(nseed)*min(pmax,1.5)
        els=98.0*exp(-2.8*tt)/ep
     &          +52.0*exp(-9.2*tt)*(1.0+rr*exp(-4.6*(bb-1.0)*tt))**2
        IF(atl_ran(nseed).GT.els/els0) go to 20
        phi=2.0*hipr1(40)*atl_ran(nseed)
C
        dbx=pcm1/ecm
        dby=pcm2/ecm
        dbz=pcm3/ecm
        db=sqrt(dbx**2+dby**2+dbz**2)
        IF(db.GT.0.99999999d0) THEN 
          dbx=dbx*(0.99999999d0/db) 
          dby=dby*(0.99999999d0/db) 
          dbz=dbz*(0.99999999d0/db) 
          db=0.99999999d0   
          WRITE(6,*) ' (HIJELS) boost vector too large' 
C               ********Rescale boost vector if too close to unity. 
        ENDIF   
        dga=1d0/sqrt(1d0-db**2)      
C
        dp1=sqrt(tt)*sin(phi)
        dp2=sqrt(tt)*cos(phi)
        dp3=sqrt(pmax-tt)
        dp4=sqrt(pmax+am1)
        dbp=dbx*dp1+dby*dp2+dbz*dp3   
        dgabp=dga*(dga*dbp/(1d0+dga)+dp4) 
        psc1(1)=dp1+dgabp*dbx
        psc1(2)=dp2+dgabp*dby  
        psc1(3)=dp3+dgabp*dbz  
        psc1(4)=dga*(dp4+dbp)    
C       
        dp1=-sqrt(tt)*sin(phi)
        dp2=-sqrt(tt)*cos(phi)
        dp3=-sqrt(pmax-tt)
        dp4=sqrt(pmax+am2)
        dbp=dbx*dp1+dby*dp2+dbz*dp3   
        dgabp=dga*(dga*dbp/(1d0+dga)+dp4) 
        psc2(1)=dp1+dgabp*dbx
        psc2(2)=dp2+dgabp*dby  
        psc2(3)=dp3+dgabp*dbz  
        psc2(4)=dga*(dp4+dbp)
        RETURN
        END