[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / hijing.F

* $Id$
C     Version 1.36
C     Nothing important has been changed here. A few 'garbage' has been
C     cleaned up here, like common block HIJJET3 for the sea quark strings
C     which were originally created to implement the DPM scheme which
C     later was abadoned in the final version. The lines which operate
C     on these data are also deleted in the program.
C
C
C     Version 1.35
C     There are some changes in the program: subroutine HARDJET is now
C     consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
C     for the triggered hard processes. Now that is done  altogether
C     with other normal hard processes in modified JETINI. In the new
C     version one calls JETINI every time one calls HIJHRD. In the new
C     version the effect of the isospin of the nucleon on hard processes,
C     especially direct photons is correctly considered.
C     For A+A collisions, one has to initilize pythia
C     separately for each type of collisions, pp, pn,np and nn,
C     or hp and hn for hA collisions. In JETINI we use the following
C     catalogue for different types of collisions:
C     h+h: h+h (I_TYPE=1)
C     h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
C     A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
C     A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
C*****************************************************************
c
C
C     Version 1.34
C     Last modification on January 5, 1998. Two misstakes are corrected in
C     function G. A Misstake in the subroutine Parton is also corrected.
C     (These are pointed out by Ysushi Nara).
C
C
C       Last modifcation on April 10, 1996. To conduct final
C       state radiation, PYTHIA reorganize the two scattered
C       partons and their final momenta will be a little
C       different. The summed total momenta of the partons
C       from the final state radiation are stored in HINT1(26-29)
C       and HINT1(36-39) which are little different from 
C       HINT1(21-24) and HINT1(41-44).
C
C       Version 1.33
C
C       Last modfication  on September 11, 1995. When HIJING and
C       PYTHIA are initialized, the shadowing is evaluated at
C       b=0 which is the maximum. This will cause overestimate
C       of shadowing for peripheral interactions. To correct this
C       problem, shadowing is set to zero when initializing. Then
C       use these maximum  cross section without shadowing as a
C       normalization of the Monte Carlo. This however increase
C       the computing time. IHNT2(16) is used to indicate whether
C       the sturcture function is called for (IHNT2(16)=1) initialization
C       or for (IHNT2(16)=0)normal collisions simulation
C
C       Last modification on Aagust 28, 1994. Two bugs associate
C       with the impact parameter dependence of the shadowing is
C       corrected.
C
C
c       Last modification on October 14, 1994. One bug is corrected
c       in the direct photon production option in subroutine
C       HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
C       Another bug associated with keeping the decay history
C       in the particle information is also corrected.(this problem
C       was reported by Matt Bloomer)
C
C
C       Last modification on July 15, 1994. The option to trig on
C       heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1) 
C       is added. To do this, set IHPR2(3)=3. For inclusive production,
C       one should reset HIPR1(10)=0.0. One can also trig larger pt
C       QQbar production by giving HIPR1(10) a nonvanishing value.
C       The mass of the heavy quark in the calculation of the cross
C       section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
C       default is the charm mass D=1.5). We also include a separate
C       K-factor for heavy quark and direct photon production by
C       HIPR1(23)(D=2.0).
C
C       Last modification on May 24, 1994.  The option to
C       retain the information of all particles including those
C       who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is 
C       added to contain the line number of the parent particle 
C       of the current line which is produced via a decay. 
C       KATT(I,4) is the status number of the particle: 11=particle
C       which has decayed; 1=finally produced particle.
C
C
C       Last modification on May 24, 1994( in HIJSFT when valence quark
C       is quenched, the following error is corrected. 1.2*IHNT2(1) --> 
C       1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
C
C
C       Last modification on March 16, 1994 (heavy flavor production
C       processes MSUB(81)=1 MSUB(82)=1 have been switched on,
C       charm production is the default, B-quark option is
C       IHPR2(18), when it is switched on, charm quark is 
C       automatically off)
C
C
C       Last modification on March 23, 1994 (an error is corrected
C       in the impact parameter dependence of the jet cross section)
C
C       Last modification Oct. 1993 to comply with non-vax
C       machines' compiler 
C
C*********************************************
C	LAST MODIFICATION April 5, 1991
CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B 
C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18), 
C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
C
C	gluon jets can form a single string system
C
C	initial state radiation is included
C	
C	all QCD subprocesses are included
c
c	direct particles production is included(currently only direct
C		photon)
c
C	Effect of high P_T trigger bias on multiple jets distribution
c
C******************************************************************
C	                        HIJING.10                         *
C	          Heavy Ion Jet INteraction Generator        	  *
C	                           by                       	  *
C		   X. N. Wang      and   M. Gyulassy           	  *
C	 	      Lawrence Berkeley Laboratory		  *
C								  *
C******************************************************************
C
C******************************************************************
C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of  *
C proj, NFP(K,4) original ID of proj.  NFP(K,5)=colli status(0=no,*
C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it   *
C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
C (0=has not been,1=is going to be, -1=has already been scattered *
C NFP(k,11) total number of interactions this proj has suffered   *
C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant  *
C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and     *
C diquark,PP(K,10)=PT of the hard scattering between the valence  *
C quarks; PP(K,14,15)=the mass of quark,diquark.       		  * 
C******************************************************************
C
C****************************************************************
C
C	SUBROUTINE HIJING
C
C****************************************************************
	SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
	CHARACTER FRAME*8
	DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
     &			IPCOL(90000),ITCOL(90000)
#include "hiparnt.inc"
C
#include "hijcrdn.inc"
#include "himain1.inc"
#include "himain2.inc"
#include "histrng.inc"
#include "hijjet1.inc"
#include "hijjet2.inc"
#include "hijjet4.inc"
C
#include "lujets_hijing.inc"
#include "ludat1_hijing.inc"
        SAVE

	BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
	BMIN=MIN(BMIN0,BMAX)
	IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
		BMIN=0.0
		BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
	ENDIF
C			********HIPR1(31) is in mb =0.1fm**2
C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS 
C       BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
C
	YP(1,1)=0.0
	YP(2,1)=0.0
	YP(3,1)=0.0
	IF(IHNT2(1).LE.1) GO TO 14
	DO 10 KP=1,IHNT2(1)
5	R=HIRND(1)
	X=RLU_HIJING(0)
	CX=2.0*X-1.0
	SX=SQRT(1.0-CX*CX)
C		********choose theta from uniform cos(theta) distr
	PHI=RLU_HIJING(0)*2.0*HIPR1(40)
C		********choose phi form uniform phi distr 0 to 2*pi
	YP(1,KP)=R*SX*COS(PHI)
	YP(2,KP)=R*SX*SIN(PHI)
	YP(3,KP)=R*CX
	IF(HIPR1(29).EQ.0.0) GO TO 10
	DO 8  KP2=1,KP-1
		DNBP1=(YP(1,KP)-YP(1,KP2))**2
		DNBP2=(YP(2,KP)-YP(2,KP2))**2
		DNBP3=(YP(3,KP)-YP(3,KP2))**2
		DNBP=DNBP1+DNBP2+DNBP3
		IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
C			********two neighbors cannot be closer than 
C				HIPR1(29)
8	CONTINUE
10	CONTINUE
	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
	Y1=YP(1,I)
	Y2=YP(2,I)
	Y3=YP(3,I)
	YP(1,I)=YP(1,J)
	YP(2,I)=YP(2,J)
	YP(3,I)=YP(3,J)
	YP(1,J)=Y1
	YP(2,J)=Y2
	YP(3,J)=Y3
12	CONTINUE
C
C******************************
14	YT(1,1)=0.0
	YT(2,1)=0.0
	YT(3,1)=0.0
	IF(IHNT2(3).LE.1) GO TO 24
	DO 20 KT=1,IHNT2(3)
15	R=HIRND(2)
	X=RLU_HIJING(0)
	CX=2.0*X-1.0
	SX=SQRT(1.0-CX*CX)
C		********choose theta from uniform cos(theta) distr
	PHI=RLU_HIJING(0)*2.0*HIPR1(40)
C		********chose phi form uniform phi distr 0 to 2*pi
	YT(1,KT)=R*SX*COS(PHI)
	YT(2,KT)=R*SX*SIN(PHI)
	YT(3,KT)=R*CX
	IF(HIPR1(29).EQ.0.0) GO TO 20
	DO 18  KT2=1,KT-1
		DNBT1=(YT(1,KT)-YT(1,KT2))**2
		DNBT2=(YT(2,KT)-YT(2,KT2))**2
		DNBT3=(YT(3,KT)-YT(3,KT2))**2
		DNBT=DNBT1+DNBT2+DNBT3
		IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
C			********two neighbors cannot be closer than 
C				HIPR1(29)
18	CONTINUE
20	CONTINUE
	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
	Y1=YT(1,I)
	Y2=YT(2,I)
	Y3=YT(3,I)
	YT(1,I)=YT(1,J)
	YT(2,I)=YT(2,J)
	YT(3,I)=YT(3,J)
	YT(1,J)=Y1
	YT(2,J)=Y2
	YT(3,J)=Y3
22	CONTINUE

C********************
24	MISS=-1

50	MISS=MISS+1
	IF(MISS.GT.50) THEN
	   WRITE(6,*) 'infinite loop happened in  HIJING'
	   STOP
	ENDIF

	NATT=0
	JATT=0
	EATT=0.0
	CALL HIJINI
        NLOP=0
C			********Initialize for a new event
60	NT=0
	NP=0
	N0=0
	N01=0
	N10=0
	N11=0
	NSG=0
	NCOLT=0

C****	BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS 
C       RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET 
C       BY THE ANGLE PHI  FOR EACH COLLISION.******************
C
	BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
	PHI=2.0*HIPR1(40)*RLU_HIJING(0)
	BBX=BB*COS(PHI)
	BBY=BB*SIN(PHI)
	HINT1(19)=BB
	HINT1(20)=PHI
C
	DO 70 JP=1,IHNT2(1)
	DO 70 JT=1,IHNT2(3)
	   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
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)
	   RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
     &          /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
	   APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
     &           *SQRT(1.0-RRB1)
	   APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
     &           *SQRT(1.0-RRB2)
	   HINT1(18)=HINT1(14)-APHX1*HINT1(15)
     &			-APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
	   IF(IHPR2(14).EQ.0.OR.
     &          (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
	      GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
	      RANTOT=RLU_HIJING(0)
	      IF(RANTOT.GT.GS) GO TO 70
	      GO TO 65
	   ENDIF
	   GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
     &             /HIPR1(31)/2.0*ROMG(0.0)))
	   R2=R2/GSTOT_0
	   GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
	   GSTOT=2.0*(1.0-SQRT(1.0-GS))
	   RANTOT=RLU_HIJING(0)*GSTOT_0
	   IF(RANTOT.GT.GSTOT) GO TO 70
	   IF(RANTOT.GT.GS) THEN
	      CALL HIJCSC(JP,JT)
	      GO TO 70
C			********perform elastic collisions
	   ENDIF
 65	   SCIP(JP,JT)=R2
	   RNIP(JP,JT)=RANTOT
	   SJIP(JP,JT)=HINT1(18)
	   NCOLT=NCOLT+1
	   IPCOL(NCOLT)=JP
	   ITCOL(NCOLT)=JT
70	CONTINUE
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
C                       interaction at all. For NN collision
C                       repeat the event until interaction happens
C
	IF(IHPR2(3).NE.0) THEN
	   NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
	   NHARD=MIN(NHARD,NCOLT)
	   JPHARD=IPCOL(NHARD)
	   JTHARD=ITCOL(NHARD)
	ENDIF
C
	IF(IHPR2(9).EQ.1) THEN
		NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
		NMINI=MIN(NMINI,NCOLT)
		JPMINI=IPCOL(NMINI)
		JTMINI=ITCOL(NMINI)
	ENDIF
C		********Specifying the location of the hard and
C			minijet if they are enforced by user
C
	DO 200 JP=1,IHNT2(1)
	DO 200 JT=1,IHNT2(3)
	IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
		NFP(JP,11)=NFP(JP,11)+1
		NFT(JT,11)=NFT(JT,11)+1
	IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
		NP=NP+1
		N01=N01+1
	ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
		NT=NT+1
		N10=N10+1
	ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
		NP=NP+1
		NT=NT+1
		N0=N0+1
	ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
		N11=N11+1
	ENDIF
	JOUT=0
	NFP(JP,10)=0
	NFT(JT,10)=0
C*****************************************************************
	IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
C		********When IHPR2(8)=0 no jets are produced
	IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
C		********jets can not be produced for (JP,JT)
C			because not enough energy avaible for 
C				JP or JT 
	R2=SCIP(JP,JT)
	HINT1(18)=SJIP(JP,JT)
	TT=ROMG(R2)*HINT1(18)/HIPR1(31)
	TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
	NJET=0
	IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
           CALL JETINI(JP,JT,1)
           CALL HIJHRD(JP,JT,0,JFLG,0)
           HINT1(26)=HINT1(47)
           HINT1(27)=HINT1(48)
           HINT1(28)=HINT1(49)
           HINT1(29)=HINT1(50)
           HINT1(36)=HINT1(67)
           HINT1(37)=HINT1(68)
           HINT1(38)=HINT1(69)
           HINT1(39)=HINT1(70)
C
	   IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
	   IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
	   IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
     &				JFLG.GE.3) IASG(NSG,3)=1
	   IHNT2(9)=IHNT2(14)
	   IHNT2(10)=IHNT2(15)
	   DO 105 I05=1,5
	      HINT1(20+I05)=HINT1(40+I05)
	      HINT1(30+I05)=HINT1(50+I05)
 105	   CONTINUE
	   JOUT=1
	   IF(IHPR2(8).EQ.0) GO TO 160
	   RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
     &		/REAL(IHNT2(1))**0.6666667,1.0)
	   RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
     &		/REAL(IHNT2(3))**0.6666667,1.0)
	   APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
     &           *SQRT(1.0-RRB1)
	   APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
     &           *SQRT(1.0-RRB2)
	   HINT1(65)=HINT1(61)-APHX1*HINT1(62)
     &			-APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
	   TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
	   NJET=-1
C		********subtract the trigger jet from total number
C			of jet production  to be done since it has
C				already been produced here
	   XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
 106	   NJET=NJET+1
	   XR1=XR1-ALOG(RLU_HIJING(0))
	   IF(XR1.LT.TTRIG) GO TO 106
	   XR=0.0
 107	   NJET=NJET+1
	   XR=XR-ALOG(RLU_HIJING(0))
	   IF(XR.LT.TT-TTRIG) GO TO 107
	   NJET=NJET-1
	   GO TO 112
	ENDIF
C		********create a hard interaction with specified P_T
c				 when IHPR2(3)>0
	IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
C		********create at least one pair of mini jets 
C			when IHPR2(9)=1
C
	IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
     &		(1.0-EXP(-TTS))) GO TO 160
C		********this is the probability for no jet production
110	XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
111	NJET=NJET+1
	XR=XR-ALOG(RLU_HIJING(0))
	IF(XR.LT.TT) GO TO 111
112	NJET=MIN(NJET,IHPR2(8))
	IF(IHPR2(8).LT.0)  NJET=ABS(IHPR2(8))
C		******** Determine number of mini jet production
C
	DO 150 I_JET=1,NJET
           CALL JETINI(JP,JT,0)
	   CALL HIJHRD(JP,JT,JOUT,JFLG,1)
C		********JFLG=1 jets valence quarks, JFLG=2 with 
C			gluon jet, JFLG=3 with q-qbar prod for
C			(JP,JT). If JFLG=0 jets can not be produced 
C			this time. If JFLG=-1, error occured abandon
C			this event. JOUT is the total hard scat for
C			(JP,JT) up to now.
	   IF(JFLG.EQ.0) GO TO 160
	   IF(JFLG.LT.0) THEN
	      IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
	      GO TO 50
	   ENDIF
	   JOUT=JOUT+1
	   IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
	   IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
	   IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
     &			JFLG.GE.3) IASG(NSG,3)=1
C		******** jet with PT>HIPR1(11) will be quenched
 150	CONTINUE
 160	CONTINUE
	CALL HIJSFT(JP,JT,JOUT,IERROR)
	IF(IERROR.NE.0) THEN
	   IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
	   GO TO 50
	ENDIF
C
C		********conduct soft scattering between JP and JT
	JATT=JATT+JOUT

200	CONTINUE

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.
     &			IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
		DO 271 I=1,IHNT2(1)
			IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
271		CONTINUE
		DO 272 I=1,IHNT2(3)
			IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
272		CONTINUE
		DO 273 ISG=1,NSG
			IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
273		CONTINUE
	ENDIF
C
C**************fragment all the string systems in the following*****
C
C********N_ST is where particle information starts
C********N_STR+1 is the number of strings in fragmentation
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
        IF(IHPR2(20).NE.0) THEN
	   DO 360 ISG=1,NSG
		CALL HIJFRG(ISG,3,IERROR)
        	IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
		   MSTU(24)=0
		   MSTU(28)=0
		   IF(IHPR2(10).NE.0) THEN
		      call LULIST_HIJING(1)
		      WRITE(6,*) 'error occured, repeat the event'
		   ENDIF
		   GO TO 50
		ENDIF
C			********Check errors
C
		N_ST=1
		IDSTR=92
		IF(IHPR2(21).EQ.0) THEN
		   CALL LUEDIT_HIJING(2)
		ELSE
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
		ENDIF
C
		IF(FRAME.EQ.'LAB') THEN
			CALL HIBOOST
		ENDIF
C		******** boost back to lab frame(if it was in)
C
		N_STR=0
		DO 360 I=N_ST,N
		   IF(K(I,2).EQ.IDSTR) THEN
		      N_STR=N_STR+1
		      GO TO 360
		   ENDIF
		   K(I,4)=N_STR
		   NATT=NATT+1
		   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
C       ****** identify the mother particle
		   PATT(NATT,1)=P(I,1)
		   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)
360	   CONTINUE
C		********Fragment the q-qbar jets systems *****
C
	   JTP(1)=IHNT2(1)
	   JTP(2)=IHNT2(3)
	   DO 400 NTP=1,2
	   DO 400 J_JTP=1,JTP(NTP)
		CALL HIJFRG(J_JTP,NTP,IERROR)
        	IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
		   MSTU(24)=0
		   MSTU(28)=0
		   IF(IHPR2(10).NE.0) THEN
		      call LULIST_HIJING(1)
		      WRITE(6,*) 'error occured, repeat the event'
		   ENDIF
		   GO TO 50
		ENDIF
C			********check errors
C
		N_ST=1
		IDSTR=92
		IF(IHPR2(21).EQ.0) THEN
		   CALL LUEDIT_HIJING(2)
		ELSE
381		   N_ST=N_ST+1
		   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
		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
		      N_STR=N_STR+1
		      GO TO 390
		   ENDIF
		   K(I,4)=N_STR
		   NATT=NATT+1
		   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
		      KATT(NATT,3)=0
		   ELSE
		      KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
		   ENDIF
C       ****** identify the mother particle
		   PATT(NATT,1)=P(I,1)
		   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)
390		CONTINUE 
400	   CONTINUE
C		********Fragment the q-qq related string systems
	ENDIF

	DO 450 I=1,NDR
		NATT=NATT+1
		KATT(NATT,1)=KFDR(I)
		KATT(NATT,2)=40
		KATT(NATT,3)=0
		PATT(NATT,1)=PDR(I,1)
		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,4)=V(I,4)

		EATT=EATT+PDR(I,4)
450	CONTINUE
C			********store the direct-produced particles
C
	DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
	IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
     &     .AND.IHPR2(21).EQ.0) THEN
	IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
     &	   event'
C		call LULIST_HIJING(1)
		GO TO 50
	ENDIF
	RETURN
	END
Commit	Line	Data
e74335a4	1	* $Id$
	2	C Version 1.36
	3	C Nothing important has been changed here. A few 'garbage' has been
	4	C cleaned up here, like common block HIJJET3 for the sea quark strings
	5	C which were originally created to implement the DPM scheme which
	6	C later was abadoned in the final version. The lines which operate
	7	C on these data are also deleted in the program.
	8	C
	9	C
	10	C Version 1.35
	11	C There are some changes in the program: subroutine HARDJET is now
	12	C consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
	13	C for the triggered hard processes. Now that is done altogether
	14	C with other normal hard processes in modified JETINI. In the new
	15	C version one calls JETINI every time one calls HIJHRD. In the new
	16	C version the effect of the isospin of the nucleon on hard processes,
	17	C especially direct photons is correctly considered.
	18	C For A+A collisions, one has to initilize pythia
	19	C separately for each type of collisions, pp, pn,np and nn,
	20	C or hp and hn for hA collisions. In JETINI we use the following
	21	C catalogue for different types of collisions:
	22	C h+h: h+h (I_TYPE=1)
	23	C h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
	24	C A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
	25	C A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
	26	C*****************************************************************
	27	c
	28	C
	29	C Version 1.34
	30	C Last modification on January 5, 1998. Two misstakes are corrected in
	31	C function G. A Misstake in the subroutine Parton is also corrected.
	32	C (These are pointed out by Ysushi Nara).
	33	C
	34	C
	35	C Last modifcation on April 10, 1996. To conduct final
	36	C state radiation, PYTHIA reorganize the two scattered
	37	C partons and their final momenta will be a little
	38	C different. The summed total momenta of the partons
	39	C from the final state radiation are stored in HINT1(26-29)
	40	C and HINT1(36-39) which are little different from
	41	C HINT1(21-24) and HINT1(41-44).
	42	C
	43	C Version 1.33
	44	C
	45	C Last modfication on September 11, 1995. When HIJING and
	46	C PYTHIA are initialized, the shadowing is evaluated at
	47	C b=0 which is the maximum. This will cause overestimate
	48	C of shadowing for peripheral interactions. To correct this
	49	C problem, shadowing is set to zero when initializing. Then
	50	C use these maximum cross section without shadowing as a
	51	C normalization of the Monte Carlo. This however increase
	52	C the computing time. IHNT2(16) is used to indicate whether
	53	C the sturcture function is called for (IHNT2(16)=1) initialization
	54	C or for (IHNT2(16)=0)normal collisions simulation
	55	C
	56	C Last modification on Aagust 28, 1994. Two bugs associate
	57	C with the impact parameter dependence of the shadowing is
	58	C corrected.
	59	C
	60	C
	61	c Last modification on October 14, 1994. One bug is corrected
	62	c in the direct photon production option in subroutine
	63	C HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
	64	C Another bug associated with keeping the decay history
65	C in the particle information is also corrected.(this problem
66	C was reported by Matt Bloomer)
67	C
68	C
69	C Last modification on July 15, 1994. The option to trig on
70	C heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1)
71	C is added. To do this, set IHPR2(3)=3. For inclusive production,
72	C one should reset HIPR1(10)=0.0. One can also trig larger pt
73	C QQbar production by giving HIPR1(10) a nonvanishing value.
74	C The mass of the heavy quark in the calculation of the cross
75	C section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
76	C default is the charm mass D=1.5). We also include a separate
77	C K-factor for heavy quark and direct photon production by
78	C HIPR1(23)(D=2.0).
79	C
80	C Last modification on May 24, 1994. The option to
81	C retain the information of all particles including those
82	C who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is
83	C added to contain the line number of the parent particle
84	C of the current line which is produced via a decay.
85	C KATT(I,4) is the status number of the particle: 11=particle
86	C which has decayed; 1=finally produced particle.
87	C
88	C
89	C Last modification on May 24, 1994( in HIJSFT when valence quark
90	C is quenched, the following error is corrected. 1.2*IHNT2(1) -->
91	C 1.2IHNT2(1)0.333333, 1.2IHNT2(3) -->1.2IHNT(3)*0.333333)
92	C
93	C
94	C Last modification on March 16, 1994 (heavy flavor production
95	C processes MSUB(81)=1 MSUB(82)=1 have been switched on,
96	C charm production is the default, B-quark option is
97	C IHPR2(18), when it is switched on, charm quark is
98	C automatically off)
99	C
100	C
101	C Last modification on March 23, 1994 (an error is corrected
102	C in the impact parameter dependence of the jet cross section)
103	C
104	C Last modification Oct. 1993 to comply with non-vax
105	C machines' compiler
106	C
107	C*********************************************
108	C LAST MODIFICATION April 5, 1991
109	CQUARK DISTRIBUTIOIN (1-X)A/(X2+C2/S)B
110	C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
111	C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
112	C TWO COLOR CHANGES, (1-W)*2,W(1-W),W(1-W),AND W2, W=HIPR1(18),
113	C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
114	C
115	C gluon jets can form a single string system
116	C
117	C initial state radiation is included
118	C
119	C all QCD subprocesses are included
120	c
121	c direct particles production is included(currently only direct
122	C photon)
123	c
124	C Effect of high P_T trigger bias on multiple jets distribution
125	c
126	C******************************************************************
127	C HIJING.10 *
128	C Heavy Ion Jet INteraction Generator *
129	C by *
130	C X. N. Wang and M. Gyulassy *
131	C Lawrence Berkeley Laboratory *
132	C *
133	C******************************************************************
134	C
135	C******************************************************************
136	C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of *
137	C proj, NFP(K,4) original ID of proj. NFP(K,5)=colli status(0=no,*
138	C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
139	C \|NFP(K,6)\| is the total # of jet production, if NFP(K,6)<0 it *
140	C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
141	C (0=has not been,1=is going to be, -1=has already been scattered *
142	C NFP(k,11) total number of interactions this proj has suffered *
143	C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant *
144	C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and *
145	C diquark,PP(K,10)=PT of the hard scattering between the valence *
146	C quarks; PP(K,14,15)=the mass of quark,diquark. *
147	C******************************************************************
148	C
149	C****************************************************************
150	C
151	C SUBROUTINE HIJING
152	C
153	C****************************************************************
154	SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
155	CHARACTER FRAME*8
156	DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
157	& IPCOL(90000),ITCOL(90000)
158	#include "hiparnt.inc"
159	C
160	#include "hijcrdn.inc"
161	#include "himain1.inc"
162	#include "himain2.inc"
163	#include "histrng.inc"
164	#include "hijjet1.inc"
165	#include "hijjet2.inc"
166	#include "hijjet4.inc"
167	C
168	#include "lujets_hijing.inc"
169	#include "ludat1_hijing.inc"
170	SAVE
171
172	BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
173	BMIN=MIN(BMIN0,BMAX)
174	IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
175	BMIN=0.0
176	BMAX=2.5SQRT(HIPR1(31)0.1/HIPR1(40))
177	ENDIF
178	C ******HIPR1(31) is in mb =0.1fm2
179	C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
180	C BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
181	C
182	YP(1,1)=0.0
183	YP(2,1)=0.0
184	YP(3,1)=0.0
185	IF(IHNT2(1).LE.1) GO TO 14
186	DO 10 KP=1,IHNT2(1)
187	5 R=HIRND(1)
188	X=RLU_HIJING(0)
189	CX=2.0*X-1.0
190	SX=SQRT(1.0-CX*CX)
191	C ********choose theta from uniform cos(theta) distr
192	PHI=RLU_HIJING(0)2.0HIPR1(40)
193	C *******choose phi form uniform phi distr 0 to 2pi
194	YP(1,KP)=RSXCOS(PHI)
195	YP(2,KP)=RSXSIN(PHI)
196	YP(3,KP)=R*CX
197	IF(HIPR1(29).EQ.0.0) GO TO 10
198	DO 8 KP2=1,KP-1
199	DNBP1=(YP(1,KP)-YP(1,KP2))**2
200	DNBP2=(YP(2,KP)-YP(2,KP2))**2
201	DNBP3=(YP(3,KP)-YP(3,KP2))**2
202	DNBP=DNBP1+DNBP2+DNBP3
203	IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
204	C ********two neighbors cannot be closer than
205	C HIPR1(29)
206	8 CONTINUE
207	10 CONTINUE
208	DO 12 I=1,IHNT2(1)-1
209	DO 12 J=I+1,IHNT2(1)
210	IF(YP(3,I).GT.YP(3,J)) GO TO 12
211	Y1=YP(1,I)
212	Y2=YP(2,I)
213	Y3=YP(3,I)
214	YP(1,I)=YP(1,J)
215	YP(2,I)=YP(2,J)
216	YP(3,I)=YP(3,J)
217	YP(1,J)=Y1
218	YP(2,J)=Y2
219	YP(3,J)=Y3
220	12 CONTINUE
221	C
222	C******************************
223	14 YT(1,1)=0.0
224	YT(2,1)=0.0
225	YT(3,1)=0.0
226	IF(IHNT2(3).LE.1) GO TO 24
227	DO 20 KT=1,IHNT2(3)
228	15 R=HIRND(2)
229	X=RLU_HIJING(0)
230	CX=2.0*X-1.0
231	SX=SQRT(1.0-CX*CX)
232	C ********choose theta from uniform cos(theta) distr
233	PHI=RLU_HIJING(0)2.0HIPR1(40)
234	C *******chose phi form uniform phi distr 0 to 2pi
235	YT(1,KT)=RSXCOS(PHI)
236	YT(2,KT)=RSXSIN(PHI)
237	YT(3,KT)=R*CX
238	IF(HIPR1(29).EQ.0.0) GO TO 20
239	DO 18 KT2=1,KT-1
240	DNBT1=(YT(1,KT)-YT(1,KT2))**2
241	DNBT2=(YT(2,KT)-YT(2,KT2))**2
242	DNBT3=(YT(3,KT)-YT(3,KT2))**2
243	DNBT=DNBT1+DNBT2+DNBT3
244	IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
245	C ********two neighbors cannot be closer than
246	C HIPR1(29)
247	18 CONTINUE
248	20 CONTINUE
249	DO 22 I=1,IHNT2(3)-1
250	DO 22 J=I+1,IHNT2(3)
251	IF(YT(3,I).LT.YT(3,J)) GO TO 22
252	Y1=YT(1,I)
253	Y2=YT(2,I)
254	Y3=YT(3,I)
255	YT(1,I)=YT(1,J)
256	YT(2,I)=YT(2,J)
257	YT(3,I)=YT(3,J)
258	YT(1,J)=Y1
259	YT(2,J)=Y2
260	YT(3,J)=Y3
261	22 CONTINUE
262
263	C********************
264	24 MISS=-1
265
266	50 MISS=MISS+1
267	IF(MISS.GT.50) THEN
268	WRITE(6,*) 'infinite loop happened in HIJING'
269	STOP
270	ENDIF
271
272	NATT=0
273	JATT=0
274	EATT=0.0
275	CALL HIJINI
276	NLOP=0
277	C ********Initialize for a new event
278	60 NT=0
279	NP=0
280	N0=0
281	N01=0
282	N10=0
283	N11=0
284	NSG=0
285	NCOLT=0
286
287	C** BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB2 IS
288	C RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
289	C BY THE ANGLE PHI FOR EACH COLLISION.******************
290	C
291	BB=SQRT(BMIN*2+RLU_HIJING(0)(BMAX2-BMIN2))
292	PHI=2.0HIPR1(40)RLU_HIJING(0)
293	BBX=BB*COS(PHI)
294	BBY=BB*SIN(PHI)
295	HINT1(19)=BB
296	HINT1(20)=PHI
297	C
298	DO 70 JP=1,IHNT2(1)
299	DO 70 JT=1,IHNT2(3)
300	SCIP(JP,JT)=-1.0
301	B2=(YP(1,JP)+BBX-YT(1,JT))2+(YP(2,JP)+BBY-YT(2,JT))2
302	R2=B2*HIPR1(40)/HIPR1(31)/0.1
303	C *******mb=0.1fm, YP is in fm,HIPR1(31) is in mb
304	RRB1=MIN((YP(1,JP)2+YP(2,JP)2)
305	& /1.22/REAL(IHNT2(1))0.6666667,1.0)
306	RRB2=MIN((YT(1,JT)2+YT(2,JT)2)
307	& /1.22/REAL(IHNT2(3))0.6666667,1.0)
308	APHX1=HIPR1(6)4.0/3.0(IHNT2(1)**0.3333333-1.0)
309	& *SQRT(1.0-RRB1)
310	APHX2=HIPR1(6)4.0/3.0(IHNT2(3)**0.3333333-1.0)
311	& *SQRT(1.0-RRB2)
312	HINT1(18)=HINT1(14)-APHX1*HINT1(15)
313	& -APHX2HINT1(16)+APHX1APHX2*HINT1(17)
314	IF(IHPR2(14).EQ.0.OR.
315	& (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
316	GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
317	RANTOT=RLU_HIJING(0)
318	IF(RANTOT.GT.GS) GO TO 70
319	GO TO 65
320	ENDIF
321	GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
322	& /HIPR1(31)/2.0*ROMG(0.0)))
323	R2=R2/GSTOT_0
324	GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
325	GSTOT=2.0*(1.0-SQRT(1.0-GS))
326	RANTOT=RLU_HIJING(0)*GSTOT_0
327	IF(RANTOT.GT.GSTOT) GO TO 70
328	IF(RANTOT.GT.GS) THEN
329	CALL HIJCSC(JP,JT)
330	GO TO 70
331	C ********perform elastic collisions
332	ENDIF
333	65 SCIP(JP,JT)=R2
334	RNIP(JP,JT)=RANTOT
335	SJIP(JP,JT)=HINT1(18)
336	NCOLT=NCOLT+1
337	IPCOL(NCOLT)=JP
338	ITCOL(NCOLT)=JT
339	70 CONTINUE
340	C ********total number interactions proj and targ has
341	C suffered
342	IF(NCOLT.EQ.0) THEN
343	NLOP=NLOP+1
344	IF(NLOP.LE.20.OR.
345	& (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
346	RETURN
347	ENDIF
348	C ********At large impact parameter, there maybe no
349	C interaction at all. For NN collision
350	C repeat the event until interaction happens
351	C
352	IF(IHPR2(3).NE.0) THEN
353	NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
354	NHARD=MIN(NHARD,NCOLT)
355	JPHARD=IPCOL(NHARD)
356	JTHARD=ITCOL(NHARD)
357	ENDIF
358	C
359	IF(IHPR2(9).EQ.1) THEN
360	NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
361	NMINI=MIN(NMINI,NCOLT)
362	JPMINI=IPCOL(NMINI)
363	JTMINI=ITCOL(NMINI)
364	ENDIF
365	C ********Specifying the location of the hard and
366	C minijet if they are enforced by user
367	C
368	DO 200 JP=1,IHNT2(1)
369	DO 200 JT=1,IHNT2(3)
370	IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
371	NFP(JP,11)=NFP(JP,11)+1
372	NFT(JT,11)=NFT(JT,11)+1
373	IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
374	NP=NP+1
375	N01=N01+1
376	ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
377	NT=NT+1
378	N10=N10+1
379	ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
380	NP=NP+1
381	NT=NT+1
382	N0=N0+1
383	ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
384	N11=N11+1
385	ENDIF
386	JOUT=0
387	NFP(JP,10)=0
388	NFT(JT,10)=0
389	C*****************************************************************
390	IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
391	C ********When IHPR2(8)=0 no jets are produced
392	IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
393	C ********jets can not be produced for (JP,JT)
394	C because not enough energy avaible for
395	C JP or JT
396	R2=SCIP(JP,JT)
397	HINT1(18)=SJIP(JP,JT)
398	TT=ROMG(R2)*HINT1(18)/HIPR1(31)
399	TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
400	NJET=0
401	IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
402	CALL JETINI(JP,JT,1)
403	CALL HIJHRD(JP,JT,0,JFLG,0)
404	HINT1(26)=HINT1(47)
405	HINT1(27)=HINT1(48)
406	HINT1(28)=HINT1(49)
407	HINT1(29)=HINT1(50)
408	HINT1(36)=HINT1(67)
409	HINT1(37)=HINT1(68)
410	HINT1(38)=HINT1(69)
411	HINT1(39)=HINT1(70)
412	C
413	IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
414	IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
415	IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
416	& JFLG.GE.3) IASG(NSG,3)=1
417	IHNT2(9)=IHNT2(14)
418	IHNT2(10)=IHNT2(15)
419	DO 105 I05=1,5
420	HINT1(20+I05)=HINT1(40+I05)
421	HINT1(30+I05)=HINT1(50+I05)
422	105 CONTINUE
423	JOUT=1
424	IF(IHPR2(8).EQ.0) GO TO 160
425	RRB1=MIN((YP(1,JP)2+YP(2,JP)2)/1.2**2
426	& /REAL(IHNT2(1))**0.6666667,1.0)
427	RRB2=MIN((YT(1,JT)2+YT(2,JT)2)/1.2**2
428	& /REAL(IHNT2(3))**0.6666667,1.0)
429	APHX1=HIPR1(6)4.0/3.0(IHNT2(1)**0.3333333-1.0)
430	& *SQRT(1.0-RRB1)
431	APHX2=HIPR1(6)4.0/3.0(IHNT2(3)**0.3333333-1.0)
432	& *SQRT(1.0-RRB2)
433	HINT1(65)=HINT1(61)-APHX1*HINT1(62)
434	& -APHX2HINT1(63)+APHX1APHX2*HINT1(64)
435	TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
436	NJET=-1
437	C ********subtract the trigger jet from total number
438	C of jet production to be done since it has
439	C already been produced here
440	XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
441	106 NJET=NJET+1
442	XR1=XR1-ALOG(RLU_HIJING(0))
443	IF(XR1.LT.TTRIG) GO TO 106
444	XR=0.0
445	107 NJET=NJET+1
446	XR=XR-ALOG(RLU_HIJING(0))
447	IF(XR.LT.TT-TTRIG) GO TO 107
448	NJET=NJET-1
449	GO TO 112
450	ENDIF
451	C ********create a hard interaction with specified P_T
452	c when IHPR2(3)>0
453	IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
454	C ********create at least one pair of mini jets
455	C when IHPR2(9)=1
456	C
457	IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
458	& (1.0-EXP(-TTS))) GO TO 160
459	C ********this is the probability for no jet production
460	110 XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
461	111 NJET=NJET+1
462	XR=XR-ALOG(RLU_HIJING(0))
463	IF(XR.LT.TT) GO TO 111
464	112 NJET=MIN(NJET,IHPR2(8))
465	IF(IHPR2(8).LT.0) NJET=ABS(IHPR2(8))
466	C ******** Determine number of mini jet production
467	C
468	DO 150 I_JET=1,NJET
469	CALL JETINI(JP,JT,0)
470	CALL HIJHRD(JP,JT,JOUT,JFLG,1)
471	C ********JFLG=1 jets valence quarks, JFLG=2 with
472	C gluon jet, JFLG=3 with q-qbar prod for
473	C (JP,JT). If JFLG=0 jets can not be produced
474	C this time. If JFLG=-1, error occured abandon
475	C this event. JOUT is the total hard scat for
476	C (JP,JT) up to now.
477	IF(JFLG.EQ.0) GO TO 160
478	IF(JFLG.LT.0) THEN
479	IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
480	GO TO 50
481	ENDIF
482	JOUT=JOUT+1
483	IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
484	IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
485	IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
486	& JFLG.GE.3) IASG(NSG,3)=1
487	C ******** jet with PT>HIPR1(11) will be quenched
488	150 CONTINUE
489	160 CONTINUE
490	CALL HIJSFT(JP,JT,JOUT,IERROR)
491	IF(IERROR.NE.0) THEN
492	IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
493	GO TO 50
494	ENDIF
495	C
496	C ********conduct soft scattering between JP and JT
497	JATT=JATT+JOUT
498
499	200 CONTINUE
500
501	C******perform jet quenching for jets with PT>HIPR1(11)********
502
503	IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
504	& IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
505	DO 271 I=1,IHNT2(1)
506	IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
507	271 CONTINUE
508	DO 272 I=1,IHNT2(3)
509	IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
510	272 CONTINUE
511	DO 273 ISG=1,NSG
512	IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
513	273 CONTINUE
514	ENDIF
515	C
516	C************fragment all the string systems in the following***
517	C
518	C********N_ST is where particle information starts
519	C********N_STR+1 is the number of strings in fragmentation
520	C********the number of strings before a line is stored in K(I,4)
521	C********IDSTR is id number of the string system (91,92 or 93)
522	C
523	IF(IHPR2(20).NE.0) THEN
524	DO 360 ISG=1,NSG
525	CALL HIJFRG(ISG,3,IERROR)
526	IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
527	MSTU(24)=0
528	MSTU(28)=0
529	IF(IHPR2(10).NE.0) THEN
530	call LULIST_HIJING(1)
531	WRITE(6,*) 'error occured, repeat the event'
532	ENDIF
533	GO TO 50
534	ENDIF
535	C ********Check errors
536	C
537	N_ST=1
538	IDSTR=92
539	IF(IHPR2(21).EQ.0) THEN
540	CALL LUEDIT_HIJING(2)
541	ELSE
542	351 N_ST=N_ST+1
543	IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 351
544	IDSTR=K(N_ST,2)
545	N_ST=N_ST+1
546	ENDIF
547	C
548	IF(FRAME.EQ.'LAB') THEN
549	CALL HIBOOST
550	ENDIF
551	C ******** boost back to lab frame(if it was in)
552	C
553	N_STR=0
554	DO 360 I=N_ST,N
555	IF(K(I,2).EQ.IDSTR) THEN
556	N_STR=N_STR+1
557	GO TO 360
558	ENDIF
559	K(I,4)=N_STR
560	NATT=NATT+1
561	KATT(NATT,1)=K(I,2)
562	KATT(NATT,2)=20
563	KATT(NATT,4)=K(I,1)
564	IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
565	KATT(NATT,3)=0
566	ELSE
567	KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
568	ENDIF
569	C ****** identify the mother particle
570	PATT(NATT,1)=P(I,1)
571	PATT(NATT,2)=P(I,2)
572	PATT(NATT,3)=P(I,3)
573	PATT(NATT,4)=P(I,4)
2818cab6	574	VATT(NATT,1)=V(I,1)
	575	VATT(NATT,2)=V(I,2)
	576	VATT(NATT,3)=V(I,3)
	577	VATT(NATT,4)=V(I,4)
	578
e74335a4	579	EATT=EATT+P(I,4)
	580	360 CONTINUE
	581	C ******Fragment the q-qbar jets systems ***
	582	C
	583	JTP(1)=IHNT2(1)
	584	JTP(2)=IHNT2(3)
	585	DO 400 NTP=1,2
	586	DO 400 J_JTP=1,JTP(NTP)
	587	CALL HIJFRG(J_JTP,NTP,IERROR)
	588	IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
	589	MSTU(24)=0
	590	MSTU(28)=0
	591	IF(IHPR2(10).NE.0) THEN
	592	call LULIST_HIJING(1)
	593	WRITE(6,*) 'error occured, repeat the event'
	594	ENDIF
	595	GO TO 50
	596	ENDIF
	597	C ********check errors
	598	C
	599	N_ST=1
	600	IDSTR=92
	601	IF(IHPR2(21).EQ.0) THEN
	602	CALL LUEDIT_HIJING(2)
	603	ELSE
	604	381 N_ST=N_ST+1
	605	IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO 381
	606	IDSTR=K(N_ST,2)
	607	N_ST=N_ST+1
	608	ENDIF
	609	IF(FRAME.EQ.'LAB') THEN
	610	CALL HIBOOST
	611	ENDIF
	612	C ******** boost back to lab frame(if it was in)
	613	C
	614	NFTP=NFP(J_JTP,5)
	615	IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
	616	N_STR=0
	617	DO 390 I=N_ST,N
	618	IF(K(I,2).EQ.IDSTR) THEN
	619	N_STR=N_STR+1
	620	GO TO 390
	621	ENDIF
	622	K(I,4)=N_STR
	623	NATT=NATT+1
	624	KATT(NATT,1)=K(I,2)
	625	KATT(NATT,2)=NFTP
	626	KATT(NATT,4)=K(I,1)
	627	IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
	628	KATT(NATT,3)=0
	629	ELSE
	630	KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
	631	ENDIF
	632	C ****** identify the mother particle
	633	PATT(NATT,1)=P(I,1)
	634	PATT(NATT,2)=P(I,2)
	635	PATT(NATT,3)=P(I,3)
	636	PATT(NATT,4)=P(I,4)
2818cab6	637	VATT(NATT,1)=V(I,1)
	638	VATT(NATT,2)=V(I,2)
	639	VATT(NATT,3)=V(I,3)
	640	VATT(NATT,4)=V(I,4)
	641
e74335a4	642	EATT=EATT+P(I,4)
	643	390 CONTINUE
	644	400 CONTINUE
	645	C ********Fragment the q-qq related string systems
	646	ENDIF
	647
	648	DO 450 I=1,NDR
	649	NATT=NATT+1
	650	KATT(NATT,1)=KFDR(I)
	651	KATT(NATT,2)=40
	652	KATT(NATT,3)=0
	653	PATT(NATT,1)=PDR(I,1)
	654	PATT(NATT,2)=PDR(I,2)
	655	PATT(NATT,3)=PDR(I,3)
	656	PATT(NATT,4)=PDR(I,4)
2818cab6	657	VATT(NATT,1)=V(I,1)
	658	VATT(NATT,2)=V(I,2)
	659	VATT(NATT,3)=V(I,3)
	660	VATT(NATT,4)=V(I,4)
	661
e74335a4	662	EATT=EATT+PDR(I,4)
	663	450 CONTINUE
	664	C ********store the direct-produced particles
	665	C
	666	DENGY=EATT/(IHNT2(1)HINT1(6)+IHNT2(3)HINT1(7))-1.0
	667	IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
	668	& .AND.IHPR2(21).EQ.0) THEN
	669	IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
	670	& event'
	671	C call LULIST_HIJING(1)
	672	GO TO 50
	673	ENDIF
	674	RETURN
	675	END