EIC code displayed by LXR master/​geant4/​examples/​extended/​hadronic/​Hadr02/​hijing1_383/​hijing1.383.f	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-04-04 08:05:07

 c     Version 1.383
 c     The variables I_SNG in HIJSFT and JL in ATTRAD were not initialized.
 c     The version initialize them. (as found by Fernando Marroquim)
 c
 c
 c
 c     Version 1.382
 c     Nuclear distribution for deuteron is taken as the Hulthen wave
 c     function as provided by Brian Cole (Columbia)
 c
 c
 c     Version 1.381
 c
 c     The parameters for Wood-Saxon distribution for deuteron are
 c     constrained to give the right rms ratius 2.116 fm
 c     (R=0.0, D=0.5882)
 c
 c
 c     Version 1.38
 c
 c     The following common block is added to record the number of elastic
 c     (NELT, NELP) and inelastic (NINT, NINP) participants
 c
 c        COMMON/HIJGLBR/NELT,NINT,NELP,NINP
 c        SAVE  /HIJGLBR/
 c
 c     Version 1.37
 c
 c     A bug in the quenching subroutine is corrected. When calculating the
 c     distance between two wounded nucleons, the displacement of the
 c     impact parameter was not inculded. This bug was discovered by
 c     Dr. V.Uzhinskii JINR, Dubna, Russia
 c
 c
 C     Version 1.36
 c
 c     Modification Oct. 8, 1998. In hijing, log(ran(nseed)) occasionally
 c     causes overfloat. It is modified to log(max(ran(nseed),1.0e-20)).
 c
 c
 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 mistakes are corrected in
 C     function G. A Mistake 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****************************************************************
 c       SUBROUTINE HIJING(FRAME,BMIN0,BMAX0)   !khaled
         SUBROUTINE HIJING(BMIN0,BMAX0)
 
         CHARACTER FRAME*8
         DIMENSION SCIP(300,300),RNIP(300,300),SJIP(300,300),JTP(3),
      &                  IPCOL(90000),ITCOL(90000)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
 C
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIJGLBR/NELT,NINT,NELP,NINP
         SAVE  /HIJGLBR/
         COMMON/HIMAIN1/NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
         SAVE  /HIMAIN1/
         COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
         SAVE  /HIMAIN2/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &                PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
         COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
         SAVE  /HIJJET4/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)   
         SAVE  /LUJETS/
         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
         SAVE  /LUDAT1/
 
         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))
 c        write(*,*)"HIJING bmin bmax",BMIN,"=====", BMAX
         ENDIF
 c       ----------------------------
         FRAME='CMS'      !khaled    //CMS
 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)
 c
         if(IHNT2(1).EQ.2) then
            rnd1=max(RLU(0),1.0e-20)
            rnd2=max(RLU(0),1.0e-20)
            rnd3=max(RLU(0),1.0e-20)
            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
         endif
 c
         X=RLU(0)
         CX=2.0*X-1.0
         SX=SQRT(1.0-CX*CX)
 C               ********choose theta from uniform cos(theta) distr
         PHI=RLU(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
 c*******************************
         if(IHNT2(1).EQ.2) then
            YP(1,2)=-YP(1,1)
            YP(2,2)=-YP(2,1)
            YP(3,2)=-YP(3,1)
         endif
 c********************************
         DO 12 I=1,IHNT2(1)-1
         DO 12 J=I+1,IHNT2(1)
         IF(YP(3,I).GT.YP(3,J)) GO TO 12
         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)
 c
         if(IHNT2(3).EQ.2) then
            rnd1=max(RLU(0),1.0e-20)
            rnd2=max(RLU(0),1.0e-20)
            rnd3=max(RLU(0),1.0e-20)
            R=-0.5*(log(rnd1)*4.38/2.0+log(rnd2)*0.85/2.0
      &          +4.38*0.85*log(rnd3)/(4.38+0.85))
         endif
 c
         X=RLU(0)
         CX=2.0*X-1.0
         SX=SQRT(1.0-CX*CX)
 C               ********choose theta from uniform cos(theta) distr
         PHI=RLU(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
 c**********************************
         if(IHNT2(3).EQ.2) then
            YT(1,2)=-YT(1,1)
            YT(2,2)=-YT(2,1)
            YT(3,2)=-YT(3,1)
         endif
 c*********************************
         DO 22 I=1,IHNT2(3)-1
         DO 22 J=I+1,IHNT2(3)
         IF(YT(3,I).LT.YT(3,J)) GO TO 22
         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
         NELT=0
         NINT=0
         NELP=0
         NINP=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(0)*(BMAX**2-BMIN**2))
 c        write(*,*)"HIJING bmin bmax bb", bmin, bmax,bb
         PHI=2.0*HIPR1(40)*RLU(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
 c           write(*,*)'R2 ROMG(R2)',R2, ROMG(R2), OMG0(R2)  KHal
               GS=1.0-EXP(-(HIPR1(30)+HINT1(18))*ROMG(R2)/HIPR1(31))
               RANTOT=RLU(0)
               
               IF(RANTOT.GT.GS) GO TO 70
               GO TO 65
            ENDIF
 c           write(*,*)'ROMG(0.0)'  !kha
            GSTOT_0=2.0*(1.0-EXP(-(HIPR1(30)+HINT1(18))
      &             /HIPR1(31)/2.0*ROMG(0.0)))
 c           write(*,*)GSTOT_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(0)*GSTOT_0
 c           write(*,*)'rantot, gstot',rantot, gstot
            IF(RANTOT.GT.GSTOT) GO TO 70
            IF(RANTOT.GT.GS) THEN
 c              write(*,*)"----------------------------------start HIJCSC"
               CALL HIJCSC(JP,JT)
 c              write(*,*)"----------------------------------end HIJCSC"
               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(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(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
 c
         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(0)*(1.0-EXP(-TTRIG)))
  106       NJET=NJET+1
            XR1=XR1-ALOG(max(RLU(0),1.0e-20))
            IF(XR1.LT.TTRIG) GO TO 106
            XR=0.0
  107       NJET=NJET+1
            XR=XR-ALOG(max(RLU(0),1.0e-20))
            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(0)*(1.0-EXP(-TT)))
 111     NJET=NJET+1
         XR=XR-ALOG(max(RLU(0),1.0e-20))
         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
 c**************************
 c
         DO 201 JP=1,IHNT2(1)
            IF(NFP(JP,5).GT.2) THEN
               NINP=NINP+1
            ELSE IF(NFP(JP,5).EQ.2.OR.NFP(JP,5).EQ.1) THEN
               NELP=NELP+1
            ENDIF
  201    continue
         DO 202 JT=1,IHNT2(3)
            IF(NFT(JT,5).GT.2) THEN
               NINT=NINT+1
            ELSE IF(NFT(JT,5).EQ.2.OR.NFT(JT,5).EQ.1) THEN
               NELT=NELT+1
            ENDIF
  202    continue
 c     
 c*******************************
 
 
 C********perform jet quenching for jets with PT>HIPR1(11)**********
 
         IF((IHPR2(8).NE.0.OR.IHPR2(3).NE.0).AND.IHPR2(4).GT.0.AND.
      &                  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(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(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)
                    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(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(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)
                    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)
                 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) THEN
         WRITE(6,*) 'Energy not conserved, repeat the event'
         ENDIF
 C               call lulist(1)
                 GO TO 50
         ENDIF
         RETURN
         END
 C
 C
 C
 c       SUBROUTINE HIJSET(EFRM,FRAME,PROJ,TARG,IAP,IZP,IAT,IZT)  !khaled
 c        SUBROUTINE HIJSET(EFRM,IAP,IZP,IAT,IZT)  !khaled
         SUBROUTINE hijset(efrm)
 
 c        real*4 efrm
         CHARACTER FRAME*4,PROJ*4,TARG*4,EFRAME*4
         DOUBLE PRECISION  DD1,DD2,DD3,DD4
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
         SAVE  /HIJDAT/
         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
         SAVE  /LUDAT1/
         EXTERNAL FNKICK,FNKICK2,FNSTRU,FNSTRUM,FNSTRUS
         CALL TITLE
 c       ---------->khaled
         FRAME='CMS'   
 c        FRAME='LAB'
         PROJ="A"
         TARG="B"
 c        write(*,*)"Fortran    khaled"
 c       write(*,*)EFRM, IHNT2(1),IHNT2(2),IHNT2(3),IHNT2(4)
 c        write(*,*)IHNT2(5), IHNT2(6), HINT1(8), HINT1(9)
 c        write(*,*)iap, izp, iat,izt
 c        write(*,*)hipr1(1), hipr1(2),hipr1(3),hipr1(4), ihpr2(12)
 c      ---------> closed by khaled-------------->
 c          ------------------------------------
 c       IHNT2(1)=IAP
 c       IHNT2(2)=IZP
 c       IHNT2(3)=IAT
 c       IHNT2(4)=IZT
 c       IHNT2(5)=0
 c       IHNT2(6)=0
 C
 
         
 c       HINT1(8)=MAX(ULMASS(2112),ULMASS(2212))
 c       HINT1(9)=HINT1(8)
 C
 c      ---------> closed by khaled
 c       IF(PROJ.NE.'A') THEN
 c               IF(PROJ.EQ.'P') THEN
 c                   IHNT2(5)=2212
 c               ELSE IF(PROJ.EQ.'PBAR') THEN 
 c                   IHNT2(5)=-2212
 c               ELSE IF(PROJ.EQ.'PI+') THEN
 c                   IHNT2(5)=211
 c               ELSE IF(PROJ.EQ.'PI-') THEN
 c                   IHNT2(5)=-211
 c               ELSE IF(PROJ.EQ.'K+') THEN
 c                   IHNT2(5)=321
 c               ELSE IF(PROJ.EQ.'K-') THEN
 c                   IHNT2(5)=-321
 c               ELSE IF(PROJ.EQ.'N') THEN
 c                   IHNT2(5)=2112
 c               ELSE IF(PROJ.EQ.'NBAR') THEN
 c                   IHNT2(5)=-2112
 c               ELSE
 c                   WRITE(6,*) PROJ, 'wrong or unavailable proj name'
 c                   STOP
 c               ENDIF
 c               HINT1(8)=ULMASS(IHNT2(5))
 c       ENDIF
 c       IF(TARG.NE.'A') THEN
 c               IF(TARG.EQ.'P') THEN
 c                   IHNT2(6)=2212
 c               ELSE IF(TARG.EQ.'PBAR') THEN 
 c                   IHNT2(6)=-2212
 c               ELSE IF(TARG.EQ.'PI+') THEN
 c                   IHNT2(6)=211
 c               ELSE IF(TARG.EQ.'PI-') THEN
 c                   IHNT2(6)=-211
 c               ELSE IF(TARG.EQ.'K+') THEN
 c                   IHNT2(6)=321
 c               ELSE IF(TARG.EQ.'K-') THEN
 c                   IHNT2(6)=-321
 c               ELSE IF(TARG.EQ.'N') THEN
 c                   IHNT2(6)=2112
 c               ELSE IF(TARG.EQ.'NBAR') THEN
 c                   IHNT2(6)=-2112
 c               ELSE
 c                   WRITE(6,*) TARG,'wrong or unavailable targ name'
 c                   STOP
 c               ENDIF
 c               HINT1(9)=ULMASS(IHNT2(6))
 c       ENDIF
 c         ------------------------>end close by khaled
 
 C...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
         IF(IHPR2(12).GT.0) THEN
         CALL LUGIVE('MDCY(C111,1)=0')
         CALL LUGIVE('MDCY(C310,1)=0')
            CALL LUGIVE('MDCY(C411,1)=0;MDCY(C-411,1)=0')
            CALL LUGIVE('MDCY(C421,1)=0;MDCY(C-421,1)=0')
            CALL LUGIVE('MDCY(C431,1)=0;MDCY(C-431,1)=0')
            CALL LUGIVE('MDCY(C511,1)=0;MDCY(C-511,1)=0')
            CALL LUGIVE('MDCY(C521,1)=0;MDCY(C-521,1)=0')
            CALL LUGIVE('MDCY(C531,1)=0;MDCY(C-531,1)=0')
         CALL LUGIVE('MDCY(C3122,1)=0;MDCY(C-3122,1)=0')
         CALL LUGIVE('MDCY(C3112,1)=0;MDCY(C-3112,1)=0')
         CALL LUGIVE('MDCY(C3212,1)=0;MDCY(C-3212,1)=0')
         CALL LUGIVE('MDCY(C3222,1)=0;MDCY(C-3222,1)=0')
         CALL LUGIVE('MDCY(C3312,1)=0;MDCY(C-3312,1)=0')
         CALL LUGIVE('MDCY(C3322,1)=0;MDCY(C-3322,1)=0')
         CALL LUGIVE('MDCY(C3334,1)=0;MDCY(C-3334,1)=0')
         ENDIF
         MSTU(12)=0
         MSTU(21)=1
         IF(IHPR2(10).EQ.0) THEN
                 MSTU(22)=0
                 MSTU(25)=0
                 MSTU(26)=0
         ENDIF
         MSTJ(12)=IHPR2(11)
         PARJ(21)=HIPR1(2)
         PARJ(41)=HIPR1(3)
         PARJ(42)=HIPR1(4)
 C                       ******** set up for jetset
         IF(FRAME.EQ.'LAB') THEN
            DD1=EFRM
            DD2=HINT1(8)
            DD3=HINT1(9)
            HINT1(1)=SQRT(HINT1(8)**2+2.0*HINT1(9)*EFRM+HINT1(9)**2)
            DD4=DSQRT(DD1**2-DD2**2)/(DD1+DD3)
            HINT1(2)=DD4
            HINT1(3)=0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
            DD4=DSQRT(DD1**2-DD2**2)/DD1
            HINT1(4)=0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
            HINT1(5)=0.0
            HINT1(6)=EFRM
            HINT1(7)=HINT1(9)
         ELSE IF(FRAME.EQ.'CMS') THEN
            HINT1(1)=EFRM
            HINT1(2)=0.0
            HINT1(3)=0.0
            DD1=HINT1(1)
            DD2=HINT1(8)
            DD3=HINT1(9)
            DD4=DSQRT(1.D0-4.D0*DD2**2/DD1**2)
            HINT1(4)=0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
            DD4=DSQRT(1.D0-4.D0*DD3**2/DD1**2)
            HINT1(5)=-0.5*DLOG((1.D0+DD4)/(1.D0-DD4))
            HINT1(6)=HINT1(1)/2.0
            HINT1(7)=HINT1(1)/2.0
         ENDIF
 C               ********define Lorentz transform to lab frame
 c
 C               ********calculate the cross sections involved with
 C                       nucleon collisions.
         IF(IHNT2(1).GT.1) THEN
                 CALL HIJWDS(IHNT2(1),1,RMAX)
                 HIPR1(34)=RMAX
 C                       ********set up Wood-Sax distr for proj.
         ENDIF
         IF(IHNT2(3).GT.1) THEN
                 CALL HIJWDS(IHNT2(3),2,RMAX)
                 HIPR1(35)=RMAX
 C                       ********set up Wood-Sax distr for  targ.
         ENDIF
 c        write(*,*)"radii hipr1(34--35)", hipr1(34),"<------>", hipr1(35)  !khaled
 c        write(*,*)"hint1(6)--Hint1(7)", hint1(6),"<------>", hint1(7)  !khaled
 C
 C
         I=0
 20      I=I+1
         IF(I.EQ.10) GO TO 30
         IF(HIDAT0(10,I).LE.HINT1(1)) GO TO 20
 30      IF(I.EQ.1) I=2
         DO 40 J=1,9
            HIDAT(J)=HIDAT0(J,I-1)+(HIDAT0(J,I)-HIDAT0(J,I-1))
      &     *(HINT1(1)-HIDAT0(10,I-1))/(HIDAT0(10,I)-HIDAT0(10,I-1))
 40      CONTINUE
         HIPR1(31)=HIDAT(5)
         HIPR1(30)=2.0*HIDAT(5)
 c       write(*,*)"hipr1(31--30)", hipr1(31),"<------>", hipr1(30)  !khaled 
 C
 C
         
          CALL HIJCRS
        
 C
         IF(IHPR2(5).NE.0) THEN
                 CALL HIFUN(3,0.0,36.0,FNKICK)
 C               ********booking for generating pt**2 for pt kick
         ENDIF
         CALL HIFUN(7,0.0,6.0,FNKICK2)
         CALL HIFUN(4,0.0,1.0,FNSTRU)
         CALL HIFUN(5,0.0,1.0,FNSTRUM)
         CALL HIFUN(6,0.0,1.0,FNSTRUS)
 C               ********booking for x distribution of valence quarks
         EFRAME='Ecm'
         IF(FRAME.EQ.'LAB') EFRAME='Elab'
         WRITE(6,100) EFRAME,EFRM,PROJ,IHNT2(1),IHNT2(2),
      &               TARG,IHNT2(3),IHNT2(4) 
 100     FORMAT(//10X,'****************************************
      &  **********'/
      &  10X,'*',48X,'*'/
      &  10X,'*         HIJING has been initialized at         *'/
      &  10X,'*',13X,A4,'= ',F10.2,' GeV/n',13X,'*'/
      &  10X,'*',48X,'*'/
      &  10X,'*',8X,'for ',
      &  A4,'(',I3,',',I3,')',' + ',A4,'(',I3,',',I3,')',7X,'*'/
      &  10X,'**************************************************')
         RETURN
         END
 C
 C
 C
         FUNCTION FNKICK(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         FNKICK=1.0/(X+HIPR1(19)**2)/(X+HIPR1(20)**2)
      &          /(1+EXP((SQRT(X)-HIPR1(20))/0.4))
         RETURN
         END
 C
 C
         FUNCTION FNKICK2(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         FNKICK2=X*EXP(-2.0*X/HIPR1(42))
         RETURN
         END
 C
 C
 C
         FUNCTION FNSTRU(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         FNSTRU=(1.0-X)**HIPR1(44)/
      &          (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
         RETURN
         END
 C
 C
 C
         FUNCTION FNSTRUM(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         FNSTRUM=1.0/((1.0-X)**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
      &          /(X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(46)
         RETURN
         END
 C
 C
         FUNCTION FNSTRUS(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         FNSTRUS=(1.0-X)**HIPR1(47)/
      &          (X**2+HIPR1(45)**2/HINT1(1)**2)**HIPR1(48)
         RETURN
         END
 C
 C
 C
 C
         SUBROUTINE HIBOOST
         IMPLICIT DOUBLE PRECISION(D)  
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5) 
         SAVE /LUJETS/ 
         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
         SAVE /LUDAT1/ 
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         DO 100 I=1,N
            DBETA=P(I,3)/P(I,4)
            IF(ABS(DBETA).GE.1.D0) THEN
               DB=HINT1(2)
               IF(DB.GT.0.99999999D0) THEN 
 C               ********Rescale boost vector if too close to unity. 
                  WRITE(6,*) '(HIBOOT:) boost vector too large' 
                  DB=0.99999999D0
               ENDIF 
               DGA=1D0/SQRT(1D0-DB**2)
               DP3=P(I,3)
               DP4=P(I,4)
               P(I,3)=(DP3+DB*DP4)*DGA  
               P(I,4)=(DP4+DB*DP3)*DGA  
               GO TO 100
            ENDIF
            Y=0.5*DLOG((1.D0+DBETA)/(1.D0-DBETA))
            AMT=SQRT(P(I,1)**2+P(I,2)**2+P(I,5)**2)
            P(I,3)=AMT*SINH(Y+HINT1(3))
            P(I,4)=AMT*COSH(Y+HINT1(3))
 100     CONTINUE
         RETURN
         END
 C
 C
 C
 C
         SUBROUTINE QUENCH(JPJT,NTP)
         DIMENSION RDP(300),LQP(300),RDT(300),LQT(300)
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
 C
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &                PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C
         BB=HINT1(19)                                            ! Uzhi
         PHI=HINT1(20)                                           ! Uzhi
         BBX=BB*COS(PHI)                                         ! Uzhi
         BBY=BB*SIN(PHI)                                         ! Uzhi
 c
         IF(NTP.EQ.2) GO TO 400
         IF(NTP.EQ.3) GO TO 2000 
 C*******************************************************
 C Jet interaction for proj jet in the direction PHIP
 C******************************************************
 C
         IF(NFP(JPJT,7).NE.1) RETURN
 
         JP=JPJT
         DO 290 I=1,NPJ(JP)
            PTJET0=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2)
            IF(PTJET0.LE.HIPR1(11)) GO TO 290
            PTOT=SQRT(PTJET0*PTJET0+PJPZ(JP,I)**2)
            IF(PTOT.LT.HIPR1(8)) GO TO 290
            PHIP=ULANGL(PJPX(JP,I),PJPY(JP,I))
 C******* find the wounded proj which can interact with jet***
            KP=0
            DO 100 I2=1,IHNT2(1)
               IF(NFP(I2,5).NE.3 .OR. I2.EQ.JP) GO TO 100
               DX=YP(1,I2)-YP(1,JP)
               DY=YP(2,I2)-YP(2,JP)
               PHI=ULANGL(DX,DY)
               DPHI=ABS(PHI-PHIP)
               IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI      ! Uzhi
               IF(DPHI.GE.HIPR1(40)/2.0) GO TO 100
               RD0=SQRT(DX*DX+DY*DY)
               IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 100
               KP=KP+1
               LQP(KP)=I2
               RDP(KP)=COS(DPHI)*RD0
  100       CONTINUE
 C*******        rearrange according decending rd************
            DO 110 I2=1,KP-1
               DO 110 J2=I2+1,KP
                  IF(RDP(I2).LT.RDP(J2)) GO TO 110
                  RD=RDP(I2)
                  LQ=LQP(I2)
                  RDP(I2)=RDP(J2)
                  LQP(I2)=LQP(J2)
                  RDP(J2)=RD
                  LQP(J2)=LQ
  110          CONTINUE
 C****** find wounded targ which can interact with jet********
               KT=0
               DO 120 I2=1,IHNT2(3)
                  IF(NFT(I2,5).NE.3) GO TO 120
                  DX=YT(1,I2)-YP(1,JP)-BBX
                  DY=YT(2,I2)-YP(2,JP)-BBY
                  PHI=ULANGL(DX,DY)
                  DPHI=ABS(PHI-PHIP)
                  IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI   ! Uzhi
                  IF(DPHI.GT.HIPR1(40)/2.0) GO TO 120
                  RD0=SQRT(DX*DX+DY*DY)
                  IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 120
                  KT=KT+1
                  LQT(KT)=I2
                  RDT(KT)=COS(DPHI)*RD0
  120          CONTINUE
 C*******        rearrange according decending rd************
               DO 130 I2=1,KT-1
                  DO 130 J2=I2+1,KT
                     IF(RDT(I2).LT.RDT(J2)) GO TO 130
                     RD=RDT(I2)
                     LQ=LQT(I2)
                     RDT(I2)=RDT(J2)
                     LQT(I2)=LQT(J2)
                     RDT(J2)=RD
                     LQT(J2)=LQ
  130             CONTINUE
                 
                  MP=0
                  MT=0
                  R0=0.0
                  NQ=0
                  DP=0.0
                  PTOT=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2+PJPZ(JP,I)**2)
                  V1=PJPX(JP,I)/PTOT
                  V2=PJPY(JP,I)/PTOT
                  V3=PJPZ(JP,I)/PTOT
 
  200             RN=RLU(0)
  210             IF(MT.GE.KT .AND. MP.GE.KP) GO TO 290
                  IF(MT.GE.KT) GO TO 220
                  IF(MP.GE.KP) GO TO 240
                  IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 240
  220             MP=MP+1
                  DRR=RDP(MP)-R0
                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
                  DP=DRR*HIPR1(14)
                  IF(KFPJ(JP,I).NE.21) DP=0.5*DP
 C       ********string tension of quark jet is 0.5 of gluon's 
                  IF(DP.LE.0.2) GO TO 210
                  IF(PTOT.LE.0.4) GO TO 290
                  IF(PTOT.LE.DP) DP=PTOT-0.2
                  DE=DP
 
                  IF(KFPJ(JP,I).NE.21) THEN
                     PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
      &                   +PP(LQP(MP),3)**2
                     DE=SQRT(PJPM(JP,I)**2+PTOT**2)
      &                  -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
                     ERSHU=(PP(LQP(MP),4)+DE-DP)**2
                     AMSHU=ERSHU-PRSHU
                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
                     PP(LQP(MP),4)=SQRT(ERSHU)
                     PP(LQP(MP),5)=SQRT(AMSHU)
                  ENDIF
 C               ********reshuffle the energy when jet has mass
                  R0=RDP(MP)
                  DP1=DP*V1
                  DP2=DP*V2
                  DP3=DP*V3
 C               ********momentum and energy transfer from jet
                  
                  NPJ(LQP(MP))=NPJ(LQP(MP))+1
                  KFPJ(LQP(MP),NPJ(LQP(MP)))=21
                  PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
                  PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
                  PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
                  PJPE(LQP(MP),NPJ(LQP(MP)))=DP
                  PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
                  GO TO 260
 
  240             MT=MT+1
                  DRR=RDT(MT)-R0
                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
                  DP=DRR*HIPR1(14)
                  IF(DP.LE.0.2) GO TO 210
                  IF(PTOT.LE.0.4) GO TO 290
                  IF(PTOT.LE.DP) DP=PTOT-0.2
                  DE=DP
 
                  IF(KFPJ(JP,I).NE.21) THEN
                     PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
      &                   +PT(LQT(MT),3)**2
                     DE=SQRT(PJPM(JP,I)**2+PTOT**2)
      &                  -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
                     ERSHU=(PT(LQT(MT),4)+DE-DP)**2
                     AMSHU=ERSHU-PRSHU
                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
                     PT(LQT(MT),4)=SQRT(ERSHU)
                     PT(LQT(MT),5)=SQRT(AMSHU)
                  ENDIF
 C               ********reshuffle the energy when jet has mass
 
                  R0=RDT(MT)
                  DP1=DP*V1
                  DP2=DP*V2
                  DP3=DP*V3
 C               ********momentum and energy transfer from jet
                  NTJ(LQT(MT))=NTJ(LQT(MT))+1
                  KFTJ(LQT(MT),NTJ(LQT(MT)))=21
                  PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
                  PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
                  PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
                  PJTE(LQT(MT),NTJ(LQT(MT)))=DP
                  PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
 
  260             PJPX(JP,I)=(PTOT-DP)*V1
                  PJPY(JP,I)=(PTOT-DP)*V2
                  PJPZ(JP,I)=(PTOT-DP)*V3
                  PJPE(JP,I)=PJPE(JP,I)-DE
 
                  PTOT=PTOT-DP
                  NQ=NQ+1
                  GO TO 200
  290          CONTINUE
 
               RETURN
 
 C*******************************************************
 C Jet interaction for target jet in the direction PHIT
 C******************************************************
 C
 C******* find the wounded proj which can interact with jet***
 
  400          IF(NFT(JPJT,7).NE.1) RETURN
               JT=JPJT
               DO 690 I=1,NTJ(JT)
                  PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
                  IF(PTJET0.LE.HIPR1(11)) GO TO 690
                  PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
                  IF(PTOT.LT.HIPR1(8)) GO TO 690
                  PHIT=ULANGL(PJTX(JT,I),PJTY(JT,I))
                  KP=0
                  DO 500 I2=1,IHNT2(1)
                     IF(NFP(I2,5).NE.3) GO TO 500
                     DX=YP(1,I2)+BBX-YT(1,JT)
                     DY=YP(2,I2)+BBY-YT(2,JT)
                     PHI=ULANGL(DX,DY)
                     DPHI=ABS(PHI-PHIT)
                     IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
                     IF(DPHI.GT.HIPR1(40)/2.0) GO TO 500
                     RD0=SQRT(DX*DX+DY*DY)
                     IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 500
                     KP=KP+1
                     LQP(KP)=I2
                     RDP(KP)=COS(DPHI)*RD0
  500             CONTINUE
 C*******        rearrange according to decending rd************
                  DO 510 I2=1,KP-1
                     DO 510 J2=I2+1,KP
                        IF(RDP(I2).LT.RDP(J2)) GO TO 510
                        RD=RDP(I2)
                        LQ=LQP(I2)
                        RDP(I2)=RDP(J2)
                        LQP(I2)=LQP(J2)
                        RDP(J2)=RD
                        LQP(J2)=LQ
  510                CONTINUE
 C****** find wounded targ which can interact with jet********
                     KT=0
                     DO 520 I2=1,IHNT2(3)
                        IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 520
                        DX=YT(1,I2)-YT(1,JT)
                        DY=YT(2,I2)-YT(2,JT)
                        PHI=ULANGL(DX,DY)
                        DPHI=ABS(PHI-PHIT)
                        IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
                        IF(DPHI.GT.HIPR1(40)/2.0) GO TO 520
                        RD0=SQRT(DX*DX+DY*DY)
                        IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 520
                        KT=KT+1
                        LQT(KT)=I2
                        RDT(KT)=COS(DPHI)*RD0
  520                CONTINUE
 C*******        rearrange according to decending rd************
                     DO 530 I2=1,KT-1
                        DO 530 J2=I2+1,KT
                           IF(RDT(I2).LT.RDT(J2)) GO TO 530
                           RD=RDT(I2)
                           LQ=LQT(I2)
                           RDT(I2)=RDT(J2)
                           LQT(I2)=LQT(J2)
                           RDT(J2)=RD
                           LQT(J2)=LQ
  530                   CONTINUE
                        
                        MP=0
                        MT=0
                        NQ=0
                        DP=0.0
                        R0=0.0
                 PTOT=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2+PJTZ(JT,I)**2)
                 V1=PJTX(JT,I)/PTOT
                 V2=PJTY(JT,I)/PTOT
                 V3=PJTZ(JT,I)/PTOT
 
  600            RN=RLU(0)
  610            IF(MT.GE.KT .AND. MP.GE.KP) GO TO 690
                 IF(MT.GE.KT) GO TO 620
                 IF(MP.GE.KP) GO TO 640
                 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 640
 620             MP=MP+1
                 DRR=RDP(MP)-R0
                 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
                 DP=DRR*HIPR1(14)
                 IF(KFTJ(JT,I).NE.21) DP=0.5*DP
 C       ********string tension of quark jet is 0.5 of gluon's 
                 IF(DP.LE.0.2) GO TO 610
                 IF(PTOT.LE.0.4) GO TO 690
                 IF(PTOT.LE.DP) DP=PTOT-0.2
                 DE=DP
 C
                 IF(KFTJ(JT,I).NE.21) THEN
                    PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
      &                   +PP(LQP(MP),3)**2
                    DE=SQRT(PJTM(JT,I)**2+PTOT**2)
      &               -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
                    ERSHU=(PP(LQP(MP),4)+DE-DP)**2
                    AMSHU=ERSHU-PRSHU
                    IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
                    PP(LQP(MP),4)=SQRT(ERSHU)
                    PP(LQP(MP),5)=SQRT(AMSHU)
                 ENDIF
 C               ********reshuffle the energy when jet has mass
 C
                 R0=RDP(MP)
                 DP1=DP*V1
                 DP2=DP*V2
                 DP3=DP*V3
 C               ********momentum and energy transfer from jet
                 NPJ(LQP(MP))=NPJ(LQP(MP))+1
                 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
                 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
                 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
                 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
                 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
                 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
 
                 GO TO 660
 
 640             MT=MT+1
                 DRR=RDT(MT)-R0
                 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
                 DP=DRR*HIPR1(14)
                 IF(DP.LE.0.2) GO TO 610
                 IF(PTOT.LE.0.4) GO TO 690
                 IF(PTOT.LE.DP) DP=PTOT-0.2
                 DE=DP
 
                 IF(KFTJ(JT,I).NE.21) THEN
                    PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
      &                   +PT(LQT(MT),3)**2
                    DE=SQRT(PJTM(JT,I)**2+PTOT**2)
      &               -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
                    ERSHU=(PT(LQT(MT),4)+DE-DP)**2
                    AMSHU=ERSHU-PRSHU
                    IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
                    PT(LQT(MT),4)=SQRT(ERSHU)
                    PT(LQT(MT),5)=SQRT(AMSHU)
                 ENDIF
 C               ********reshuffle the energy when jet has mass
 
                 R0=RDT(MT)
                 DP1=DP*V1
                 DP2=DP*V2
                 DP3=DP*V3
 C               ********momentum and energy transfer from jet
                 NTJ(LQT(MT))=NTJ(LQT(MT))+1
                 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
                 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
                 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
                 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
                 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
                 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
 
 660             PJTX(JT,I)=(PTOT-DP)*V1
                 PJTY(JT,I)=(PTOT-DP)*V2
                 PJTZ(JT,I)=(PTOT-DP)*V3
                 PJTE(JT,I)=PJTE(JT,I)-DE
 
                 PTOT=PTOT-DP
                 NQ=NQ+1
                 GO TO 600
 690     CONTINUE
         RETURN
 C********************************************************
 C       Q-QBAR jet interaction
 C********************************************************
 2000    ISG=JPJT
         IF(IASG(ISG,3).NE.1) RETURN
 C
         JP=IASG(ISG,1)
         JT=IASG(ISG,2)
         XJ=(YP(1,JP)+BBX+YT(1,JT))/2.0
         YJ=(YP(2,JP)+BBY+YT(2,JT))/2.0
         DO 2690 I=1,NJSG(ISG)
            PTJET0=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2)
            IF(PTJET0.LE.HIPR1(11).OR.PESG(ISG,I).LT.HIPR1(1))
      &            GO TO 2690
            PTOT=SQRT(PTJET0*PTJET0+PZSG(ISG,I)**2)
            IF(PTOT.LT.MAX(HIPR1(1),HIPR1(8))) GO TO 2690
            PHIQ=ULANGL(PXSG(ISG,I),PYSG(ISG,I))
            KP=0
            DO 2500 I2=1,IHNT2(1)
               IF(NFP(I2,5).NE.3.OR.I2.EQ.JP) GO TO 2500
               DX=YP(1,I2)+BBX-XJ
               DY=YP(2,I2)+BBY-YJ
               PHI=ULANGL(DX,DY)
               DPHI=ABS(PHI-PHIQ)
               IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI      ! Uzhi
               IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2500
               RD0=SQRT(DX*DX+DY*DY)
               IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2500
               KP=KP+1
               LQP(KP)=I2
               RDP(KP)=COS(DPHI)*RD0
  2500      CONTINUE
 C*******        rearrange according to decending rd************
            DO 2510 I2=1,KP-1
               DO 2510 J2=I2+1,KP
                  IF(RDP(I2).LT.RDP(J2)) GO TO 2510
                  RD=RDP(I2)
                  LQ=LQP(I2)
                  RDP(I2)=RDP(J2)
                  LQP(I2)=LQP(J2)
                  RDP(J2)=RD
                  LQP(J2)=LQ
  2510         CONTINUE
 C****** find wounded targ which can interact with jet********
               KT=0
               DO 2520 I2=1,IHNT2(3)
                  IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 2520
                  DX=YT(1,I2)-XJ
                  DY=YT(2,I2)-YJ
                  PHI=ULANGL(DX,DY)
                  DPHI=ABS(PHI-PHIQ)
                  IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
                  IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2520
                  RD0=SQRT(DX*DX+DY*DY)
                  IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2520
                  KT=KT+1
                  LQT(KT)=I2
                  RDT(KT)=COS(DPHI)*RD0
  2520         CONTINUE
 C*******        rearrange according to decending rd************
               DO 2530 I2=1,KT-1
                  DO 2530 J2=I2+1,KT
                     IF(RDT(I2).LT.RDT(J2)) GO TO 2530
                     RD=RDT(I2)
                     LQ=LQT(I2)
                     RDT(I2)=RDT(J2)
                     LQT(I2)=LQT(J2)
                     RDT(J2)=RD
                     LQT(J2)=LQ
  2530            CONTINUE
                 
                  MP=0
                  MT=0
                  NQ=0
                  DP=0.0
                  R0=0.0
                  PTOT=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2
      &                +PZSG(ISG,I)**2)
                  V1=PXSG(ISG,I)/PTOT
                  V2=PYSG(ISG,I)/PTOT
                  V3=PZSG(ISG,I)/PTOT
 
  2600            RN=RLU(0)
  2610            IF(MT.GE.KT .AND. MP.GE.KP) GO TO 2690
                  IF(MT.GE.KT) GO TO 2620
                  IF(MP.GE.KP) GO TO 2640
                  IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 2640
  2620            MP=MP+1
                  DRR=RDP(MP)-R0
                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
                  DP=DRR*HIPR1(14)/2.0
                  IF(DP.LE.0.2) GO TO 2610
                  IF(PTOT.LE.0.4) GO TO 2690
                  IF(PTOT.LE.DP) DP=PTOT-0.2
                  DE=DP
 C
                  IF(K2SG(ISG,I).NE.21) THEN
                     IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
                     PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
      &                    +PP(LQP(MP),3)**2
                     DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
      &                 -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
                     ERSHU=(PP(LQP(MP),4)+DE-DP)**2
                     AMSHU=ERSHU-PRSHU
                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
                     PP(LQP(MP),4)=SQRT(ERSHU)
                     PP(LQP(MP),5)=SQRT(AMSHU)
                  ENDIF
 C               ********reshuffle the energy when jet has mass
 C
                  R0=RDP(MP)
                  DP1=DP*V1
                  DP2=DP*V2
                  DP3=DP*V3
 C               ********momentum and energy transfer from jet
                  NPJ(LQP(MP))=NPJ(LQP(MP))+1
                  KFPJ(LQP(MP),NPJ(LQP(MP)))=21
                  PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
                  PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
                  PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
                  PJPE(LQP(MP),NPJ(LQP(MP)))=DP
                  PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
 
                  GO TO 2660
 
  2640            MT=MT+1
                  DRR=RDT(MT)-R0
                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
                  DP=DRR*HIPR1(14)
                  IF(DP.LE.0.2) GO TO 2610
                  IF(PTOT.LE.0.4) GO TO 2690
                  IF(PTOT.LE.DP) DP=PTOT-0.2
                  DE=DP
 
                  IF(K2SG(ISG,I).NE.21) THEN
                     IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
                     PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
      &                    +PT(LQT(MT),3)**2
                     DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
      &                 -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
                     ERSHU=(PT(LQT(MT),4)+DE-DP)**2
                     AMSHU=ERSHU-PRSHU
                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
                     PT(LQT(MT),4)=SQRT(ERSHU)
                     PT(LQT(MT),5)=SQRT(AMSHU)
                  ENDIF
 C               ********reshuffle the energy when jet has mass
 
                  R0=RDT(MT)
                  DP1=DP*V1
                  DP2=DP*V2
                  DP3=DP*V3
 C               ********momentum and energy transfer from jet
                  NTJ(LQT(MT))=NTJ(LQT(MT))+1
                  KFTJ(LQT(MT),NTJ(LQT(MT)))=21
                  PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
                  PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
                  PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
                  PJTE(LQT(MT),NTJ(LQT(MT)))=DP
                  PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
 
  2660            PXSG(ISG,I)=(PTOT-DP)*V1
                  PYSG(ISG,I)=(PTOT-DP)*V2
                  PZSG(ISG,I)=(PTOT-DP)*V3
                  PESG(ISG,I)=PESG(ISG,I)-DE
 
                  PTOT=PTOT-DP
                  NQ=NQ+1
                  GO TO 2600
  2690   CONTINUE
         RETURN
         END
 
 C
 C
 C
 C
         SUBROUTINE HIJFRG(JTP,NTP,IERROR)
 C       NTP=1, fragment proj string, NTP=2, targ string, 
 C       NTP=3, independent 
 C       strings from jets.  JTP is the line number of the string
 C*******Fragment all leading strings of proj and targ**************
 C       IHNT2(1)=atomic #, IHNT2(2)=proton #(=-1 if anti-proton)  *
 C******************************************************************
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
         SAVE  /HIJDAT/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &                PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
 C
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
         SAVE  /LUJETS/
         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
         SAVE  /LUDAT1/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         
         IERROR=0
         CALL LUEDIT(0)
         N=0
 C                       ********initialize the document lines
         IF(NTP.EQ.3) THEN
                 ISG=JTP
                 N=NJSG(ISG)
                 DO 100 I=1,NJSG(ISG)
                         K(I,1)=K1SG(ISG,I)
                         K(I,2)=K2SG(ISG,I)
                         P(I,1)=PXSG(ISG,I)
                         P(I,2)=PYSG(ISG,I)
                         P(I,3)=PZSG(ISG,I)
                         P(I,4)=PESG(ISG,I)
                         P(I,5)=PMSG(ISG,I)
 100             CONTINUE
 C               IF(IHPR2(1).GT.0) CALL ATTRAD(IERROR)
 c               IF(IERROR.NE.0) RETURN
 C               CALL LULIST(1)
                 CALL LUEXEC
                 RETURN
         ENDIF
 C
         IF(NTP.EQ.2) GO TO 200
         IF(JTP.GT.IHNT2(1))   RETURN
         IF(NFP(JTP,5).NE.3.AND.NFP(JTP,3).NE.0
      &      .AND.NPJ(JTP).EQ.0.AND.NFP(JTP,10).EQ.0) GO TO 1000
         IF(NFP(JTP,15).EQ.-1) THEN
                 KF1=NFP(JTP,2)
                 KF2=NFP(JTP,1)
                 PQ21=PP(JTP,6)
                 PQ22=PP(JTP,7)
                 PQ11=PP(JTP,8)
                 PQ12=PP(JTP,9)
                 AM1=PP(JTP,15)
                 AM2=PP(JTP,14)
         ELSE
                 KF1=NFP(JTP,1)
                 KF2=NFP(JTP,2)
                 PQ21=PP(JTP,8)
                 PQ22=PP(JTP,9)
                 PQ11=PP(JTP,6)
                 PQ12=PP(JTP,7)
                 AM1=PP(JTP,14)
                 AM2=PP(JTP,15)  
         ENDIF
 C       ********for NFP(JTP,15)=-1 NFP(JTP,1) IS IN -Z DIRECTION
         PB1=PQ11+PQ21
         PB2=PQ12+PQ22
         PB3=PP(JTP,3)
         PECM=PP(JTP,5)
         BTZ=PB3/PP(JTP,4)
         IF((ABS(PB1-PP(JTP,1)).GT.0.01.OR.
      &    ABS(PB2-PP(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
      &    WRITE(6,*) '  Pt of Q and QQ do not sum to the total'
 
         GO TO 300
 
 200     IF(JTP.GT.IHNT2(3))  RETURN
         IF(NFT(JTP,5).NE.3.AND.NFT(JTP,3).NE.0
      &     .AND.NTJ(JTP).EQ.0.AND.NFT(JTP,10).EQ.0) GO TO 1200
         IF(NFT(JTP,15).EQ.1) THEN
                 KF1=NFT(JTP,1)
                 KF2=NFT(JTP,2)
                 PQ11=PT(JTP,6)
                 PQ12=PT(JTP,7)
                 PQ21=PT(JTP,8)
                 PQ22=PT(JTP,9)
                 AM1=PT(JTP,14)
                 AM2=PT(JTP,15)
         ELSE
                 KF1=NFT(JTP,2)
                 KF2=NFT(JTP,1)
                 PQ11=PT(JTP,8)
                 PQ12=PT(JTP,9)
                 PQ21=PT(JTP,6)
                 PQ22=PT(JTP,7)
                 AM1=PT(JTP,15)
                 AM2=PT(JTP,14)
         ENDIF   
 C       ********for NFT(JTP,15)=1 NFT(JTP,1) IS IN +Z DIRECTION
         PB1=PQ11+PQ21
         PB2=PQ12+PQ22
         PB3=PT(JTP,3)
         PECM=PT(JTP,5)
         BTZ=PB3/PT(JTP,4)
 
         IF((ABS(PB1-PT(JTP,1)).GT.0.01.OR.
      &     ABS(PB2-PT(JTP,2)).GT.0.01).AND.IHPR2(10).NE.0)
      &     WRITE(6,*) '  Pt of Q and QQ do not sum to the total'
 
 300     IF(PECM.LT.HIPR1(1)) THEN
            IERROR=1
            IF(IHPR2(10).EQ.0) RETURN
            WRITE(6,*) ' ECM=',PECM,' energy of the string is too small'
            RETURN
         ENDIF
         AMT=PECM**2+PB1**2+PB2**2
         AMT1=AM1**2+PQ11**2+PQ12**2
         AMT2=AM2**2+PQ21**2+PQ22**2
         PZCM=SQRT(ABS(AMT**2+AMT1**2+AMT2**2-2.0*AMT*AMT1
      &       -2.0*AMT*AMT2-2.0*AMT1*AMT2))/2.0/SQRT(AMT)
 C               *******PZ of end-partons in c.m. frame of the string
         K(1,1)=2
         K(1,2)=KF1
         P(1,1)=PQ11
         P(1,2)=PQ12
         P(1,3)=PZCM
         P(1,4)=SQRT(AMT1+PZCM**2)
         P(1,5)=AM1
         K(2,1)=1
         K(2,2)=KF2
         P(2,1)=PQ21
         P(2,2)=PQ22
         P(2,3)=-PZCM
         P(2,4)=SQRT(AMT2+PZCM**2)
         P(2,5)=AM2
         N=2
 C*****
         CALL HIROBO(0.0,0.0,0.0,0.0,BTZ)
         JETOT=0
         IF((PQ21**2+PQ22**2).GT.(PQ11**2+PQ12**2)) THEN
                 PMAX1=P(2,1)
                 PMAX2=P(2,2)
                 PMAX3=P(2,3)
         ELSE
                 PMAX1=P(1,1)
                 PMAX2=P(1,2)
                 PMAX3=P(1,3)
         ENDIF
         IF(NTP.EQ.1) THEN
                 PP(JTP,10)=PMAX1
                 PP(JTP,11)=PMAX2
                 PP(JTP,12)=PMAX3
         ELSE IF(NTP.EQ.2) THEN
                 PT(JTP,10)=PMAX1
                 PT(JTP,11)=PMAX2
                 PT(JTP,12)=PMAX3
         ENDIF
 C*******************attach produced jets to the leading partons****
         IF(NTP.EQ.1.AND.NPJ(JTP).NE.0) THEN
                 JETOT=NPJ(JTP)
 C               IF(NPJ(JTP).GE.2) CALL HIJSRT(JTP,1)
 C                       ********sort jets in order of y
                 IEX=0
                 IF((ABS(KF1).GT.1000.AND.KF1.LT.0)
      &                  .OR.(ABS(KF1).LT.1000.AND.KF1.GT.0)) IEX=1
                 DO 520 I=N,2,-1
                 DO 520 J=1,5
                         II=NPJ(JTP)+I
                         K(II,J)=K(I,J)
                         P(II,J)=P(I,J)
                         V(II,J)=V(I,J)
 520             CONTINUE
                 DO 540 I=1,NPJ(JTP)
                         DO 542 J=1,5
                                 K(I+1,J)=0
                                 V(I+1,J)=0
 542                     CONTINUE                                
                         I0=I
                         IF(IEX.EQ.1) I0=NPJ(JTP)-I+1
 C                               ********reverse the order of jets
                         KK1=KFPJ(JTP,I0)
                         K(I+1,1)=2
                         K(I+1,2)=KK1
                         IF(KK1.NE.21 .AND. KK1.NE.0)  K(I+1,1)=
      &                    1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
                         P(I+1,1)=PJPX(JTP,I0)
                         P(I+1,2)=PJPY(JTP,I0)
                         P(I+1,3)=PJPZ(JTP,I0)
                         P(I+1,4)=PJPE(JTP,I0)
                         P(I+1,5)=PJPM(JTP,I0)
 540             CONTINUE
                 N=N+NPJ(JTP)
         ELSE IF(NTP.EQ.2.AND.NTJ(JTP).NE.0) THEN
                 JETOT=NTJ(JTP)
 c               IF(NTJ(JTP).GE.2)  CALL HIJSRT(JTP,2)
 C                       ********sort jets in order of y
                 IEX=1
                 IF((ABS(KF2).GT.1000.AND.KF2.LT.0)
      &                  .OR.(ABS(KF2).LT.1000.AND.KF2.GT.0)) IEX=0
                 DO 560 I=N,2,-1
                 DO 560 J=1,5
                         II=NTJ(JTP)+I
                         K(II,J)=K(I,J)
                         P(II,J)=P(I,J)
                         V(II,J)=V(I,J)
 560             CONTINUE
                 DO 580 I=1,NTJ(JTP)
                         DO 582 J=1,5
                                 K(I+1,J)=0
                                 V(I+1,J)=0
 582                     CONTINUE                                
                         I0=I
                         IF(IEX.EQ.1) I0=NTJ(JTP)-I+1
 C                               ********reverse the order of jets
                         KK1=KFTJ(JTP,I0)
                         K(I+1,1)=2
                         K(I+1,2)=KK1
                         IF(KK1.NE.21 .AND. KK1.NE.0) K(I+1,1)=
      &                     1+(ABS(KK1)+(2*IEX-1)*KK1)/2/ABS(KK1)
                         P(I+1,1)=PJTX(JTP,I0)
                         P(I+1,2)=PJTY(JTP,I0)
                         P(I+1,3)=PJTZ(JTP,I0)
                         P(I+1,4)=PJTE(JTP,I0)
                         P(I+1,5)=PJTM(JTP,I0)
 580             CONTINUE
                 N=N+NTJ(JTP)
         ENDIF
         IF(IHPR2(1).GT.0.AND.RLU(0).LE.HIDAT(3)) THEN
              HIDAT20=HIDAT(2)
              HIPR150=HIPR1(5)
              IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
      &                  HIDAT(2)=2.0
              IF(HINT1(1).GE.1000.0.AND.JETOT.EQ.0)THEN
                 HIDAT(2)=3.0
                 HIPR1(5)=5.0
              ENDIF
              CALL ATTRAD(IERROR)
              HIDAT(2)=HIDAT20
              HIPR1(5)=HIPR150
         ELSE IF(JETOT.EQ.0.AND.IHPR2(1).GT.0.AND.
      &                       HINT1(1).GE.1000.0.AND.
      &          RLU(0).LE.0.8) THEN
                 HIDAT20=HIDAT(2)
                 HIPR150=HIPR1(5)
                 HIDAT(2)=3.0
                 HIPR1(5)=5.0
              IF(IHPR2(8).EQ.0.AND.IHPR2(3).EQ.0.AND.IHPR2(9).EQ.0)
      &                  HIDAT(2)=2.0
                 CALL ATTRAD(IERROR)
                 HIDAT(2)=HIDAT20
                 HIPR1(5)=HIPR150
         ENDIF
         IF(IERROR.NE.0) RETURN
 C               ******** conduct soft radiations
 C****************************
 C
 C
 C       CALL LULIST(1)
         CALL LUEXEC
         RETURN
 
 1000    N=1
         K(1,1)=1
         K(1,2)=NFP(JTP,3)
         DO 1100 JJ=1,5
                 P(1,JJ)=PP(JTP,JJ)
 1100            CONTINUE
 C                       ********proj remain as a nucleon or delta
         CALL LUEXEC
 C       call lulist(1)
         RETURN
 C
 1200    N=1
         K(1,1)=1
         K(1,2)=NFT(JTP,3)
         DO 1300 JJ=1,5
                 P(1,JJ)=PT(JTP,JJ)
 1300    CONTINUE
 C                       ********targ remain as a nucleon or delta
         CALL LUEXEC
 C       call lulist(1)
         RETURN
         END
 C
 C
 C
 C********************************************************************
 C       Sort the jets associated with a nucleon in order of their
 C       rapdities
 C********************************************************************
         SUBROUTINE HIJSRT(JPJT,NPT)
         DIMENSION KF(100),PX(100),PY(100),PZ(100),PE(100),PM(100)
         DIMENSION Y(100),IP(100,2)
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &                PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         IF(NPT.EQ.2) GO TO 500
         JP=JPJT
         IQ=0
         I=1
 100     KF(I)=KFPJ(JP,I)
         PX(I)=PJPX(JP,I)
         PY(I)=PJPY(JP,I)
         PZ(I)=PJPZ(JP,I)
         PE(I)=PJPE(JP,I)
         PM(I)=PJPM(JP,I)
         Y(I-IQ)=0.5*ALOG((ABS(PE(I)+PZ(I))+1.E-5)
      &          /(ABS(PE(I)-PZ(I))+1.E-5))
         IP(I-IQ,1)=I
         IP(I-IQ,2)=0
         IF(KF(I).NE.21) THEN
                 IP(I-IQ,2)=1
                 IQ=IQ+1
                 I=I+1
                 KF(I)=KFPJ(JP,I)
                 PX(I)=PJPX(JP,I)
                 PY(I)=PJPY(JP,I)
                 PZ(I)=PJPZ(JP,I)
                 PE(I)=PJPE(JP,I)
                 PM(I)=PJPM(JP,I)
         ENDIF
         I=I+1
         IF(I.LE.NPJ(JP)) GO TO 100
                         
         DO 200 I=1,NPJ(JP)-IQ
         DO 200 J=I+1,NPJ(JP)-IQ
                 IF(Y(I).GT.Y(J)) GO TO 200
                 IP1=IP(I,1)
                 IP2=IP(I,2)
                 IP(I,1)=IP(J,1)
                 IP(I,2)=IP(J,2)
                 IP(J,1)=IP1
                 IP(J,2)=IP2
 200     CONTINUE
 C                       ********sort in decending y
         IQQ=0
         I=1
 300     KFPJ(JP,I)=KF(IP(I-IQQ,1))
         PJPX(JP,I)=PX(IP(I-IQQ,1))
         PJPY(JP,I)=PY(IP(I-IQQ,1))
         PJPZ(JP,I)=PZ(IP(I-IQQ,1))
         PJPE(JP,I)=PE(IP(I-IQQ,1))
         PJPM(JP,I)=PM(IP(I-IQQ,1))
         IF(IP(I-IQQ,2).EQ.1) THEN
                 KFPJ(JP,I+1)=KF(IP(I-IQQ,1)+1)
                 PJPX(JP,I+1)=PX(IP(I-IQQ,1)+1)
                 PJPY(JP,I+1)=PY(IP(I-IQQ,1)+1)
                 PJPZ(JP,I+1)=PZ(IP(I-IQQ,1)+1)
                 PJPE(JP,I+1)=PE(IP(I-IQQ,1)+1)
                 PJPM(JP,I+1)=PM(IP(I-IQQ,1)+1)
                 I=I+1
                 IQQ=IQQ+1
         ENDIF
         I=I+1
         IF(I.LE.NPJ(JP)) GO TO 300
 
         RETURN
 
 500     JT=JPJT
         IQ=0
         I=1
 600     KF(I)=KFTJ(JT,I)
         PX(I)=PJTX(JT,I)
         PY(I)=PJTY(JT,I)
         PZ(I)=PJTZ(JT,I)
         PE(I)=PJTE(JT,I)
         PM(I)=PJTM(JT,I)
         Y(I-IQ)=0.5*ALOG((ABS(PE(I)+PZ(I))+1.E-5)
      &          /(ABS(PE(I)-PZ(I))+1.E-5))
         IP(I-IQ,1)=I
         IP(I-IQ,2)=0
         IF(KF(I).NE.21) THEN
                 IP(I-IQ,2)=1
                 IQ=IQ+1
                 I=I+1
                 KF(I)=KFTJ(JT,I)
                 PX(I)=PJTX(JT,I)
                 PY(I)=PJTY(JT,I)
                 PZ(I)=PJTZ(JT,I)
                 PE(I)=PJTE(JT,I)
                 PM(I)=PJTM(JT,I)
         ENDIF
         I=I+1
         IF(I.LE.NTJ(JT)) GO TO 600
                         
         DO 700 I=1,NTJ(JT)-IQ
         DO 700 J=I+1,NTJ(JT)-IQ
                 IF(Y(I).LT.Y(J)) GO TO 700
                 IP1=IP(I,1)
                 IP2=IP(I,2)
                 IP(I,1)=IP(J,1)
                 IP(I,2)=IP(J,2)
                 IP(J,1)=IP1
                 IP(J,2)=IP2
 700     CONTINUE
 C                       ********sort in acending y
         IQQ=0
         I=1
 800     KFTJ(JT,I)=KF(IP(I-IQQ,1))
         PJTX(JT,I)=PX(IP(I-IQQ,1))
         PJTY(JT,I)=PY(IP(I-IQQ,1))
         PJTZ(JT,I)=PZ(IP(I-IQQ,1))
         PJTE(JT,I)=PE(IP(I-IQQ,1))
         PJTM(JT,I)=PM(IP(I-IQQ,1))
         IF(IP(I-IQQ,2).EQ.1) THEN
                 KFTJ(JT,I+1)=KF(IP(I-IQQ,1)+1)
                 PJTX(JT,I+1)=PX(IP(I-IQQ,1)+1)
                 PJTY(JT,I+1)=PY(IP(I-IQQ,1)+1)
                 PJTZ(JT,I+1)=PZ(IP(I-IQQ,1)+1)
                 PJTE(JT,I+1)=PE(IP(I-IQQ,1)+1)
                 PJTM(JT,I+1)=PM(IP(I-IQQ,1)+1)
                 I=I+1
                 IQQ=IQQ+1
         ENDIF
         I=I+1
         IF(I.LE.NTJ(JT)) GO TO 800
         RETURN
         END     
 
 C
 C
 C
 C****************************************************************
 C       conduct soft radiation according to dipole approxiamtion
 C****************************************************************
         SUBROUTINE ATTRAD(IERROR)
 C
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
         SAVE  /HIJDAT/
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
         SAVE  /LUJETS/
         IERROR=0
 
 C.....S INVARIANT MASS-SQUARED BETWEEN PARTONS I AND I+1......
 C.....SM IS THE LARGEST MASS-SQUARED....
 
 40      SM=0.
         JL=1
         DO 30 I=1,N-1
            S=2.*(P(I,4)*P(I+1,4)-P(I,1)*P(I+1,1)-P(I,2)*P(I+1,2)
      &          -P(I,3)*P(I+1,3))+P(I,5)**2+P(I+1,5)**2
            IF(S.LT.0.) S=0.
            WP=SQRT(S)-1.5*(P(I,5)+P(I+1,5))
            IF(WP.GT.SM) THEN
               PBT1=P(I,1)+P(I+1,1)
               PBT2=P(I,2)+P(I+1,2)
               PBT3=P(I,3)+P(I+1,3)
               PBT4=P(I,4)+P(I+1,4)
               BTT=(PBT1**2+PBT2**2+PBT3**2)/PBT4**2
               IF(BTT.GE.1.0-1.0E-10) GO TO 30
               IF((I.NE.1.OR.I.NE.N-1).AND.
      &             (K(I,2).NE.21.AND.K(I+1,2).NE.21)) GO TO 30
               JL=I
               SM=WP
            ENDIF
 30      CONTINUE
         S=(SM+1.5*(P(JL,5)+P(JL+1,5)))**2
         IF(SM.LT.HIPR1(5)) GOTO 2
      
 C.....MAKE PLACE FOR ONE GLUON.....
         IF(JL+1.EQ.N) GOTO 190
         DO 160 J=N,JL+2,-1
                 K(J+1,1)=K(J,1)
                 K(J+1,2)=K(J,2)
                 DO 150 M=1,5
 150                     P(J+1,M)=P(J,M)
 160             CONTINUE
 190     N=N+1
      
 C.....BOOST TO REST SYSTEM FOR PARTICLES JL AND JL+1.....
         P1=P(JL,1)+P(JL+1,1)
         P2=P(JL,2)+P(JL+1,2)
         P3=P(JL,3)+P(JL+1,3)
         P4=P(JL,4)+P(JL+1,4)
         BEX=-P1/P4
         BEY=-P2/P4
         BEZ=-P3/P4
         IMIN=JL
         IMAX=JL+1
         CALL ATROBO(0.,0.,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
         IF(IERROR.NE.0) RETURN
 C.....ROTATE TO Z-AXIS....
         CTH=P(JL,3)/SQRT(P(JL,4)**2-P(JL,5)**2)
         IF(ABS(CTH).GT.1.0)  CTH=MAX(-1.,MIN(1.,CTH))
         THETA=ACOS(CTH)
         PHI=ULANGL(P(JL,1),P(JL,2))
         CALL ATROBO(0.,-PHI,0.,0.,0.,IMIN,IMAX,IERROR)
         CALL ATROBO(-THETA,0.,0.,0.,0.,IMIN,IMAX,IERROR)
      
 C.....CREATE ONE GLUON AND ORIENTATE.....
      
 1       CALL AR3JET(S,X1,X3,JL)
         CALL ARORIE(S,X1,X3,JL)         
         IF(HIDAT(2).GT.0.0) THEN
            PTG1=SQRT(P(JL,1)**2+P(JL,2)**2)
            PTG2=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
            PTG3=SQRT(P(JL+2,1)**2+P(JL+2,2)**2)
            PTG=MAX(PTG1,PTG2,PTG3)
            IF(PTG.GT.HIDAT(2)) THEN
               FMFACT=EXP(-(PTG**2-HIDAT(2)**2)/HIPR1(2)**2)
               IF(RLU(0).GT.FMFACT) GO TO 1
            ENDIF
         ENDIF
 C.....ROTATE AND BOOST BACK.....
         IMIN=JL
         IMAX=JL+2
         CALL ATROBO(THETA,PHI,-BEX,-BEY,-BEZ,IMIN,IMAX,IERROR)
         IF(IERROR.NE.0) RETURN
 C.....ENUMERATE THE GLUONS.....
         K(JL+2,1)=K(JL+1,1)
         K(JL+2,2)=K(JL+1,2)
         K(JL+2,3)=K(JL+1,3)
         K(JL+2,4)=K(JL+1,4)
         K(JL+2,5)=K(JL+1,5)
         P(JL+2,5)=P(JL+1,5)
         K(JL+1,1)=2
         K(JL+1,2)=21
         K(JL+1,3)=0
         K(JL+1,4)=0
         K(JL+1,5)=0
         P(JL+1,5)=0.
 C----THETA FUNCTION DAMPING OF THE EMITTED GLUONS. FOR HADRON-HADRON.
 C----R0=VFR(2)
 C       IF(VFR(2).GT.0.) THEN
 C       PTG=SQRT(P(JL+1,1)**2+P(JL+1,2)**2)
 C       PTGMAX=WSTRI/2.
 C       DOPT=SQRT((4.*PAR(71)*VFR(2))/WSTRI)
 C       PTOPT=(DOPT*WSTRI)/(2.*VFR(2))
 C       IF(PTG.GT.PTOPT) IORDER=IORDER-1
 C       IF(PTG.GT.PTOPT) GOTO 1
 C       ENDIF
 C-----
         IF(SM.GE.HIPR1(5)) GOTO 40
 
 2       K(1,1)=2
         K(1,3)=0
         K(1,4)=0
         K(1,5)=0
         K(N,1)=1
         K(N,3)=0
         K(N,4)=0
         K(N,5)=0
 
         RETURN
         END
 
 
         SUBROUTINE AR3JET(S,X1,X3,JL)
 C     
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
         SAVE  /LUJETS/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C     
         C=1./3.
         IF(K(JL,2).NE.21 .AND. K(JL+1,2).NE.21) C=8./27.
         EXP1=3
         EXP3=3
         IF(K(JL,2).NE.21) EXP1=2
         IF(K(JL+1,2).NE.21) EXP3=2
         A=0.24**2/S
         YMA=ALOG(.5/SQRT(A)+SQRT(.25/A-1))
         D=4.*C*YMA
         SM1=P(JL,5)**2/S
         SM3=P(JL+1,5)**2/S
         XT2M=(1.-2.*SQRT(SM1)+SM1-SM3)*(1.-2.*SQRT(SM3)-SM1+SM3)
         XT2M=MIN(.25,XT2M)
         NTRY=0
 1       IF(NTRY.EQ.5000) THEN
                 X1=.5*(2.*SQRT(SM1)+1.+SM1-SM3)
                 X3=.5*(2.*SQRT(SM3)+1.-SM1+SM3)
                 RETURN
         ENDIF
         NTRY=NTRY+1
      
         XT2=A*(XT2M/A)**(RLU(0)**(1./D))
      
         YMAX=ALOG(.5/SQRT(XT2)+SQRT(.25/XT2-1.))
         Y=(2.*RLU(0)-1.)*YMAX
         X1=1.-SQRT(XT2)*EXP(Y)
         X3=1.-SQRT(XT2)*EXP(-Y)
         X2=2.-X1-X3
         NEG=0
         IF(K(JL,2).NE.21 .OR. K(JL+1,2).NE.21) THEN
         IF((1.-X1)*(1.-X2)*(1.-X3)-X2*SM1*(1.-X1)-X2*SM3*(1.-X3).
      &  LE.0..OR.X1.LE.2.*SQRT(SM1)-SM1+SM3.OR.X3.LE.2.*SQRT(SM3)
      &  -SM3+SM1) NEG=1
         X1=X1+SM1-SM3
         X3=X3-SM1+SM3
         ENDIF
         IF(NEG.EQ.1) GOTO 1
      
         FG=2.*YMAX*C*(X1**EXP1+X3**EXP3)/D
         XT2M=XT2
         IF(FG.LT.RLU(0)) GOTO 1
      
         RETURN
         END
 C*************************************************************
 
 
         SUBROUTINE ARORIE(S,X1,X3,JL)
 C     
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
         SAVE  /LUJETS/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C     
         W=SQRT(S)
         X2=2.-X1-X3
         E1=.5*X1*W
         E3=.5*X3*W
         P1=SQRT(E1**2-P(JL,5)**2)
         P3=SQRT(E3**2-P(JL+1,5)**2)
         CBET=1.
         IF(P1.GT.0..AND.P3.GT.0.) CBET=(P(JL,5)**2
      &           +P(JL+1,5)**2+2.*E1*E3-S*(1.-X2))/(2.*P1*P3)
         IF(ABS(CBET).GT.1.0) CBET=MAX(-1.,MIN(1.,CBET))
         BET=ACOS(CBET)
      
 C.....MINIMIZE PT1-SQUARED PLUS PT3-SQUARED.....
         IF(P1.GE.P3) THEN
            PSI=.5*ULANGL(P1**2+P3**2*COS(2.*BET),-P3**2*SIN(2.*BET))
            PT1=P1*SIN(PSI)
            PZ1=P1*COS(PSI)
            PT3=P3*SIN(PSI+BET)
            PZ3=P3*COS(PSI+BET)
         ELSE IF(P3.GT.P1) THEN
            PSI=.5*ULANGL(P3**2+P1**2*COS(2.*BET),-P1**2*SIN(2.*BET))
            PT1=P1*SIN(BET+PSI)
            PZ1=-P1*COS(BET+PSI)
            PT3=P3*SIN(PSI)
            PZ3=-P3*COS(PSI)
         ENDIF
      
         DEL=2.0*HIPR1(40)*RLU(0)
         P(JL,4)=E1
         P(JL,1)=PT1*SIN(DEL)
         P(JL,2)=-PT1*COS(DEL)
         P(JL,3)=PZ1
         P(JL+2,4)=E3
         P(JL+2,1)=PT3*SIN(DEL)
         P(JL+2,2)=-PT3*COS(DEL)
         P(JL+2,3)=PZ3
         P(JL+1,4)=W-E1-E3
         P(JL+1,1)=-P(JL,1)-P(JL+2,1)
         P(JL+1,2)=-P(JL,2)-P(JL+2,2)
         P(JL+1,3)=-P(JL,3)-P(JL+2,3)
         RETURN
         END
 
 
 C
 C*******************************************************************
 C       make  boost and rotation to entries from IMIN to IMAX
 C*******************************************************************
         SUBROUTINE ATROBO(THE,PHI,BEX,BEY,BEZ,IMIN,IMAX,IERROR)
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
         SAVE  /LUJETS/
         DIMENSION ROT(3,3),PV(3)
         DOUBLE PRECISION DP(4),DBEX,DBEY,DBEZ,DGA,DGA2,DBEP,DGABEP
         IERROR=0
      
         IF(IMIN.LE.0 .OR. IMAX.GT.N .OR. IMIN.GT.IMAX) RETURN
 
         IF(THE**2+PHI**2.GT.1E-20) THEN
 C...ROTATE (TYPICALLY FROM Z AXIS TO DIRECTION THETA,PHI)
            ROT(1,1)=COS(THE)*COS(PHI)
            ROT(1,2)=-SIN(PHI)
            ROT(1,3)=SIN(THE)*COS(PHI)
            ROT(2,1)=COS(THE)*SIN(PHI)
            ROT(2,2)=COS(PHI)
            ROT(2,3)=SIN(THE)*SIN(PHI)
            ROT(3,1)=-SIN(THE)
            ROT(3,2)=0.
            ROT(3,3)=COS(THE)
            DO 120 I=IMIN,IMAX
 C**************    IF(MOD(K(I,1)/10000,10).GE.6) GOTO 120
               DO 100 J=1,3
  100             PV(J)=P(I,J)
                  DO 110 J=1,3
  110                P(I,J)=ROT(J,1)*PV(1)+ROT(J,2)*PV(2)
      &                     +ROT(J,3)*PV(3)
  120             CONTINUE
         ENDIF
      
         IF(BEX**2+BEY**2+BEZ**2.GT.1E-20) THEN
 C...LORENTZ BOOST (TYPICALLY FROM REST TO MOMENTUM/ENERGY=BETA)
                 DBEX=BEX
                 DBEY=BEY
                 DBEZ=BEZ
                 DGA2=1D0-DBEX**2-DBEY**2-DBEZ**2
                 IF(DGA2.LE.0D0) THEN
                         IERROR=1
                         RETURN
                 ENDIF
                 DGA=1D0/DSQRT(DGA2)
                 DO 140 I=IMIN,IMAX
 C*************     IF(MOD(K(I,1)/10000,10).GE.6) GOTO 140
                    DO 130 J=1,4
 130                DP(J)=P(I,J)
                    DBEP=DBEX*DP(1)+DBEY*DP(2)+DBEZ*DP(3)
                    DGABEP=DGA*(DGA*DBEP/(1D0+DGA)+DP(4))
                    P(I,1)=DP(1)+DGABEP*DBEX
                    P(I,2)=DP(2)+DGABEP*DBEY
                    P(I,3)=DP(3)+DGABEP*DBEZ
                    P(I,4)=DGA*(DP(4)+DBEP)
 140             CONTINUE
         ENDIF
      
         RETURN
         END
 C
 C
 C
         SUBROUTINE HIJHRD(JP,JT,JOUT,JFLG,IOPJET0)
 C
 C       IOPTJET=1, ALL JET WILL FORM SINGLE STRING SYSTEM
 C               0, ONLY Q-QBAR JET FORM SINGLE STRING SYSTEM
 C*******Perform jets production and fragmentation when JP JT *******
 C     scatter. JOUT-> number of hard scatterings precede this one  *
 C     for the the same pair(JP,JT). JFLG->a flag to show whether   *
 C     jets can be produced (with valence quark=1,gluon=2, q-qbar=3)*
 C     or not(0). Information of jets are in  COMMON/ATTJET and     *
 C     /MINJET. ABS(NFP(JP,6)) is the total number jets produced by *
 C    JP. If NFP(JP,6)<0 JP can not produce jet anymore.            *
 C*******************************************************************
         DIMENSION IP(100,2),IPQ(50),IPB(50),IT(100,2),ITQ(50),ITB(50)
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
         SAVE  /HIJDAT/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &                PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
         COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
         SAVE  /HIJJET4/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C************************************ HIJING common block
         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
         SAVE  /LUJETS/
         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
         SAVE  /LUDAT1/
         COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
         SAVE  /PYSUBS/
         COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
         SAVE  /PYPARS/
         COMMON/PYINT1/MINT(400),VINT(400)
         SAVE  /PYINT1/
         COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
         SAVE  /PYINT2/
         COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
         SAVE  /PYINT5/
         COMMON/HIPYINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
         SAVE  /HIPYINT/
 C*********************************** LU common block
         MXJT=500
 C               SIZE OF COMMON BLOCK FOR # OF PARTON PER STRING
         MXSG=900
 C               SIZE OF COMMON BLOCK FOR # OF SINGLE STRINGS
         MXSJ=100
 C               SIZE OF COMMON BLOCK FOR # OF PARTON PER SINGLE
 C               STRING
         JFLG=0
         IHNT2(11)=JP
         IHNT2(12)=JT
 C
         IOPJET=IOPJET0
         IF(IOPJET.EQ.1.AND.(NFP(JP,6).NE.0.OR.NFT(JT,6).NE.0))
      &                   IOPJET=0
         IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
         IF(NFP(JP,6).LT.0 .OR. NFT(JT,6).LT.0) RETURN
 C               ******** JP or JT can not produce jet anymore
 C
         IF(JOUT.EQ.0) THEN
                 EPP=PP(JP,4)+PP(JP,3)
                 EPM=PP(JP,4)-PP(JP,3)
                 ETP=PT(JT,4)+PT(JT,3)
                 ETM=PT(JT,4)-PT(JT,3)
                 IF(EPP.LT.0.0) GO TO 1000
                 IF(EPM.LT.0.0) GO TO 1000
                 IF(ETP.LT.0.0) GO TO 1000
                 IF(ETM.LT.0.0) GO TO 1000
                 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
         ENDIF
 C               ********for the first hard scattering of (JP,JT)
 C                       have collision only when Ycm(JP)>Ycm(JT)
 
         ECUT1=HIPR1(1)+HIPR1(8)+PP(JP,14)+PP(JP,15)
         ECUT2=HIPR1(1)+HIPR1(8)+PT(JT,14)+PT(JT,15)
         IF(PP(JP,4).LE.ECUT1) THEN
                 NFP(JP,6)=-ABS(NFP(JP,6))
                 RETURN
         ENDIF
         IF(PT(JT,4).LE.ECUT2) THEN
                 NFT(JT,6)=-ABS(NFT(JT,6))
                 RETURN
         ENDIF
 C               *********must have enough energy to produce jets
 
         MISS=0
         MISP=0
         MIST=0
 C
         IF(NFP(JP,10).EQ.0 .AND. NFT(JT,10).EQ.0) THEN
                 MINT(44)=MINT4
                 MINT(45)=MINT5
                 XSEC(0,1)=ATXS(0)
                 XSEC(11,1)=ATXS(11)
                 XSEC(12,1)=ATXS(12)
                 XSEC(28,1)=ATXS(28)
                 DO 120 I=1,20
                 COEF(11,I)=ATCO(11,I)
                 COEF(12,I)=ATCO(12,I)
                 COEF(28,I)=ATCO(28,I)
 120             CONTINUE
         ELSE
                 ISUB11=0
                 ISUB12=0
                 ISUB28=0
                 IF(XSEC(11,1).NE.0) ISUB11=1
                 IF(XSEC(12,1).NE.0) ISUB12=1
                 IF(XSEC(28,1).NE.0) ISUB28=1            
                 MINT(44)=MINT4-ISUB11-ISUB12-ISUB28
                 MINT(45)=MINT5-ISUB11-ISUB12-ISUB28
                 XSEC(0,1)=ATXS(0)-ATXS(11)-ATXS(12)-ATXS(28)
                 XSEC(11,1)=0.0
                 XSEC(12,1)=0.0
                 XSEC(28,1)=0.0  
                 DO 110 I=1,20
                 COEF(11,I)=0.0
                 COEF(12,I)=0.0
                 COEF(28,I)=0.0
 110             CONTINUE
         ENDIF           
 C       ********Scatter the valence quarks only once per NN 
 C       collision,
 C               afterwards only gluon can have hard scattering.
  155    CALL PYTHIA
         JJ=MINT(31)
         IF(JJ.NE.1) GO TO 155
 C               *********one hard collision at a time
         IF(K(7,2).EQ.-K(8,2)) THEN
                 QMASS2=(P(7,4)+P(8,4))**2-(P(7,1)+P(8,1))**2
      &                  -(P(7,2)+P(8,2))**2-(P(7,3)+P(8,3))**2
                 QM=ULMASS(K(7,2))
                 IF(QMASS2.LT.(2.0*QM+HIPR1(1))**2) GO TO 155
         ENDIF
 C               ********q-qbar jets must has minimum mass HIPR1(1)
         PXP=PP(JP,1)-P(3,1)
         PYP=PP(JP,2)-P(3,2)
         PZP=PP(JP,3)-P(3,3)
         PEP=PP(JP,4)-P(3,4)
         PXT=PT(JT,1)-P(4,1)
         PYT=PT(JT,2)-P(4,2)
         PZT=PT(JT,3)-P(4,3)
         PET=PT(JT,4)-P(4,4)
 
         IF(PEP.LE.ECUT1) THEN
                 MISP=MISP+1
                 IF(MISP.LT.50) GO TO 155
                 NFP(JP,6)=-ABS(NFP(JP,6))
                 RETURN
         ENDIF
         IF(PET.LE.ECUT2) THEN
                 MIST=MIST+1
                 IF(MIST.LT.50) GO TO 155
                 NFT(JT,6)=-ABS(NFT(JT,6))
                 RETURN
         ENDIF
 C               ******** if the remain energy<ECUT the proj or targ
 C                        can not produce jet anymore
 
         WP=PEP+PZP+PET+PZT
         WM=PEP-PZP+PET-PZT
         IF(WP.LT.0.0 .OR. WM.LT.0.0) THEN
                 MISS=MISS+1
                 IF(MISS.LT.50) GO TO 155
                 RETURN
         ENDIF
 C               ********the total W+, W- must be positive
         SW=WP*WM
         AMPX=SQRT((ECUT1-HIPR1(8))**2+PXP**2+PYP**2+0.01)
         AMTX=SQRT((ECUT2-HIPR1(8))**2+PXT**2+PYT**2+0.01)
         SXX=(AMPX+AMTX)**2
         IF(SW.LT.SXX.OR.VINT(43).LT.HIPR1(1)) THEN
                 MISS=MISS+1
                 IF(MISS.LT.50) GO TO 155
                 RETURN
         ENDIF  
 C               ********the proj and targ remnants must have at least
 C                       a CM energy that can produce two strings
 C                       with minimum mass HIPR1(1)(see HIJSFT HIJFRG)
 C
         HINT1(41)=P(7,1)
         HINT1(42)=P(7,2)
         HINT1(43)=P(7,3)
         HINT1(44)=P(7,4)
         HINT1(45)=P(7,5)
         HINT1(46)=SQRT(P(7,1)**2+P(7,2)**2)
         HINT1(51)=P(8,1)
         HINT1(52)=P(8,2)
         HINT1(53)=P(8,3)
         HINT1(54)=P(8,4)
         HINT1(55)=P(8,5)
         HINT1(56)=SQRT(P(8,1)**2+P(8,2)**2) 
         IHNT2(14)=K(7,2)
         IHNT2(15)=K(8,2)
 C
         PINIRAD=(1.0-EXP(-2.0*(VINT(47)-HIDAT(1))))
      &          /(1.0+EXP(-2.0*(VINT(47)-HIDAT(1))))
         I_INIRAD=0
         IF(RLU(0).LE.PINIRAD) I_INIRAD=1
         IF(K(7,2).EQ.-K(8,2)) GO TO 190
         IF(K(7,2).EQ.21.AND.K(8,2).EQ.21.AND.IOPJET.EQ.1) GO TO 190
 C*******************************************************************
 C       gluon  jets are going to be connectd with
 C       the final leading string of quark-aintquark
 C*******************************************************************
         JFLG=2
         JPP=0
         LPQ=0
         LPB=0
         JTT=0
         LTQ=0
         LTB=0
         IS7=0
         IS8=0
         HINT1(47)=0.0
         HINT1(48)=0.0
         HINT1(49)=0.0
         HINT1(50)=0.0
         HINT1(67)=0.0
         HINT1(68)=0.0
         HINT1(69)=0.0
         HINT1(70)=0.0
         DO 180 I=9,N
            IF(K(I,3).EQ.1 .OR. K(I,3).EQ.2.OR.
      &                   ABS(K(I,2)).GT.30) GO TO 180
 C************************************************************
            IF(K(I,3).EQ.7) THEN
               HINT1(47)=HINT1(47)+P(I,1)
               HINT1(48)=HINT1(48)+P(I,2)
               HINT1(49)=HINT1(49)+P(I,3)
               HINT1(50)=HINT1(50)+P(I,4)
            ENDIF
            IF(K(I,3).EQ.8) THEN
               HINT1(67)=HINT1(67)+P(I,1)
               HINT1(68)=HINT1(68)+P(I,2)
               HINT1(69)=HINT1(69)+P(I,3)
               HINT1(70)=HINT1(70)+P(I,4)
            ENDIF
 C************************modifcation made on Apr 10. 1996*****
            IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
               NDR=NDR+1
               IADR(NDR,1)=JP
               IADR(NDR,2)=JT
               KFDR(NDR)=K(I,2)
               PDR(NDR,1)=P(I,1)
               PDR(NDR,2)=P(I,2)
               PDR(NDR,3)=P(I,3)
               PDR(NDR,4)=P(I,4)
               PDR(NDR,5)=P(I,5)
 C************************************************************
               GO TO 180
 C************************correction made on Oct. 14,1994*****
            ENDIF
            IF(K(I,3).EQ.7.OR.K(I,3).EQ.3) THEN
               IF(K(I,3).EQ.7.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(7,2)
      &               .AND.IS7.EQ.0) THEN
                  PP(JP,10)=P(I,1)
                  PP(JP,11)=P(I,2)
                  PP(JP,12)=P(I,3)
                  PZP=PZP+P(I,3)
                  PEP=PEP+P(I,4)
                  NFP(JP,10)=1
                  IS7=1
                  GO TO 180
               ENDIF
               IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
      &                               I_INIRAD.EQ.0)) THEN
                  PXP=PXP+P(I,1)
                  PYP=PYP+P(I,2)
                  PZP=PZP+P(I,3)
                  PEP=PEP+P(I,4)
                  GO TO 180 
               ENDIF
               JPP=JPP+1
               IP(JPP,1)=I
               IP(JPP,2)=0
               IF(K(I,2).NE.21) THEN
                  IF(K(I,2).GT.0) THEN
                     LPQ=LPQ+1
                     IPQ(LPQ)=JPP
                     IP(JPP,2)=LPQ
                  ELSE IF(K(I,2).LT.0) THEN
                     LPB=LPB+1
                     IPB(LPB)=JPP
                     IP(JPP,2)=-LPB
                  ENDIF
               ENDIF
            ELSE IF(K(I,3).EQ.8.OR.K(I,3).EQ.4) THEN
               IF(K(I,3).EQ.8.AND.K(I,2).NE.21.AND.K(I,2).EQ.K(8,2)
      &                          .AND.IS8.EQ.0) THEN
                  PT(JT,10)=P(I,1)
                  PT(JT,11)=P(I,2)
                  PT(JT,12)=P(I,3)
                  PZT=PZT+P(I,3)
                  PET=PET+P(I,4)
                  NFT(JT,10)=1
                  IS8=1
                  GO TO 180
               ENDIF                     
               IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
      &                             I_INIRAD.EQ.0)) THEN
                  PXT=PXT+P(I,1)
                  PYT=PYT+P(I,2)
                  PZT=PZT+P(I,3)
                  PET=PET+P(I,4)
                  GO TO 180
               ENDIF
               JTT=JTT+1
               IT(JTT,1)=I
               IT(JTT,2)=0
               IF(K(I,2).NE.21) THEN
                  IF(K(I,2).GT.0) THEN
                     LTQ=LTQ+1
                     ITQ(LTQ)=JTT
                     IT(JTT,2)=LTQ
                  ELSE IF(K(I,2).LT.0) THEN
                     LTB=LTB+1
                     ITB(LTB)=JTT
                     IT(JTT,2)=-LTB
                  ENDIF
               ENDIF
            ENDIF
  180    CONTINUE
 c
 c
         IF(LPQ.NE.LPB .OR. LTQ.NE.LTB) THEN
                 MISS=MISS+1
                 IF(MISS.LE.50) GO TO 155
                 WRITE(6,*) ' Q -QBAR NOT MATCHED IN HIJHRD'
                 JFLG=0
                 RETURN
         ENDIF
 C****The following will rearrange the partons so that a quark is***
 C****allways followed by an anti-quark ****************************
 
         J=0
 181     J=J+1
         IF(J.GT.JPP) GO TO 182
         IF(IP(J,2).EQ.0) THEN
                 GO TO 181
         ELSE IF(IP(J,2).NE.0) THEN
                 LP=ABS(IP(J,2))
                 IP1=IP(J,1)
                 IP2=IP(J,2)
                 IP(J,1)=IP(IPQ(LP),1)
                 IP(J,2)=IP(IPQ(LP),2)
                 IP(IPQ(LP),1)=IP1
                 IP(IPQ(LP),2)=IP2
                 IF(IP2.GT.0) THEN
                         IPQ(IP2)=IPQ(LP)
                 ELSE IF(IP2.LT.0) THEN
                         IPB(-IP2)=IPQ(LP)
                 ENDIF
 C               ********replace J with a quark
                 IP1=IP(J+1,1)
                 IP2=IP(J+1,2)
                 IP(J+1,1)=IP(IPB(LP),1)
                 IP(J+1,2)=IP(IPB(LP),2)
                 IP(IPB(LP),1)=IP1
                 IP(IPB(LP),2)=IP2
                 IF(IP2.GT.0) THEN
                         IPQ(IP2)=IPB(LP)
                 ELSE IF(IP2.LT.0) THEN
                         IPB(-IP2)=IPB(LP)
                 ENDIF
 C               ******** replace J+1 with anti-quark
                 J=J+1
                 GO TO 181
         ENDIF
 
 182     J=0
 183     J=J+1
         IF(J.GT.JTT) GO TO 184
         IF(IT(J,2).EQ.0) THEN
                 GO TO 183
         ELSE IF(IT(J,2).NE.0) THEN
                 LT=ABS(IT(J,2))
                 IT1=IT(J,1)
                 IT2=IT(J,2)
                 IT(J,1)=IT(ITQ(LT),1)
                 IT(J,2)=IT(ITQ(LT),2)
                 IT(ITQ(LT),1)=IT1
                 IT(ITQ(LT),2)=IT2
                 IF(IT2.GT.0) THEN
                         ITQ(IT2)=ITQ(LT)
                 ELSE IF(IT2.LT.0) THEN
                         ITB(-IT2)=ITQ(LT)
                 ENDIF
 C               ********replace J with a quark
                 IT1=IT(J+1,1)
                 IT2=IT(J+1,2)
                 IT(J+1,1)=IT(ITB(LT),1)
                 IT(J+1,2)=IT(ITB(LT),2)
                 IT(ITB(LT),1)=IT1
                 IT(ITB(LT),2)=IT2
                 IF(IT2.GT.0) THEN
                         ITQ(IT2)=ITB(LT)
                 ELSE IF(IT2.LT.0) THEN
                         ITB(-IT2)=ITB(LT)
                 ENDIF
 C               ******** replace J+1 with anti-quark
                 J=J+1
                 GO TO 183
 
         ENDIF
 
 184     CONTINUE
         IF(NPJ(JP)+JPP.GT.MXJT.OR.NTJ(JT)+JTT.GT.MXJT) THEN
                 JFLG=0
                 WRITE(6,*) 'number of partons per string exceeds'
                 WRITE(6,*) 'the common block size'
                 RETURN
         ENDIF
 C                       ********check the bounds of common blocks
         DO 186 J=1,JPP
                 KFPJ(JP,NPJ(JP)+J)=K(IP(J,1),2)
                 PJPX(JP,NPJ(JP)+J)=P(IP(J,1),1)
                 PJPY(JP,NPJ(JP)+J)=P(IP(J,1),2)
                 PJPZ(JP,NPJ(JP)+J)=P(IP(J,1),3)
                 PJPE(JP,NPJ(JP)+J)=P(IP(J,1),4)
                 PJPM(JP,NPJ(JP)+J)=P(IP(J,1),5)
 186     CONTINUE
         NPJ(JP)=NPJ(JP)+JPP
         DO 188 J=1,JTT
                 KFTJ(JT,NTJ(JT)+J)=K(IT(J,1),2)
                 PJTX(JT,NTJ(JT)+J)=P(IT(J,1),1)
                 PJTY(JT,NTJ(JT)+J)=P(IT(J,1),2)
                 PJTZ(JT,NTJ(JT)+J)=P(IT(J,1),3)
                 PJTE(JT,NTJ(JT)+J)=P(IT(J,1),4)
                 PJTM(JT,NTJ(JT)+J)=P(IT(J,1),5)
 188     CONTINUE
         NTJ(JT)=NTJ(JT)+JTT
         GO TO 900
 C*****************************************************************
 CThis is the case of a quark-antiquark jet it will fragment alone
 C****************************************************************
 190     JFLG=3
         IF(K(7,2).NE.21.AND.K(8,2).NE.21.AND.
      &                   K(7,2)*K(8,2).GT.0) GO TO 155
         JPP=0
         LPQ=0
         LPB=0
         DO 200 I=9,N
            IF(K(I,3).EQ.1.OR.K(I,3).EQ.2.OR.
      &                  ABS(K(I,2)).GT.30) GO TO 200
                 IF(K(I,2).GT.21.AND.K(I,2).LE.30) THEN
                         NDR=NDR+1
                         IADR(NDR,1)=JP
                         IADR(NDR,2)=JT
                         KFDR(NDR)=K(I,2)
                         PDR(NDR,1)=P(I,1)
                         PDR(NDR,2)=P(I,2)
                         PDR(NDR,3)=P(I,3)
                         PDR(NDR,4)=P(I,4)
                         PDR(NDR,5)=P(I,5)
 C************************************************************
                         GO TO 200
 C************************correction made on Oct. 14,1994*****
                 ENDIF
                 IF(K(I,3).EQ.3.AND.(K(I,2).NE.21.OR.
      &                              I_INIRAD.EQ.0)) THEN
                         PXP=PXP+P(I,1)
                         PYP=PYP+P(I,2)
                         PZP=PZP+P(I,3)
                         PEP=PEP+P(I,4)
                         GO TO 200
                 ENDIF
                 IF(K(I,3).EQ.4.AND.(K(I,2).NE.21.OR.
      &                                I_INIRAD.EQ.0)) THEN
                         PXT=PXT+P(I,1)
                         PYT=PYT+P(I,2)
                         PZT=PZT+P(I,3)
                         PET=PET+P(I,4)
                         GO TO 200
                 ENDIF
                 JPP=JPP+1
                 IP(JPP,1)=I
                 IP(JPP,2)=0
                 IF(K(I,2).NE.21) THEN
                         IF(K(I,2).GT.0) THEN
                                 LPQ=LPQ+1
                                 IPQ(LPQ)=JPP
                                 IP(JPP,2)=LPQ
                         ELSE IF(K(I,2).LT.0) THEN
                                 LPB=LPB+1
                                 IPB(LPB)=JPP
                                 IP(JPP,2)=-LPB
                         ENDIF
                 ENDIF
 200     CONTINUE
         IF(LPQ.NE.LPB) THEN
            MISS=MISS+1
            IF(MISS.LE.50) GO TO 155
            WRITE(6,*) LPQ,LPB, 'Q-QBAR NOT CONSERVED OR NOT MATCHED'
            JFLG=0
            RETURN
         ENDIF
 
 C**** The following will rearrange the partons so that a quark is***
 C**** allways followed by an anti-quark ****************************
         J=0
 220     J=J+1
         IF(J.GT.JPP) GO TO 222
         IF(IP(J,2).EQ.0) GO TO 220
                 LP=ABS(IP(J,2))
                 IP1=IP(J,1)
                 IP2=IP(J,2)
                 IP(J,1)=IP(IPQ(LP),1)
                 IP(J,2)=IP(IPQ(LP),2)
                 IP(IPQ(LP),1)=IP1
                 IP(IPQ(LP),2)=IP2
                 IF(IP2.GT.0) THEN
                         IPQ(IP2)=IPQ(LP)
                 ELSE IF(IP2.LT.0) THEN
                         IPB(-IP2)=IPQ(LP)
                 ENDIF
                 IPQ(LP)=J
 C               ********replace J with a quark
                 IP1=IP(J+1,1)
                 IP2=IP(J+1,2)
                 IP(J+1,1)=IP(IPB(LP),1)
                 IP(J+1,2)=IP(IPB(LP),2)
                 IP(IPB(LP),1)=IP1
                 IP(IPB(LP),2)=IP2
                 IF(IP2.GT.0) THEN
                         IPQ(IP2)=IPB(LP)
                 ELSE IF(IP2.LT.0) THEN
                         IPB(-IP2)=IPB(LP)
                 ENDIF
 C               ******** replace J+1 with an anti-quark
                 IPB(LP)=J+1
                 J=J+1
                 GO TO 220
 
 222     CONTINUE
         IF(LPQ.GE.1) THEN
                 DO 240 L0=2,LPQ
                         IP1=IP(2*L0-3,1)
                         IP2=IP(2*L0-3,2)
                         IP(2*L0-3,1)=IP(IPQ(L0),1)
                         IP(2*L0-3,2)=IP(IPQ(L0),2)
                         IP(IPQ(L0),1)=IP1
                         IP(IPQ(L0),2)=IP2
                         IF(IP2.GT.0) THEN
                                 IPQ(IP2)=IPQ(L0)
                         ELSE IF(IP2.LT.0) THEN
                                 IPB(-IP2)=IPQ(L0)
                         ENDIF
                         IPQ(L0)=2*L0-3
 C
                         IP1=IP(2*L0-2,1)
                         IP2=IP(2*L0-2,2)
                         IP(2*L0-2,1)=IP(IPB(L0),1)
                         IP(2*L0-2,2)=IP(IPB(L0),2)
                         IP(IPB(L0),1)=IP1
                         IP(IPB(L0),2)=IP2
                         IF(IP2.GT.0) THEN
                                 IPQ(IP2)=IPB(L0)
                         ELSE IF(IP2.LT.0) THEN
                                 IPB(-IP2)=IPB(L0)
                         ENDIF
                         IPB(L0)=2*L0-2
 240             CONTINUE
 C               ********move all the qqbar pair to the front of 
 C                               the list, except the first pair
                 IP1=IP(2*LPQ-1,1)
                 IP2=IP(2*LPQ-1,2)
                 IP(2*LPQ-1,1)=IP(IPQ(1),1)
                 IP(2*LPQ-1,2)=IP(IPQ(1),2)
                 IP(IPQ(1),1)=IP1
                 IP(IPQ(1),2)=IP2
                 IF(IP2.GT.0) THEN
                         IPQ(IP2)=IPQ(1)
                 ELSE IF(IP2.LT.0) THEN
                         IPB(-IP2)=IPQ(1)
                 ENDIF
                 IPQ(1)=2*LPQ-1
 C               ********move the first quark to the beginning of
 C                               the last string system
                 IP1=IP(JPP,1)
                 IP2=IP(JPP,2)
                 IP(JPP,1)=IP(IPB(1),1)
                 IP(JPP,2)=IP(IPB(1),2)
                 IP(IPB(1),1)=IP1
                 IP(IPB(1),2)=IP2
                 IF(IP2.GT.0) THEN
                         IPQ(IP2)=IPB(1)
                 ELSE IF(IP2.LT.0) THEN
                         IPB(-IP2)=IPB(1)
                 ENDIF
                 IPB(1)=JPP
 C               ********move the first anti-quark to the end of the 
 C                       last string system
         ENDIF
         IF(NSG.GE.MXSG) THEN
            JFLG=0
            WRITE(6,*) 'number of jets forming single strings exceeds'
            WRITE(6,*) 'the common block size'
            RETURN
         ENDIF
         IF(JPP.GT.MXSJ) THEN
            JFLG=0
            WRITE(6,*) 'number of partons per single jet system'
            WRITE(6,*) 'exceeds the common block size'
            RETURN
         ENDIF
 C               ********check the bounds of common block size
         NSG=NSG+1
         NJSG(NSG)=JPP
         IASG(NSG,1)=JP
         IASG(NSG,2)=JT
         IASG(NSG,3)=0
         DO 300 I=1,JPP
                 K1SG(NSG,I)=2
                 K2SG(NSG,I)=K(IP(I,1),2)
                 IF(K2SG(NSG,I).LT.0) K1SG(NSG,I)=1
                 PXSG(NSG,I)=P(IP(I,1),1)
                 PYSG(NSG,I)=P(IP(I,1),2)
                 PZSG(NSG,I)=P(IP(I,1),3)
                 PESG(NSG,I)=P(IP(I,1),4)
                 PMSG(NSG,I)=P(IP(I,1),5)
 300     CONTINUE
         K1SG(NSG,1)=2
         K1SG(NSG,JPP)=1
 C******* reset the energy-momentum of incoming particles ********
 900     PP(JP,1)=PXP
         PP(JP,2)=PYP
         PP(JP,3)=PZP
         PP(JP,4)=PEP
         PP(JP,5)=0.0
         PT(JT,1)=PXT
         PT(JT,2)=PYT
         PT(JT,3)=PZT
         PT(JT,4)=PET
         PT(JT,5)=0.0
 
         NFP(JP,6)=NFP(JP,6)+1
         NFT(JT,6)=NFT(JT,6)+1
         RETURN
 C
 1000    JFLG=-1
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) 'Fatal HIJHRD error'
         WRITE(6,*) JP, ' proj E+,E-',EPP,EPM,' status',NFP(JP,5)
         WRITE(6,*) JT, ' targ E+,E_',ETP,ETM,' status',NFT(JT,5)
         RETURN
         END
 C
 C
 C
 C
 C
         SUBROUTINE JETINI(JP,JT,I_TRIG)
 C*******Initialize PYTHIA for jet production**********************
 C       I_TRIG=0: for normal processes
 C       I_TRIG=1: for triggered processes
 C       JP: sequence number of the projectile
 C       JT: sequence number of the target
 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 this subroutine we use the following
 C     catalogue for different type 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*****************************************************************
         CHARACTER BEAM*16,TARG*16
         DIMENSION XSEC0(8,0:200),COEF0(8,200,20),INI(8),
      &          MINT44(8),MINT45(8)
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIPYINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
         SAVE  /HIPYINT/
 C
         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
         SAVE  /LUDAT1/
         COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)    
         SAVE  /LUDAT3/
         COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
         SAVE  /PYSUBS/
         COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
         SAVE  /PYPARS/
         COMMON/PYINT1/MINT(400),VINT(400)
         SAVE  /PYINT1/
         COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
         SAVE  /PYINT2/
         COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
         SAVE  /PYINT5/
         DATA INI/8*0/I_LAST/-1/
 
         save INI, I_LAST    !uzhi
 
         IHNT2(11)=JP
         IHNT2(12)=JT
         IF(IHNT2(5).NE.0 .AND. IHNT2(6).NE.0) THEN
            I_TYPE=1
         ELSE IF(IHNT2(5).NE.0 .AND. IHNT2(6).EQ.0) THEN
            I_TYPE=1
            IF(NFT(JT,4).EQ.2112) I_TYPE=2
         ELSE IF(IHNT2(5).EQ.0 .AND. IHNT2(6).NE.0) THEN
            I_TYPE=1
            IF(NFP(JP,4).EQ.2112) I_TYPE=2
         ELSE
            IF(NFP(JP,4).EQ.2212 .AND. NFT(JT,4).EQ.2212) THEN
               I_TYPE=1
            ELSE IF(NFP(JP,4).EQ.2212 .AND. NFT(JT,4).EQ.2112) THEN
               I_TYPE=2
            ELSE IF(NFP(JP,4).EQ.2112 .AND. NFT(JT,4).EQ.2212) THEN
               I_TYPE=3
            ELSE
               I_TYPE=4
            ENDIF
         ENDIF
 c
         IF(I_TRIG.NE.0) GO TO 160
         IF(I_TRIG.EQ.I_LAST) GO TO 150
         MSTP(2)=2
 c                       ********second order running alpha_strong
         MSTP(33)=1
         PARP(31)=HIPR1(17)
 C                       ********inclusion of K factor
         MSTP(51)=3
 C                       ********Duke-Owens set 1 structure functions
         MSTP(61)=1
 C                       ********INITIAL STATE RADIATION
         MSTP(71)=1
 C                       ********FINAL STATE RADIATION
         IF(IHPR2(2).EQ.0.OR.IHPR2(2).EQ.2) MSTP(61)=0
         IF(IHPR2(2).EQ.0.OR.IHPR2(2).EQ.1) MSTP(71)=0
 c
         MSTP(81)=0
 C                       ******** NO MULTIPLE INTERACTION
         MSTP(82)=1
 C                       *******STRUCTURE OF MUTLIPLE INTERACTION
         MSTP(111)=0
 C               ********frag off(have to be done by local call)
         IF(IHPR2(10).EQ.0) MSTP(122)=0
 C               ********No printout of initialization information
         PARP(81)=HIPR1(8)
         CKIN(5)=HIPR1(8)
         CKIN(3)=HIPR1(8)
         CKIN(4)=HIPR1(9)
         IF(HIPR1(9).LE.HIPR1(8)) CKIN(4)=-1.0
         CKIN(9)=-10.0
         CKIN(10)=10.0
         MSEL=0
         DO 100 ISUB=1,200
            MSUB(ISUB)=0
  100    CONTINUE
         MSUB(11)=1
         MSUB(12)=1
         MSUB(13)=1
         MSUB(28)=1
         MSUB(53)=1
         MSUB(68)=1
         MSUB(81)=1
         MSUB(82)=1
         DO 110 J=1,MIN(8,MDCY(21,3))
  110    MDME(MDCY(21,2)+J-1,1)=0
         ISEL=4
         IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
         MDME(MDCY(21,2)+ISEL-1,1)=1
 C                       ********QCD subprocesses
         MSUB(14)=1
         MSUB(18)=1
         MSUB(29)=1
 C                       ******* direct photon production
  150    IF(INI(I_TYPE).NE.0) GO TO 800
         GO TO 400
 C
 C       *****triggered subprocesses, jet, photon, heavy quark and DY
 C
  160    I_TYPE=4+I_TYPE
         IF(I_TRIG.EQ.I_LAST) GO TO 260
         PARP(81)=ABS(HIPR1(10))-0.25
         CKIN(5)=ABS(HIPR1(10))-0.25
         CKIN(3)=ABS(HIPR1(10))-0.25
         CKIN(4)=ABS(HIPR1(10))+0.25
         IF(HIPR1(10).LT.HIPR1(8)) CKIN(4)=-1.0
 c
         MSEL=0
         DO 101 ISUB=1,200
            MSUB(ISUB)=0
  101    CONTINUE
         IF(IHPR2(3).EQ.1) THEN
            MSUB(11)=1
            MSUB(12)=1
            MSUB(13)=1
            MSUB(28)=1
            MSUB(53)=1
            MSUB(68)=1
            MSUB(81)=1
            MSUB(82)=1
            MSUB(14)=1
            MSUB(18)=1
            MSUB(29)=1
            DO 102 J=1,MIN(8,MDCY(21,3))
  102       MDME(MDCY(21,2)+J-1,1)=0
            ISEL=4
            IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
            MDME(MDCY(21,2)+ISEL-1,1)=1
 C                       ********QCD subprocesses
         ELSE IF(IHPR2(3).EQ.2) THEN
            MSUB(14)=1
            MSUB(18)=1
            MSUB(29)=1
 C               ********Direct photon production
 c               q+qbar->g+gamma,q+qbar->gamma+gamma, q+g->q+gamma
         ELSE IF(IHPR2(3).EQ.3) THEN
            CKIN(3)=MAX(0.0,HIPR1(10))
            CKIN(5)=HIPR1(8)
            PARP(81)=HIPR1(8)
            MSUB(81)=1
            MSUB(82)=1
            DO 105 J=1,MIN(8,MDCY(21,3))
  105       MDME(MDCY(21,2)+J-1,1)=0
            ISEL=4
            IF(HINT1(1).GE.20.0 .and. IHPR2(18).EQ.1) ISEL=5
            MDME(MDCY(21,2)+ISEL-1,1)=1
 C             **********Heavy quark production
         ENDIF
 260     IF(INI(I_TYPE).NE.0) GO TO 800
 C
 C
 400     INI(I_TYPE)=1
         IF(IHPR2(10).EQ.0) MSTP(122)=0
         IF(NFP(JP,4).EQ.2212) THEN
                 BEAM='P'
         ELSE IF(NFP(JP,4).EQ.-2212) THEN
                 BEAM='P~'
         ELSE IF(NFP(JP,4).EQ.2112) THEN
                 BEAM='N'
         ELSE IF(NFP(JP,4).EQ.-2112) THEN
                 BEAM='N~'
         ELSE IF(NFP(JP,4).EQ.211) THEN
                 BEAM='PI+'
         ELSE IF(NFP(JP,4).EQ.-211) THEN
                 BEAM='PI-'
         ELSE IF(NFP(JP,4).EQ.321) THEN
                 BEAM='PI+'
         ELSE IF(NFP(JP,4).EQ.-321) THEN
                 BEAM='PI-'
         ELSE
                 WRITE(6,*) 'unavailable beam type', NFP(JP,4)
         ENDIF
         IF(NFT(JT,4).EQ.2212) THEN
                 TARG='P'
         ELSE IF(NFT(JT,4).EQ.-2212) THEN
                 TARG='P~'
         ELSE IF(NFT(JT,4).EQ.2112) THEN
                 TARG='N'
         ELSE IF(NFT(JT,4).EQ.-2112) THEN
                 TARG='N~'
         ELSE IF(NFT(JT,4).EQ.211) THEN
                 TARG='PI+'
         ELSE IF(NFT(JT,4).EQ.-211) THEN
                 TARG='PI-'
         ELSE IF(NFT(JT,4).EQ.321) THEN
                 TARG='PI+'
         ELSE IF(NFT(JT,4).EQ.-321) THEN
                 TARG='PI-'
         ELSE
                 WRITE(6,*) 'unavailable target type', NFT(JT,4)
         ENDIF
 C
         IHNT2(16)=1
 C       ******************indicate for initialization use when
 C                         structure functions are called in PYTHIA
 C
         CALL PYINIT('CMS',BEAM,TARG,HINT1(1))
         MINT4=MINT(44)
         MINT5=MINT(45)
         MINT44(I_TYPE)=MINT(44)
         MINT45(I_TYPE)=MINT(45)
         ATXS(0)=XSEC(0,1)
         XSEC0(I_TYPE,0)=XSEC(0,1)
         DO 500 I=1,200
                 ATXS(I)=XSEC(I,1)
                 XSEC0(I_TYPE,I)=XSEC(I,1)
                 DO 500 J=1,20
                         ATCO(I,J)=COEF(I,J)
                         COEF0(I_TYPE,I,J)=COEF(I,J)
 500     CONTINUE
 C
         IHNT2(16)=0
 C
         RETURN
 C               ********Store the initialization information for
 C                               late use
 C
 C
 800     MINT(44)=MINT44(I_TYPE)
         MINT(45)=MINT45(I_TYPE)
         MINT4=MINT(44)
         MINT5=MINT(45)
         XSEC(0,1)=XSEC0(I_TYPE,0)
         ATXS(0)=XSEC(0,1)
         DO 900 I=1,200
                 XSEC(I,1)=XSEC0(I_TYPE,I)
                 ATXS(I)=XSEC(I,1)
         DO 900 J=1,20
                 COEF(I,J)=COEF0(I_TYPE,I,J)
                 ATCO(I,J)=COEF(I,J)
 900     CONTINUE
         I_LAST=I_TRIG
         MINT(11)=NFP(JP,4)
         MINT(12)=NFT(JT,4)
         RETURN
         END
 C            
 C
 C
         SUBROUTINE HIJINI
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &                PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
         COMMON/HIJJET4/NDR,IADR(900,2),KFDR(900),PDR(900,5)
         SAVE  /HIJJET4/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C****************Reset the momentum of initial particles************
 C             and assign flavors to the proj and targ string       *
 C*******************************************************************
         NSG=0
         NDR=0
         IPP=2212
         IPT=2212
         IF(IHNT2(5).NE.0) IPP=IHNT2(5)
         IF(IHNT2(6).NE.0) IPT=IHNT2(6)
 C               ********in case the proj or targ is a hadron.
 C
         DO 100 I=1,IHNT2(1)
         PP(I,1)=0.0
         PP(I,2)=0.0
         PP(I,3)=SQRT(HINT1(1)**2/4.0-HINT1(8)**2)
         PP(I,4)=HINT1(1)/2
         PP(I,5)=HINT1(8)
         PP(I,6)=0.0
         PP(I,7)=0.0
         PP(I,8)=0.0
         PP(I,9)=0.0
         PP(I,10)=0.0
         NFP(I,3)=IPP
         NFP(I,4)=IPP
         NFP(I,5)=0
         NFP(I,6)=0
         NFP(I,7)=0
         NFP(I,8)=0
         NFP(I,9)=0
         NFP(I,10)=0
         NFP(I,11)=0
         NPJ(I)=0
         IF(I.GT.ABS(IHNT2(2))) NFP(I,3)=2112
         CALL ATTFLV(NFP(I,3),IDQ,IDQQ)
         NFP(I,1)=IDQ
         NFP(I,2)=IDQQ
         NFP(I,15)=-1
         IF(ABS(IDQ).GT.1000.OR.(ABS(IDQ*IDQQ).LT.100.AND.
      &          RLU(0).LT.0.5)) NFP(I,15)=1
         PP(I,14)=ULMASS(IDQ)
         PP(I,15)=ULMASS(IDQQ)
 100     CONTINUE
 C
         DO 200 I=1,IHNT2(3)
         PT(I,1)=0.0
         PT(I,2)=0.0
         PT(I,3)=-SQRT(HINT1(1)**2/4.0-HINT1(9)**2)
         PT(I,4)=HINT1(1)/2.0
         PT(I,5)=HINT1(9)
         PT(I,6)=0.0
         PT(I,7)=0.0
         PT(I,8)=0.0
         PT(I,9)=0.0
         PT(I,10)=0.0
         NFT(I,3)=IPT
         NFT(I,4)=IPT
         NFT(I,5)=0
         NFT(I,6)=0
         NFT(I,7)=0
         NFT(I,8)=0
         NFT(I,9)=0
         NFT(I,10)=0
         NFT(I,11)=0
         NTJ(I)=0
         IF(I.GT.ABS(IHNT2(4))) NFT(I,3)=2112
         CALL ATTFLV(NFT(I,3),IDQ,IDQQ)
         NFT(I,1)=IDQ
         NFT(I,2)=IDQQ
         NFT(I,15)=1
         IF(ABS(IDQ).GT.1000.OR.(ABS(IDQ*IDQQ).LT.100.AND.
      &                  RLU(0).LT.0.5)) NFT(I,15)=-1
         PT(I,14)=ULMASS(IDQ)
         PT(I,15)=ULMASS(IDQQ)
 200     CONTINUE
         RETURN
         END
 C
 C
 C
         SUBROUTINE ATTFLV(ID,IDQ,IDQQ)
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C
         IF(ABS(ID).LT.100) THEN
                 NSIGN=1
                 IDQ=ID/100
                 IDQQ=-ID/10+IDQ*10
                 IF(ABS(IDQ).EQ.3) NSIGN=-1
                 IDQ=NSIGN*IDQ
                 IDQQ=NSIGN*IDQQ
                 IF(IDQ.LT.0) THEN
                         ID0=IDQ
                         IDQ=IDQQ
                         IDQQ=ID0
                 ENDIF
                 RETURN
         ENDIF
 C               ********return ID of quark(IDQ) and anti-quark(IDQQ)
 C                       for pions and kaons
 c
 C       Return LU ID for quarks and diquarks for proton(ID=2212) 
 C       anti-proton(ID=-2212) and nuetron(ID=2112)
 C       LU ID for d=1,u=2, (ud)0=2101, (ud)1=2103, 
 C       (dd)1=1103,(uu)1=2203.
 C       Use SU(6)  weight  proton=1/3d(uu)1 + 1/6u(ud)1 + 1/2u(ud)0
 C                         nurtron=1/3u(dd)1 + 1/6d(ud)1 + 1/2d(ud)0
 C 
         IDQ=2
         IF(ABS(ID).EQ.2112) IDQ=1
         IDQQ=2101
         X=RLU(0)
         IF(X.LE.0.5) GO TO 30
         IF(X.GT.0.666667) GO TO 10
         IDQQ=2103
         GO TO 30
 10      IDQ=1
         IDQQ=2203
         IF(ABS(ID).EQ.2112) THEN
                 IDQ=2
                 IDQQ=1103
         ENDIF
 30      IF(ID.LT.0) THEN
                 ID00=IDQQ
                 IDQQ=-IDQ
                 IDQ=-ID00
         ENDIF
         RETURN
         END     
 C
 C*******************************************************************
 C       This subroutine performs elastic scatterings and possible 
 C       elastic cascading within their own nuclei
 c*******************************************************************
         SUBROUTINE HIJCSC(JP,JT)
         DIMENSION PSC1(5),PSC2(5)
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         IF(JP.EQ.0 .OR. JT.EQ.0) GO TO 25
         DO 10 I=1,5
         PSC1(I)=PP(JP,I)
         PSC2(I)=PT(JT,I)
 10      CONTINUE
         CALL HIJELS(PSC1,PSC2)
         DPP1=PSC1(1)-PP(JP,1)
         DPP2=PSC1(2)-PP(JP,2)
         DPT1=PSC2(1)-PT(JT,1)
         DPT2=PSC2(2)-PT(JT,2)
         PP(JP,6)=PP(JP,6)+DPP1/2.0
         PP(JP,7)=PP(JP,7)+DPP2/2.0
         PP(JP,8)=PP(JP,8)+DPP1/2.0
         PP(JP,9)=PP(JP,9)+DPP2/2.0
         PT(JT,6)=PT(JT,6)+DPT1/2.0
         PT(JT,7)=PT(JT,7)+DPT2/2.0
         PT(JT,8)=PT(JT,8)+DPT1/2.0
         PT(JT,9)=PT(JT,9)+DPT2/2.0
         DO 20 I=1,4
         PP(JP,I)=PSC1(I)
         PT(JT,I)=PSC2(I)
 20      CONTINUE
         NFP(JP,5)=MAX(1,NFP(JP,5))
         NFT(JT,5)=MAX(1,NFT(JT,5))
 C               ********Perform elastic scattering between JP and JT
         RETURN
 C               ********The following is for possible elastic cascade
 c
 25      IF(JP.EQ.0) GO TO 45
         PABS=SQRT(PP(JP,1)**2+PP(JP,2)**2+PP(JP,3)**2)
         BX=PP(JP,1)/PABS
         BY=PP(JP,2)/PABS
         BZ=PP(JP,3)/PABS
         DO 40 I=1,IHNT2(1)
                 IF(I.EQ.JP) GO TO 40
                 DX=YP(1,I)-YP(1,JP)
                 DY=YP(2,I)-YP(2,JP)
                 DZ=YP(3,I)-YP(3,JP)
                 DIS=DX*BX+DY*BY+DZ*BZ
                 IF(DIS.LE.0) GO TO 40
                 BB=DX**2+DY**2+DZ**2-DIS**2
                 R2=BB*HIPR1(40)/HIPR1(31)/0.1
 C               ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
                 GS=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
      &                  *ROMG(R2))**2
                 GS0=1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
      &                  *ROMG(0.0))**2
                 IF(RLU(0).GT.GS/GS0) GO TO 40
                 DO 30 K=1,5
                         PSC1(K)=PP(JP,K)
                         PSC2(K)=PP(I,K)
 30              CONTINUE
                 CALL HIJELS(PSC1,PSC2)
                 DPP1=PSC1(1)-PP(JP,1)
                 DPP2=PSC1(2)-PP(JP,2)
                 DPT1=PSC2(1)-PP(I,1)
                 DPT2=PSC2(2)-PP(I,2)
                 PP(JP,6)=PP(JP,6)+DPP1/2.0
                 PP(JP,7)=PP(JP,7)+DPP2/2.0
                 PP(JP,8)=PP(JP,8)+DPP1/2.0
                 PP(JP,9)=PP(JP,9)+DPP2/2.0
                 PP(I,6)=PP(I,6)+DPT1/2.0
                 PP(I,7)=PP(I,7)+DPT2/2.0
                 PP(I,8)=PP(I,8)+DPT1/2.0
                 PP(I,9)=PP(I,9)+DPT2/2.0
                 DO 35 K=1,5
                         PP(JP,K)=PSC1(K)
                         PP(I,K)=PSC2(K)
 35              CONTINUE
                 NFP(I,5)=MAX(1,NFP(I,5))
                 GO TO 45
 40      CONTINUE
 45      IF(JT.EQ.0) GO TO 80
 50      PABS=SQRT(PT(JT,1)**2+PT(JT,2)**2+PT(JT,3)**2)
         BX=PT(JT,1)/PABS
         BY=PT(JT,2)/PABS
         BZ=PT(JT,3)/PABS
         DO 70 I=1,IHNT2(3)
                 IF(I.EQ.JT) GO TO 70
                 DX=YT(1,I)-YT(1,JT)
                 DY=YT(2,I)-YT(2,JT)
                 DZ=YT(3,I)-YT(3,JT)
                 DIS=DX*BX+DY*BY+DZ*BZ
                 IF(DIS.LE.0) GO TO 70
                 BB=DX**2+DY**2+DZ**2-DIS**2
                 R2=BB*HIPR1(40)/HIPR1(31)/0.1
 C               ********mb=0.1*fm, YP is in fm,HIPR1(31) is in mb
                 GS=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
      &                  *ROMG(R2)))**2
                 GS0=(1.0-EXP(-(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
      &                  *ROMG(0.0)))**2
                 IF(RLU(0).GT.GS/GS0) GO TO 70
                 DO 60 K=1,5
                         PSC1(K)=PT(JT,K)
                         PSC2(K)=PT(I,K)
 60              CONTINUE
                 CALL HIJELS(PSC1,PSC2)
                 DPP1=PSC1(1)-PT(JT,1)
                 DPP2=PSC1(2)-PT(JT,2)
                 DPT1=PSC2(1)-PT(I,1)
                 DPT2=PSC2(2)-PT(I,2)
                 PT(JT,6)=PT(JT,6)+DPP1/2.0
                 PT(JT,7)=PT(JT,7)+DPP2/2.0
                 PT(JT,8)=PT(JT,8)+DPP1/2.0
                 PT(JT,9)=PT(JT,9)+DPP2/2.0
                 PT(I,6)=PT(I,6)+DPT1/2.0
                 PT(I,7)=PT(I,7)+DPT2/2.0
                 PT(I,8)=PT(I,8)+DPT1/2.0
                 PT(I,9)=PT(I,9)+DPT2/2.0
                 DO 65 K=1,5
                         PT(JT,K)=PSC1(K)
                         PT(I,K)=PSC2(K)
 65              CONTINUE
                 NFT(I,5)=MAX(1,NFT(I,5))
                 GO TO 80
 70      CONTINUE
 80      RETURN
         END
 C
 C
 C*******************************************************************
 CThis subroutine performs elastic scattering between two nucleons
 C
 C*******************************************************************
         SUBROUTINE HIJELS(PSC1,PSC2)
         IMPLICIT DOUBLE PRECISION(D)
         DIMENSION PSC1(5),PSC2(5)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
 C
         CC=1.0-HINT1(12)/HINT1(13)
         RR=(1.0-CC)*HINT1(13)/HINT1(12)/(1.0-HIPR1(33))-1.0
         BB=0.5*(3.0+RR+SQRT(9.0+10.0*RR+RR**2))
         EP=SQRT((PSC1(1)-PSC2(1))**2+(PSC1(2)-PSC2(2))**2
      &          +(PSC1(3)-PSC2(3))**2)
         IF(EP.LE.0.1) RETURN
         ELS0=98.0/EP+52.0*(1.0+RR)**2
         PCM1=PSC1(1)+PSC2(1)
         PCM2=PSC1(2)+PSC2(2)
         PCM3=PSC1(3)+PSC2(3)
         ECM=PSC1(4)+PSC2(4)
         AM1=PSC1(5)**2
         AM2=PSC2(5)**2
         AMM=ECM**2-PCM1**2-PCM2**2-PCM3**2
         IF(AMM.LE.PSC1(5)+PSC2(5)) RETURN
 C               ********elastic scattering only when approaching
 C                               to each other
         PMAX=(AMM**2+AM1**2+AM2**2-2.0*AMM*AM1-2.0*AMM*AM2
      &                  -2.0*AM1*AM2)/4.0/AMM
         PMAX=ABS(PMAX)
 20      TT=RLU(0)*MIN(PMAX,1.5)
         ELS=98.0*EXP(-2.8*TT)/EP
      &          +52.0*EXP(-9.2*TT)*(1.0+RR*EXP(-4.6*(BB-1.0)*TT))**2
         IF(RLU(0).GT.ELS/ELS0) GO TO 20
         PHI=2.0*HIPR1(40)*RLU(0)
 C
         DBX=PCM1/ECM
         DBY=PCM2/ECM
         DBZ=PCM3/ECM
         DB=SQRT(DBX**2+DBY**2+DBZ**2)
         IF(DB.GT.0.99999999D0) THEN 
           DBX=DBX*(0.99999999D0/DB) 
           DBY=DBY*(0.99999999D0/DB) 
           DBZ=DBZ*(0.99999999D0/DB) 
           DB=0.99999999D0   
           WRITE(6,*) ' (HIJELS) boost vector too large' 
 C               ********Rescale boost vector if too close to unity. 
         ENDIF   
         DGA=1D0/SQRT(1D0-DB**2)      
 C
         DP1=SQRT(TT)*SIN(PHI)
         DP2=SQRT(TT)*COS(PHI)
         DP3=SQRT(PMAX-TT)
         DP4=SQRT(PMAX+AM1)
         DBP=DBX*DP1+DBY*DP2+DBZ*DP3   
         DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4) 
         PSC1(1)=DP1+DGABP*DBX
         PSC1(2)=DP2+DGABP*DBY  
         PSC1(3)=DP3+DGABP*DBZ  
         PSC1(4)=DGA*(DP4+DBP)    
 C       
         DP1=-SQRT(TT)*SIN(PHI)
         DP2=-SQRT(TT)*COS(PHI)
         DP3=-SQRT(PMAX-TT)
         DP4=SQRT(PMAX+AM2)
         DBP=DBX*DP1+DBY*DP2+DBZ*DP3   
         DGABP=DGA*(DGA*DBP/(1D0+DGA)+DP4) 
         PSC2(1)=DP1+DGABP*DBX
         PSC2(2)=DP2+DGABP*DBY  
         PSC2(3)=DP3+DGABP*DBZ  
         PSC2(4)=DGA*(DP4+DBP)
         RETURN
         END
 C
 C       
 C*******************************************************************
 C                                                                  *
 C               Subroutine HIJSFT                                  *
 C                                                                  *
 C  Scatter two excited strings, JP from proj and JT from target    *
 C*******************************************************************
         SUBROUTINE HIJSFT(JP,JT,JOUT,IERROR)
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
         SAVE  /HIJDAT/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &               PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &               PJPM(300,500),NTJ(300),KFTJ(300,500),
      &               PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &               PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/DPMCOM1/JJP,JJT,AMP,AMT,APX0,ATX0,AMPN,AMTN,AMP0,AMT0,
      &          NFDP,NFDT,WP,WM,SW,XREMP,XREMT,DPKC1,DPKC2,PP11,PP12,
      &          PT11,PT12,PTP2,PTT2
         SAVE  /DPMCOM1/
         COMMON/DPMCOM2/NDPM,KDPM(20,2),PDPM1(20,5),PDPM2(20,5)
         SAVE  /DPMCOM2/
 C*******************************************************************
 C       JOUT-> the number
 C       of hard scatterings preceding this soft collision. 
 C       IHNT2(13)-> 1=
 C       double diffrac 2=single diffrac, 3=non-single diffrac.
 C*******************************************************************
         IERROR=0
         JJP=JP
         JJT=JT
         NDPM=0
         IOPMAIN=0
         IF(JP.GT.IHNT2(1) .OR. JT.GT.IHNT2(3)) RETURN
 
         EPP=PP(JP,4)+PP(JP,3)
         EPM=PP(JP,4)-PP(JP,3)
         ETP=PT(JT,4)+PT(JT,3)
         ETM=PT(JT,4)-PT(JT,3)
 
         WP=EPP+ETP
         WM=EPM+ETM
         SW=WP*WM
 C               ********total W+,W- and center-of-mass energy
 
         IF(WP.LT.0.0 .OR. WM.LT.0.0) GO TO 1000
 
         IF(JOUT.EQ.0) THEN
                 IF(EPP.LT.0.0) GO TO 1000
                 IF(EPM.LT.0.0) GO TO 1000
                 IF(ETP.LT.0.0) GO TO 1000
                 IF(ETM.LT.0.0) GO TO 1000    
                 IF(EPP/(EPM+0.01).LE.ETP/(ETM+0.01)) RETURN
         ENDIF
 C               ********For strings which does not follow a jet-prod,
 C                       scatter only if Ycm(JP)>Ycm(JT). When jets
 C                       are produced just before this collision
 C                       this requirement has already be enforced
 C                       (see SUBROUTINE HIJHRD)
         IHNT2(11)=JP
         IHNT2(12)=JT
 C
 C
 C
         MISS=0
         PKC1=0.0
         PKC2=0.0
         PKC11=0.0
         PKC12=0.0
         PKC21=0.0
         PKC22=0.0
         DPKC11=0.0
         DPKC12=0.0
         DPKC21=0.0
         DPKC22=0.0
         IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
            IF(NFP(JP,10).EQ.1) THEN
               PHI1=ULANGL(PP(JP,10),PP(JP,11))
               PPJET=SQRT(PP(JP,10)**2+PP(JP,11)**2)
               PKC1=PPJET
               PKC11=PP(JP,10)
               PKC12=PP(JP,11)
            ENDIF
            IF(NFT(JT,10).EQ.1) THEN
               PHI2=ULANGL(PT(JT,10),PT(JT,11))
               PTJET=SQRT(PT(JT,10)**2+PT(JT,11)**2)
               PKC2=PTJET
               PKC21=PT(JT,10)
               PKC22=PT(JT,11)
            ENDIF
            IF(IHPR2(4).GT.0.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
               IF(NFP(JP,10).EQ.0) THEN
                  PHI=-PHI2
               ELSE IF(NFT(JT,10).EQ.0) THEN
                  PHI=PHI1
               ELSE
                  PHI=(PHI1+PHI2-HIPR1(40))/2.0
               ENDIF
               BX=HINT1(19)*COS(HINT1(20))
               BY=HINT1(19)*SIN(HINT1(20))
               XP0=YP(1,JP)
               YP0=YP(2,JP)
               XT0=YT(1,JT)+BX
               YT0=YT(2,JT)+BY
               R1=MAX(1.2*IHNT2(1)**0.3333333,
      &               SQRT(XP0**2+YP0**2))
               R2=MAX(1.2*IHNT2(3)**0.3333333,
      &               SQRT((XT0-BX)**2+(YT0-BY)**2))
               IF(ABS(COS(PHI)).LT.1.0E-5) THEN
                  DD1=R1
                  DD2=R1
                  DD3=ABS(BY+SQRT(R2**2-(XP0-BX)**2)-YP0)
                  DD4=ABS(BY-SQRT(R2**2-(XP0-BX)**2)-YP0)
                  GO TO 5
               ENDIF
               BB=2.0*SIN(PHI)*(COS(PHI)*YP0-SIN(PHI)*XP0)
               CC=(YP0**2-R1**2)*COS(PHI)**2+XP0*SIN(PHI)*(
      &                          XP0*SIN(PHI)-2.0*YP0*COS(PHI))
               DD=BB**2-4.0*CC
               IF(DD.LT.0.0) GO TO 10
               XX1=(-BB+SQRT(DD))/2.0
               XX2=(-BB-SQRT(DD))/2.0
               DD1=ABS((XX1-XP0)/COS(PHI))
               DD2=ABS((XX2-XP0)/COS(PHI))
 C                       
               BB=2.0*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)-2.0*BX
               CC=(BX**2+(YT0-BY)**2-R2**2)*COS(PHI)**2+XT0*SIN(PHI)
      &           *(XT0*SIN(PHI)-2.0*COS(PHI)*(YT0-BY))
      &           -2.0*BX*SIN(PHI)*(COS(PHI)*(YT0-BY)-SIN(PHI)*XT0)
               DD=BB**2-4.0*CC
               IF(DD.LT.0.0) GO TO 10
               XX1=(-BB+SQRT(DD))/2.0
               XX2=(-BB-SQRT(DD))/2.0
               DD3=ABS((XX1-XT0)/COS(PHI))
               DD4=ABS((XX2-XT0)/COS(PHI))
 C
  5            DD1=MIN(DD1,DD3)
               DD2=MIN(DD2,DD4)
               IF(DD1.LT.HIPR1(13)) DD1=0.0
               IF(DD2.LT.HIPR1(13)) DD2=0.0
               IF(NFP(JP,10).EQ.1.AND.PPJET.GT.HIPR1(11)) THEN
                  DP1=DD1*HIPR1(14)/2.0
                  DP1=MIN(DP1,PPJET-HIPR1(11))
                  PKC1=PPJET-DP1
                  DPX1=COS(PHI1)*DP1
                  DPY1=SIN(PHI1)*DP1
                  PKC11=PP(JP,10)-DPX1
                  PKC12=PP(JP,11)-DPY1
                  IF(DP1.GT.0.0) THEN
                     CTHEP=PP(JP,12)/SQRT(PP(JP,12)**2+PPJET**2)
                     DPZ1=DP1*CTHEP/SQRT(1.0-CTHEP**2)
                     DPE1=SQRT(DPX1**2+DPY1**2+DPZ1**2)
                     EPPPRM=PP(JP,4)+PP(JP,3)-DPE1-DPZ1
                     EPMPRM=PP(JP,4)-PP(JP,3)-DPE1+DPZ1
                     IF(EPPPRM.LE.0.0.OR.EPMPRM.LE.0.0) GO TO 15
                     EPP=EPPPRM
                     EPM=EPMPRM
                     PP(JP,10)=PKC11
                     PP(JP,11)=PKC12
                     NPJ(JP)=NPJ(JP)+1
                     KFPJ(JP,NPJ(JP))=21
                     PJPX(JP,NPJ(JP))=DPX1
                     PJPY(JP,NPJ(JP))=DPY1
                     PJPZ(JP,NPJ(JP))=DPZ1
                     PJPE(JP,NPJ(JP))=DPE1
                     PJPM(JP,NPJ(JP))=0.0
                     PP(JP,3)=PP(JP,3)-DPZ1
                     PP(JP,4)=PP(JP,4)-DPE1
                  ENDIF
               ENDIF
  15           IF(NFT(JT,10).EQ.1.AND.PTJET.GT.HIPR1(11)) THEN
                  DP2=DD2*HIPR1(14)/2.0
                  DP2=MIN(DP2,PTJET-HIPR1(11))
                  PKC2=PTJET-DP2
                  DPX2=COS(PHI2)*DP2
                  DPY2=SIN(PHI2)*DP2
                  PKC21=PT(JT,10)-DPX2
                  PKC22=PT(JT,11)-DPY2
                  IF(DP2.GT.0.0) THEN
                     CTHET=PT(JT,12)/SQRT(PT(JT,12)**2+PTJET**2)
                     DPZ2=DP2*CTHET/SQRT(1.0-CTHET**2)
                     DPE2=SQRT(DPX2**2+DPY2**2+DPZ2**2)
                     ETPPRM=PT(JT,4)+PT(JT,3)-DPE2-DPZ2
                     ETMPRM=PT(JT,4)-PT(JT,3)-DPE2+DPZ2
                     IF(ETPPRM.LE.0.0.OR.ETMPRM.LE.0.0) GO TO 16
                     ETP=ETPPRM
                     ETM=ETMPRM
                     PT(JT,10)=PKC21
                     PT(JT,11)=PKC22
                     NTJ(JT)=NTJ(JT)+1
                     KFTJ(JT,NTJ(JT))=21
                     PJTX(JT,NTJ(JT))=DPX2
                     PJTY(JT,NTJ(JT))=DPY2
                     PJTZ(JT,NTJ(JT))=DPZ2
                     PJTE(JT,NTJ(JT))=DPE2
                     PJTM(JT,NTJ(JT))=0.0
                     PT(JT,3)=PT(JT,3)-DPZ2
                     PT(JT,4)=PT(JT,4)-DPE2
                  ENDIF
               ENDIF
  16           DPKC11=-(PP(JP,10)-PKC11)/2.0
               DPKC12=-(PP(JP,11)-PKC12)/2.0
               DPKC21=-(PT(JT,10)-PKC21)/2.0
               DPKC22=-(PT(JT,11)-PKC22)/2.0
               WP=EPP+ETP
               WM=EPM+ETM
               SW=WP*WM
            ENDIF
         ENDIF
 C               ********If jet is quenched the pt from valence quark
 C                       hard scattering has to reduced by d*kapa
 C
 C   
 10      PTP02=PP(JP,1)**2+PP(JP,2)**2
         PTT02=PT(JT,1)**2+PT(JT,2)**2
 C       
         AMQ=MAX(PP(JP,14)+PP(JP,15),PT(JT,14)+PT(JT,15))
         AMX=HIPR1(1)+AMQ
 C               ********consider mass cut-off for strings which
 C                       must also include quark's mass
         AMP0=AMX
         DPM0=AMX
         NFDP=0
         IF(NFP(JP,5).LE.2.AND.NFP(JP,3).NE.0) THEN
                 AMP0=ULMASS(NFP(JP,3))
                 NFDP=NFP(JP,3)+2*NFP(JP,3)/ABS(NFP(JP,3))
                 DPM0=ULMASS(NFDP)
                 IF(DPM0.LE.0.0) THEN
                         NFDP=NFDP-2*NFDP/ABS(NFDP)
                         DPM0=ULMASS(NFDP)
                 ENDIF
         ENDIF
         AMT0=AMX
         DTM0=AMX
         NFDT=0
         IF(NFT(JT,5).LE.2.AND.NFT(JT,3).NE.0) THEN
                 AMT0=ULMASS(NFT(JT,3))
                 NFDT=NFT(JT,3)+2*NFT(JT,3)/ABS(NFT(JT,3))
                 DTM0=ULMASS(NFDT)
                 IF(DTM0.LE.0.0) THEN
                         NFDT=NFDT-2*NFDT/ABS(NFDT)
                         DTM0=ULMASS(NFDT)
                 ENDIF
         ENDIF
 C       
         AMPN=SQRT(AMP0**2+PTP02)
         AMTN=SQRT(AMT0**2+PTT02)
         SNN=(AMPN+AMTN)**2+0.001
 C
         IF(SW.LT.SNN+0.001) GO TO 4000
 C               ********Scatter only if SW>SNN
 C*****give some PT kick to the two exited strings******************
 20      SWPTN=4.0*(MAX(AMP0,AMT0)**2+MAX(PTP02,PTT02))
         SWPTD=4.0*(MAX(DPM0,DTM0)**2+MAX(PTP02,PTT02))
         SWPTX=4.0*(AMX**2+MAX(PTP02,PTT02))
         IF(SW.LE.SWPTN) THEN
                 PKCMX=0.0
         ELSE IF(SW.GT.SWPTN .AND. SW.LE.SWPTD
      &          .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
            PKCMX=SQRT(SW/4.0-MAX(AMP0,AMT0)**2)
      &           -SQRT(MAX(PTP02,PTT02))
         ELSE IF(SW.GT.SWPTD .AND. SW.LE.SWPTX
      &          .AND.NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0) THEN
            PKCMX=SQRT(SW/4.0-MAX(DPM0,DTM0)**2)
      &           -SQRT(MAX(PTP02,PTT02))
         ELSE IF(SW.GT.SWPTX) THEN
            PKCMX=SQRT(SW/4.0-AMX**2)-SQRT(MAX(PTP02,PTT02))
         ENDIF
 C               ********maximun PT kick
 C*********************************************************
 C
         IF(NFP(JP,10).EQ.1.OR.NFT(JT,10).EQ.1) THEN
                 IF(PKC1.GT.PKCMX) THEN
                         PKC1=PKCMX
                         PKC11=PKC1*COS(PHI1)
                         PKC12=PKC1*SIN(PHI1)
                         DPKC11=-(PP(JP,10)-PKC11)/2.0
                         DPKC12=-(PP(JP,11)-PKC12)/2.0
                 ENDIF
                 IF(PKC2.GT.PKCMX) THEN
                         PKC2=PKCMX
                         PKC21=PKC2*COS(PHI2)
                         PKC22=PKC2*SIN(PHI2)
                         DPKC21=-(PT(JT,10)-PKC21)/2.0
                         DPKC22=-(PT(JT,11)-PKC22)/2.0
                 ENDIF
                 DPKC1=DPKC11+DPKC21
                 DPKC2=DPKC12+DPKC22
                 NFP(JP,10)=-NFP(JP,10)
                 NFT(JT,10)=-NFT(JT,10)
                 GO TO 40
         ENDIF
 C               ********If the valence quarks had a hard-collision
 C                       the pt kick is the pt from hard-collision.
         I_SNG=0
         IF(IHPR2(13).NE.0 .AND. RLU(0).LE.HIDAT(4)) I_SNG=1
         IF((NFP(JP,5).EQ.3 .OR.NFT(JT,5).EQ.3).OR.
      &          (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
      &          (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) I_SNG=0
 C
 C               ********decite whether to have single-diffractive
         IF(IHPR2(5).EQ.0) THEN
                 PKC=HIPR1(2)*SQRT(-ALOG(1.0-RLU(0)
      &                  *(1.0-EXP(-PKCMX**2/HIPR1(2)**2))))
                 GO TO 30
         ENDIF
         PKC=HIRND2(3,0.0,PKCMX**2)
         PKC=SQRT(PKC)
         IF(PKC.GT.HIPR1(20)) 
      &     PKC=HIPR1(2)*SQRT(-ALOG(EXP(-HIPR1(20)**2/HIPR1(2)**2)
      &         -RLU(0)*(EXP(-HIPR1(20)**2/HIPR1(2)**2)-
      &         EXP(-PKCMX**2/HIPR1(2)**2))))
 C
         IF(I_SNG.EQ.1) PKC=0.65*SQRT(
      &          -ALOG(1.0-RLU(0)*(1.0-EXP(-PKCMX**2/0.65**2))))
 C                       ********select PT kick
 30      PHI0=2.0*HIPR1(40)*RLU(0)
         PKC11=PKC*SIN(PHI0)
         PKC12=PKC*COS(PHI0)
         PKC21=-PKC11
         PKC22=-PKC12
         DPKC1=0.0
         DPKC2=0.0
 40      PP11=PP(JP,1)+PKC11-DPKC1
         PP12=PP(JP,2)+PKC12-DPKC2
         PT11=PT(JT,1)+PKC21-DPKC1
         PT12=PT(JT,2)+PKC22-DPKC2
         PTP2=PP11**2+PP12**2
         PTT2=PT11**2+PT12**2
 C
         AMPN=SQRT(AMP0**2+PTP2)
         AMTN=SQRT(AMT0**2+PTT2)
         SNN=(AMPN+AMTN)**2+0.001
 C***************************************
         WP=EPP+ETP
         WM=EPM+ETM
         SW=WP*WM
 C****************************************
         IF(SW.LT.SNN) THEN
            MISS=MISS+1
            IF(MISS.LE.100) then
               PKC=0.0
               GO TO 30
            ENDIF
            IF(IHPR2(10).NE.0) 
      &       WRITE(6,*) 'Error occured in Pt kick section of HIJSFT'
            GO TO 4000
         ENDIF
 C******************************************************************
         AMPD=SQRT(DPM0**2+PTP2)
         AMTD=SQRT(DTM0**2+PTT2)
 
         AMPX=SQRT(AMX**2+PTP2)
         AMTX=SQRT(AMX**2+PTT2)
 
         DPN=AMPN**2/SW
         DTN=AMTN**2/SW
         DPD=AMPD**2/SW
         DTD=AMTD**2/SW
         DPX=AMPX**2/SW
         DTX=AMTX**2/SW
 C
         SPNTD=(AMPN+AMTD)**2
         SPNTX=(AMPN+AMTX)**2
 C                       ********CM energy if proj=N,targ=N*
         SPDTN=(AMPD+AMTN)**2
         SPXTN=(AMPX+AMTN)**2
 C                       ********CM energy if proj=N*,targ=N
         SPDTX=(AMPD+AMTX)**2
         SPXTD=(AMPX+AMTD)**2
         SDD=(AMPD+AMTD)**2
         SXX=(AMPX+AMTX)**2
 
 C
 C       
 C               ********CM energy if proj=delta, targ=delta
 C****************There are many different cases**********
 c       IF(IHPR2(15).EQ.1) GO TO 500
 C
 C               ********to have DPM type soft interactions
 C
  45     CONTINUE
         IF(SW.GT.SXX+0.001) THEN
            IF(I_SNG.EQ.0) THEN
               D1=DPX
               D2=DTX
               NFP3=0
               NFT3=0
               GO TO 400
            ELSE
 c**** 5/30/1998 this is identical to the above statement. Added to
 c**** avoid questional branching to block.
               IF((NFP(JP,5).EQ.3 .AND.NFT(JT,5).EQ.3).OR.
      &           (NPJ(JP).NE.0.OR.NFP(JP,10).NE.0).OR.
      &           (NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
                  D1=DPX
                  D2=DTX
                  NFP3=0
                  NFT3=0
                  GO TO 400
               ENDIF
 C               ********do not allow excited strings to have 
 C                       single-diffr 
               IF(RLU(0).GT.0.5.OR.(NFT(JT,5).GT.2.OR.
      &                NTJ(JT).NE.0.OR.NFT(JT,10).NE.0)) THEN
                  D1=DPN
                  D2=DTX
                  NFP3=NFP(JP,3)
                  NFT3=0
                  GO TO 220
               ELSE
                  D1=DPX
                  D2=DTN
                  NFP3=0
                  NFT3=NFT(JT,3)
                  GO TO 240
               ENDIF
 C               ********have single diffractive collision
            ENDIF
         ELSE IF(SW.GT.MAX(SPDTX,SPXTD)+0.001 .AND.
      &                          SW.LE.SXX+0.001) THEN
            IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0.AND.
      &         RLU(0).GT.0.5).OR.(NPJ(JP).EQ.0
      &         .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
               D1=DPD
               D2=DTX
               NFP3=NFDP
               NFT3=0
               GO TO 220
            ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
               D1=DPX
               D2=DTD
               NFP3=0
               NFT3=NFDT
               GO TO 240
            ENDIF
            GO TO 4000
         ELSE IF(SW.GT.MIN(SPDTX,SPXTD)+0.001.AND.
      &                  SW.LE.MAX(SPDTX,SPXTD)+0.001) THEN
            IF(SPDTX.LE.SPXTD.AND.NPJ(JP).EQ.0
      &                       .AND.NFP(JP,5).LE.2) THEN
               D1=DPD
               D2=DTX
               NFP3=NFDP
               NFT3=0
               GO TO 220
            ELSE IF(SPDTX.GT.SPXTD.AND.NTJ(JT).EQ.0
      &                       .AND.NFT(JT,5).LE.2) THEN
               D1=DPX
               D2=DTD
               NFP3=0
               NFT3=NFDT
               GO TO 240
            ENDIF
 c*** 5/30/1998 added to avoid questional branching to another block
 c*** this is identical to the statement following the next ELSE IF
            IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
      &       .AND.RLU(0).GT.0.5).OR.(NPJ(JP).EQ.0
      &        .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
               D1=DPN
               D2=DTX
               NFP3=NFP(JP,3)
               NFT3=0
               GO TO 220
            ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
               D1=DPX
               D2=DTN
               NFP3=0
               NFT3=NFT(JT,3)
               GO TO 240
            ENDIF
            GO TO 4000
         ELSE IF(SW.GT.MAX(SPNTX,SPXTN)+0.001 .AND.
      &                  SW.LE.MIN(SPDTX,SPXTD)+0.001) THEN
            IF(((NPJ(JP).EQ.0.AND.NTJ(JT).EQ.0
      &       .AND.RLU(0).GT.0.5).OR.(NPJ(JP).EQ.0
      &        .AND.NTJ(JT).NE.0)).AND.NFP(JP,5).LE.2) THEN
               D1=DPN
               D2=DTX
               NFP3=NFP(JP,3)
               NFT3=0
               GO TO 220
            ELSE IF(NTJ(JT).EQ.0.AND.NFT(JT,5).LE.2) THEN
               D1=DPX
               D2=DTN
               NFP3=0
               NFT3=NFT(JT,3)
               GO TO 240
            ENDIF
            GO TO 4000
         ELSE IF(SW.GT.MIN(SPNTX,SPXTN)+0.001 .AND.
      &                  SW.LE.MAX(SPNTX,SPXTN)+0.001) THEN
            IF(SPNTX.LE.SPXTN.AND.NPJ(JP).EQ.0
      &                           .AND.NFP(JP,5).LE.2) THEN
               D1=DPN
               D2=DTX
               NFP3=NFP(JP,3)
               NFT3=0
               GO TO 220
            ELSEIF(SPNTX.GT.SPXTN.AND.NTJ(JT).EQ.0
      &                           .AND.NFT(JT,5).LE.2) THEN
               D1=DPX
               D2=DTN
               NFP3=0
               NFT3=NFT(JT,3)
               GO TO 240
            ENDIF
            GO TO 4000
         ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
      &                  (NPJ(JP).NE.0 .OR.NTJ(JT).NE.0)) THEN
            GO TO 4000
         ELSE IF(SW.LE.MIN(SPNTX,SPXTN)+0.001 .AND.
      &          NFP(JP,5).GT.2.AND.NFT(JT,5).GT.2) THEN
            GO TO 4000
         ELSE IF(SW.GT.SDD+0.001.AND.SW.LE.
      &                     MIN(SPNTX,SPXTN)+0.001) THEN
            D1=DPD
            D2=DTD
            NFP3=NFDP
            NFT3=NFDT
            GO TO 100
         ELSE IF(SW.GT.MAX(SPNTD,SPDTN)+0.001 
      &                      .AND. SW.LE.SDD+0.001) THEN
            IF(RLU(0).GT.0.5) THEN
               D1=DPD
               D2=DTN
               NFP3=NFDP
               NFT3=NFT(JT,3)
               GO TO 100
            ELSE
               D1=DPN
               D2=DTD
               NFP3=NFP(JP,3)
               NFT3=NFDT
               GO TO 100
            ENDIF
         ELSE IF(SW.GT.MIN(SPNTD,SPDTN)+0.001
      &          .AND. SW.LE.MAX(SPNTD,SPDTN)+0.001) THEN
            IF(SPNTD.GT.SPDTN) THEN
               D1=DPD
               D2=DTN
               NFP3=NFDP
               NFT3=NFT(JT,3)
               GO TO 100
            ELSE
               D1=DPN
               D2=DTD
               NFP3=NFP(JP,3)
               NFT3=NFDT
               GO TO 100
            ENDIF
         ELSE IF(SW.LE.MIN(SPNTD,SPDTN)+0.001) THEN
            D1=DPN
            D2=DTN
            NFP3=NFP(JP,3)
            NFT3=NFT(JT,3)
            GO TO 100
         ENDIF
         WRITE(6,*) ' Error in HIJSFT: There is no path to here'
         RETURN
 C
 C***************  elastic scattering ***************
 C       this is like elastic, both proj and targ mass
 C       must be fixed
 C***************************************************
 100     NFP5=MAX(2,NFP(JP,5))
         NFT5=MAX(2,NFT(JT,5))
         BB1=1.0+D1-D2
         BB2=1.0+D2-D1
         IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
                 MISS=MISS+1
                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
                 PKC=PKC*0.5
                 GO TO 30
         ENDIF
         IF(RLU(0).LT.0.5) THEN
                 X1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
                 X2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
         ELSE
                 X1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
                 X2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
         ENDIF
         IHNT2(13)=2
         GO TO 600
 C
 C********** Single diffractive ***********************
 C either proj or targ's mass is fixed
 C*****************************************************
 220     NFP5=MAX(2,NFP(JP,5))
         NFT5=3
         IF(NFP3.EQ.0) NFP5=3
         BB2=1.0+D2-D1
         IF(BB2**2.LT.4.0*D2) THEN
                 MISS=MISS+1
                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
                 PKC=PKC*0.5
                 GO TO 30
         ENDIF
         XMIN=(BB2-SQRT(BB2**2-4.0*D2))/2.0
         XMAX=(BB2+SQRT(BB2**2-4.0*D2))/2.0
         MISS4=0
 222     X2=HIRND2(6,XMIN,XMAX)
         X1=D1/(1.0-X2)
         IF(X2*(1.0-X1).LT.(D2+1.E-4/SW)) THEN
                 MISS4=MISS4+1
                 IF(MISS4.LE.1000) GO TO 222
                 GO TO 5000
         ENDIF
         IHNT2(13)=2
         GO TO 600
 C                       ********Fix proj mass*********
 240     NFP5=3
         NFT5=MAX(2,NFT(JT,5))
         IF(NFT3.EQ.0) NFT5=3
         BB1=1.0+D1-D2
         IF(BB1**2.LT.4.0*D1) THEN
                 MISS=MISS+1
                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
                 PKC=PKC*0.5
                 GO TO 30
         ENDIF
         XMIN=(BB1-SQRT(BB1**2-4.0*D1))/2.0
         XMAX=(BB1+SQRT(BB1**2-4.0*D1))/2.0
         MISS4=0
 242     X1=HIRND2(6,XMIN,XMAX)
         X2=D2/(1.0-X1)
         IF(X1*(1.0-X2).LT.(D1+1.E-4/SW)) THEN
                 MISS4=MISS4+1
                 IF(MISS4.LE.1000) GO TO 242
                 GO TO 5000
         ENDIF
         IHNT2(13)=2
         GO TO 600
 C                       ********Fix targ mass*********
 C
 C*************non-single diffractive**********************
 C       both proj and targ may not be fixed in mass 
 C*********************************************************
 C
 400     NFP5=3
         NFT5=3
         BB1=1.0+D1-D2
         BB2=1.0+D2-D1
         IF(BB1**2.LT.4.0*D1 .OR. BB2**2.LT.4.0*D2) THEN
                 MISS=MISS+1
                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 3000
                 PKC=PKC*0.5
                 GO TO 30
         ENDIF
         XMIN1=(BB1-SQRT(BB1**2-4.0*D1))/2.0
         XMAX1=(BB1+SQRT(BB1**2-4.0*D1))/2.0
         XMIN2=(BB2-SQRT(BB2**2-4.0*D2))/2.0
         XMAX2=(BB2+SQRT(BB2**2-4.0*D2))/2.0
         MISS4=0 
 410     X1=HIRND2(4,XMIN1,XMAX1)
         X2=HIRND2(4,XMIN2,XMAX2)
         IF(NFP(JP,5).EQ.3.OR.NFT(JT,5).EQ.3) THEN
                 X1=HIRND2(6,XMIN1,XMAX1)
                 X2=HIRND2(6,XMIN2,XMAX2)
         ENDIF
 C                       ********
         IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
      &                  ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
                 X1=HIRND2(5,XMIN1,XMAX1)
         ENDIF
         IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
      &                  ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
                 X2=HIRND2(5,XMIN2,XMAX2)
         ENDIF
 c       IF(IOPMAIN.EQ.3) X1=HIRND2(6,XMIN1,XMAX1)
 c       IF(IOPMAIN.EQ.2) X2=HIRND2(6,XMIN2,XMAX2) 
 C       ********For q-qbar or (qq)-(qq)bar system use symetric
 C               distribution, for q-(qq) or qbar-(qq)bar use
 C               unsymetrical distribution
 C
         IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000) X1=1.0-X1
         XXP=X1*(1.0-X2)
         XXT=X2*(1.0-X1)
         IF(XXP.LT.(D1+1.E-4/SW) .OR. XXT.LT.(D2+1.E-4/SW)) THEN
                 MISS4=MISS4+1
                 IF(MISS4.LE.1000) GO TO 410
                 GO TO 5000
         ENDIF
         IHNT2(13)=3
 C***************************************************
 C***************************************************
 600     CONTINUE
         IF(X1*(1.0-X2).LT.(AMPN**2-1.E-4)/SW.OR.
      &                  X2*(1.0-X1).LT.(AMTN**2-1.E-4)/SW) THEN
                 MISS=MISS+1
                 IF(MISS.GT.100.OR.PKC.EQ.0.0) GO TO 2000
                 PKC=0.0
                 GO TO 30
         ENDIF
 C
         EPP=(1.0-X2)*WP
         EPM=X1*WM
         ETP=X2*WP
         ETM=(1.0-X1)*WM
         PP(JP,3)=(EPP-EPM)/2.0
         PP(JP,4)=(EPP+EPM)/2.0
         IF(EPP*EPM-PTP2.LT.0.0) GO TO 6000
         PP(JP,5)=SQRT(EPP*EPM-PTP2)
         NFP(JP,3)=NFP3
         NFP(JP,5)=NFP5
 
         PT(JT,3)=(ETP-ETM)/2.0
         PT(JT,4)=(ETP+ETM)/2.0
         IF(ETP*ETM-PTT2.LT.0.0) GO TO 6000
         PT(JT,5)=SQRT(ETP*ETM-PTT2)
         NFT(JT,3)=NFT3
         NFT(JT,5)=NFT5
 C*****recoil PT from hard-inter is shared by two end-partons 
 C       so that pt=p1+p2
         PP(JP,1)=PP11-PKC11
         PP(JP,2)=PP12-PKC12
 
         KICKDIP=1
         KICKDIT=1
         IF(ABS(NFP(JP,1)*NFP(JP,2)).GT.1000000.OR.
      &                  ABS(NFP(JP,1)*NFP(JP,2)).LT.100) THEN
                 KICKDIP=0
         ENDIF
         IF(ABS(NFT(JT,1)*NFT(JT,2)).GT.1000000.OR.
      &                  ABS(NFT(JT,1)*NFT(JT,2)).LT.100) THEN
                 KICKDIT=0
         ENDIF
         IF((KICKDIP.EQ.0.AND.RLU(0).LT.0.5)
      &     .OR.(KICKDIP.NE.0.AND.RLU(0)
      &     .LT.0.5/(1.0+(PKC11**2+PKC12**2)/HIPR1(22)**2))) THEN
            PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,6)
            PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,7)
            PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
      &              +PP(JP,8)+PKC11
            PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
      &              +PP(JP,9)+PKC12
         ELSE
            PP(JP,8)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0+PP(JP,8)
            PP(JP,9)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0+PP(JP,9)
            PP(JP,6)=(PP(JP,1)-PP(JP,6)-PP(JP,8)-DPKC1)/2.0
      &              +PP(JP,6)+PKC11
            PP(JP,7)=(PP(JP,2)-PP(JP,7)-PP(JP,9)-DPKC2)/2.0
      &              +PP(JP,7)+PKC12
         ENDIF
         PP(JP,1)=PP(JP,6)+PP(JP,8)
         PP(JP,2)=PP(JP,7)+PP(JP,9)
 C                               ********pt kick for proj
         PT(JT,1)=PT11-PKC21
         PT(JT,2)=PT12-PKC22
         IF((KICKDIT.EQ.0.AND.RLU(0).LT.0.5)
      &     .OR.(KICKDIT.NE.0.AND.RLU(0)
      &     .LT.0.5/(1.0+(PKC21**2+PKC22**2)/HIPR1(22)**2))) THEN
            PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,6)
            PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,7)
            PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
      &              +PT(JT,8)+PKC21
            PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
      &              +PT(JT,9)+PKC22
         ELSE
            PT(JT,8)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0+PT(JT,8)
            PT(JT,9)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0+PT(JT,9)
            PT(JT,6)=(PT(JT,1)-PT(JT,6)-PT(JT,8)-DPKC1)/2.0
      &              +PT(JT,6)+PKC21
            PT(JT,7)=(PT(JT,2)-PT(JT,7)-PT(JT,9)-DPKC2)/2.0
      &              +PT(JT,7)+PKC22
         ENDIF
         PT(JT,1)=PT(JT,6)+PT(JT,8)
         PT(JT,2)=PT(JT,7)+PT(JT,9)
 C                       ********pt kick for targ
 
         IF(NPJ(JP).NE.0) NFP(JP,5)=3
         IF(NTJ(JT).NE.0) NFT(JT,5)=3
 C                       ********jets must be connected to string
         IF(EPP/(EPM+0.0001).LT.ETP/(ETM+0.0001).AND.
      &                  ABS(NFP(JP,1)*NFP(JP,2)).LT.1000000)THEN
                 DO 620 JSB=1,15
                 PSB=PP(JP,JSB)
                 PP(JP,JSB)=PT(JT,JSB)
                 PT(JT,JSB)=PSB
                 NSB=NFP(JP,JSB)
                 NFP(JP,JSB)=NFT(JT,JSB)
                 NFT(JT,JSB)=NSB
 620             CONTINUE
 C               ********when Ycm(JP)<Ycm(JT) after the collision
 C                       exchange the positions of the two   
         ENDIF
 C
         RETURN
 C**************************************************
 C**************************************************
 1000    IERROR=1
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) '    Fatal HIJSFT start error,abandon this event'
         WRITE(6,*) '    PROJ E+,E-,W+',EPP,EPM,WP
         WRITE(6,*) '    TARG E+,E-,W-',ETP,ETM,WM
         WRITE(6,*) '    W+*W-, (APN+ATN)^2',SW,SNN
         RETURN
 2000    IERROR=0
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) '    (2)energy partition fail,'
         WRITE(6,*) '    HIJSFT not performed, but continue'
         WRITE(6,*) '    MP1,MPN',X1*(1.0-X2)*SW,AMPN**2
         WRITE(6,*) '    MT2,MTN',X2*(1.0-X1)*SW,AMTN**2
         RETURN
 3000    IERROR=0
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) '    (3)something is wrong with the pt kick, '
         WRITE(6,*) '    HIJSFT not performed, but continue'
         WRITE(6,*) '    D1=',D1,' D2=',D2,' SW=',SW
         WRITE(6,*) '    HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
         WRITE(6,*) '    THIS COLLISON NFP5=',NFP5, ' NFT5=',NFT5
         WRITE(6,*) '    # OF JET IN PROJ',NPJ(JP),' IN TARG',NTJ(JT)
         RETURN
 4000    IERROR=0
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) '    (4)unable to choose process, but not harmful'
         WRITE(6,*) '    HIJSFT not performed, but continue'
         WRITE(6,*) '    PTP=',SQRT(PTP2),' PTT=',SQRT(PTT2),' SW=',SW
         WRITE(6,*) '    AMCUT=',AMX,' JP=',JP,' JT=',JT
         WRITE(6,*) '    HISTORY NFP5=',NFP(JP,5),' NFT5=',NFT(JT,5)
         RETURN
 5000    IERROR=0
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) '    energy partition failed(5),for limited try'
         WRITE(6,*) '    HIJSFT not performed, but continue'
         WRITE(6,*) '    NFP5=',NFP5,' NFT5=',NFT5
         WRITE(6,*) '    D1',D1,' X1(1-X2)',X1*(1.0-X2)
         WRITE(6,*) '    D2',D2,' X2(1-X1)',X2*(1.0-X1)
         RETURN
 6000    PKC=0.0
         MISS=MISS+1
         IF(MISS.LT.100) GO TO 30
         IERROR=1
         IF(IHPR2(10).EQ.0) RETURN
         WRITE(6,*) ' ERROR OCCURED, HIJSFT NOT PERFORMED'
         WRITE(6,*) ' Abort this event'
         WRITE(6,*) 'MTP,PTP2',EPP*EPM,PTP2,'  MTT,PTT2',ETP*ETM,PTT2 
         RETURN
         END
 C
 C       
 C***************************************
         SUBROUTINE HIJFLV(ID)
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         ID=1
         RNID=RLU(0)
         IF(RNID.GT.0.43478) THEN
                 ID=2
                 IF(RNID.GT.0.86956) ID=3
         ENDIF
         RETURN
         END
 C
 C
 C
         SUBROUTINE HIPTDI(PT,PTMAX,IOPT)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         IF(IOPT.EQ.2) THEN
                 PT=HIRND2(7,0.0,PTMAX)
                 IF(PT.GT.HIPR1(8)) 
      &          PT=HIPR1(2)*SQRT(-ALOG(EXP(-HIPR1(8)**2/HIPR1(2)**2)
      &                  -RLU(0)*(EXP(-HIPR1(8)**2/HIPR1(2)**2)-
      &                  EXP(-PTMAX**2/HIPR1(2)**2))))
 
         ELSE
                 PT=HIPR1(2)*SQRT(-ALOG(1.0-RLU(0)*
      &                  (1.0-EXP(-PTMAX**2/HIPR1(2)**2))))
         ENDIF
         PTMAX0=MAX(PTMAX,0.01)
         PT=MIN(PTMAX0-0.01,PT)
         RETURN
         END
 C*************************
 C
 C
 C
 C
 C ********************************************************
 C ************************              WOOD-SAX
         SUBROUTINE HIJWDS(IA,IDH,XHIGH)
 C     SETS UP HISTOGRAM IDH WITH RADII FOR
 C     NUCLEUS IA DISTRIBUTED ACCORDING TO THREE PARAM WOOD SAXON
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/WOOD/R,D,FNORM,W
         SAVE  /WOOD/
         DIMENSION IAA(20),RR(20),DD(20),WW(20),RMS(20)
         EXTERNAL RWDSAX,WDSAX
 C
 C   PARAMETERS OF SPECIAL NUCLEI FROM ATOMIC DATA AND NUC DATA TABLES
 C     VOL 14, 5-6 1974
         DATA IAA/2,4,12,16,27,32,40,56,63,93,184,197,208,7*0./
         DATA RR/0.01,.964,2.355,2.608,2.84,3.458,3.766,3.971,4.214,
      1        4.87,6.51,6.38,6.624,7*0./
         DATA DD/0.5882,.322,.522,.513,.569,.61,.586,.5935,.586,.573,
      1        .535,.535,.549,7*0./
         DATA WW/0.0,.517,-0.149,-0.051,0.,-0.208,-0.161,13*0./
         DATA RMS/2.11,1.71,2.46,2.73,3.05,3.247,3.482,3.737,3.925,4.31,
      1        5.42,5.33,5.521,7*0./
       
         SAVE IAA, RR, DD, WW, RMS
 C
         A=IA
 C
 C               ********SET WOOD-SAX PARAMS FIRST  AS IN DATE ET AL
         D=0.54
 C                       ********D IS WOOD SAX DIFFUSE PARAM IN FM
         R=1.19*A**(1./3.) - 1.61*A**(-1./3.)
 C                       ********R IS RADIUS PARAM
         W=0.
 C               ********W IS The third of three WOOD-SAX PARAM
 C
 C               ********CHECK TABLE FOR SPECIAL CASES
         DO 10 I=1,13
                 IF (IA.EQ.IAA(I)) THEN
                         R=RR(I)
                         D=DD(I)
                         W=WW(I)
                         RS=RMS(I)
                 END IF
 10      CONTINUE
 C                       ********FNORM is the normalize factor
         FNORM=1.0
         XLOW=0.
         XHIGH=R+ 12.*D
         IF (W.LT.-0.01)  THEN
                 IF (XHIGH.GT.R/SQRT(ABS(W))) XHIGH=R/SQRT(ABS(W))
         END IF
         FGAUS=GAUSS1(RWDSAX,XLOW,XHIGH,0.001)
         FNORM=1./FGAUS
 C
         IF (IDH.EQ.1) THEN
            HINT1(72)=R
            HINT1(73)=D
            HINT1(74)=W
            HINT1(75)=FNORM/4.0/HIPR1(40)
         ELSE IF (IDH.EQ.2) THEN
            HINT1(76)=R
            HINT1(77)=D
            HINT1(78)=W
            HINT1(79)=FNORM/4.0/HIPR1(40)
         ENDIF
 C
 C       NOW SET UP HBOOK FUNCTIONS IDH FOR  R**2*RHO(R)
 C       THESE HISTOGRAMS ARE USED TO GENERATE RANDOM RADII
         CALL HIFUN(IDH,XLOW,XHIGH,RWDSAX)
         RETURN
         END
 C
 C
         FUNCTION WDSAX(X)
 C                       ********THREE PARAMETER WOOD SAXON
         COMMON/WOOD/R,D,FNORM,W
         SAVE  /WOOD/
         WDSAX=FNORM*(1.+W*(X/R)**2)/(1+EXP((X-R)/D))
         IF (W.LT.0.) THEN
                 IF (X.GE.R/SQRT(ABS(W))) WDSAX=0.
         ENDIF
         RETURN
         END
 C
 C
         FUNCTION RWDSAX(X)
         RWDSAX=X*X*WDSAX(X)
         RETURN
         END
 C
 C
 C
 C
         FUNCTION WDSAX1(X)
 C                       ********THREE PARAMETER WOOD SAXON 
 C                               FOR  PROJECTILE
 C       HINT1(72)=R, HINT1(73)=D, HINT1(74)=W, HINT1(75)=FNORM
 C
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         WDSAX1=HINT1(75)*(1.+HINT1(74)*(X/HINT1(72))**2)/
      &       (1+EXP((X-HINT1(72))/HINT1(73)))
         IF (HINT1(74).LT.0.) THEN
                 IF (X.GE.HINT1(72)/SQRT(ABS(HINT1(74)))) WDSAX1=0.
         ENDIF
         RETURN
         END
 C
 C
         FUNCTION WDSAX2(X)
 C                       ********THREE PARAMETER WOOD SAXON 
 C                               FOR  TARGET
 C       HINT1(76)=R,HINT1(77)=D, HINT1(78)=W, HINT1(79)=FNORM
 C
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         WDSAX2=HINT1(79)*(1.+HINT1(78)*(X/HINT1(76))**2)/
      &       (1+EXP((X-HINT1(76))/HINT1(77)))
         IF (HINT1(78).LT.0.) THEN
                 IF (X.GE.HINT1(76)/SQRT(ABS(HINT1(78)))) WDSAX2=0.
         ENDIF
         RETURN
         END
 C
 C
 C       THIS FUNCTION IS TO CALCULATE THE NUCLEAR PROFILE FUNCTION
 C       OF THE  COLLIDERING SYSTEM (IN UNITS OF 1/mb)
 C
         FUNCTION  PROFILE(XB)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         EXTERNAL FLAP, FLAP1, FLAP2
 C
         BB=XB
         PROFILE=1.0
         IF(IHNT2(1).GT.1 .AND. IHNT2(3).GT.1) THEN
            PROFILE=float(IHNT2(1))*float(IHNT2(3))*0.1*
      &          GAUSS1(FLAP,0.0,HIPR1(34),0.01)
         ELSE IF(IHNT2(1).EQ.1 .AND. IHNT2(3).GT.1) THEN
            PROFILE=0.2*float(IHNT2(3))*
      &          GAUSS1(FLAP2,0.0,HIPR1(35),0.001)
         ELSE IF(IHNT2(1).GT.1 .AND. IHNT2(3).EQ.1) THEN
            PROFILE=0.2*float(IHNT2(1))*
      &          GAUSS1(FLAP1,0.0,HIPR1(34),0.001)
         ENDIF
         RETURN
         END
 C
 C
         FUNCTION FLAP(X)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         EXTERNAL FGP1
         B1=X
         FLAP=GAUSS2(FGP1,0.0,2.0*HIPR1(40),0.01)
         RETURN
         END
 C
         FUNCTION FGP1(X)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         EXTERNAL FGP2
         PHI=X
         FGP1=2.0*GAUSS3(FGP2,0.0,HIPR1(34),0.01)
         RETURN
         END
 C
         FUNCTION FGP2(X)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         EXTERNAL FGP3
         Z1=X
         FGP2=2.0*GAUSS4(FGP3,0.0,HIPR1(35),0.01)
         RETURN
         END
 C
         FUNCTION FGP3(X)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         R1=SQRT(B1**2+Z1**2)
         R2=SQRT(BB**2+B1**2-2.0*B1*BB*COS(PHI)+X**2)
         FGP3=B1*WDSAX1(R1)*WDSAX2(R2)
         RETURN
         END
 C
 C
         FUNCTION FLAP1(X)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         R=SQRT(BB**2+X**2)
         FLAP1=WDSAX1(R)
         RETURN
         END
 C
 C
         FUNCTION FLAP2(X)
         COMMON/PACT/BB,B1,PHI,Z1
         SAVE  /PACT/
         R=SQRT(BB**2+X**2)
         FLAP2=WDSAX2(R)
         RETURN
         END
 C
 C
 C The next three subroutines are for Monte Carlo generation 
 C according to a given function FHB. One calls first HIFUN 
 C with assigned channel number I, low and up limits. Then to 
 C generate the distribution one can call HIRND(I) which gives 
 C you a random number generated according to the given function.
 C 
         SUBROUTINE HIFUN(I,XMIN,XMAX,FHB)
         COMMON/HIJHB/RR(10,201),XX(10,201)
         SAVE  /HIJHB/
         EXTERNAL FHB
         FNORM=GAUSS1(FHB,XMIN,XMAX,0.001)
         DO 100 J=1,201
                 XX(I,J)=XMIN+(XMAX-XMIN)*(J-1)/200.0
                 XDD=XX(I,J)
                 RR(I,J)=GAUSS1(FHB,XMIN,XDD,0.001)/FNORM
 100     CONTINUE
         RETURN
         END
 C
 C
 C
         FUNCTION HIRND(I)
         COMMON/HIJHB/RR(10,201),XX(10,201)
         SAVE  /HIJHB/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         RX=RLU(0)
         JL=0
         JU=202
 10      IF(JU-JL.GT.1) THEN
            JM=(JU+JL)/2
            IF((RR(I,201).GT.RR(I,1)).EQV.(RX.GT.RR(I,JM))) THEN
               JL=JM
            ELSE
               JU=JM
            ENDIF
         GO TO 10
         ENDIF
         J=JL
         IF(J.LT.1) J=1
         IF(J.GE.201) J=200
         HIRND=(XX(I,J)+XX(I,J+1))/2.0
         RETURN
         END     
 C
 C
 C
 C
 C       This generate random number between XMIN and XMAX
         FUNCTION HIRND2(I,XMIN,XMAX)
         COMMON/HIJHB/RR(10,201),XX(10,201)
         SAVE  /HIJHB/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         IF(XMIN.LT.XX(I,1)) XMIN=XX(I,1)
         IF(XMAX.GT.XX(I,201)) XMAX=XX(I,201)
         JMIN=1+200*(XMIN-XX(I,1))/(XX(I,201)-XX(I,1))
         JMAX=1+200*(XMAX-XX(I,1))/(XX(I,201)-XX(I,1))
         RX=RR(I,JMIN)+(RR(I,JMAX)-RR(I,JMIN))*RLU(0)
         JL=0
         JU=202
 10      IF(JU-JL.GT.1) THEN
            JM=(JU+JL)/2
            IF((RR(I,201).GT.RR(I,1)).EQV.(RX.GT.RR(I,JM))) THEN
               JL=JM
            ELSE
               JU=JM
            ENDIF
         GO TO 10
         ENDIF
         J=JL
         IF(J.LT.1) J=1
         IF(J.GE.201) J=200
         HIRND2=(XX(I,J)+XX(I,J+1))/2.0
         RETURN
         END     
 C
 C
 C
 C
         SUBROUTINE HIJCRS
 C       THIS IS TO CALCULATE THE CROSS SECTIONS OF JET PRODUCTION AND
 C       THE TOTAL INELASTIC CROSS SECTIONS.
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/NJET/N,IP_CRS
         SAVE  /NJET/
         EXTERNAL FHIN,FTOT,FNJET,FTOTJET,FTOTRIG
         IF(HINT1(1).GE.10.0) CALL CRSJET
 C                       ********calculate jet cross section(in mb)
 C
         APHX1=HIPR1(6)*(IHNT2(1)**0.3333333-1.0)
         APHX2=HIPR1(6)*(IHNT2(3)**0.3333333-1.0)
         HINT1(11)=HINT1(14)-APHX1*HINT1(15)
      &                  -APHX2*HINT1(16)+APHX1*APHX2*HINT1(17)
         HINT1(10)=GAUSS1(FTOTJET,0.0,20.0,0.01)
         HINT1(12)=GAUSS1(FHIN,0.0,20.0,0.01)
         HINT1(13)=GAUSS1(FTOT,0.0,20.0,0.01)
         HINT1(60)=HINT1(61)-APHX1*HINT1(62)
      &                  -APHX2*HINT1(63)+APHX1*APHX2*HINT1(64)
         HINT1(59)=GAUSS1(FTOTRIG,0.0,20.0,0.01)
         IF(HINT1(59).EQ.0.0) HINT1(59)=HINT1(60)
         IF(HINT1(1).GE.10.0) Then
            DO 20 I=0,20
               N=I
               HINT1(80+I)=GAUSS1(FNJET,0.0,20.0,0.01)/HINT1(12)
  20        CONTINUE
         ENDIF
         HINT1(10)=HINT1(10)*HIPR1(31)
         HINT1(12)=HINT1(12)*HIPR1(31)
         HINT1(13)=HINT1(13)*HIPR1(31)
         HINT1(59)=HINT1(59)*HIPR1(31)
 C               ********Total and Inel cross section are calculated
 C                       by Gaussian integration.
         IF(IHPR2(13).NE.0) THEN
         HIPR1(33)=1.36*(1.0+36.0/HINT1(1)**2)
      &             *ALOG(0.6+0.1*HINT1(1)**2)
         HIPR1(33)=HIPR1(33)/HINT1(12)
         ENDIF
 C               ********Parametrized cross section for single
 C                       diffractive reaction(Goulianos)
         RETURN
         END
 C
 C
 C
 C
         FUNCTION FTOT(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         OMG=OMG0(X)*(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
         FTOT=2.0*(1.0-EXP(-OMG))
         RETURN
         END
 C
 C
 C
         FUNCTION FHIN(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         OMG=OMG0(X)*(HIPR1(30)+HINT1(11))/HIPR1(31)/2.0
         FHIN=1.0-EXP(-2.0*OMG)
         RETURN
         END
 C
 C
 C
         FUNCTION FTOTJET(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         OMG=OMG0(X)*HINT1(11)/HIPR1(31)/2.0
         FTOTJET=1.0-EXP(-2.0*OMG)
         RETURN
         END
 C
 C
 C
         FUNCTION FTOTRIG(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         OMG=OMG0(X)*HINT1(60)/HIPR1(31)/2.0
         FTOTRIG=1.0-EXP(-2.0*OMG)
         RETURN
         END
 C
 C
 C
 C
         FUNCTION FNJET(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/NJET/N,IP_CRS
         SAVE  /NJET/
         OMG1=OMG0(X)*HINT1(11)/HIPR1(31)
         C0=EXP(N*ALOG(OMG1)-SGMIN(N+1))
         IF(N.EQ.0) C0=1.0-EXP(-2.0*OMG0(X)*HIPR1(30)/HIPR1(31)/2.0)
         FNJET=C0*EXP(-OMG1)
         RETURN
         END
 C
 C
 C
 C
 C
         FUNCTION SGMIN(N)
         GA=0.
         IF(N.LE.2) GO TO 20
         DO 10 I=1,N-1
         Z=I
         GA=GA+ALOG(Z)
 10      CONTINUE
 20      SGMIN=GA
         RETURN
         END
 C
 C
 C
         FUNCTION OMG0(X)
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON /BESEL/X4
         SAVE   /BESEL/
         EXTERNAL BK
         X4=HIPR1(32)*SQRT(X)
         OMG0=HIPR1(32)**2*GAUSS2(BK,X4,X4+20.0,0.01)/96.0
 c        write(*,*)'OMG0 start',x,x4,OMG0
         RETURN
         END
 C
 C
 C
         FUNCTION ROMG(X)
 C               ********This gives the eikonal function from a table
 C                       calculated in the first call
         DIMENSION FR(0:1000)
 c       --------------Khaled
 c       we need to save FR because it is calculated only once
         COMMON/EIKONAL/FR             
 c       --------------->end Khaled
 
         DATA I0/0/
         save I0    !Uzhi
         IF(I0.NE.0) GO TO 100
         DO 50 I=1,1001
         XR=(I-1)*0.01
         FR(I-1)=OMG0(XR)
     
 50      CONTINUE
 100     I0=1                             
 
         IF(X.GE.10.0) THEN
                 ROMG=0.0
                 RETURN
         ENDIF
         IX=INT(X*100)
 
         ROMG=(FR(IX)*((IX+1)*0.01-X)+FR(IX+1)*(X-IX*0.01))/0.01
 
         RETURN
         END
 C
 C
 C
 
         FUNCTION BK(X)
         COMMON /BESEL/X4
         SAVE   /BESEL/
         BK=EXP(-X)*(X**2-X4**2)**2.50/15.0
         RETURN
         END
 C
 C
 C       THIS PROGRAM IS TO CALCULATE THE JET CROSS SECTION
 C       THE INTEGRATION IS DONE BY USING VEGAS
 C
         SUBROUTINE CRSJET
         IMPLICIT REAL*8(A-H,O-Z)
         REAL HIPR1(100),HINT1(100)
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/NJET/N,IP_CRS
         SAVE  /NJET/
         COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
         SAVE  /BVEG1/
         COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
         SAVE  /BVEG2/
         COMMON/BVEG3/F,TI,TSI
         SAVE  /BVEG3/
         COMMON/SEEDVAX/NUM1
         SAVE  /SEEDVAX/
         EXTERNAL FJET,FJETRIG
 C
 c************************
 c       NCALL give the number of inner-iteration, ITMX 
 C       gives the limit of out-iteration. Nprn is an option
 C       ( 1: print the integration process. 0: do not print)
 C
         NDIM=3
         IP_CRS=0
         CALL VEGAS(FJET,AVGI,SD,CHI2A)
         HINT1(14)=AVGI/2.5682
         IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
                 IP_CRS=1
                 CALL VEGAS(FJET,AVGI,SD,CHI2A)
                 HINT1(15)=AVGI/2.5682
         ENDIF
         IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
                 IP_CRS=2
                 CALL VEGAS(FJET,AVGI,SD,CHI2A)
                 HINT1(16)=AVGI/2.5682
         ENDIF
         IF(IHPR2(6).EQ.1.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
                 IP_CRS=3
                 CALL VEGAS(FJET,AVGI,SD,CHI2A)
                 HINT1(17)=AVGI/2.5682
         ENDIF
 C               ********Total inclusive jet cross section(Pt>P0) 
 C
         IF(IHPR2(3).NE.0) THEN
            IP_CRS=0
            CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
            HINT1(61)=AVGI/2.5682
            IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
               IP_CRS=1
               CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
               HINT1(62)=AVGI/2.5682
            ENDIF
            IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
               IP_CRS=2
               CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
               HINT1(63)=AVGI/2.5682
            ENDIF
            IF(IHPR2(6).EQ.1.AND.IHNT2(1).GT.1.AND.IHNT2(3).GT.1) THEN
               IP_CRS=3
               CALL VEGAS(FJETRIG,AVGI,SD,CHI2A)
               HINT1(64)=AVGI/2.5682
            ENDIF
         ENDIF
 C                       ********cross section of trigger jet
 C
         RETURN
         END
 C
 C
 C
         FUNCTION FJET(X,WGT)
         IMPLICIT REAL*8(A-H,O-Z)
         REAL HIPR1(100),HINT1(100)
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         DIMENSION X(10)
         PT2=(HINT1(1)**2/4.0-HIPR1(8)**2)*X(1)+HIPR1(8)**2
         XT=2.0*DSQRT(PT2)/HINT1(1)
         YMX1=DLOG(1.0/XT+DSQRT(1.0/XT**2-1.0))
         Y1=2.0*YMX1*X(2)-YMX1
         YMX2=DLOG(2.0/XT-DEXP(Y1))
         YMN2=DLOG(2.0/XT-DEXP(-Y1))
         Y2=(YMX2+YMN2)*X(3)-YMN2
         FJET=2.0*YMX1*(YMX2+YMN2)*(HINT1(1)**2/4.0-HIPR1(8)**2)
      &          *G(Y1,Y2,PT2)/2.0
         RETURN
         END
 C
 C
 C
         FUNCTION FJETRIG(X,WGT)
         IMPLICIT REAL*8(A-H,O-Z)
         REAL HIPR1(100),HINT1(100),PTMAX,PTMIN
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         DIMENSION X(10)
         PTMIN=ABS(HIPR1(10))-0.25
         PTMIN=MAX(PTMIN,HIPR1(8))
         AM2=0.D0
         IF(IHPR2(3).EQ.3) THEN
            AM2=HIPR1(7)**2
            PTMIN=MAX(0.0,HIPR1(10))
         ENDIF
         PTMAX=ABS(HIPR1(10))+0.25
         IF(HIPR1(10).LE.0.0) PTMAX=HINT1(1)/2.0-AM2
         IF(PTMAX.LE.PTMIN) PTMAX=PTMIN+0.25
         PT2=(PTMAX**2-PTMIN**2)*X(1)+PTMIN**2
         AMT2=PT2+AM2
         XT=2.0*DSQRT(AMT2)/HINT1(1)
         YMX1=DLOG(1.0/XT+DSQRT(1.0/XT**2-1.0))
         Y1=2.0*YMX1*X(2)-YMX1
         YMX2=DLOG(2.0/XT-DEXP(Y1))
         YMN2=DLOG(2.0/XT-DEXP(-Y1))
         Y2=(YMX2+YMN2)*X(3)-YMN2
         IF(IHPR2(3).EQ.3) THEN
            GTRIG=2.0*GHVQ(Y1,Y2,AMT2)
         ELSE IF(IHPR2(3).EQ.2) THEN
            GTRIG=2.0*GPHOTON(Y1,Y2,PT2)
         ELSE
            GTRIG=G(Y1,Y2,PT2)
         ENDIF
         FJETRIG=2.0*YMX1*(YMX2+YMN2)*(PTMAX**2-PTMIN**2)
      &          *GTRIG/2.0
         RETURN
         END
 C
 C
 C
         FUNCTION GHVQ(Y1,Y2,AMT2)
         IMPLICIT REAL*8 (A-H,O-Z)
         REAL HIPR1(100),HINT1(100)
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         DIMENSION F(2,7)
         XT=2.0*DSQRT(AMT2)/HINT1(1)
         X1=0.50*XT*(DEXP(Y1)+DEXP(Y2))
         X2=0.50*XT*(DEXP(-Y1)+DEXP(-Y2))
         SS=X1*X2*HINT1(1)**2
         AF=4.0
         IF(IHPR2(18).NE.0) AF=5.0
         DLAM=HIPR1(15)
         APH=12.0*3.1415926/(33.0-2.0*AF)/DLOG(AMT2/DLAM**2)
 C
         CALL PARTON(F,X1,X2,AMT2)
 C
         Gqq=4.0*(COSH(Y1-Y2)+HIPR1(7)**2/AMT2)/(1.D0+COSH(Y1-Y2))/9.0
      &      *(F(1,1)*F(2,2)+F(1,2)*F(2,1)+F(1,3)*F(2,4)
      &        +F(1,4)*F(2,3)+F(1,5)*F(2,6)+F(1,6)*F(2,5))
         Ggg=(8.D0*COSH(Y1-Y2)-1.D0)*(COSH(Y1-Y2)+2.0*HIPR1(7)**2/AMT2
      &      -2.0*HIPR1(7)**4/AMT2**2)/(1.0+COSH(Y1-Y2))/24.D0
      &      *F(1,7)*F(2,7)
 C
         GHVQ=(Gqq+Ggg)*HIPR1(23)*3.14159*APH**2/SS**2
         RETURN
         END
 C
 C
 C
         FUNCTION GPHOTON(Y1,Y2,PT2)
         IMPLICIT REAL*8 (A-H,O-Z)
         REAL HIPR1(100),HINT1(100)
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         DIMENSION F(2,7)
         XT=2.0*DSQRT(PT2)/HINT1(1)
         X1=0.50*XT*(DEXP(Y1)+DEXP(Y2))
         X2=0.50*XT*(DEXP(-Y1)+DEXP(-Y2))
         Z=DSQRT(1.D0-XT**2/X1/X2)
         SS=X1*X2*HINT1(1)**2
         T=-(1.0-Z)/2.0
         U=-(1.0+Z)/2.0
         AF=3.0
         DLAM=HIPR1(15)
         APH=12.0*3.1415926/(33.0-2.0*AF)/DLOG(PT2/DLAM**2)
         APHEM=1.0/137.0
 C
         CALL PARTON(F,X1,X2,PT2)
 C
         G11=-(U**2+1.0)/U/3.0*F(1,7)*(4.0*F(2,1)+4.0*F(2,2)
      &      +F(2,3)+F(2,4)+F(2,5)+F(2,6))/9.0
         G12=-(T**2+1.0)/T/3.0*F(2,7)*(4.0*F(1,1)+4.0*F(1,2)
      &      +F(1,3)+F(1,4)+F(1,5)+F(1,6))/9.0
         G2=8.0*(U**2+T**2)/U/T/9.0*(4.0*F(1,1)*F(2,2)
      &     +4.0*F(1,2)*F(2,1)+F(1,3)*F(2,4)+F(1,4)*F(2,3)
      &     +F(1,5)*F(2,6)+F(1,6)*F(2,5))/9.0
 C
         GPHOTON=(G11+G12+G2)*HIPR1(23)*3.14159*APH*APHEM/SS**2
         RETURN
         END
 C
 C
 C
 C
         FUNCTION G(Y1,Y2,PT2)
         IMPLICIT REAL*8 (A-H,O-Z)
         REAL HIPR1(100),HINT1(100)
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         DIMENSION F(2,7)
         XT=2.0*DSQRT(PT2)/HINT1(1)
         X1=0.50*XT*(DEXP(Y1)+DEXP(Y2))
         X2=0.50*XT*(DEXP(-Y1)+DEXP(-Y2))
         Z=DSQRT(1.D0-XT**2/X1/X2)
         SS=X1*X2*HINT1(1)**2
         T=-(1.0-Z)/2.0
         U=-(1.0+Z)/2.0
         AF=3.0
         DLAM=HIPR1(15)
         APH=12.0*3.1415926/(33.0-2.0*AF)/DLOG(PT2/DLAM**2)
 C
         CALL PARTON(F,X1,X2,PT2)
 C
         G11=( (F(1,1)+F(1,2))*(F(2,3)+F(2,4)+F(2,5)+F(2,6))
      &      +(F(1,3)+F(1,4))*(F(2,5)+F(2,6)) )*SUBCRS1(T,U)
 C
         G12=( (F(2,1)+F(2,2))*(F(1,3)+F(1,4)+F(1,5)+F(1,6))
      &      +(F(2,3)+F(2,4))*(F(1,5)+F(1,6)) )*SUBCRS1(U,T)
 C
         G13=(F(1,1)*F(2,1)+F(1,2)*F(2,2)+F(1,3)*F(2,3)+F(1,4)*F(2,4)
      &      +F(1,5)*F(2,5)+F(1,6)*F(2,6))*(SUBCRS1(U,T)
      &      +SUBCRS1(T,U)-8.D0/T/U/27.D0)
 C
         G2=(AF-1)*(F(1,1)*F(2,2)+F(2,1)*F(1,2)+F(1,3)*F(2,4)
      &     +F(2,3)*F(1,4)+F(1,5)*F(2,6)+F(2,5)*F(1,6))*SUBCRS2(T,U)
 C
         G31=(F(1,1)*F(2,2)+F(1,3)*F(2,4)+F(1,5)*F(2,6))*SUBCRS3(T,U)
         G32=(F(2,1)*F(1,2)+F(2,3)*F(1,4)+F(2,5)*F(1,6))*SUBCRS3(U,T)
 C
         G4=(F(1,1)*F(2,2)+F(2,1)*F(1,2)+F(1,3)*F(2,4)+F(2,3)*F(1,4)+
      1  F(1,5)*F(2,6)+F(2,5)*F(1,6))*SUBCRS4(T,U)
 C
         G5=AF*F(1,7)*F(2,7)*SUBCRS5(T,U)
 C
         G61=F(1,7)*(F(2,1)+F(2,2)+F(2,3)+F(2,4)+F(2,5)
      &      +F(2,6))*SUBCRS6(T,U)
         G62=F(2,7)*(F(1,1)+F(1,2)+F(1,3)+F(1,4)+F(1,5)
      &      +F(1,6))*SUBCRS6(U,T)
 C
         G7=F(1,7)*F(2,7)*SUBCRS7(T,U)
 C
         G=(G11+G12+G13+G2+G31+G32+G4+G5+G61+G62+G7)*HIPR1(17)*
      1  3.14159D0*APH**2/SS**2
         RETURN
         END
 C
 C
 C
         FUNCTION SUBCRS1(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS1=4.D0/9.D0*(1.D0+U**2)/T**2
         RETURN
         END
 C
 C
         FUNCTION SUBCRS2(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS2=4.D0/9.D0*(T**2+U**2)
         RETURN
         END
 C
 C
         FUNCTION SUBCRS3(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS3=4.D0/9.D0*(T**2+U**2+(1.D0+U**2)/T**2
      1  -2.D0*U**2/3.D0/T)
         RETURN
         END
 C
 C
         FUNCTION SUBCRS4(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS4=8.D0/3.D0*(T**2+U**2)*(4.D0/9.D0/T/U-1.D0)
         RETURN
         END
 C
 C
 C
         FUNCTION SUBCRS5(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS5=3.D0/8.D0*(T**2+U**2)*(4.D0/9.D0/T/U-1.D0)
         RETURN
         END
 C
 C
         FUNCTION SUBCRS6(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS6=(1.D0+U**2)*(1.D0/T**2-4.D0/U/9.D0)
         RETURN
         END
 C
 C
         FUNCTION SUBCRS7(T,U)
         IMPLICIT REAL*8(A-H,O-Z)
         SUBCRS7=9.D0/2.D0*(3.D0-T*U-U/T**2-T/U**2)
         RETURN
         END
 C
 C
 C
         SUBROUTINE PARTON(F,X1,X2,QQ)
         IMPLICIT REAL*8(A-H,O-Z)
         REAL HIPR1(100),HINT1(100)
         COMMON/HIPARNT/HIPR1,IHPR2(50),HINT1,IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/NJET/N,IP_CRS
         SAVE  /NJET/
         DIMENSION F(2,7) 
         DLAM=HIPR1(15)
         Q0=HIPR1(16)
         S=DLOG(DLOG(QQ/DLAM**2)/DLOG(Q0**2/DLAM**2))
         IF(IHPR2(7).EQ.2) GO TO 200
 C*******************************************************
         AT1=0.419+0.004*S-0.007*S**2
         AT2=3.460+0.724*S-0.066*S**2
         GMUD=4.40-4.86*S+1.33*S**2
         AT3=0.763-0.237*S+0.026*S**2
         AT4=4.00+0.627*S-0.019*S**2
         GMD=-0.421*S+0.033*S**2
 C*******************************************************
         CAS=1.265-1.132*S+0.293*S**2
         AS=-0.372*S-0.029*S**2
         BS=8.05+1.59*S-0.153*S**2
         APHS=6.31*S-0.273*S**2
         BTAS=-10.5*S-3.17*S**2
         GMS=14.7*S+9.80*S**2
 C********************************************************
 C       CAC=0.135*S-0.075*S**2
 C       AC=-0.036-0.222*S-0.058*S**2
 C       BC=6.35+3.26*S-0.909*S**2
 C       APHC=-3.03*S+1.50*S**2
 C       BTAC=17.4*S-11.3*S**2
 C       GMC=-17.9*S+15.6*S**2
 C***********************************************************
         CAG=1.56-1.71*S+0.638*S**2
         AG=-0.949*S+0.325*S**2
         BG=6.0+1.44*S-1.05*S**2
         APHG=9.0-7.19*S+0.255*S**2
         BTAG=-16.5*S+10.9*S**2
         GMG=15.3*S-10.1*S**2
         GO TO 300
 C********************************************************
 200     AT1=0.374+0.014*S
         AT2=3.33+0.753*S-0.076*S**2
         GMUD=6.03-6.22*S+1.56*S**2
         AT3=0.761-0.232*S+0.023*S**2
         AT4=3.83+0.627*S-0.019*S**2
         GMD=-0.418*S+0.036*S**2
 C************************************
         CAS=1.67-1.92*S+0.582*S**2
         AS=-0.273*S-0.164*S**2
         BS=9.15+0.530*S-0.763*S**2
         APHS=15.7*S-2.83*S**2
         BTAS=-101.0*S+44.7*S**2
         GMS=223.0*S-117.0*S**2
 C*********************************
 C       CAC=0.067*S-0.031*S**2
 C       AC=-0.120-0.233*S-0.023*S**2
 C       BC=3.51+3.66*S-0.453*S**2
 C       APHC=-0.474*S+0.358*S**2
 C       BTAC=9.50*S-5.43*S**2
 C       GMC=-16.6*S+15.5*S**2
 C**********************************
         CAG=0.879-0.971*S+0.434*S**2
         AG=-1.16*S+0.476*S**2
         BG=4.0+1.23*S-0.254*S**2
         APHG=9.0-5.64*S-0.817*S**2
         BTAG=-7.54*S+5.50*S**2
         GMG=-0.596*S+1.26*S**2
 C*********************************
 300     B12=DEXP(GMRE(AT1)+GMRE(AT2+1.D0)-GMRE(AT1+AT2+1.D0))
         B34=DEXP(GMRE(AT3)+GMRE(AT4+1.D0)-GMRE(AT3+AT4+1.D0))
         CNUD=3.D0/B12/(1.D0+GMUD*AT1/(AT1+AT2+1.D0))
         CND=1.D0/B34/(1.D0+GMD*AT3/(AT3+AT4+1.D0))
 C********************************************************
 C       FUD=X*(U+D)
 C       FS=X*2(UBAR+DBAR+SBAR)  AND UBAR=DBAR=SBAR
 C*******************************************************
         FUD1=CNUD*X1**AT1*(1.D0-X1)**AT2*(1.D0+GMUD*X1)
         FS1=CAS*X1**AS*(1.D0-X1)**BS*(1.D0+APHS*X1
      &      +BTAS*X1**2+GMS*X1**3)
         F(1,3)=CND*X1**AT3*(1.D0-X1)**AT4*(1.D0+GMD*X1)+FS1/6.D0
         F(1,1)=FUD1-F(1,3)+FS1/3.D0
         F(1,2)=FS1/6.D0
         F(1,4)=FS1/6.D0
         F(1,5)=FS1/6.D0
         F(1,6)=FS1/6.D0
         F(1,7)=CAG*X1**AG*(1.D0-X1)**BG*(1.D0+APHG*X1
      &         +BTAG*X1**2+GMG*X1**3)
 C
         FUD2=CNUD*X2**AT1*(1.D0-X2)**AT2*(1.D0+GMUD*X2)
         FS2=CAS*X2**AS*(1.D0-X2)**BS*(1.D0+APHS*X2
      &      +BTAS*X2**2+GMS*X2**3)
         F(2,3)=CND*X2**AT3*(1.D0-X2)**AT4*(1.D0+GMD*X2)+FS2/6.D0
         F(2,1)=FUD2-F(2,3)+FS2/3.D0
         F(2,2)=FS2/6.D0
         F(2,4)=FS2/6.D0
         F(2,5)=FS2/6.D0
         F(2,6)=FS2/6.D0
         F(2,7)=CAG*X2**AG*(1.D0-X2)**BG*(1.D0+APHG*X2
      &         +BTAG*X2**2+GMG*X2**3)
 C***********Nuclear effect on the structure function****************
 C
         IF(IHPR2(6).EQ.1 .AND. IHNT2(1).GT.1) THEN
            AAX=1.193*ALOG(FLOAT(IHNT2(1)))**0.16666666
            RRX=AAX*(X1**3-1.2*X1**2+0.21*X1)+1.0
      &         +1.079*(FLOAT(IHNT2(1))**0.33333333-1.0)
      &         /DLOG(IHNT2(1)+1.0D0)*DSQRT(X1)*DEXP(-X1**2/0.01)
            IF(IP_CRS.EQ.1 .OR.IP_CRS.EQ.3) RRX=DEXP(-X1**2/0.01)
            DO 400 I=1,7
               F(1,I)=RRX*F(1,I)
  400       CONTINUE
         ENDIF
         IF(IHPR2(6).EQ.1 .AND. IHNT2(3).GT.1) THEN
            AAX=1.193*ALOG(FLOAT(IHNT2(3)))**0.16666666
            RRX=AAX*(X2**3-1.2*X2**2+0.21*X2)+1.0
      &         +1.079*(FLOAT(IHNT2(3))**0.33333-1.0)
      &         /DLOG(IHNT2(3)+1.0D0)*DSQRT(X2)*DEXP(-X2**2/0.01)
            IF(IP_CRS.EQ.2 .OR. IP_CRS.EQ.3) RRX=DEXP(-X2**2/0.01)
            DO 500 I=1,7
               F(2,I)=RRX*F(2,I)
  500       CONTINUE
         ENDIF
 c
         RETURN
         END
 C
 C
 C
         FUNCTION GMRE(X)
         IMPLICIT REAL*8(A-H,O-Z)
         Z=X
         IF(X.GT.3.0D0) GO TO 10
         Z=X+3.D0
 10      GMRE=0.5D0*DLOG(2.D0*3.14159265D0/Z)+Z*DLOG(Z)-Z+DLOG(1.D0
      1  +1.D0/12.D0/Z+1.D0/288.D0/Z**2-139.D0/51840.D0/Z**3
      1  -571.D0/2488320.D0/Z**4)
         IF(Z.EQ.X) GO TO 20
         GMRE=GMRE-DLOG(Z-1.D0)-DLOG(Z-2.D0)-DLOG(Z-3.D0)
 20      CONTINUE
         RETURN
         END
 c
 C
 C
         FUNCTION GMIN(N)
         IMPLICIT REAL*8(A-H,O-Z)
         GA=0.
         IF(N.LE.2) GO TO 20
         DO 10 I=1,N-1
         Z=I
         GA=GA+DLOG(Z)
 10      CONTINUE
 20      GMIN=GA
         RETURN
         END
 C
 C
 C***************************************************************
 
         BLOCK DATA HIDATA
         REAL*8 XL(10),XU(10),ACC
         COMMON/BVEG1/XL,XU,ACC,NDIM,NCALL,ITMX,NPRN
         SAVE  /BVEG1/
         COMMON/SEEDVAX/NUM1
         SAVE  /SEEDVAX/
         COMMON/HIPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
         SAVE  /HIPARNT/
         COMMON/RANSEED/NSEED
         SAVE  /RANSEED/
         COMMON/HIMAIN1/ NATT,EATT,JATT,NT,NP,N0,N01,N10,N11
         SAVE  /HIMAIN1/
         COMMON/HIMAIN2/KATT(130000,4),PATT(130000,4)
         SAVE  /HIMAIN2/
         COMMON/HISTRNG/NFP(300,15),PP(300,15),NFT(300,15),PT(300,15)
         SAVE  /HISTRNG/
         COMMON/HIJCRDN/YP(3,300),YT(3,300)
         SAVE  /HIJCRDN/
         COMMON/HIJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
      &               PJPY(300,500),PJPZ(300,500),PJPE(300,500),
      &               PJPM(300,500),NTJ(300),KFTJ(300,500),
      &               PJTX(300,500),PJTY(300,500),PJTZ(300,500),
      &               PJTE(300,500),PJTM(300,500)
         SAVE  /HIJJET1/
         COMMON/HIJJET2/NSG,NJSG(900),IASG(900,3),K1SG(900,100),
      &          K2SG(900,100),PXSG(900,100),PYSG(900,100),
      &          PZSG(900,100),PESG(900,100),PMSG(900,100)
         SAVE  /HIJJET2/
         COMMON/HIJDAT/HIDAT0(10,10),HIDAT(10)
         SAVE  /HIJDAT/
         COMMON/HIPYINT/MINT4,MINT5,ATCO(200,20),ATXS(0:200)
         SAVE  /HIPYINT/
         DATA NUM1/30123984/,XL/10*0.D0/,XU/10*1.D0/
         DATA NCALL/1000/ITMX/100/ACC/0.01/NPRN/0/
 C...give all the switchs and parameters the default values
 
         DATA NSEED/74769375/
         DATA HIPR1/
      &  1.5,  0.35, 0.5,  0.9,  2.0,  0.1,  1.5,  2.0, -1.0, -2.25,
      &  2.0,  0.5,  1.0,  2.0,  0.2,  2.0,  2.5,  0.3,  0.1,  1.4,
      &  1.6,  1.0,  2.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.4,  57.0,
      &  28.5, 3.9,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  
      &  3.141592654,
      &  0.0,  0.4,  0.1,  1.5,  0.1, 0.25, 0.0,  0.5,  0.0,  0.0,
      &  50*0.0/
 
         DATA IHPR2/
      &  1,    3,    0,    1,    1,    1,    1,    10,    0,    0,
      &  1,    1,    1,    1,    0,    0,    1,     0,    0,    1,
      &  30*0/
 
         DATA HINT1/100*0/
         DATA IHNT2/50*0/
 
 C...initialize all the data common blocks
         DATA NATT/0/EATT/0.0/JATT/0/NT/0/NP/0/N0/0/N01/0/N10/0/N11/0/
         DATA KATT/520000*0/PATT/520000*0.0/
 
         DATA NFP/4500*0/PP/4500*0.0/NFT/4500*0/PT/4500*0.0/
 
         DATA YP/900*0.0/YT/900*0.0/
 
         DATA NPJ/300*0/KFPJ/150000*0/PJPX/150000*0.0/PJPY/150000*0.0/
      &  PJPZ/150000*0.0/PJPE/150000*0.0/PJPM/150000*0.0/
         DATA NTJ/300*0/KFTJ/150000*0/PJTX/150000*0.0/PJTY/150000*0.0/
      &  PJTZ/150000*0.0/PJTE/150000*0.0/PJTM/150000*0.0/
 
         DATA NSG/0/NJSG/900*0/IASG/2700*0/K1SG/90000*0/K2SG/90000*0/
      &  PXSG/90000*0.0/PYSG/90000*0.0/PZSG/90000*0.0/PESG/90000*0.0/
      &  PMSG/90000*0.0/
         DATA MINT4/0/MINT5/0/ATCO/4000*0.0/ATXS/201*0.0/
         DATA (HIDAT0(1,I),I=1,10)/0.0,0.0,0.0,0.0,0.0,0.0,2.25,
      &          2.5,4.0,4.1/
         DATA (HIDAT0(2,I),I=1,10)/2.0,3.0,5.0,6.0,7.0,8.0,8.0,10.0,
      &          10.0,10.0/
         DATA (HIDAT0(3,I),I=1,10)/1.0,0.8,0.8,0.7,0.45,0.215,
      &          0.21,0.19,0.19,0.19/
         DATA (HIDAT0(4,I),I=1,10)/0.35,0.35,0.3,0.3,0.3,0.3,
      &          0.5,0.6,0.6,0.6/
         DATA (HIDAT0(5,I),I=1,10)/23.8,24.0,26.0,26.2,27.0,28.5,28.5,
      &          28.5,28.5,28.5/
         DATA ((HIDAT0(J,I),I=1,10),J=6,9)/40*0.0/
         DATA (HIDAT0(10,I),I=1,10)/5.0,20.0,53.0,62.0,100.0,200.0,
      &          546.0,900.0,1800.0,4000.0/
         DATA HIDAT/10*0.0/
         END
 C*******************************************************************
 C
 C
 C
 C
 C*******************************************************************
 C   SUBROUTINE PERFORMS N-DIMENSIONAL MONTE CARLO INTEG'N
 C      - BY G.P. LEPAGE   SEPT 1976/(REV)APR 1978
 C*******************************************************************
 C
       SUBROUTINE VEGAS(FXN,AVGI,SD,CHI2A)
       IMPLICIT REAL*8(A-H,O-Z)
       COMMON/BVEG1/XL(10),XU(10),ACC,NDIM,NCALL,ITMX,NPRN
       SAVE  /BVEG1/
       COMMON/BVEG2/XI(50,10),SI,SI2,SWGT,SCHI,NDO,IT
       SAVE  /BVEG2/
       COMMON/BVEG3/F,TI,TSI   
       SAVE  /BVEG3/
       EXTERNAL FXN
       DIMENSION D(50,10),DI(50,10),XIN(50),R(50),DX(10),DT(10),X(10)
      1   ,KG(10),IA(10)
       REAL*4 QRAN(10)
       DATA NDMX/50/,ALPH/1.5D0/,ONE/1.D0/,MDS/-1/
 
        save NDMX, ALPH, ONE, MDS    !uzhi
 C
       NDO=1
       DO 1 J=1,NDIM
 1     XI(1,J)=ONE
 C
       ENTRY VEGAS1(FXN,AVGI,SD,CHI2A)
 C         - INITIALIZES CUMMULATIVE VARIABLES, BUT NOT GRID
       IT=0
       SI=0.
       SI2=SI
       SWGT=SI
       SCHI=SI
 C
       ENTRY VEGAS2(FXN,AVGI,SD,CHI2A)
 C         - NO INITIALIZATION
       ND=NDMX
       NG=1
       IF(MDS.EQ.0) GO TO 2
       NG=(NCALL/2.)**(1./NDIM)
       MDS=1
       IF((2*NG-NDMX).LT.0) GO TO 2
       MDS=-1
       NPG=NG/NDMX+1
       ND=NG/NPG
       NG=NPG*ND
 2     K=NG**NDIM
       NPG=NCALL/K
       IF(NPG.LT.2) NPG=2
       CALLS=NPG*K
       DXG=ONE/NG
       DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
       XND=ND
       NDM=ND-1
       DXG=DXG*XND
       XJAC=ONE/CALLS
       DO 3 J=1,NDIM
 c***this is the line 50
       DX(J)=XU(J)-XL(J)
 3     XJAC=XJAC*DX(J)
 C
 C   REBIN PRESERVING BIN DENSITY
 C
       IF(ND.EQ.NDO) GO TO 8
       RC=NDO/XND
       DO 7 J=1,NDIM
       K=0
       XN=0.
       DR=XN
       I=K
 4     K=K+1
       DR=DR+ONE
       XO=XN
       XN=XI(K,J)
 5     IF(RC.GT.DR) GO TO 4
       I=I+1
       DR=DR-RC
       XIN(I)=XN-(XN-XO)*DR
       IF(I.LT.NDM) GO TO 5
       DO 6 I=1,NDM
 6     XI(I,J)=XIN(I)
 7     XI(ND,J)=ONE
       NDO=ND
 C
 8     CONTINUE
       IF(NPRN.NE.0) WRITE(16,200) NDIM,CALLS,IT,ITMX,ACC,MDS,ND
      1                           ,(XL(J),XU(J),J=1,NDIM)
 C
       ENTRY VEGAS3(FXN,AVGI,SD,CHI2A)
 C         - MAIN INTEGRATION LOOP
 9     IT=IT+1
       TI=0.
       TSI=TI
       DO 10 J=1,NDIM
       KG(J)=1
       DO 10 I=1,ND
       D(I,J)=TI
 10    DI(I,J)=TI
 C
 11    FB=0.
       F2B=FB
       K=0
 12    K=K+1
       CALL ARAN9(QRAN,NDIM)
       WGT=XJAC
       DO 15 J=1,NDIM
       XN=(KG(J)-QRAN(J))*DXG+ONE
 c*****this is the line 100
       IA(J)=XN
       IF(IA(J).GT.1) GO TO 13
       XO=XI(IA(J),J)
       RC=(XN-IA(J))*XO
       GO TO 14
 13    XO=XI(IA(J),J)-XI(IA(J)-1,J)
       RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
 14    X(J)=XL(J)+RC*DX(J)
       WGT=WGT*XO*XND
 15    CONTINUE
 C
       F=WGT
       F=F*FXN(X,WGT)
       F2=F*F
       FB=FB+F
       F2B=F2B+F2
       DO 16 J=1,NDIM
       DI(IA(J),J)=DI(IA(J),J)+F
 16    IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
       IF(K.LT.NPG) GO TO 12
 C
       F2B=DSQRT(F2B*NPG)
       F2B=(F2B-FB)*(F2B+FB)
       TI=TI+FB
       TSI=TSI+F2B
       IF(MDS.GE.0) GO TO 18
       DO 17 J=1,NDIM
 17    D(IA(J),J)=D(IA(J),J)+F2B
 18    K=NDIM
 19    KG(K)=MOD(KG(K),NG)+1
       IF(KG(K).NE.1) GO TO 11
       K=K-1
       IF(K.GT.0) GO TO 19
 C
 C   FINAL RESULTS FOR THIS ITERATION
 C
       TSI=TSI*DV2G
       TI2=TI*TI
       WGT=TI2/(TSI+1.0d-37)
       SI=SI+TI*WGT
       SI2=SI2+TI2
       SWGT=SWGT+WGT
       SWGT=SWGT+1.0D-37
       SI2=SI2+1.0D-37
       SCHI=SCHI+TI2*WGT
       AVGI=SI/(SWGT)
       SD=SWGT*IT/(SI2)
       CHI2A=SD*(SCHI/SWGT-AVGI*AVGI)/(IT-.999)
       SD=DSQRT(ONE/SD)
 C****this is the line 150
       IF(NPRN.EQ.0) GO TO 21
       TSI=DSQRT(TSI)
       WRITE(16,201) IT,TI,TSI,AVGI,SD,CHI2A
       IF(NPRN.GE.0) GO TO 21
       DO 20 J=1,NDIM
 20    WRITE(16,202) J,(XI(I,J),DI(I,J),D(I,J),I=1,ND)
 C
 C   REFINE GRID
 C
 21    DO 23 J=1,NDIM
       XO=D(1,J)
       XN=D(2,J)
       D(1,J)=(XO+XN)/2.
       DT(J)=D(1,J)
       DO 22 I=2,NDM
       D(I,J)=XO+XN
       XO=XN
       XN=D(I+1,J)
       D(I,J)=(D(I,J)+XN)/3.
 22    DT(J)=DT(J)+D(I,J)
       D(ND,J)=(XN+XO)/2.
 23    DT(J)=DT(J)+D(ND,J)
 C
       DO 28 J=1,NDIM
       RC=0.
       DO 24 I=1,ND
       R(I)=0.
       IF (DT(J).GE.1.0D18) THEN
        WRITE(6,*) '************** A SINGULARITY >1.0D18'
 C      WRITE(5,1111)
 C1111  FORMAT(1X,'**************IMPORTANT NOTICE***************')
 C      WRITE(5,1112)
 C1112  FORMAT(1X,'THE INTEGRAND GIVES RISE A SINGULARITY >1.0D18')
 C      WRITE(5,1113)
 C1113  FORMAT(1X,'PLEASE CHECK THE INTEGRAND AND THE LIMITS')
 C      WRITE(5,1114)
 C1114  FORMAT(1X,'**************END NOTICE*************')
       END IF    
       IF(D(I,J).LE.1.0D-18) GO TO 24
       XO=DT(J)/D(I,J)
       R(I)=((XO-ONE)/XO/DLOG(XO))**ALPH
 24    RC=RC+R(I)
       RC=RC/XND
       K=0
       XN=0.
       DR=XN
       I=K
 25    K=K+1
       DR=DR+R(K)
       XO=XN
 c****this is the line 200
       XN=XI(K,J)
 26    IF(RC.GT.DR) GO TO 25
       I=I+1
       DR=DR-RC
       XIN(I)=XN-(XN-XO)*DR/(R(K)+1.0d-30)
       IF(I.LT.NDM) GO TO 26
       DO 27 I=1,NDM
 27    XI(I,J)=XIN(I)
 28    XI(ND,J)=ONE
 C
       IF(IT.LT.ITMX.AND.ACC*DABS(AVGI).LT.SD) GO TO 9
 200   FORMAT('0INPUT PARAMETERS FOR VEGAS:  NDIM=',I3,'  NCALL=',F8.0
      1    /28X,'  IT=',I5,'  ITMX=',I5/28X,'  ACC=',G9.3
      2    /28X,'  MDS=',I3,'   ND=',I4/28X,'  (XL,XU)=',
      3    (T40,'( ',G12.6,' , ',G12.6,' )'))
 201   FORMAT(///' INTEGRATION BY VEGAS' / '0ITERATION NO.',I3,
      1    ':   INTEGRAL =',G14.8/21X,'STD DEV  =',G10.4 /
      2    ' ACCUMULATED RESULTS:   INTEGRAL =',G14.8 /
      3    24X,'STD DEV  =',G10.4 / 24X,'CHI**2 PER IT''N =',G10.4)
 202   FORMAT('0DATA FOR AXIS',I2 / ' ',6X,'X',7X,'  DELT I  ',
      1    2X,' CONV''CE  ',11X,'X',7X,'  DELT I  ',2X,' CONV''CE  '
      2   ,11X,'X',7X,'  DELT I  ',2X,' CONV''CE  ' /
      2    (' ',3G12.4,5X,3G12.4,5X,3G12.4))
       RETURN
       END
 C
 C
       SUBROUTINE ARAN9(QRAN,NDIM)
       DIMENSION QRAN(10)
       COMMON/SEEDVAX/NUM1
       DO 1 I=1,NDIM
     1 QRAN(I)=RLU(0)
       RETURN
       END
 
 C
 C
 C*********GAUSSIAN ONE-DIMENSIONAL INTEGRATION PROGRAM*************
 C
         FUNCTION GAUSS1(F,A,B,EPS)
         EXTERNAL F
         DIMENSION W(12),X(12)
         DATA CONST/1.0E-12/
         DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
      &         .0622535,0.0951585,.1246290,.1495960,.1691565,
      &         .1826034,.1894506/
         DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
      &         .9445750,0.8656312,.7554044,.6178762,.4580168,
      &         .2816036,.0950125/
         DELTA=CONST*ABS(A-B)
         GAUSS1=0.0
         AA=A
 5       Y=B-AA
         IF(ABS(Y).LE.DELTA) RETURN
 2       BB=AA+Y
         C1=0.5*(AA+BB)
         C2=C1-AA
         S8=0.0
         S16=0.0
         DO 1 I=1,4
         U=X(I)*C2
 1       S8=S8+W(I)*(F(C1+U)+F(C1-U))
         DO 3 I=5,12
         U=X(I)*C2
 3       S16=S16+W(I)*(F(C1+U)+F(C1-U))
         S8=S8*C2
         S16=S16*C2
         IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
         GAUSS1=GAUSS1+S16
         AA=BB
         GOTO 5
 4       Y=0.5*Y
         IF(ABS(Y).GT.DELTA) GOTO 2
         WRITE(6,7)
         GAUSS1=0.0
         RETURN
 7       FORMAT(1X,'GAUSS1....TOO HIGH ACURACY REQUIRED')
         END
 C
 C
 C
         FUNCTION GAUSS2(F,A,B,EPS)
         EXTERNAL F
         DIMENSION W(12),X(12)
         DATA CONST/1.0E-12/
         DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
      &         .0622535,0.0951585,.1246290,.1495960,.1691565,
      &         .1826034,.1894506/
         DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
      &         .9445750,0.8656312,.7554044,.6178762,.4580168,
      &         .2816036,.0950125/
         DELTA=CONST*ABS(A-B)
         GAUSS2=0.0
         AA=A
 5       Y=B-AA
         IF(ABS(Y).LE.DELTA) RETURN
 2       BB=AA+Y
         C1=0.5*(AA+BB)
         C2=C1-AA
         S8=0.0
         S16=0.0
         DO 1 I=1,4
         U=X(I)*C2
 1       S8=S8+W(I)*(F(C1+U)+F(C1-U))
         DO 3 I=5,12
         U=X(I)*C2
 3       S16=S16+W(I)*(F(C1+U)+F(C1-U))
         S8=S8*C2
         S16=S16*C2
         IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
         GAUSS2=GAUSS2+S16
         AA=BB
         GOTO 5
 4       Y=0.5*Y
         IF(ABS(Y).GT.DELTA) GOTO 2
         WRITE(6,7)
         GAUSS2=0.0
         RETURN
 7       FORMAT(1X,'GAUSS2....TOO HIGH ACURACY REQUIRED')
         END
 C
 C
 C
         FUNCTION GAUSS3(F,A,B,EPS)
         EXTERNAL F
         DIMENSION W(12),X(12)
         DATA CONST/1.0E-12/
         DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
      &         .0622535,0.0951585,.1246290,.1495960,.1691565,
      &         .1826034,.1894506/
         DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
      &         .9445750,0.8656312,.7554044,.6178762,.4580168,
      &         .2816036,.0950125/
         DELTA=CONST*ABS(A-B)
         GAUSS3=0.0
         AA=A
 5       Y=B-AA
         IF(ABS(Y).LE.DELTA) RETURN
 2       BB=AA+Y
         C1=0.5*(AA+BB)
         C2=C1-AA
         S8=0.0
         S16=0.0
         DO 1 I=1,4
         U=X(I)*C2
 1       S8=S8+W(I)*(F(C1+U)+F(C1-U))
         DO 3 I=5,12
         U=X(I)*C2
 3       S16=S16+W(I)*(F(C1+U)+F(C1-U))
         S8=S8*C2
         S16=S16*C2
         IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
         GAUSS3=GAUSS3+S16
         AA=BB
         GOTO 5
 4       Y=0.5*Y
         IF(ABS(Y).GT.DELTA) GOTO 2
         WRITE(6,7)
         GAUSS3=0.0
         RETURN
 7       FORMAT(1X,'GAUSS3....TOO HIGH ACURACY REQUIRED')
         END
 C
 C
 C
 C
         FUNCTION GAUSS4(F,A,B,EPS)
         EXTERNAL F
         DIMENSION W(12),X(12)
         DATA CONST/1.0E-12/
         DATA W/0.1012285,.2223810,.3137067,.3623838,.0271525,
      &         .0622535,0.0951585,.1246290,.1495960,.1691565,
      &         .1826034,.1894506/
         DATA X/0.9602899,.7966665,.5255324,.1834346,.9894009,
      &         .9445750,0.8656312,.7554044,.6178762,.4580168,
      &         .2816036,.0950125/
         DELTA=CONST*ABS(A-B)
         GAUSS4=0.0
         AA=A
 5       Y=B-AA
         IF(ABS(Y).LE.DELTA) RETURN
 2       BB=AA+Y
         C1=0.5*(AA+BB)
         C2=C1-AA
         S8=0.0
         S16=0.0
         DO 1 I=1,4
         U=X(I)*C2
 1       S8=S8+W(I)*(F(C1+U)+F(C1-U))
         DO 3 I=5,12
         U=X(I)*C2
 3       S16=S16+W(I)*(F(C1+U)+F(C1-U))
         S8=S8*C2
         S16=S16*C2
         IF(ABS(S16-S8).GT.EPS*(1.+ABS(S16))) GOTO 4
         GAUSS4=GAUSS4+S16
         AA=BB
         GOTO 5
 4       Y=0.5*Y
         IF(ABS(Y).GT.DELTA) GOTO 2
         WRITE(6,7)
         GAUSS4=0.0
         RETURN
 7       FORMAT(1X,'GAUSS4....TOO HIGH ACURACY REQUIRED')
         END
 C
 C
 C
 C
 C
         SUBROUTINE TITLE
         WRITE(6,200)
 200     FORMAT(//10X,
      &  '**************************************************'/10X,
      &  '*     |      \\       _______      /  ------/     *'/10X,
      &  '*   ----- ------     |_____|     /_/     /       *'/10X,
      &  '*    ||\\    /        |_____|      /    / \\       *'/10X,
      &  '*    /| \\  /_/       /_______    /_  /    \\_     *'/10X,
      &  '*   / |     / /     /  /  / |        -------     *'/10X,
      &  '*     |    / /\\       /  /  |     /     |        *'/10X,
      &  '*     |   / /  \\     /  / \\_|    /   -------     *'/10X,
      &  '*                                                *'/10X,
      &  '**************************************************'/10X,
      &  '                      HIJING                      '/10X,
      &  '       Heavy Ion Jet INteraction Generator        '/10X,
      &  '                        by                        '/10X,
      &  '           X. N. Wang  and  M. Gyulassy           '/10X,
      &  '           Lawrence Berkeley Laboratory           '//) 
         RETURN
         END