hijing/src/hijels.f

0001 C
0002 C
0003 C*******************************************************************
0004 CThis subroutine performs elastic scattering between two nucleons
0005 C
0006 C*******************************************************************
0007         SUBROUTINE HIJELS(PSC1,PSC2)
0008         IMPLICIT DOUBLE PRECISION(D)
0009         DIMENSION PSC1(5),PSC2(5)
0010         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0011         SAVE  /HIPARNT/
0012         COMMON/RANSEED/NSEED
0013         SAVE  /RANSEED/
0014 C
0015         CC=1.0-HINT1(12)/HINT1(13)
0016         RR=(1.0-CC)*HINT1(13)/HINT1(12)/(1.0-HIPR1(33))-1.0
0017         BB=0.5*(3.0+RR+SQRT(9.0+10.0*RR+RR**2))
0018         EP=SQRT((PSC1(1)-PSC2(1))**2+(PSC1(2)-PSC2(2))**2
0019      &          +(PSC1(3)-PSC2(3))**2)
0020         IF(EP.LE.0.1) RETURN
0021         ELS0=98.0/EP+52.0*(1.0+RR)**2
0022         PCM1=PSC1(1)+PSC2(1)
0023         PCM2=PSC1(2)+PSC2(2)
0024         PCM3=PSC1(3)+PSC2(3)
0025         ECM=PSC1(4)+PSC2(4)
0026         AM1=PSC1(5)**2
0027         AM2=PSC2(5)**2
0028         AMM=ECM**2-PCM1**2-PCM2**2-PCM3**2
0029         IF(AMM.LE.PSC1(5)+PSC2(5)) RETURN
0030 C               ********elastic scattering only when approaching
0031 C                               to each other
0032         PMAX=(AMM**2+AM1**2+AM2**2-2.0*AMM*AM1-2.0*AMM*AM2
0033      &                  -2.0*AM1*AM2)/4.0/AMM
0034         PMAX=ABS(PMAX)
0035 20      TT=ATL_RAN(NSEED)*MIN(PMAX,1.5)
0036         ELS=98.0*EXP(-2.8*TT)/EP
0037      &          +52.0*EXP(-9.2*TT)*(1.0+RR*EXP(-4.6*(BB-1.0)*TT))**2
0038         IF(ATL_RAN(NSEED).GT.ELS/ELS0) GO TO 20
0039         PHI=2.0*HIPR1(40)*ATL_RAN(NSEED)
0040 C
0041         DBX=PCM1/ECM
0042         DBY=PCM2/ECM
0043         DBZ=PCM3/ECM
0044         DB=SQRT(DBX**2+DBY**2+DBZ**2)
0045         IF(DB.GT.0.99999999D0) THEN
0046           DBX=DBX*(0.99999999D0/DB)
0047           DBY=DBY*(0.99999999D0/DB)
0048           DBZ=DBZ*(0.99999999D0/DB)
0049           DB=0.99999999D0
0050           WRITE(6,*) ' (HIJELS) boost vector too large'
0051 C               ********Rescale boost vector if too close to unity.
0052         ENDIF
0053         DGA=1D0/SQRT(1D0-DB**2)
0054 C
0055         DP1=SQRT(TT)*SIN(PHI)
0056         DP2=SQRT(TT)*COS(PHI)
0057         DP3=SQRT(PMAX-TT)
0058         DP4=SQRT(PMAX+AM1)
0059         DBP=DBX*DP1+DBY*DP2+DBZ*DP3
0060         DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4)
0061         PSC1(1)=DP1+DGABP*DBX
0062         PSC1(2)=DP2+DGABP*DBY
0063         PSC1(3)=DP3+DGABP*DBZ
0064         PSC1(4)=DGA*(DP4+DBP)
0065 C
0066         DP1=-SQRT(TT)*SIN(PHI)
0067         DP2=-SQRT(TT)*COS(PHI)
0068         DP3=-SQRT(PMAX-TT)
0069         DP4=SQRT(PMAX+AM2)
0070         DBP=DBX*DP1+DBY*DP2+DBZ*DP3
0071         DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4)
0072         PSC2(1)=DP1+DGABP*DBX
0073         PSC2(2)=DP2+DGABP*DBY
0074         PSC2(3)=DP3+DGABP*DBZ
0075         PSC2(4)=DGA*(DP4+DBP)
0076         RETURN
0077         END