hijing/src/hijing.f

0001 c     Version 1.383
0002 c     The variables I_SNG in HIJSFT and JL in ATTRAD were not initialized.
0003 c     The version initialize them. (as found by Fernando Marroquim)
0004 c
0005 c
0006 c
0007 c     Version 1.382
0008 c     Nuclear distribution for deuteron is taken as the Hulthen wave
0009 c     function as provided by Brian Cole (Columbia)
0010 c
0011 c
0012 c     Version 1.381
0013 c
0014 c     The parameters for Wood-Saxon distribution for deuteron are
0015 c     constrained to give the right rms ratius 2.116 fm
0016 c     (R=0.0, D=0.5882)
0017 c
0018 c
0019 c     Version 1.38
0020 c
0021 c     The following common block is added to record the number of elastic
0022 c     (NELT, NELP) and inelastic (NINT, NINP) participants
0023 c
0024 c        COMMON/HIJGLBR/NELT,NINT,NELP,NINP
0025 c        SAVE  /HIJGLBR/
0026 c
0027 c     Version 1.37
0028 c
0029 c     A bug in the quenching subroutine is corrected. When calculating the
0030 c     distance between two wounded nucleons, the displacement of the
0031 c     impact parameter was not inculded. This bug was discovered by
0032 c     Dr. V.Uzhinskii JINR, Dubna, Russia
0033 c
0034 c
0035 C     Version 1.36
0036 c
0037 c     Modification Oct. 8, 1998. In hijing, log(ran(nseed)) occasionally
0038 c     causes overfloat. It is modified to log(max(ran(nseed),1.0e-20)).
0039 c
0040 c
0041 C     Nothing important has been changed here. A few 'garbage' has been
0042 C     cleaned up here, like common block HIJJET3 for the sea quark strings
0043 C     which were originally created to implement the DPM scheme which
0044 C     later was abadoned in the final version. The lines which operate
0045 C     on these data are also deleted in the program.
0046 C
0047 C
0048 C     Version 1.35
0049 C     There are some changes in the program: subroutine HARDJET is now
0050 C     consolidated with HIJHRD. HARDJET is used to re-initiate PYTHIA
0051 C     for the triggered hard processes. Now that is done  altogether
0052 C     with other normal hard processes in modified JETINI. In the new
0053 C     version one calls JETINI every time one calls HIJHRD. In the new
0054 C     version the effect of the isospin of the nucleon on hard processes,
0055 C     especially direct photons is correctly considered.
0056 C     For A+A collisions, one has to initilize pythia
0057 C     separately for each type of collisions, pp, pn,np and nn,
0058 C     or hp and hn for hA collisions. In JETINI we use the following
0059 C     catalogue for different types of collisions:
0060 C     h+h: h+h (I_TYPE=1)
0061 C     h+A: h+p (I_TYPE=1), h+n (I_TYPE=2)
0062 C     A+h: p+h (I_TYPE=1), n+h (I_TYPE=2)
0063 C     A+A: p+p (I_TYPE=1), p+n (I_TYPE=2), n+p (I_TYPE=3), n+n (I_TYPE=4)
0064 C*****************************************************************
0065 c
0066 C
0067 C     Version 1.34
0068 C     Last modification on January 5, 1998. Two mistakes are corrected in
0069 C     function G. A Mistake in the subroutine Parton is also corrected.
0070 C     (These are pointed out by Ysushi Nara).
0071 C
0072 C
0073 C       Last modifcation on April 10, 1996. To conduct final
0074 C       state radiation, PYTHIA reorganize the two scattered
0075 C       partons and their final momenta will be a little
0076 C       different. The summed total momenta of the partons
0077 C       from the final state radiation are stored in HINT1(26-29)
0078 C       and HINT1(36-39) which are little different from
0079 C       HINT1(21-24) and HINT1(41-44).
0080 C
0081 C       Version 1.33
0082 C
0083 C       Last modfication  on September 11, 1995. When HIJING and
0084 C       PYTHIA are initialized, the shadowing is evaluated at
0085 C       b=0 which is the maximum. This will cause overestimate
0086 C       of shadowing for peripheral interactions. To correct this
0087 C       problem, shadowing is set to zero when initializing. Then
0088 C       use these maximum  cross section without shadowing as a
0089 C       normalization of the Monte Carlo. This however increase
0090 C       the computing time. IHNT2(16) is used to indicate whether
0091 C       the sturcture function is called for (IHNT2(16)=1) initialization
0092 C       or for (IHNT2(16)=0)normal collisions simulation
0093 C
0094 C       Last modification on Aagust 28, 1994. Two bugs associate
0095 C       with the impact parameter dependence of the shadowing is
0096 C       corrected.
0097 C
0098 C
0099 c       Last modification on October 14, 1994. One bug is corrected
0100 c       in the direct photon production option in subroutine
0101 C       HIJHRD.( this problem was reported by Jim Carroll and Mike Beddo).
0102 C       Another bug associated with keeping the decay history
0103 C       in the particle information is also corrected.(this problem
0104 C       was reported by Matt Bloomer)
0105 C
0106 C
0107 C       Last modification on July 15, 1994. The option to trig on
0108 C       heavy quark production (charm IHPR2(18)=0 or beauty IHPR2(18)=1)
0109 C       is added. To do this, set IHPR2(3)=3. For inclusive production,
0110 C       one should reset HIPR1(10)=0.0. One can also trig larger pt
0111 C       QQbar production by giving HIPR1(10) a nonvanishing value.
0112 C       The mass of the heavy quark in the calculation of the cross
0113 C       section (HINT1(59)--HINT1(65)) is given by HIPR1(7) (the
0114 C       default is the charm mass D=1.5). We also include a separate
0115 C       K-factor for heavy quark and direct photon production by
0116 C       HIPR1(23)(D=2.0).
0117 C
0118 C       Last modification on May 24, 1994.  The option to
0119 C       retain the information of all particles including those
0120 C       who have decayed is IHPR(21)=1 (default=0). KATT(I,3) is
0121 C       added to contain the line number of the parent particle
0122 C       of the current line which is produced via a decay.
0123 C       KATT(I,4) is the status number of the particle: 11=particle
0124 C       which has decayed; 1=finally produced particle.
0125 C
0126 C
0127 C       Last modification on May 24, 1994( in HIJSFT when valence quark
0128 C       is quenched, the following error is corrected. 1.2*IHNT2(1) -->
0129 C       1.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(3)**0.333333)
0130 C
0131 C
0132 C       Last modification on March 16, 1994 (heavy flavor production
0133 C       processes MSUB(81)=1 MSUB(82)=1 have been switched on,
0134 C       charm production is the default, B-quark option is
0135 C       IHPR2(18), when it is switched on, charm quark is
0136 C       automatically off)
0137 C
0138 C
0139 C       Last modification on March 23, 1994 (an error is corrected
0140 C       in the impact parameter dependence of the jet cross section)
0141 C
0142 C       Last modification Oct. 1993 to comply with non-vax
0143 C       machines' compiler
0144 C
0145 C*********************************************
0146 C       LAST MODIFICATION April 5, 1991
0147 CQUARK DISTRIBUTIOIN (1-X)**A/(X**2+C**2/S)**B
0148 C(A=HIPR1(44),B=HIPR1(46),C=HIPR1(45))
0149 C STRING FLIP, VENUS OPTION IHPR2(15)=1,IN WHICH ONE CAN HAVE ONE AND
0150 C TWO COLOR CHANGES, (1-W)**2,W*(1-W),W*(1-W),AND W*2, W=HIPR1(18),
0151 C AMONG PT DISTRIBUTION OF SEA QUARKS IS CONTROLLED BY HIPR1(42)
0152 C
0153 C       gluon jets can form a single string system
0154 C
0155 C       initial state radiation is included
0156 C
0157 C       all QCD subprocesses are included
0158 c
0159 c       direct particles production is included(currently only direct
0160 C               photon)
0161 c
0162 C       Effect of high P_T trigger bias on multiple jets distribution
0163 c
0164 C******************************************************************
0165 C                               HIJING.10                         *
0166 C                 Heavy Ion Jet INteraction Generator             *
0167 C                                  by                             *
0168 C                  X. N. Wang      and   M. Gyulassy              *
0169 C                     Lawrence Berkeley Laboratory                *
0170 C                                                                 *
0171 C******************************************************************
0172 C
0173 C******************************************************************
0174 C NFP(K,1),NFP(K,2)=flavor of q and di-q, NFP(K,3)=present ID of  *
0175 C proj, NFP(K,4) original ID of proj.  NFP(K,5)=colli status(0=no,*
0176 C 1=elastic,2=the diffrac one in single-diffrac,3= excited string.*
0177 C |NFP(K,6)| is the total # of jet production, if NFP(K,6)<0 it   *
0178 C can not produce jet anymore. NFP(K,10)=valence quarks scattering*
0179 C (0=has not been,1=is going to be, -1=has already been scattered *
0180 C NFP(k,11) total number of interactions this proj has suffered   *
0181 C PP(K,1)=PX,PP(K,2)=PY,PP(K,3)=PZ,PP(K,4)=E,PP(K,5)=M(invariant  *
0182 C mass), PP(K,6,7),PP(K,8,9)=transverse momentum of quark and     *
0183 C diquark,PP(K,10)=PT of the hard scattering between the valence  *
0184 C quarks; PP(K,14,15)=the mass of quark,diquark.                  *
0185 C******************************************************************
0186 C
0187 C****************************************************************
0188 C
0189 C       SUBROUTINE HIJING
0190 C
0191 C****************************************************************
0192         SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)
0193         CHARACTER FRAME*8
0194         DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
0195      &                  IPCOL(90000),ITCOL(90000)
0196         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0197         SAVE  /HIPARNT/
0198 C
0199         COMMON/HIJCRDN/YP(3,300),YT(3,300)
0200         SAVE  /HIJCRDN/
0201         COMMON/HIJGLBR/NELT,NINT,NELP,NINP
0202         SAVE  /HIJGLBR/
0203
0204 c+++BAC
0205 c
0206 c     Add an error entry (IERRSTAT) to HIMAIN1 common block
0207 c
0208 c       COMMON/HIMAIN1/NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
0209 c
0210         COMMON/HIMAIN1/NATT,EATT,JATT,NT,NP,N0,N01,N10,N11,IERRSTAT
0211 c---BAC
0212         SAVE  /HIMAIN1/
0213
0214 C+++BAC
0215 C       COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
0216 C       SAVE  /HIMAIN2/
0217 C
0218         COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4), VATT(130000,4)
0219         SAVE  /HIMAIN2/
0220 C---BAC
0221
0222         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
0223         SAVE  /HISTRNG/
0224         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
0225      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
0226      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
0227      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
0228      &                PJTE(300,500),PJTM(300,500)
0229         SAVE  /HIJJET1/
0230         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
0231      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
0232      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
0233         SAVE  /HIJJET2/
0234 C+++BAC
0235 C
0236 C       COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
0237 C       SAVE  /HIJJET4/
0238 C
0239         COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5), VDR(900,5)
0240         SAVE  /HIJJET4/
0241
0242 C---BAC
0243
0244         COMMON/RANSEED/NSEED
0245         SAVE  /RANSEED/
0246 C
0247         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0248         SAVE  /LUJETS/
0249         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0250         SAVE  /LUDAT1/
0251
0252 c+++BAC
0253 c
0254 c  Initialize error return code
0255 c
0256 c---BAC
0257         IERRSTAT = 0
0258
0259         BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
0260         BMIN=MIN(BMIN0,BMAX)
0261         IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
0262                 BMIN=0.0
0263                 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
0264         ENDIF
0265 C                       ********HIPR1(31) is in mb =0.1fm**2
0266 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
0267 C       BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
0268 C
0269         YP(1,1)=0.0
0270         YP(2,1)=0.0
0271         YP(3,1)=0.0
0272         IF(IHNT2(1).LE.1) GO TO 14
0273         DO 10 KP=1,IHNT2(1)
0274 5       R=HIRND(1)
0275 c
0276         if(IHNT2(1).EQ.2) then
0277            rnd1=max(ATL_RAN(NSEED),1.0e-20)
0278            rnd2=max(ATL_RAN(NSEED),1.0e-20)
0279            rnd3=max(ATL_RAN(NSEED),1.0e-20)
0280            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
0281      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
0282         endif
0283 c
0284         X=ATL_RAN(NSEED)
0285         CX=2.0*X-1.0
0286         SX=SQRT(1.0-CX*CX)
0287 C               ********choose theta from uniform cos(theta) distr
0288         PHI=ATL_RAN(NSEED)*2.0*HIPR1(40)
0289 C               ********choose phi form uniform phi distr 0 to 2*pi
0290         YP(1,KP)=R*SX*COS(PHI)
0291         YP(2,KP)=R*SX*SIN(PHI)
0292         YP(3,KP)=R*CX
0293         IF(HIPR1(29).EQ.0.0) GO TO 10
0294         DO 8  KP2=1,KP-1
0295                 DNBP1=(YP(1,KP)-YP(1,KP2))**2
0296                 DNBP2=(YP(2,KP)-YP(2,KP2))**2
0297                 DNBP3=(YP(3,KP)-YP(3,KP2))**2
0298                 DNBP=DNBP1+DNBP2+DNBP3
0299                 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
0300 C                       ********two neighbors cannot be closer than
0301 C                               HIPR1(29)
0302 8       CONTINUE
0303 10      CONTINUE
0304 c*******************************
0305         if(IHNT2(1).EQ.2) then
0306            YP(1,2)=-YP(1,1)
0307            YP(2,2)=-YP(2,1)
0308            YP(3,2)=-YP(3,1)
0309         endif
0310 c********************************
0311         DO 12 I=1,IHNT2(1)-1
0312         DO 12 J=I+1,IHNT2(1)
0313         IF(YP(3,I).GT.YP(3,J)) GO TO 12
0314         Y1=YP(1,I)
0315         Y2=YP(2,I)
0316         Y3=YP(3,I)
0317         YP(1,I)=YP(1,J)
0318         YP(2,I)=YP(2,J)
0319         YP(3,I)=YP(3,J)
0320         YP(1,J)=Y1
0321         YP(2,J)=Y2
0322         YP(3,J)=Y3
0323 12      CONTINUE
0324 C
0325 C******************************
0326 14      YT(1,1)=0.0
0327         YT(2,1)=0.0
0328         YT(3,1)=0.0
0329         IF(IHNT2(3).LE.1) GO TO 24
0330         DO 20 KT=1,IHNT2(3)
0331 15      R=HIRND(2)
0332 c
0333         if(IHNT2(3).EQ.2) then
0334            rnd1=max(ATL_RAN(NSEED),1.0e-20)
0335            rnd2=max(ATL_RAN(NSEED),1.0e-20)
0336            rnd3=max(ATL_RAN(NSEED),1.0e-20)
0337            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
0338      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
0339         endif
0340 c
0341         X=ATL_RAN(NSEED)
0342         CX=2.0*X-1.0
0343         SX=SQRT(1.0-CX*CX)
0344 C               ********choose theta from uniform cos(theta) distr
0345         PHI=ATL_RAN(NSEED)*2.0*HIPR1(40)
0346 C               ********chose phi form uniform phi distr 0 to 2*pi
0347         YT(1,KT)=R*SX*COS(PHI)
0348         YT(2,KT)=R*SX*SIN(PHI)
0349         YT(3,KT)=R*CX
0350         IF(HIPR1(29).EQ.0.0) GO TO 20
0351         DO 18  KT2=1,KT-1
0352                 DNBT1=(YT(1,KT)-YT(1,KT2))**2
0353                 DNBT2=(YT(2,KT)-YT(2,KT2))**2
0354                 DNBT3=(YT(3,KT)-YT(3,KT2))**2
0355                 DNBT=DNBT1+DNBT2+DNBT3
0356                 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
0357 C                       ********two neighbors cannot be closer than
0358 C                               HIPR1(29)
0359 18      CONTINUE
0360 20      CONTINUE
0361 c**********************************
0362         if(IHNT2(3).EQ.2) then
0363            YT(1,2)=-YT(1,1)
0364            YT(2,2)=-YT(2,1)
0365            YT(3,2)=-YT(3,1)
0366         endif
0367 c*********************************
0368         DO 22 I=1,IHNT2(3)-1
0369         DO 22 J=I+1,IHNT2(3)
0370         IF(YT(3,I).LT.YT(3,J)) GO TO 22
0371         Y1=YT(1,I)
0372         Y2=YT(2,I)
0373         Y3=YT(3,I)
0374         YT(1,I)=YT(1,J)
0375         YT(2,I)=YT(2,J)
0376         YT(3,I)=YT(3,J)
0377         YT(1,J)=Y1
0378         YT(2,J)=Y2
0379         YT(3,J)=Y3
0380 22      CONTINUE
0381 C********************
0382 24      MISS=-1
0383
0384 50      MISS=MISS+1
0385         IF(MISS.GT.50) THEN
0386            WRITE(6,*) 'infinite loop happened in  HIJING'
0387            STOP
0388         ENDIF
0389
0390         NATT=0
0391         JATT=0
0392         EATT=0.0
0393         CALL HIJINI
0394         NLOP=0
0395 C                       ********Initialize for a new event
0396 60      NT=0
0397         NP=0
0398         N0=0
0399         N01=0
0400         N10=0
0401         N11=0
0402         NELT=0
0403         NINT=0
0404         NELP=0
0405         NINP=0
0406         NSG=0
0407         NCOLT=0
0408
0409 C****   BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
0410 C       RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
0411 C       BY THE ANGLE PHI  FOR EACH COLLISION.******************
0412 C
0413         BB=SQRT(BMIN**2+ATL_RAN(NSEED)*(BMAX**2-BMIN**2))
0414         PHI=2.0*HIPR1(40)*ATL_RAN(NSEED)
0415         BBX=BB*COS(PHI)
0416         BBY=BB*SIN(PHI)
0417         HINT1(19)=BB
0418         HINT1(20)=PHI
0419 C
0420         DO 70 JP=1,IHNT2(1)
0421         DO 70 JT=1,IHNT2(3)
0422            SCIP(JP,JT)=-1.0
0423            B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
0424            R2=B2*HIPR1(40)/HIPR1(31)/0.1
0425 C               ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
0426            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
0427      &          /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
0428            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
0429      &          /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
0430            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
0431      &           *SQRT(1.0-RRB1)
0432            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
0433      &           *SQRT(1.0-RRB2)
0434            HINT1(18)=HINT1(14)-APHX1*HINT1(15)
0435      &                  -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
0436            IF(IHPR2(14).EQ.0.OR.
0437      &          (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
0438               GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
0439               RANTOT=ATL_RAN(NSEED)
0440               IF(RANTOT.GT.GS) GO TO 70
0441               GO TO 65
0442            ENDIF
0443            GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
0444      &             /HIPR1(31)/2.0*ROMG(0.0)))
0445            R2=R2/GSTOT_0
0446            GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
0447            GSTOT=2.0*(1.0-SQRT(1.0-GS))
0448            RANTOT=ATL_RAN(NSEED)*GSTOT_0
0449            IF(RANTOT.GT.GSTOT) GO TO 70
0450            IF(RANTOT.GT.GS) THEN
0451               CALL HIJCSC(JP,JT)
0452               GO TO 70
0453 C                       ********perform elastic collisions
0454            ENDIF
0455  65        SCIP(JP,JT)=R2
0456            RNIP(JP,JT)=RANTOT
0457            SJIP(JP,JT)=HINT1(18)
0458            NCOLT=NCOLT+1
0459            IPCOL(NCOLT)=JP
0460            ITCOL(NCOLT)=JT
0461 70      CONTINUE
0462 C               ********total number interactions proj and targ has
0463 C                               suffered
0464         IF(NCOLT.EQ.0) THEN
0465            NLOP=NLOP+1
0466            IF(NLOP.LE.20.OR.
0467      &           (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
0468            RETURN
0469         ENDIF
0470 C               ********At large impact parameter, there maybe no
0471 C                       interaction at all. For NN collision
0472 C                       repeat the event until interaction happens
0473 C
0474         IF(IHPR2(3).NE.0) THEN
0475            NHARD=1+INT(ATL_RAN(NSEED)*(NCOLT-1)+0.5)
0476            NHARD=MIN(NHARD,NCOLT)
0477            JPHARD=IPCOL(NHARD)
0478            JTHARD=ITCOL(NHARD)
0479         ENDIF
0480 C
0481         IF(IHPR2(9).EQ.1) THEN
0482                 NMINI=1+INT(ATL_RAN(NSEED)*(NCOLT-1)+0.5)
0483                 NMINI=MIN(NMINI,NCOLT)
0484                 JPMINI=IPCOL(NMINI)
0485                 JTMINI=ITCOL(NMINI)
0486         ENDIF
0487 C               ********Specifying the location of the hard and
0488 C                       minijet if they are enforced by user
0489 C
0490         DO 200 JP=1,IHNT2(1)
0491         DO 200 JT=1,IHNT2(3)
0492         IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
0493                 NFP(JP,11)=NFP(JP,11)+1
0494                 NFT(JT,11)=NFT(JT,11)+1
0495         IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
0496                 NP=NP+1
0497                 N01=N01+1
0498         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
0499                 NT=NT+1
0500                 N10=N10+1
0501         ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
0502                 NP=NP+1
0503                 NT=NT+1
0504                 N0=N0+1
0505         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
0506                 N11=N11+1
0507         ENDIF
0508 c
0509         JOUT=0
0510         NFP(JP,10)=0
0511         NFT(JT,10)=0
0512 C*****************************************************************
0513         IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
0514 C               ********When IHPR2(8)=0 no jets are produced
0515         IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
0516 C               ********jets can not be produced for (JP,JT)
0517 C                       because not enough energy avaible for
0518 C                               JP or JT
0519         R2=SCIP(JP,JT)
0520         HINT1(18)=SJIP(JP,JT)
0521         TT=ROMG(R2)*HINT1(18)/HIPR1(31)
0522         TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
0523         NJET=0
0524         IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
0525            CALL JETINI(JP,JT,1)
0526            CALL HIJHRD(JP,JT,0,JFLG,0)
0527            HINT1(26)=HINT1(47)
0528            HINT1(27)=HINT1(48)
0529            HINT1(28)=HINT1(49)
0530            HINT1(29)=HINT1(50)
0531            HINT1(36)=HINT1(67)
0532            HINT1(37)=HINT1(68)
0533            HINT1(38)=HINT1(69)
0534            HINT1(39)=HINT1(70)
0535 C
0536            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
0537            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
0538            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
0539      &                          JFLG.GE.3) IASG(NSG,3)=1
0540            IHNT2(9)=IHNT2(14)
0541            IHNT2(10)=IHNT2(15)
0542            DO 105 I05=1,5
0543               HINT1(20+I05)=HINT1(40+I05)
0544               HINT1(30+I05)=HINT1(50+I05)
0545  105       CONTINUE
0546            JOUT=1
0547            IF(IHPR2(8).EQ.0) GO TO 160
0548            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
0549      &          /REAL(IHNT2(1))**0.6666667,1.0)
0550            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
0551      &          /REAL(IHNT2(3))**0.6666667,1.0)
0552            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
0553      &           *SQRT(1.0-RRB1)
0554            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
0555      &           *SQRT(1.0-RRB2)
0556            HINT1(65)=HINT1(61)-APHX1*HINT1(62)
0557      &                  -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
0558            TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
0559            NJET=-1
0560 C               ********subtract the trigger jet from total number
0561 C                       of jet production  to be done since it has
0562 C                               already been produced here
0563            XR1=-ALOG(EXP(-TTRIG)+ATL_RAN(NSEED)*(1.0-EXP(-TTRIG)))
0564  106       NJET=NJET+1
0565            XR1=XR1-ALOG(max(ATL_RAN(NSEED),1.0e-20))
0566            IF(XR1.LT.TTRIG) GO TO 106
0567            XR=0.0
0568  107       NJET=NJET+1
0569            XR=XR-ALOG(max(ATL_RAN(NSEED),1.0e-20))
0570            IF(XR.LT.TT-TTRIG) GO TO 107
0571            NJET=NJET-1
0572            GO TO 112
0573         ENDIF
0574 C               ********create a hard interaction with specified P_T
0575 c                                when IHPR2(3)>0
0576         IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
0577 C               ********create at least one pair of mini jets
0578 C                       when IHPR2(9)=1
0579 C
0580         IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
0581      &          (1.0-EXP(-TTS))) GO TO 160
0582 C               ********this is the probability for no jet production
0583 110     XR=-ALOG(EXP(-TT)+ATL_RAN(NSEED)*(1.0-EXP(-TT)))
0584 111     NJET=NJET+1
0585         XR=XR-ALOG(max(ATL_RAN(NSEED),1.0e-20))
0586         IF(XR.LT.TT) GO TO 111
0587 112     NJET=MIN(NJET,IHPR2(8))
0588         IF(IHPR2(8).LT.0)  NJET=ABS(IHPR2(8))
0589 C               ******** Determine number of mini jet production
0590 C
0591         DO 150 I_JET=1,NJET
0592            CALL JETINI(JP,JT,0)
0593            CALL HIJHRD(JP,JT,JOUT,JFLG,1)
0594 C               ********JFLG=1 jets valence quarks, JFLG=2 with
0595 C                       gluon jet, JFLG=3 with q-qbar prod for
0596 C                       (JP,JT). If JFLG=0 jets can not be produced
0597 C                       this time. If JFLG=-1, error occured abandon
0598 C                       this event. JOUT is the total hard scat for
0599 C                       (JP,JT) up to now.
0600            IF(JFLG.EQ.0) GO TO 160
0601            IF(JFLG.LT.0) THEN
0602               IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
0603               GO TO 50
0604            ENDIF
0605            JOUT=JOUT+1
0606            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
0607            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
0608            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
0609      &                  JFLG.GE.3) IASG(NSG,3)=1
0610 C               ******** jet with PT>HIPR1(11) will be quenched
0611  150    CONTINUE
0612  160    CONTINUE
0613         CALL HIJSFT(JP,JT,JOUT,IERROR)
0614         IF(IERROR.NE.0) THEN
0615            IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
0616            GO TO 50
0617         ENDIF
0618 C
0619 C               ********conduct soft scattering between JP and JT
0620         JATT=JATT+JOUT
0621
0622 200     CONTINUE
0623 c
0624 c**************************
0625 c
0626         DO 201 JP=1,IHNT2(1)
0627            IF(NFP(JP,5).GT.2) THEN
0628               NINP=NINP+1
0629            ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
0630               NELP=NELP+1
0631            ENDIF
0632  201    continue
0633         DO 202 JT=1,IHNT2(3)
0634            IF(NFT(JT,5).GT.2) THEN
0635               NINT=NINT+1
0636            ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
0637               NELT=NELT+1
0638            ENDIF
0639  202    continue
0640 c
0641 c*******************************
0642
0643
0644 C********perform jet quenching for jets with PT>HIPR1(11)**********
0645
0646         IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
0647      &                  IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
0648                 DO 271 I=1,IHNT2(1)
0649                         IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
0650 271             CONTINUE
0651                 DO 272 I=1,IHNT2(3)
0652                         IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
0653 272             CONTINUE
0654                 DO 273 ISG=1,NSG
0655                         IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
0656 273             CONTINUE
0657         ENDIF
0658 C
0659 C**************fragment all the string systems in the following*****
0660 C
0661 C********N_ST is where particle information starts
0662 C********N_STR+1 is the number of strings in fragmentation
0663 C********the number of strings before a line is stored in K(I,4)
0664 C********IDSTR is id number of the string system (91,92 or 93)
0665 C
0666         IF(IHPR2(20).NE.0) THEN
0667            DO 360 ISG=1,NSG
0668                 CALL HIJFRG(ISG,3,IERROR)
0669                 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
0670                    MSTU(24)=0
0671                    MSTU(28)=0
0672                    IF(IHPR2(10).NE.0) THEN
0673                       call lulist(1)
0674                       WRITE(6,*) 'error occured, repeat the event'
0675                    ENDIF
0676                    GO TO 50
0677                 ENDIF
0678 C                       ********Check errors
0679 C
0680 C
0681 C DPM: call fastjet and record "jets" in LUJETS
0682 C
0683 c                NSAVE = N
0684                 CALL HIJFST(N,9000,K,P,V)
0685 c                WRITE(*,*)
0686 c                WRITE(*,*) N-NSAVE, ' "jets" produced'
0687 c                WRITE(*,*) ' id     px      py      pz      E       m  '
0688 c                WRITE(*,*) '==========================================='
0689 c                DO I = NSAVE+1, N
0690 c                   WRITE(*,2718) K(I,1), (P(I,J), J=1,5)
0691 c                ENDDO
0692 c 2718           FORMAT(I4,5(2X,F6.3))
0693 C
0694 C
0695 C
0696
0697                 N_ST=1
0698                 IDSTR=92
0699                 IF(IHPR2(21).EQ.0) THEN
0700                    CALL LUEDIT(2)
0701 c+++BAC
0702 c
0703 c   Don't perform the search for string if N=1 because it will search beyond the end of the valid particle list
0704 c
0705 c               ELSE
0706 c
0707                 ELSE IF (N .GT. 1) THEN
0708 c---BAC
0709
0710 351                N_ST=N_ST+1
0711
0712 c+++BAC
0713 c
0714 c  Check for inconsistency -- no string line found
0715 c
0716 c---BAC
0717                    if (N_ST .GT. N) then
0718                       IERRSTAT = 2
0719                       RETURN
0720                    ENDIF
0721
0722                    IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  351
0723                    IDSTR=K(N_ST,2)
0724                    N_ST=N_ST+1
0725                 ENDIF
0726 C
0727                 IF(FRAME.EQ.'LAB') THEN
0728                         CALL HIBOOST
0729                 ENDIF
0730 C               ******** boost back to lab frame(if it was in)
0731 C
0732                 N_STR=0
0733                 DO 361 I=N_ST,N
0734                    IF(K(I,2).EQ.IDSTR) THEN
0735 C+++BAC
0736                       IF (K(I,3) .LT. N_ST) THEN
0737 C---BAC
0738                          N_STR=N_STR+1
0739                          GO TO 360
0740                       ENDIF
0741                    ENDIF
0742
0743                    K(I,4)=N_STR
0744                    NATT=NATT+1
0745 c BAC+++
0746 c
0747 c   Add a check on array overflow
0748 c
0749 c BAC---
0750                    if (NATT .GT. 130000) THEN
0751                       IERRSTAT = 1
0752                       RETURN
0753                    ENDIF
0754
0755                    KATT(NATT,1)=K(I,2)
0756                    KATT(NATT,2)=20
0757                    KATT(NATT,4)=K(I,1)
0758 C+++BAC
0759 C                  IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
0760                    IF(K(I,3).EQ.0 .OR.
0761      &                K(I,3).LT.N_ST .OR.
0762      &                (K(K(I,3),2) .EQ. IDSTR .AND.
0763      &                 K(K(I,3),3) .LT. N_ST)) THEN
0764 C---BAC
0765                       KATT(NATT,3)=0
0766                    ELSE
0767                       KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
0768                    ENDIF
0769 C       ****** identify the mother particle
0770                    PATT(NATT,1)=P(I,1)
0771                    PATT(NATT,2)=P(I,2)
0772                    PATT(NATT,3)=P(I,3)
0773                    PATT(NATT,4)=P(I,4)
0774                    EATT=EATT+P(I,4)
0775
0776                    VATT(NATT,1)=V(I,1)
0777                    VATT(NATT,2)=V(I,2)
0778                    VATT(NATT,3)=V(I,3)
0779
0780                    IF ((ABS(VATT(NATT,3)) .GT. 0.00001) .and.
0781      &                 (KATT(NATT,3) .eq. 0 )) THEN
0782                       CALL LULIST(3)
0783                    ENDIF
0784
0785 C+++BAC
0786                    KPARENT =  KATT(NATT,3)
0787                    if (KPARENT .ne. 0) then
0788                       R = sqrt (VATT(NATT,1)**2 + VATT(NATT,2)**2 +
0789      &                     VATT(NATT,3)**2)
0790
0791                       RPARENT = sqrt (VATT(KPARENT,1)**2 +
0792      &                     VATT(KPARENT,2)**2 +
0793      &                     VATT(KPARENT,3)**2)
0794                       IF (R/RPARENT .LT. 0.85) THEN
0795                          CALL LULIST(3)
0796                       ENDIF
0797                    ENDIF
0798 C---BAC
0799
0800                    VATT(NATT,4)=V(I,4)
0801  361            CONTINUE
0802  360         CONTINUE
0803 C               ********Fragment the q-qbar jets systems *****
0804 C
0805            JTP(1)=IHNT2(1)
0806            JTP(2)=IHNT2(3)
0807            DO 400 NTP=1,2
0808            DO 400 J_JTP=1,JTP(NTP)
0809                 CALL HIJFRG(J_JTP,NTP,IERROR)
0810                 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
0811                    MSTU(24)=0
0812                    MSTU(28)=0
0813                    IF(IHPR2(10).NE.0) THEN
0814                       call lulist(1)
0815                       WRITE(6,*) 'error occured, repeat the event'
0816                    ENDIF
0817                    GO TO 50
0818                 ENDIF
0819 C                       ********check errors
0820 C
0821 C
0822 C DPM: call fastjet
0823 C
0824                 CALL HIJFST(N,9000,K,P,V)
0825 C
0826 C
0827 C
0828
0829                 N_ST=1
0830                 IDSTR=92
0831                 IF(IHPR2(21).EQ.0) THEN
0832                    CALL LUEDIT(2)
0833 c+++BAC
0834 c
0835 c   Don't perform the search for string if N=1 because it will search beyond the end of the valid particle list
0836 c
0837 c               ELSE
0838 c
0839                 ELSE IF (N .GT. 1) THEN
0840 c---BAC
0841 381                N_ST=N_ST+1
0842
0843 c+++BAC
0844 c
0845 c  Check for inconsistency -- no string line found
0846 c
0847 c---BAC
0848                    if (N_ST .GT. N) then
0849                       print *, 'inconsistency, n = ', N, ', n_st = ', N_ST
0850                       IERRSTAT = 2
0851                       RETURN
0852                    ENDIF
0853
0854                    IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  381
0855                    IDSTR=K(N_ST,2)
0856                    N_ST=N_ST+1
0857                 ENDIF
0858
0859                 IF(FRAME.EQ.'LAB') THEN
0860                         CALL HIBOOST
0861                 ENDIF
0862 C               ******** boost back to lab frame(if it was in)
0863 C
0864                 NFTP=NFP(J_JTP,5)
0865                 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
0866                 N_STR=0
0867                 DO 390 I=N_ST,N
0868                    IF(K(I,2).EQ.IDSTR) THEN
0869 C+++BAC
0870                       IF (K(I,3) .LT. N_ST) THEN
0871 C---BAC
0872                          N_STR=N_STR+1
0873                          GO TO 390
0874                       ENDIF
0875                    ENDIF
0876                    K(I,4)=N_STR
0877                    NATT=NATT+1
0878
0879 c BAC+++
0880 c
0881 c   Add a check on array overflow
0882 c
0883 c BAC---
0884                    if (NATT .GT. 130000) THEN
0885                       IERRSTAT = 1
0886                       RETURN
0887                    ENDIF
0888
0889                    KATT(NATT,1)=K(I,2)
0890                    KATT(NATT,2)=NFTP
0891                    KATT(NATT,4)=K(I,1)
0892 C+++BAC
0893 C                  IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
0894                    IF(K(I,3).EQ.0 .OR.
0895      &                K(I,3).LT.N_ST .OR.
0896      &                (K(K(I,3),2) .EQ. IDSTR .AND.
0897      &                 K(K(I,3),3) .LT. N_ST)) THEN
0898 C---BAC
0899                       KATT(NATT,3)=0
0900                    ELSE
0901                       KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
0902                    ENDIF
0903 C       ****** identify the mother particle
0904                    PATT(NATT,1)=P(I,1)
0905                    PATT(NATT,2)=P(I,2)
0906                    PATT(NATT,3)=P(I,3)
0907                    PATT(NATT,4)=P(I,4)
0908                    EATT=EATT+P(I,4)
0909
0910                    VATT(NATT,1)=V(I,1)
0911                    VATT(NATT,2)=V(I,2)
0912                    VATT(NATT,3)=V(I,3)
0913
0914                    if ((abs(VATT(NATT,3)) .gt. 0.00001) .and.
0915      &                 (KATT(NATT,3) .eq. 0 )) then
0916                       call lulist(3)
0917                    endif
0918
0919 C+++BAC
0920                    KPARENT =  KATT(NATT,3)
0921                    if (KPARENT .ne. 0) then
0922                       R = sqrt (VATT(NATT,1)**2 + VATT(NATT,2)**2 +
0923      &                     VATT(NATT,3)**2)
0924
0925                       RPARENT = sqrt (VATT(KPARENT,1)**2 +
0926      &                     VATT(KPARENT,2)**2 +
0927      &                     VATT(KPARENT,3)**2)
0928                       IF (R/RPARENT .LT. 0.85) THEN
0929                          CALL LULIST(3)
0930                       ENDIF
0931                    ENDIF
0932 C---BAC
0933
0934                    VATT(NATT,4)=V(I,4)
0935 390             CONTINUE
0936 400        CONTINUE
0937 C               ********Fragment the q-qq related string systems
0938         ENDIF
0939
0940         DO 450 I=1,NDR
0941                 NATT=NATT+1
0942
0943 c BAC+++
0944 c
0945 c   Add a check on array overflow
0946 c
0947 c BAC---
0948                 if (NATT .GT. 130000) THEN
0949                    IERRSTAT = 1
0950                    RETURN
0951                 ENDIF
0952
0953                 KATT(NATT,1)=KFDR(I)
0954                 KATT(NATT,2)=40
0955                 KATT(NATT,3)=0
0956                 PATT(NATT,1)=PDR(I,1)
0957                 PATT(NATT,2)=PDR(I,2)
0958                 PATT(NATT,3)=PDR(I,3)
0959                 PATT(NATT,4)=PDR(I,4)
0960                 EATT=EATT+PDR(I,4)
0961
0962                 VATT(NATT,1)=VDR(I,1)
0963                 VATT(NATT,2)=VDR(I,2)
0964                 VATT(NATT,3)=VDR(I,3)
0965                 VATT(NATT,4)=VDR(I,4)
0966 450     CONTINUE
0967 C                       ********store the direct-produced particles
0968 C
0969         DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
0970         IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
0971      &     .AND.IHPR2(21).EQ.0) THEN
0972         IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat event'
0973 C               call lulist(1)
0974                 GO TO 50
0975         ENDIF
0976         RETURN
0977         END