SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
CHARACTER FRAME*8
DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
- & IPCOL(90000),ITCOL(90000)
+ & IPCOL(90000),ITCOL(90000),SCIP2(300,300)
+#define BLANKET_SAVE
#include "hiparnt.inc"
C
#include "hijcrdn.inc"
+#include "hijglbr.inc"
#include "himain1.inc"
#include "himain2.inc"
#include "histrng.inc"
IF(IHNT2(1).LE.1) GO TO 14
DO 10 KP=1,IHNT2(1)
5 R=HIRND(1)
+c
+ if(IHNT2(1).EQ.2) then
+ rnd1=max(RLU_HIJING(NSEED),1.0e-20)
+ rnd2=max(RLU_HIJING(NSEED),1.0e-20)
+ rnd3=max(RLU_HIJING(NSEED),1.0e-20)
+ R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
+ & +4.38*0.85*log(rnd3)/(4.38+0.85))
+ endif
+c
X=RLU_HIJING(0)
CX=2.0*X-1.0
SX=SQRT(1.0-CX*CX)
C HIPR1(29)
8 CONTINUE
10 CONTINUE
+c*******************************
+ if(IHNT2(1).EQ.2) then
+ YP(1,2)=-YP(1,1)
+ YP(2,2)=-YP(2,1)
+ YP(3,2)=-YP(3,1)
+ endif
+c********************************
DO 12 I=1,IHNT2(1)-1
DO 12 J=I+1,IHNT2(1)
IF(YP(3,I).GT.YP(3,J)) GO TO 12
IF(IHNT2(3).LE.1) GO TO 24
DO 20 KT=1,IHNT2(3)
15 R=HIRND(2)
+c
+ if(IHNT2(3).EQ.2) then
+ rnd1=max(RLU_HIJING(NSEED),1.0e-20)
+ rnd2=max(RLU_HIJING(NSEED),1.0e-20)
+ rnd3=max(RLU_HIJING(NSEED),1.0e-20)
+ R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
+ & +4.38*0.85*log(rnd3)/(4.38+0.85))
+ endif
+c
X=RLU_HIJING(0)
CX=2.0*X-1.0
SX=SQRT(1.0-CX*CX)
C HIPR1(29)
18 CONTINUE
20 CONTINUE
+c**********************************
+ if(IHNT2(3).EQ.2) then
+ YT(1,2)=-YT(1,1)
+ YT(2,2)=-YT(2,1)
+ YT(3,2)=-YT(3,1)
+ endif
+c*********************************
DO 22 I=1,IHNT2(3)-1
DO 22 J=I+1,IHNT2(3)
IF(YT(3,I).LT.YT(3,J)) GO TO 22
N01=0
N10=0
N11=0
+ NELT=0
+ NINT=0
+ NELP=0
+ NINP=0
NSG=0
NCOLT=0
+ NPSPECP=0
+ NNSPECP=0
+ NPSPECT=0
+ NNSPECT=0
C**** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
SCIP(JP,JT)=-1.0
B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
R2=B2*HIPR1(40)/HIPR1(31)/0.1
+ SCIP2(JP,JT)=R2
C ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
& /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
RNIP(JP,JT)=RANTOT
SJIP(JP,JT)=HINT1(18)
NCOLT=NCOLT+1
+ if (R2.GT.2.D0) THEN
+ write (8,*) R2
+ ENDIF
IPCOL(NCOLT)=JP
ITCOL(NCOLT)=JT
70 CONTINUE
+
+c *** cl glauber ***
+ npart=0
+ xmeana=0D0
+ ymeana=0D0
+ xmeanb=0D0
+ ymeanb=0D0
+ xmeanp=0D0
+ ymeanp=0D0
+ xm2=0D0
+ ym2=0D0
+ xym=0D0
+
+ IF(NCOLT>0) THEN
+ DO 1110 JP=1,IHNT2(1)
+ YP(1,JP)=YP(1,JP)+BBX
+ YP(2,JP)=YP(2,JP)+BBY
+ xmeana=xmeana+YP(1,JP)
+ ymeana=ymeana+YP(2,JP)
+ DO 1120 JT=1,IHNT2(3)
+ IF(SCIP2(JP,JT).LT.2.0D0) THEN
+ npart=npart+1
+ xmeanp=xmeanp+YP(1,JP)
+ ymeanp=ymeanp+YP(2,JP)
+ xm2=xm2+YP(1,JP)*YP(1,JP)
+ ym2=ym2+YP(2,JP)*YP(2,JP)
+ xym=xym+YP(1,JP)*YP(2,JP)
+ goto 1110
+ end if
+ 1120 continue
+ 1110 continue
+
+ DO 1130 JT=1,IHNT2(3)
+ xmeanb=xmeanb+YT(1,JT)
+ ymeanb=ymeanb+YT(2,JT)
+ DO 1140 JP=1,IHNT2(1)
+ IF(SCIP2(JP,JT).LT.2.0D0) THEN
+ npart=npart+1
+ xmeanp=xmeanp+YT(1,JT)
+ ymeanp=ymeanp+YT(2,JT)
+ xm2=xm2+YT(1,JT)*YT(1,JT)
+ ym2=ym2+YT(2,JT)*YT(2,JT)
+ xym=xym+YT(1,JT)*YT(2,JT)
+ goto 1130
+ end if
+ 1140 continue
+ 1130 continue
+
+ IF (npart.GT.0) THEN
+ xmeana=xmeana/IHNT2(1)
+ ymeana=ymeana/IHNT2(1)
+ xmeanb=xmeanb/IHNT2(3)
+ ymeanb=ymeanb/IHNT2(3)
+ xmeanp=xmeanp/npart
+ ymeanp=ymeanp/npart
+ xm2=xm2/npart
+ ym2=ym2/npart
+ xym=xym/npart
+
+ sx2=xm2-xmeanp*xmeanp
+ sy2=ym2-ymeanp*ymeanp
+ sxy=xym-xmeanp*ymeanp
+
+ delx=xmeanb-xmeana
+ dely=ymeanb-ymeana
+ dtmp=delx**2+dely**2
+ bbtrue=sqrt(dtmp)
+ dnumt=(sy2-sx2)*(delx**2-dely**2)-4D0*sxy*delx*dely
+ ddent=(sy2+sx2)*bbtrue**2
+ eccrp=dnumt/ddent
+ dtmp=(sy2-sx2)*(sy2-sx2)+4D0*sxy*sxy
+ eccpart=sqrt(dtmp)/(sx2+sy2)
+ eccmc=(sy2-sx2)/(sy2+sx2)
+ write(*,*),'HOUT: ',bb,' ',bbtrue,' ',ncolt,' ',npart,
+ 1 ' ',eccrp,' ',eccpart
+ end if
+ end if
+
C ********total number interactions proj and targ has
C suffered
+
IF(NCOLT.EQ.0) THEN
NLOP=NLOP+1
IF(NLOP.LE.20.OR.
& (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
+
+
RETURN
ENDIF
C ********At large impact parameter, there maybe no
200 CONTINUE
+c
+c**************************
+c
+ DO 201 JP=1,IHNT2(1)
+c write(6,*) JP, NFP(JP,3), NFP(JP,4), NFP(JP,5)
+ IF(NFP(JP,5).GT.2) THEN
+ NINP=NINP+1
+ ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
+ NELP=NELP+1
+ ENDIF
+
+ IF(NFP(JP,5).LE.2) THEN
+ IF (NFP(JP,3) .EQ. 2212) THEN
+ NPSPECP = NPSPECP + 1
+ ELSE IF (NFP(JP,3) .EQ. 2112) THEN
+ NNSPECP = NNSPECP + 1
+ ENDIF
+ ENDIF
+ 201 continue
+ DO 202 JT=1,IHNT2(3)
+ IF(NFT(JT,5).GT.2) THEN
+ NINT=NINT+1
+ ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
+ NELT=NELT+1
+ ENDIF
+
+ IF(NFT(JT,5).LE.2) THEN
+ IF (NFT(JT,3) .EQ. 2212) THEN
+ NPSPECT = NPSPECT + 1
+ ELSE IF (NFT(JT,3) .EQ. 2112) THEN
+ NNSPECT = NNSPECT + 1
+ ENDIF
+ ENDIF
+ 202 continue
+c
+c*******************************
+
C********perform jet quenching for jets with PT>HIPR1(11)**********
IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
C********the number of strings before a line is stored in K(I,4)
C********IDSTR is id number of the string system (91,92 or 93)
C
+c
IF(IHPR2(20).NE.0) THEN
DO 360 ISG=1,NSG
CALL HIJFRG(ISG,3,IERROR)
IF(IHPR2(21).EQ.0) THEN
CALL LUEDIT_HIJING(2)
ELSE
-351 N_ST=N_ST+1
+351 N_ST=N_ST+1
IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351
- IDSTR=K(N_ST,2)
- N_ST=N_ST+1
+ IDSTR=K(N_ST,2)
+ N_ST=N_ST+1
ENDIF
C
IF(FRAME.EQ.'LAB') THEN
KATT(NATT,1)=K(I,2)
KATT(NATT,2)=20
KATT(NATT,4)=K(I,1)
- IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
- KATT(NATT,3)=0
- ELSE
- KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
- ENDIF
+ IF(K(I,3).EQ.0 ) THEN
+ KATT(NATT,3)=0
+ ELSE
+ IF(K(K(I,3),2).EQ.IDSTR) THEN
+ KATT(NATT,3)=0
+ ELSE
+ KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
+ ENDIF
+ ENDIF
C ****** identify the mother particle
PATT(NATT,1)=P(I,1)
PATT(NATT,2)=P(I,2)
VATT(NATT,2)=V(I,2)
VATT(NATT,3)=V(I,3)
VATT(NATT,4)=V(I,4)
-
+
EATT=EATT+P(I,4)
360 CONTINUE
C ********Fragment the q-qbar jets systems *****
C
N_ST=1
IDSTR=92
+
+ NFTP=NFP(J_JTP,5)
+ IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
+
IF(IHPR2(21).EQ.0) THEN
CALL LUEDIT_HIJING(2)
- ELSE
+ ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN
381 N_ST=N_ST+1
-C PH Problem if N_ST > 9000, some additional protection added
-C PH 9000 is the first dimension of LUJETS/K(...)
- IF (N_ST.GT.9000) THEN
- call LULIST_HIJING(1)
- WRITE(6,*) 'error occured, repeat the event'
- GOTO 50
- ENDIF
IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
IDSTR=K(N_ST,2)
N_ST=N_ST+1
- ENDIF
+ ENDIF
IF(FRAME.EQ.'LAB') THEN
CALL HIBOOST
ENDIF
C ******** boost back to lab frame(if it was in)
C
- NFTP=NFP(J_JTP,5)
- IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
+
N_STR=0
DO 390 I=N_ST,N
IF(K(I,2).EQ.IDSTR) THEN
KATT(NATT,1)=K(I,2)
KATT(NATT,2)=NFTP
KATT(NATT,4)=K(I,1)
- IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
+ IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
KATT(NATT,3)=0
ELSE
KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
PATT(NATT,2)=P(I,2)
PATT(NATT,3)=P(I,3)
PATT(NATT,4)=P(I,4)
- VATT(NATT,1)=V(I,1)
- VATT(NATT,2)=V(I,2)
- VATT(NATT,3)=V(I,3)
- VATT(NATT,4)=V(I,4)
-
EATT=EATT+P(I,4)
+ VATT(NATT,1)=V(I,1)
+ VATT(NATT,2)=V(I,2)
+ VATT(NATT,3)=V(I,3)
+ VATT(NATT,4)=V(I,4)
390 CONTINUE
400 CONTINUE
C ********Fragment the q-qq related string systems
PATT(NATT,2)=PDR(I,2)
PATT(NATT,3)=PDR(I,3)
PATT(NATT,4)=PDR(I,4)
- VATT(NATT,1)=V(I,1)
- VATT(NATT,2)=V(I,2)
- VATT(NATT,3)=V(I,3)
+ VATT(NATT,1)=V(I,1)
+ VATT(NATT,2)=V(I,2)
+ VATT(NATT,3)=V(I,3)
VATT(NATT,4)=V(I,4)
-
EATT=EATT+PDR(I,4)
450 CONTINUE
C ********store the direct-produced particles
C call LULIST_HIJING(1)
GO TO 50
ENDIF
+
RETURN
END