hijing1.383.f

Go to the documentation of this file.

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