HIJING/hijing1_36/hijing.F

   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.2*IHNT2(1)**0.333333, 1.2*IHNT2(3) -->1.2*IHNT(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/(X**2+C**2/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 W*2, 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),SCIP2(300,300)
 158 #define BLANKET_SAVE
 159 #include "hiparnt.inc"
 160 C
 161 #include "hijcrdn.inc"
 162 #include "hijglbr.inc"
 163 #include "himain1.inc"
 164 #include "himain2.inc"
 165 #include "histrng.inc"
 166 #include "hijjet1.inc"
 167 #include "hijjet2.inc"
 168 #include "hijjet4.inc"
 169 C
 170 #include "lujets_hijing.inc"
 171 #include "ludat1_hijing.inc"
 172         SAVE
 173
 174         BMAX=MIN(BMAX0,HIPR1(34)+HIPR1(35))
 175         BMIN=MIN(BMIN0,BMAX)
 176         IF(IHNT2(1).LE.1 .AND. IHNT2(3).LE.1) THEN
 177                 BMIN=0.0
 178                 BMAX=2.5*SQRT(HIPR1(31)*0.1/HIPR1(40))
 179         ENDIF
 180 C                       ********HIPR1(31) is in mb =0.1fm**2
 181 C*******THE FOLLOWING IS TO SELECT THE COORDINATIONS OF NUCLEONS
 182 C       BOTH IN PROJECTILE AND TARGET NUCLEAR( in fm)
 183 C
 184         YP(1,1)=0.0
 185         YP(2,1)=0.0
 186         YP(3,1)=0.0
 187         IF(IHNT2(1).LE.1) GO TO 14
 188         DO 10 KP=1,IHNT2(1)
 189 5       R=HIRND(1)
 190 c
 191         if(IHNT2(1).EQ.2) then
 192            rnd1=max(RLU_HIJING(NSEED),1.0e-20)
 193            rnd2=max(RLU_HIJING(NSEED),1.0e-20)
 194            rnd3=max(RLU_HIJING(NSEED),1.0e-20)
 195            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
 196      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
 197         endif
 198 c
 199         X=RLU_HIJING(0)
 200         CX=2.0*X-1.0
 201         SX=SQRT(1.0-CX*CX)
 202 C               ********choose theta from uniform cos(theta) distr
 203         PHI=RLU_HIJING(0)*2.0*HIPR1(40)
 204 C               ********choose phi form uniform phi distr 0 to 2*pi
 205         YP(1,KP)=R*SX*COS(PHI)
 206         YP(2,KP)=R*SX*SIN(PHI)
 207         YP(3,KP)=R*CX
 208         IF(HIPR1(29).EQ.0.0) GO TO 10
 209         DO 8  KP2=1,KP-1
 210                 DNBP1=(YP(1,KP)-YP(1,KP2))**2
 211                 DNBP2=(YP(2,KP)-YP(2,KP2))**2
 212                 DNBP3=(YP(3,KP)-YP(3,KP2))**2
 213                 DNBP=DNBP1+DNBP2+DNBP3
 214                 IF(DNBP.LT.HIPR1(29)*HIPR1(29)) GO TO 5
 215 C                       ********two neighbors cannot be closer than
 216 C                               HIPR1(29)
 217 8       CONTINUE
 218 10      CONTINUE
 219 c*******************************
 220         if(IHNT2(1).EQ.2) then
 221            YP(1,2)=-YP(1,1)
 222            YP(2,2)=-YP(2,1)
 223            YP(3,2)=-YP(3,1)
 224         endif
 225 c********************************
 226         DO 12 I=1,IHNT2(1)-1
 227         DO 12 J=I+1,IHNT2(1)
 228         IF(YP(3,I).GT.YP(3,J)) GO TO 12
 229         Y1=YP(1,I)
 230         Y2=YP(2,I)
 231         Y3=YP(3,I)
 232         YP(1,I)=YP(1,J)
 233         YP(2,I)=YP(2,J)
 234         YP(3,I)=YP(3,J)
 235         YP(1,J)=Y1
 236         YP(2,J)=Y2
 237         YP(3,J)=Y3
 238 12      CONTINUE
 239 C
 240 C******************************
 241 14      YT(1,1)=0.0
 242         YT(2,1)=0.0
 243         YT(3,1)=0.0
 244         IF(IHNT2(3).LE.1) GO TO 24
 245         DO 20 KT=1,IHNT2(3)
 246 15      R=HIRND(2)
 247 c
 248          if(IHNT2(3).EQ.2) then
 249             rnd1=max(RLU_HIJING(NSEED),1.0e-20)
 250             rnd2=max(RLU_HIJING(NSEED),1.0e-20)
 251             rnd3=max(RLU_HIJING(NSEED),1.0e-20)
 252             R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
 253      &           +4.38*0.85*log(rnd3)/(4.38+0.85))
 254          endif
 255 c
 256         X=RLU_HIJING(0)
 257         CX=2.0*X-1.0
 258         SX=SQRT(1.0-CX*CX)
 259 C               ********choose theta from uniform cos(theta) distr
 260         PHI=RLU_HIJING(0)*2.0*HIPR1(40)
 261 C               ********chose phi form uniform phi distr 0 to 2*pi
 262         YT(1,KT)=R*SX*COS(PHI)
 263         YT(2,KT)=R*SX*SIN(PHI)
 264         YT(3,KT)=R*CX
 265         IF(HIPR1(29).EQ.0.0) GO TO 20
 266         DO 18  KT2=1,KT-1
 267                 DNBT1=(YT(1,KT)-YT(1,KT2))**2
 268                 DNBT2=(YT(2,KT)-YT(2,KT2))**2
 269                 DNBT3=(YT(3,KT)-YT(3,KT2))**2
 270                 DNBT=DNBT1+DNBT2+DNBT3
 271                 IF(DNBT.LT.HIPR1(29)*HIPR1(29)) GO TO 15
 272 C                       ********two neighbors cannot be closer than
 273 C                               HIPR1(29)
 274 18      CONTINUE
 275 20      CONTINUE
 276 c**********************************
 277          if(IHNT2(3).EQ.2) then
 278             YT(1,2)=-YT(1,1)
 279             YT(2,2)=-YT(2,1)
 280             YT(3,2)=-YT(3,1)
 281          endif
 282 c*********************************
 283         DO 22 I=1,IHNT2(3)-1
 284         DO 22 J=I+1,IHNT2(3)
 285         IF(YT(3,I).LT.YT(3,J)) GO TO 22
 286         Y1=YT(1,I)
 287         Y2=YT(2,I)
 288         Y3=YT(3,I)
 289         YT(1,I)=YT(1,J)
 290         YT(2,I)=YT(2,J)
 291         YT(3,I)=YT(3,J)
 292         YT(1,J)=Y1
 293         YT(2,J)=Y2
 294         YT(3,J)=Y3
 295 22      CONTINUE
 296
 297 C********************
 298 24      MISS=-1
 299
 300 50      MISS=MISS+1
 301         IF(MISS.GT.50) THEN
 302            WRITE(6,*) 'infinite loop happened in  HIJING'
 303            STOP
 304         ENDIF
 305
 306         NATT=0
 307         JATT=0
 308         EATT=0.0
 309         CALL HIJINI
 310         NLOP=0
 311 C                       ********Initialize for a new event
 312 60      NT=0
 313         NP=0
 314         N0=0
 315         N01=0
 316         N10=0
 317         N11=0
 318         NELT=0
 319         NINT=0
 320         NELP=0
 321         NINP=0
 322         NSG=0
 323         NCOLT=0
 324         NPSPECP=0
 325         NNSPECP=0
 326         NPSPECT=0
 327         NNSPECT=0
 328
 329 C****   BB IS THE ABSOLUTE VALUE OF IMPACT PARAMETER,BB**2 IS
 330 C       RANDOMLY GENERATED AND ITS ORIENTATION IS RANDOMLY SET
 331 C       BY THE ANGLE PHI  FOR EACH COLLISION.******************
 332 C
 333         BB=SQRT(BMIN**2+RLU_HIJING(0)*(BMAX**2-BMIN**2))
 334         PHI=2.0*HIPR1(40)*RLU_HIJING(0)
 335         BBX=BB*COS(PHI)
 336         BBY=BB*SIN(PHI)
 337         HINT1(19)=BB
 338         HINT1(20)=PHI
 339 C
 340         DO 70 JP=1,IHNT2(1)
 341         DO 70 JT=1,IHNT2(3)
 342            SCIP(JP,JT)=-1.0
 343            B2=(YP(1,JP)+BBX-YT(1,JT))**2+(YP(2,JP)+BBY-YT(2,JT))**2
 344            R2=B2*HIPR1(40)/HIPR1(31)/0.1
 345            SCIP2(JP,JT)=R2
 346 C               ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
 347            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)
 348      &          /1.2**2/REAL(IHNT2(1))**0.6666667,1.0)
 349            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)
 350      &          /1.2**2/REAL(IHNT2(3))**0.6666667,1.0)
 351            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
 352      &           *SQRT(1.0-RRB1)
 353            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
 354      &           *SQRT(1.0-RRB2)
 355            HINT1(18)=HINT1(14)-APHX1*HINT1(15)
 356      &                  -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
 357            IF(IHPR2(14).EQ.0.OR.
 358      &          (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) THEN
 359               GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
 360               RANTOT=RLU_HIJING(0)
 361               IF(RANTOT.GT.GS) GO TO 70
 362               GO TO 65
 363            ENDIF
 364            GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
 365      &             /HIPR1(31)/2.0*ROMG(0.0)))
 366            R2=R2/GSTOT_0
 367            GS=1.0-EXP(-(HIPR1(30)+HINT1(18))/HIPR1(31)*ROMG(R2))
 368            GSTOT=2.0*(1.0-SQRT(1.0-GS))
 369            RANTOT=RLU_HIJING(0)*GSTOT_0
 370            IF(RANTOT.GT.GSTOT) GO TO 70
 371            IF(RANTOT.GT.GS) THEN
 372               CALL HIJCSC(JP,JT)
 373               GO TO 70
 374 C                       ********perform elastic collisions
 375            ENDIF
 376  65        SCIP(JP,JT)=R2
 377            RNIP(JP,JT)=RANTOT
 378            SJIP(JP,JT)=HINT1(18)
 379            NCOLT=NCOLT+1
 380            if (R2.GT.2.D0) THEN
 381               write (8,*) R2
 382            ENDIF
 383            IPCOL(NCOLT)=JP
 384            ITCOL(NCOLT)=JT
 385 70      CONTINUE
 386
 387 c *** cl glauber ***
 388         npart=0
 389         xmeana=0D0
 390         ymeana=0D0
 391         xmeanb=0D0
 392         ymeanb=0D0
 393         xmeanp=0D0
 394         ymeanp=0D0
 395         xm2=0D0
 396         ym2=0D0
 397         xym=0D0
 398
 399         IF(NCOLT>0) THEN
 400            DO 1110 JP=1,IHNT2(1)
 401               YP(1,JP)=YP(1,JP)+BBX
 402               YP(2,JP)=YP(2,JP)+BBY
 403               xmeana=xmeana+YP(1,JP)
 404               ymeana=ymeana+YP(2,JP)
 405               DO 1120 JT=1,IHNT2(3)
 406                  IF(SCIP2(JP,JT).LT.2.0D0) THEN
 407                     npart=npart+1
 408                     xmeanp=xmeanp+YP(1,JP)
 409                     ymeanp=ymeanp+YP(2,JP)
 410                     xm2=xm2+YP(1,JP)*YP(1,JP)
 411                     ym2=ym2+YP(2,JP)*YP(2,JP)
 412                     xym=xym+YP(1,JP)*YP(2,JP)
 413                     goto 1110
 414                  end if
 415  1120         continue
 416  1110      continue
 417
 418            DO 1130 JT=1,IHNT2(3)
 419               xmeanb=xmeanb+YT(1,JT)
 420               ymeanb=ymeanb+YT(2,JT)
 421               DO 1140 JP=1,IHNT2(1)
 422                  IF(SCIP2(JP,JT).LT.2.0D0) THEN
 423                     npart=npart+1
 424                     xmeanp=xmeanp+YT(1,JT)
 425                     ymeanp=ymeanp+YT(2,JT)
 426                     xm2=xm2+YT(1,JT)*YT(1,JT)
 427                     ym2=ym2+YT(2,JT)*YT(2,JT)
 428                     xym=xym+YT(1,JT)*YT(2,JT)
 429                     goto 1130
 430                  end if
 431  1140         continue
 432  1130      continue
 433
 434            xmeana=xmeana/IHNT2(1)
 435            ymeana=ymeana/IHNT2(1)
 436            xmeanb=xmeanb/IHNT2(3)
 437            ymeanb=ymeanb/IHNT2(3)
 438            xmeanp=xmeanp/npart
 439            ymeanp=ymeanp/npart
 440            xm2=xm2/npart
 441            ym2=ym2/npart
 442            xym=xym/npart
 443
 444            sx2=xm2-xmeanp*xmeanp
 445            sy2=ym2-ymeanp*ymeanp
 446            sxy=xym-xmeanp*ymeanp
 447
 448            delx=xmeanb-xmeana
 449            dely=ymeanb-ymeana
 450            dtmp=delx**2+dely**2
 451            bbtrue=sqrt(dtmp)
 452            dnumt=(sy2-sx2)*(delx**2-dely**2)-4D0*sxy*delx*dely
 453            ddent=(sy2+sx2)*bbtrue**2
 454            eccrp=dnumt/ddent
 455            dtmp=(sy2-sx2)*(sy2-sx2)+4D0*sxy*sxy
 456            eccpart=sqrt(dtmp)/(sx2+sy2)
 457            eccmc=(sy2-sx2)/(sy2+sx2)
 458            write(*,*),'HOUT: ',bb,' ',bbtrue,' ',ncolt,' ',npart,
 459      1          ' ',eccrp,' ',eccpart
 460         end if
 461
 462 C               ********total number interactions proj and targ has
 463 C                               suffered
 464
 465         IF(NCOLT.EQ.0) THEN
 466            NLOP=NLOP+1
 467            IF(NLOP.LE.20.OR.
 468      &           (IHNT2(1).EQ.1.AND.IHNT2(3).EQ.1)) GO TO 60
 469
 470
 471            RETURN
 472         ENDIF
 473 C               ********At large impact parameter, there maybe no
 474 C                       interaction at all. For NN collision
 475 C                       repeat the event until interaction happens
 476 C
 477         IF(IHPR2(3).NE.0) THEN
 478            NHARD=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
 479            NHARD=MIN(NHARD,NCOLT)
 480            JPHARD=IPCOL(NHARD)
 481            JTHARD=ITCOL(NHARD)
 482         ENDIF
 483 C
 484         IF(IHPR2(9).EQ.1) THEN
 485                 NMINI=1+INT(RLU_HIJING(0)*(NCOLT-1)+0.5)
 486                 NMINI=MIN(NMINI,NCOLT)
 487                 JPMINI=IPCOL(NMINI)
 488                 JTMINI=ITCOL(NMINI)
 489         ENDIF
 490 C               ********Specifying the location of the hard and
 491 C                       minijet if they are enforced by user
 492 C
 493         DO 200 JP=1,IHNT2(1)
 494         DO 200 JT=1,IHNT2(3)
 495         IF(SCIP(JP,JT).EQ.-1.0) GO TO 200
 496                 NFP(JP,11)=NFP(JP,11)+1
 497                 NFT(JT,11)=NFT(JT,11)+1
 498         IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).GT.1) THEN
 499                 NP=NP+1
 500                 N01=N01+1
 501         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).LE.1) THEN
 502                 NT=NT+1
 503                 N10=N10+1
 504         ELSE IF(NFP(JP,5).LE.1 .AND. NFT(JT,5).LE.1) THEN
 505                 NP=NP+1
 506                 NT=NT+1
 507                 N0=N0+1
 508         ELSE IF(NFP(JP,5).GT.1 .AND. NFT(JT,5).GT.1) THEN
 509                 N11=N11+1
 510         ENDIF
 511         JOUT=0
 512         NFP(JP,10)=0
 513         NFT(JT,10)=0
 514 C*****************************************************************
 515         IF(IHPR2(8).EQ.0 .AND. IHPR2(3).EQ.0) GO TO 160
 516 C               ********When IHPR2(8)=0 no jets are produced
 517         IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) GO TO 160
 518 C               ********jets can not be produced for (JP,JT)
 519 C                       because not enough energy avaible for
 520 C                               JP or JT
 521         R2=SCIP(JP,JT)
 522         HINT1(18)=SJIP(JP,JT)
 523         TT=ROMG(R2)*HINT1(18)/HIPR1(31)
 524         TTS=HIPR1(30)*ROMG(R2)/HIPR1(31)
 525         NJET=0
 526         IF(IHPR2(3).NE.0 .AND. JP.EQ.JPHARD .AND. JT.EQ.JTHARD) THEN
 527            CALL JETINI(JP,JT,1)
 528            CALL HIJHRD(JP,JT,0,JFLG,0)
 529            HINT1(26)=HINT1(47)
 530            HINT1(27)=HINT1(48)
 531            HINT1(28)=HINT1(49)
 532            HINT1(29)=HINT1(50)
 533            HINT1(36)=HINT1(67)
 534            HINT1(37)=HINT1(68)
 535            HINT1(38)=HINT1(69)
 536            HINT1(39)=HINT1(70)
 537 C
 538            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
 539            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
 540            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
 541      &                          JFLG.GE.3) IASG(NSG,3)=1
 542            IHNT2(9)=IHNT2(14)
 543            IHNT2(10)=IHNT2(15)
 544            DO 105 I05=1,5
 545               HINT1(20+I05)=HINT1(40+I05)
 546               HINT1(30+I05)=HINT1(50+I05)
 547  105       CONTINUE
 548            JOUT=1
 549            IF(IHPR2(8).EQ.0) GO TO 160
 550            RRB1=MIN((YP(1,JP)**2+YP(2,JP)**2)/1.2**2
 551      &          /REAL(IHNT2(1))**0.6666667,1.0)
 552            RRB2=MIN((YT(1,JT)**2+YT(2,JT)**2)/1.2**2
 553      &          /REAL(IHNT2(3))**0.6666667,1.0)
 554            APHX1=HIPR1(6)*4.0/3.0*(IHNT2(1)**0.3333333-1.0)
 555      &           *SQRT(1.0-RRB1)
 556            APHX2=HIPR1(6)*4.0/3.0*(IHNT2(3)**0.3333333-1.0)
 557      &           *SQRT(1.0-RRB2)
 558            HINT1(65)=HINT1(61)-APHX1*HINT1(62)
 559      &                  -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
 560            TTRIG=ROMG(R2)*HINT1(65)/HIPR1(31)
 561            NJET=-1
 562 C               ********subtract the trigger jet from total number
 563 C                       of jet production  to be done since it has
 564 C                               already been produced here
 565            XR1=-ALOG(EXP(-TTRIG)+RLU_HIJING(0)*(1.0-EXP(-TTRIG)))
 566  106       NJET=NJET+1
 567            XR1=XR1-ALOG(RLU_HIJING(0))
 568            IF(XR1.LT.TTRIG) GO TO 106
 569            XR=0.0
 570  107       NJET=NJET+1
 571            XR=XR-ALOG(RLU_HIJING(0))
 572            IF(XR.LT.TT-TTRIG) GO TO 107
 573            NJET=NJET-1
 574            GO TO 112
 575         ENDIF
 576 C               ********create a hard interaction with specified P_T
 577 c                                when IHPR2(3)>0
 578         IF(IHPR2(9).EQ.1.AND.JP.EQ.JPMINI.AND.JT.EQ.JTMINI) GO TO 110
 579 C               ********create at least one pair of mini jets
 580 C                       when IHPR2(9)=1
 581 C
 582         IF(IHPR2(8).GT.0 .AND.RNIP(JP,JT).LT.EXP(-TT)*
 583      &          (1.0-EXP(-TTS))) GO TO 160
 584 C               ********this is the probability for no jet production
 585 110     XR=-ALOG(EXP(-TT)+RLU_HIJING(0)*(1.0-EXP(-TT)))
 586 111     NJET=NJET+1
 587         XR=XR-ALOG(RLU_HIJING(0))
 588         IF(XR.LT.TT) GO TO 111
 589 112     NJET=MIN(NJET,IHPR2(8))
 590         IF(IHPR2(8).LT.0)  NJET=ABS(IHPR2(8))
 591 C               ******** Determine number of mini jet production
 592 C
 593         DO 150 I_JET=1,NJET
 594            CALL JETINI(JP,JT,0)
 595            CALL HIJHRD(JP,JT,JOUT,JFLG,1)
 596 C               ********JFLG=1 jets valence quarks, JFLG=2 with
 597 C                       gluon jet, JFLG=3 with q-qbar prod for
 598 C                       (JP,JT). If JFLG=0 jets can not be produced
 599 C                       this time. If JFLG=-1, error occured abandon
 600 C                       this event. JOUT is the total hard scat for
 601 C                       (JP,JT) up to now.
 602            IF(JFLG.EQ.0) GO TO 160
 603            IF(JFLG.LT.0) THEN
 604               IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJHRD'
 605               GO TO 50
 606            ENDIF
 607            JOUT=JOUT+1
 608            IF(ABS(HINT1(46)).GT.HIPR1(11).AND.JFLG.EQ.2) NFP(JP,7)=1
 609            IF(ABS(HINT1(56)).GT.HIPR1(11).AND.JFLG.EQ.2) NFT(JT,7)=1
 610            IF(MAX(ABS(HINT1(46)),ABS(HINT1(56))).GT.HIPR1(11).AND.
 611      &                  JFLG.GE.3) IASG(NSG,3)=1
 612 C               ******** jet with PT>HIPR1(11) will be quenched
 613  150    CONTINUE
 614  160    CONTINUE
 615         CALL HIJSFT(JP,JT,JOUT,IERROR)
 616         IF(IERROR.NE.0) THEN
 617            IF(IHPR2(10).NE.0) WRITE(6,*) 'error occured in HIJSFT'
 618            GO TO 50
 619         ENDIF
 620 C
 621 C               ********conduct soft scattering between JP and JT
 622         JATT=JATT+JOUT
 623
 624 200     CONTINUE
 625
 626 c
 627 c**************************
 628 c
 629         DO 201 JP=1,IHNT2(1)
 630 c           write(6,*) JP, NFP(JP,3), NFP(JP,4), NFP(JP,5)
 631            IF(NFP(JP,5).GT.2) THEN
 632               NINP=NINP+1
 633            ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
 634               NELP=NELP+1
 635            ENDIF
 636
 637            IF(NFP(JP,5).LE.2) THEN
 638               IF (NFP(JP,3) .EQ. 2212) THEN
 639                  NPSPECP = NPSPECP + 1
 640               ELSE IF (NFP(JP,3) .EQ. 2112) THEN
 641                  NNSPECP = NNSPECP + 1
 642               ENDIF
 643            ENDIF
 644  201    continue
 645         DO 202 JT=1,IHNT2(3)
 646            IF(NFT(JT,5).GT.2) THEN
 647               NINT=NINT+1
 648            ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
 649               NELT=NELT+1
 650            ENDIF
 651
 652            IF(NFT(JT,5).LE.2) THEN
 653               IF (NFT(JT,3) .EQ. 2212) THEN
 654                  NPSPECT = NPSPECT + 1
 655               ELSE IF (NFT(JT,3) .EQ. 2112) THEN
 656                  NNSPECT = NNSPECT + 1
 657               ENDIF
 658            ENDIF
 659  202    continue
 660 c
 661 c*******************************
 662
 663 C********perform jet quenching for jets with PT>HIPR1(11)**********
 664
 665         IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
 666      &                  IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
 667                 DO 271 I=1,IHNT2(1)
 668                         IF(NFP(I,7).EQ.1) CALL QUENCH(I,1)
 669 271             CONTINUE
 670                 DO 272 I=1,IHNT2(3)
 671                         IF(NFT(I,7).EQ.1) CALL QUENCH(I,2)
 672 272             CONTINUE
 673                 DO 273 ISG=1,NSG
 674                         IF(IASG(ISG,3).EQ.1) CALL QUENCH(ISG,3)
 675 273             CONTINUE
 676         ENDIF
 677 C
 678 C**************fragment all the string systems in the following*****
 679 C
 680 C********N_ST is where particle information starts
 681 C********N_STR+1 is the number of strings in fragmentation
 682 C********the number of strings before a line is stored in K(I,4)
 683 C********IDSTR is id number of the string system (91,92 or 93)
 684 C
 685 c
 686         IF(IHPR2(20).NE.0) THEN
 687            DO 360 ISG=1,NSG
 688                 CALL HIJFRG(ISG,3,IERROR)
 689                 IF(MSTU(24).NE.0 .OR.IERROR.GT.0) THEN
 690                    MSTU(24)=0
 691                    MSTU(28)=0
 692                    IF(IHPR2(10).NE.0) THEN
 693                       call LULIST_HIJING(1)
 694                       WRITE(6,*) 'error occured, repeat the event'
 695                    ENDIF
 696                    GO TO 50
 697                 ENDIF
 698 C                       ********Check errors
 699 C
 700                 N_ST=1
 701                 IDSTR=92
 702                 IF(IHPR2(21).EQ.0) THEN
 703                    CALL LUEDIT_HIJING(2)
 704                 ELSE
 705 351                N_ST=N_ST+1
 706                    IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  351
 707                    IDSTR=K(N_ST,2)
 708                    N_ST=N_ST+1
 709                 ENDIF
 710 C
 711                 IF(FRAME.EQ.'LAB') THEN
 712                         CALL HIBOOST
 713                 ENDIF
 714 C               ******** boost back to lab frame(if it was in)
 715 C
 716                 N_STR=0
 717                 DO 360 I=N_ST,N
 718                    IF(K(I,2).EQ.IDSTR) THEN
 719                       N_STR=N_STR+1
 720                       GO TO 360
 721                    ENDIF
 722                    K(I,4)=N_STR
 723                    NATT=NATT+1
 724                    KATT(NATT,1)=K(I,2)
 725                    KATT(NATT,2)=20
 726                    KATT(NATT,4)=K(I,1)
 727                    IF(K(I,3).EQ.0 ) THEN
 728                       KATT(NATT,3)=0
 729                    ELSE
 730                       IF(K(K(I,3),2).EQ.IDSTR) THEN
 731                          KATT(NATT,3)=0
 732                       ELSE
 733                          KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
 734                       ENDIF
 735                    ENDIF
 736 C       ****** identify the mother particle
 737                    PATT(NATT,1)=P(I,1)
 738                    PATT(NATT,2)=P(I,2)
 739                    PATT(NATT,3)=P(I,3)
 740                    PATT(NATT,4)=P(I,4)
 741                    VATT(NATT,1)=V(I,1)
 742                    VATT(NATT,2)=V(I,2)
 743                    VATT(NATT,3)=V(I,3)
 744                    VATT(NATT,4)=V(I,4)
 745
 746                    EATT=EATT+P(I,4)
 747 360        CONTINUE
 748 C               ********Fragment the q-qbar jets systems *****
 749 C
 750            JTP(1)=IHNT2(1)
 751            JTP(2)=IHNT2(3)
 752            DO 400 NTP=1,2
 753            DO 400 J_JTP=1,JTP(NTP)
 754                 CALL HIJFRG(J_JTP,NTP,IERROR)
 755                 IF(MSTU(24).NE.0 .OR. IERROR.GT.0) THEN
 756                    MSTU(24)=0
 757                    MSTU(28)=0
 758                    IF(IHPR2(10).NE.0) THEN
 759                       call LULIST_HIJING(1)
 760                       WRITE(6,*) 'error occured, repeat the event'
 761                    ENDIF
 762                    GO TO 50
 763                 ENDIF
 764 C                       ********check errors
 765 C
 766                 N_ST=1
 767                 IDSTR=92
 768
 769                 NFTP=NFP(J_JTP,5)
 770                 IF(NTP.EQ.2) NFTP=10+NFT(J_JTP,5)
 771
 772                 IF(IHPR2(21).EQ.0) THEN
 773                    CALL LUEDIT_HIJING(2)
 774                 ELSE IF (NFTP.EQ. 3 .OR. NFTP .EQ. 13) THEN
 775 381                N_ST=N_ST+1
 776                    IF(K(N_ST,2).LT.91.OR.K(N_ST,2).GT.93) GO TO  381
 777                    IDSTR=K(N_ST,2)
 778                    N_ST=N_ST+1
 779                 ENDIF
 780                 IF(FRAME.EQ.'LAB') THEN
 781                         CALL HIBOOST
 782                 ENDIF
 783 C               ******** boost back to lab frame(if it was in)
 784 C
 785
 786                 N_STR=0
 787                 DO 390 I=N_ST,N
 788                    IF(K(I,2).EQ.IDSTR) THEN
 789                       N_STR=N_STR+1
 790                       GO TO 390
 791                    ENDIF
 792                    K(I,4)=N_STR
 793                    NATT=NATT+1
 794                    KATT(NATT,1)=K(I,2)
 795                    KATT(NATT,2)=NFTP
 796                    KATT(NATT,4)=K(I,1)
 797                    IF(K(I,3).EQ.0 .OR. K(K(I,3),2).EQ.IDSTR) THEN
 798                       KATT(NATT,3)=0
 799                    ELSE
 800                       KATT(NATT,3)=NATT-I+K(I,3)+N_STR-K(K(I,3),4)
 801                    ENDIF
 802 C       ****** identify the mother particle
 803                    PATT(NATT,1)=P(I,1)
 804                    PATT(NATT,2)=P(I,2)
 805                    PATT(NATT,3)=P(I,3)
 806                    PATT(NATT,4)=P(I,4)
 807                    EATT=EATT+P(I,4)
 808                    VATT(NATT,1)=V(I,1)
 809                    VATT(NATT,2)=V(I,2)
 810                    VATT(NATT,3)=V(I,3)
 811                    VATT(NATT,4)=V(I,4)
 812 390             CONTINUE
 813 400        CONTINUE
 814 C               ********Fragment the q-qq related string systems
 815         ENDIF
 816
 817         DO 450 I=1,NDR
 818                 NATT=NATT+1
 819                 KATT(NATT,1)=KFDR(I)
 820                 KATT(NATT,2)=40
 821                 KATT(NATT,3)=0
 822                 PATT(NATT,1)=PDR(I,1)
 823                 PATT(NATT,2)=PDR(I,2)
 824                 PATT(NATT,3)=PDR(I,3)
 825                 PATT(NATT,4)=PDR(I,4)
 826                 VATT(NATT,1)=V(I,1)
 827                 VATT(NATT,2)=V(I,2)
 828                 VATT(NATT,3)=V(I,3)
 829                 VATT(NATT,4)=V(I,4)
 830                 EATT=EATT+PDR(I,4)
 831 450     CONTINUE
 832 C                       ********store the direct-produced particles
 833 C
 834         DENGY=EATT/(IHNT2(1)*HINT1(6)+IHNT2(3)*HINT1(7))-1.0
 835         IF(ABS(DENGY).GT.HIPR1(43).AND.IHPR2(20).NE.0
 836      &     .AND.IHPR2(21).EQ.0) THEN
 837         IF(IHPR2(10).NE.0) WRITE(6,*) 'Energy not conserved, repeat the
 838      &     event'
 839 C               call LULIST_HIJING(1)
 840                 GO TO 50
 841         ENDIF
 842
 843         RETURN
 844         END