hijing1.383.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	hijing (BMIN0, BMAX0)
subroutine	hijset (efrm)
function	fnkick (X)
function	fnkick2 (X)
function	fnstru (X)
function	fnstrum (X)
function	fnstrus (X)
subroutine	hiboost
subroutine	quench (JPJT, NTP)
subroutine	hijfrg (JTP, NTP, IERROR)
subroutine	hijsrt (JPJT, NPT)
subroutine	attrad (IERROR)
subroutine	ar3jet (S, X1, X3, JL)
subroutine	arorie (S, X1, X3, JL)
subroutine	atrobo (THE, PHI, BEX, BEY, BEZ, IMIN, IMAX, IERROR)
subroutine	hijhrd (JP, JT, JOUT, JFLG, IOPJET0)
subroutine	jetini (JP, JT, I_TRIG)
subroutine	hijini
subroutine	attflv (ID, IDQ, IDQQ)
subroutine	hijcsc (JP, JT)
subroutine	hijels (PSC1, PSC2)
subroutine	hijsft (JP, JT, JOUT, IERROR)
subroutine	hijflv (ID)
subroutine	hiptdi (PT, PTMAX, IOPT)
subroutine	hijwds (IA, IDH, XHIGH)
function	wdsax (X)
function	rwdsax (X)
function	wdsax1 (X)
function	wdsax2 (X)
function	profile (XB)
function	flap (X)
function	fgp1 (X)
function	fgp2 (X)
function	fgp3 (X)
function	flap1 (X)
function	flap2 (X)
subroutine	hifun (I, XMIN, XMAX, FHB)
function	hirnd (I)
function	hirnd2 (I, XMIN, XMAX)
subroutine	hijcrs
function	ftot (X)
function	fhin (X)
function	ftotjet (X)
function	ftotrig (X)
function	fnjet (X)
function	sgmin (N)
function	omg0 (X)
function	romg (X)
function	bk (X)
subroutine	crsjet
function	fjet (X, WGT)
function	fjetrig (X, WGT)
function	ghvq (Y1, Y2, AMT2)
function	gphoton (Y1, Y2, PT2)
function	g (Y1, Y2, PT2)
function	subcrs1 (T, U)
function	subcrs2 (T, U)
function	subcrs3 (T, U)
function	subcrs4 (T, U)
function	subcrs5 (T, U)
function	subcrs6 (T, U)
function	subcrs7 (T, U)
subroutine	parton (F, X1, X2, QQ)
function	gmre (X)
function	gmin (N)
subroutine	vegas (FXN, AVGI, SD, CHI2A)
subroutine	aran9 (QRAN, NDIM)
function	gauss1 (F, A, B, EPS)
function	gauss2 (F, A, B, EPS)
function	gauss3 (F, A, B, EPS)
function	gauss4 (F, A, B, EPS)
subroutine	title

Function/Subroutine Documentation

ar3jet()

subroutine ar3jet	(		S,
			X1,
			X3,
			JL
	)

Definition at line 2195 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

aran9()

subroutine aran9	(	dimension(10)	QRAN,
			NDIM
	)

Definition at line 5798 of file hijing1.383.f.

      dimension qran(10)
      common/seedvax/num1
      DO 1 i=1,ndim
    1 qran(i)=rlu(0)
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

arorie()

subroutine arorie	(		S,
			X1,
			X3,
			JL
	)

Definition at line 2251 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

atrobo()

subroutine atrobo	(		THE,
			PHI,
			BEX,
			BEY,
			BEZ,
			IMIN,
			IMAX,
			IERROR
	)

Definition at line 2308 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

attflv()

subroutine attflv	(		ID,
			IDQ,
			IDQQ
	)

Definition at line 3393 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

attrad()

subroutine attrad

(

IERROR

)

Definition at line 2071 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

bk()

function bk

(

X

)

Definition at line 5002 of file hijing1.383.f.

    COMMON /besel/x4
    SAVE   /besel/
    bk=exp(-x)*(x**2-x4**2)**2.50/15.0
    RETURN

Here is the caller graph for this function:

crsjet()

subroutine crsjet

(

)

Definition at line 5013 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fgp1()

function fgp1

(

X

)

Definition at line 4708 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

fgp2()

function fgp2

(

X

)

Definition at line 4719 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

fgp3()

function fgp3

(

X

)

Definition at line 4730 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

fhin()

function fhin

(

X

)

Definition at line 4895 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fjet()

function fjet	(	dimension(10)	X,
			WGT
	)

Definition at line 5083 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fjetrig()

function fjetrig	(	dimension(10)	X,
			WGT
	)

Definition at line 5103 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

flap()

function flap

(

X

)

Definition at line 4697 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

flap1()

function flap1

(

X

)

Definition at line 4740 of file hijing1.383.f.

        common/pact/bb,b1,phi,z1
        SAVE  /pact/
        r=sqrt(bb**2+x**2)
        flap1=wdsax1(r)
        RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

flap2()

function flap2

(

X

)

Definition at line 4749 of file hijing1.383.f.

        common/pact/bb,b1,phi,z1
        SAVE  /pact/
        r=sqrt(bb**2+x**2)
        flap2=wdsax2(r)
        RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

fnjet()

function fnjet

(

X

)

Definition at line 4926 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fnkick()

function fnkick

(

X

)

Definition at line 1006 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fnkick2()

function fnkick2

(

X

)

Definition at line 1015 of file hijing1.383.f.

    common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
    SAVE  /hiparnt/
    fnkick2=x*exp(-2.0*x/hipr1(42))
    RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

fnstru()

function fnstru

(

X

)

Definition at line 1024 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fnstrum()

function fnstrum

(

X

)

Definition at line 1034 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

fnstrus()

function fnstrus

(

X

)

Definition at line 1043 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

ftot()

function ftot

(

X

)

Definition at line 4885 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

ftotjet()

function ftotjet

(

X

)

Definition at line 4905 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

ftotrig()

function ftotrig

(

X

)

Definition at line 4915 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

g()

function g	(		Y1,
			Y2,
			PT2
	)

Definition at line 5206 of file hijing1.383.f.

    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

Here is the call graph for this function:

gauss1()

function gauss1	(	external	F,
			A,
			B,
			EPS
	)

Definition at line 5810 of file hijing1.383.f.

    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')

Here is the caller graph for this function:

gauss2()

function gauss2	(	external	F,
			A,
			B,
			EPS
	)

Definition at line 5852 of file hijing1.383.f.

    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')

Here is the caller graph for this function:

gauss3()

function gauss3	(	external	F,
			A,
			B,
			EPS
	)

Definition at line 5894 of file hijing1.383.f.

    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')

Here is the caller graph for this function:

gauss4()

function gauss4	(	external	F,
			A,
			B,
			EPS
	)

Definition at line 5937 of file hijing1.383.f.

    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')

Here is the caller graph for this function:

ghvq()

function ghvq	(		Y1,
			Y2,
			AMT2
	)

Definition at line 5141 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

gmin()

function gmin

(

N

)

Definition at line 5457 of file hijing1.383.f.

    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

Here is the call graph for this function:

gmre()

function gmre

(

X

)

Definition at line 5441 of file hijing1.383.f.

        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

Here is the caller graph for this function:

gphoton()

function gphoton	(		Y1,
			Y2,
			PT2
	)

Definition at line 5171 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hiboost()

subroutine hiboost

(

)

Definition at line 1054 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

hifun()

subroutine hifun	(		I,
			XMIN,
			XMAX,
		external	FHB
	)

Definition at line 4764 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hijcrs()

subroutine hijcrs

(

)

Definition at line 4838 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

hijcsc()

subroutine hijcsc	(		JP,
			JT
	)

Definition at line 3447 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hijels()

subroutine hijels	(	dimension(5)	PSC1,
		dimension(5)	PSC2
	)

Definition at line 3581 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hijflv()

subroutine hijflv

(

ID

)

Definition at line 4507 of file hijing1.383.f.

    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

Here is the call graph for this function:

hijfrg()

subroutine hijfrg	(		JTP,
			NTP,
			IERROR
	)

Definition at line 1630 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

hijhrd()

subroutine hijhrd	(		JP,
			JT,
			JOUT,
			JFLG,
			IOPJET0
	)

Definition at line 2367 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

hijing()

subroutine hijing	(		BMIN0,
			BMAX0
	)

Definition at line 194 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hijini()

subroutine hijini

(

)

Definition at line 3295 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hijset()

subroutine hijset

(

efrm

)

Definition at line 796 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

hijsft()

subroutine hijsft	(		JP,
			JT,
			JOUT,
			IERROR
	)

Definition at line 3660 of file hijing1.383.f.

    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
.LE..AND..NE.   IF(NFP(JP,5)2NFP(JP,3)0) THEN
        AMP0=ULMASS(NFP(JP,3))
        NFDP=NFP(JP,3)+2*NFP(JP,3)/ABS(NFP(JP,3))
        DPM0=ULMASS(NFDP)
.LE.        IF(DPM00.0) THEN
            NFDP=NFDP-2*NFDP/ABS(NFDP)
            DPM0=ULMASS(NFDP)
        ENDIF
    ENDIF
    AMT0=AMX
    DTM0=AMX
    NFDT=0
.LE..AND..NE.   IF(NFT(JT,5)2NFT(JT,3)0) THEN
        AMT0=ULMASS(NFT(JT,3))
        NFDT=NFT(JT,3)+2*NFT(JT,3)/ABS(NFT(JT,3))
        DTM0=ULMASS(NFDT)
.LE.        IF(DTM00.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
.LT.    IF(SWSNN+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))
.LE.    IF(SWSWPTN) THEN
        PKCMX=0.0
.GT..AND..LE.   ELSE IF(SWSWPTN  SWSWPTD
.AND..EQ..AND..EQ.     &        NPJ(JP)0NTJ(JT)0) THEN
       PKCMX=SQRT(SW/4.0-MAX(AMP0,AMT0)**2)
     &           -SQRT(MAX(PTP02,PTT02))
.GT..AND..LE.   ELSE IF(SWSWPTD  SWSWPTX
.AND..EQ..AND..EQ.     &        NPJ(JP)0NTJ(JT)0) THEN
       PKCMX=SQRT(SW/4.0-MAX(DPM0,DTM0)**2)
     &           -SQRT(MAX(PTP02,PTT02))
.GT.    ELSE IF(SWSWPTX) THEN
       PKCMX=SQRT(SW/4.0-AMX**2)-SQRT(MAX(PTP02,PTT02))
    ENDIF
C       ********maximun PT kick
C*********************************************************
C
.EQ..OR..EQ.    IF(NFP(JP,10)1NFT(JT,10)1) THEN
.GT.        IF(PKC1PKCMX) 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
.GT.        IF(PKC2PKCMX) 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
.NE..AND..LE.   IF(IHPR2(13)0  RLU(0)HIDAT(4)) I_SNG=1
.EQ..OR..EQ..OR.    IF((NFP(JP,5)3 NFT(JT,5)3)
.NE..OR..NE..OR.     &      (NPJ(JP)0NFP(JP,10)0)
.NE..OR..NE.     &      (NTJ(JT)0NFT(JT,10)0)) I_SNG=0
C
C               ********decite whether to have single-diffractive
.EQ.    IF(IHPR2(5)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)
.GT.    IF(PKCHIPR1(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
.EQ.    IF(I_SNG1) 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****************************************
.LT.    IF(SWSNN) THEN
       MISS=MISS+1
.LE.       IF(MISS100) then
          PKC=0.0
          GO TO 30
       ENDIF
.NE.       IF(IHPR2(10)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**********
.EQ.c   IF(IHPR2(15)1) GO TO 500
C
C       ********to have DPM type soft interactions
C
 45 CONTINUE
.GT.    IF(SWSXX+0.001) THEN
.EQ.       IF(I_SNG0) 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.
.EQ..AND..EQ..OR.         IF((NFP(JP,5)3 NFT(JT,5)3)
.NE..OR..NE..OR.     &           (NPJ(JP)0NFP(JP,10)0)
.NE..OR..NE.     &       (NTJ(JT)0NFT(JT,10)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 
.GT..OR..GT..OR.          IF(RLU(0)0.5(NFT(JT,5)2
.NE..OR..NE.     &            NTJ(JT)0NFT(JT,10)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
.GT..AND.   ELSE IF(SWMAX(SPDTX,SPXTD)+0.001 
.LE.     &              SWSXX+0.001) THEN
.EQ..AND..EQ..AND.     IF(((NPJ(JP)0NTJ(JT)0
.GT..OR..EQ.     &         RLU(0)0.5)(NPJ(JP)0
.AND..NE..AND..LE.     &         NTJ(JT)0))NFP(JP,5)2) THEN
          D1=DPD
          D2=DTX
          NFP3=NFDP
          NFT3=0
          GO TO 220
.EQ..AND..LE.      ELSE IF(NTJ(JT)0NFT(JT,5)2) THEN
          D1=DPX
          D2=DTD
          NFP3=0
          NFT3=NFDT
          GO TO 240
       ENDIF
       GO TO 4000
.GT..AND.   ELSE IF(SWMIN(SPDTX,SPXTD)+0.001
.LE.     &          SWMAX(SPDTX,SPXTD)+0.001) THEN
.LE..AND..EQ.      IF(SPDTXSPXTDNPJ(JP)0
.AND..LE.     &                       NFP(JP,5)2) THEN
          D1=DPD
          D2=DTX
          NFP3=NFDP
          NFT3=0
          GO TO 220
.GT..AND..EQ.      ELSE IF(SPDTXSPXTDNTJ(JT)0
.AND..LE.     &                       NFT(JT,5)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
.EQ..AND..EQ.      IF(((NPJ(JP)0NTJ(JT)0
.AND..GT..OR..EQ.     &       RLU(0)0.5)(NPJ(JP)0
.AND..NE..AND..LE.     &        NTJ(JT)0))NFP(JP,5)2) THEN
          D1=DPN
          D2=DTX
          NFP3=NFP(JP,3)
          NFT3=0
          GO TO 220
.EQ..AND..LE.      ELSE IF(NTJ(JT)0NFT(JT,5)2) THEN
          D1=DPX
          D2=DTN
          NFP3=0
          NFT3=NFT(JT,3)
          GO TO 240
       ENDIF
       GO TO 4000
.GT..AND.   ELSE IF(SWMAX(SPNTX,SPXTN)+0.001 
.LE.     &          SWMIN(SPDTX,SPXTD)+0.001) THEN
.EQ..AND..EQ.      IF(((NPJ(JP)0NTJ(JT)0
.AND..GT..OR..EQ.     &       RLU(0)0.5)(NPJ(JP)0
.AND..NE..AND..LE.     &        NTJ(JT)0))NFP(JP,5)2) THEN
          D1=DPN
          D2=DTX
          NFP3=NFP(JP,3)
          NFT3=0
          GO TO 220
.EQ..AND..LE.      ELSE IF(NTJ(JT)0NFT(JT,5)2) THEN
          D1=DPX
          D2=DTN
          NFP3=0
          NFT3=NFT(JT,3)
          GO TO 240
       ENDIF
       GO TO 4000
.GT..AND.   ELSE IF(SWMIN(SPNTX,SPXTN)+0.001 
.LE.     &          SWMAX(SPNTX,SPXTN)+0.001) THEN
.LE..AND..EQ.      IF(SPNTXSPXTNNPJ(JP)0
.AND..LE.     &                           NFP(JP,5)2) THEN
          D1=DPN
          D2=DTX
          NFP3=NFP(JP,3)
          NFT3=0
          GO TO 220
.GT..AND..EQ.      ELSEIF(SPNTXSPXTNNTJ(JT)0
.AND..LE.     &                           NFT(JT,5)2) THEN
          D1=DPX
          D2=DTN
          NFP3=0
          NFT3=NFT(JT,3)
          GO TO 240
       ENDIF
       GO TO 4000
.LE..AND.   ELSE IF(SWMIN(SPNTX,SPXTN)+0.001 
.NE..OR..NE.     &          (NPJ(JP)0 NTJ(JT)0)) THEN
       GO TO 4000
.LE..AND.   ELSE IF(SWMIN(SPNTX,SPXTN)+0.001 
.GT..AND..GT.     &     NFP(JP,5)2NFT(JT,5)2) THEN
       GO TO 4000
.GT..AND..LE.   ELSE IF(SWSDD+0.001SW
     &                     MIN(SPNTX,SPXTN)+0.001) THEN
       D1=DPD
       D2=DTD
       NFP3=NFDP
       NFT3=NFDT
       GO TO 100
.GT.    ELSE IF(SWMAX(SPNTD,SPDTN)+0.001 
.AND..LE.     &                       SWSDD+0.001) THEN
.GT.       IF(RLU(0)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
.GT.    ELSE IF(SWMIN(SPNTD,SPDTN)+0.001
.AND..LE.     &      SWMAX(SPNTD,SPDTN)+0.001) THEN
.GT.       IF(SPNTDSPDTN) 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
.LE.    ELSE IF(SWMIN(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
.LT..OR..LT.    IF(BB1**24.0*D1  BB2**24.0*D2) THEN
        MISS=MISS+1
.GT..OR..EQ.        IF(MISS100PKC0.0) GO TO 3000
        PKC=PKC*0.5
        GO TO 30
    ENDIF
.LT.    IF(RLU(0)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

Here is the call graph for this function:

Here is the caller graph for this function:

hijsrt()

subroutine hijsrt	(		JPJT,
			NPT
	)

Definition at line 1933 of file hijing1.383.f.

    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

Here is the call graph for this function:

hijwds()

subroutine hijwds	(		IA,
			IDH,
			XHIGH
	)

Definition at line 4548 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

hiptdi()

subroutine hiptdi	(		PT,
			PTMAX,
			IOPT
	)

Definition at line 4521 of file hijing1.383.f.

    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

Here is the call graph for this function:

hirnd()

function hirnd

(

I

)

Definition at line 4779 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

hirnd2()

function hirnd2	(		I,
			XMIN,
			XMAX
	)

Definition at line 4807 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

jetini()

subroutine jetini	(		JP,
			JT,
			I_TRIG
	)

Definition at line 3030 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

omg0()

function omg0

(

X

)

Definition at line 4955 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

parton()

subroutine parton	(	dimension(2,7)	F,
			X1,
			X2,
			QQ
	)

Definition at line 5312 of file hijing1.383.f.

    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

Here is the call graph for this function:

Here is the caller graph for this function:

profile()

function profile

(

XB

)

Definition at line 4674 of file hijing1.383.f.

        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

Here is the call graph for this function:

Here is the caller graph for this function:

quench()

subroutine quench	(		JPJT,
			NTP
	)

Definition at line 1089 of file hijing1.383.f.

    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 
.LE.         IF(DP0.2) GO TO 210
.LE.         IF(PTOT0.4) GO TO 290
.LE.         IF(PTOTDP) DP=PTOT-0.2
         DE=DP

.NE.         IF(KFPJ(JP,I)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
.LT.            IF(AMSHUHIPR1(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
.GE.         IF(RN1.0-EXP(-DRR/HIPR1(13))) GO TO 210
         DP=DRR*HIPR1(14)
.LE.         IF(DP0.2) GO TO 210
.LE.         IF(PTOT0.4) GO TO 290
.LE.         IF(PTOTDP) DP=PTOT-0.2
         DE=DP

.NE.         IF(KFPJ(JP,I)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
.LT.            IF(AMSHUHIPR1(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***

.NE. 400          IF(NFT(JPJT,7)1) RETURN
          JT=JPJT
          DO 690 I=1,NTJ(JT)
         PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
.LE.         IF(PTJET0HIPR1(11)) GO TO 690
         PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
.LT.         IF(PTOTHIPR1(8)) GO TO 690
         PHIT=ULANGL(PJTX(JT,I),PJTY(JT,I))
         KP=0
         DO 500 I2=1,IHNT2(1)
.NE.            IF(NFP(I2,5)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)
.GE.                    IF(DPHIHIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
.GT.            IF(DPHIHIPR1(40)/2.0) GO TO 500
            RD0=SQRT(DX*DX+DY*DY)
.GT.            IF(RD0*SIN(DPHI)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
.LT.               IF(RDP(I2)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)
.NE..OR..EQ.               IF(NFT(I2,5)3  I2JT) 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)
.GE.                       IF(DPHIHIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
.GT.               IF(DPHIHIPR1(40)/2.0) GO TO 520
               RD0=SQRT(DX*DX+DY*DY)
.GT.               IF(RD0*SIN(DPHI)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
.LT.              IF(RDT(I2)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)
.GE..AND..GE. 610       IF(MTKT  MPKP) GO TO 690
.GE.        IF(MTKT) GO TO 620
.GE.        IF(MPKP) GO TO 640
.GT.        IF(RDP(MP+1)RDT(MT+1)) GO TO 640
620     MP=MP+1
        DRR=RDP(MP)-R0
.GE.        IF(RN1.0-EXP(-DRR/HIPR1(13))) GO TO 610
        DP=DRR*HIPR1(14)
.NE.        IF(KFTJ(JT,I)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

Here is the call graph for this function:

Here is the caller graph for this function:

romg()

function romg

(

X

)

Definition at line 4969 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

rwdsax()

function rwdsax

(

X

)

Definition at line 4631 of file hijing1.383.f.

        rwdsax=x*x*wdsax(x)
        RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

sgmin()

function sgmin

(

N

)

Definition at line 4942 of file hijing1.383.f.

        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

Here is the caller graph for this function:

subcrs1()

function subcrs1	(		T,
			U
	)

Definition at line 5260 of file hijing1.383.f.

    IMPLICIT REAL*8(a-h,o-z)
    subcrs1=4.d0/9.d0*(1.d0+u**2)/t**2
    RETURN

Here is the caller graph for this function:

subcrs2()

function subcrs2	(		T,
			U
	)

Definition at line 5267 of file hijing1.383.f.

    IMPLICIT REAL*8(a-h,o-z)
    subcrs2=4.d0/9.d0*(t**2+u**2)
    RETURN

Here is the caller graph for this function:

subcrs3()

function subcrs3	(		T,
			U
	)

Definition at line 5274 of file hijing1.383.f.

    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

Here is the caller graph for this function:

subcrs4()

function subcrs4	(		T,
			U
	)

Definition at line 5282 of file hijing1.383.f.

    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

Here is the caller graph for this function:

subcrs5()

function subcrs5	(		T,
			U
	)

Definition at line 5290 of file hijing1.383.f.

    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

Here is the caller graph for this function:

subcrs6()

function subcrs6	(		T,
			U
	)

Definition at line 5297 of file hijing1.383.f.

    IMPLICIT REAL*8(a-h,o-z)
    subcrs6=(1.d0+u**2)*(1.d0/t**2-4.d0/u/9.d0)
    RETURN

Here is the caller graph for this function:

subcrs7()

function subcrs7	(		T,
			U
	)

Definition at line 5304 of file hijing1.383.f.

    IMPLICIT REAL*8(a-h,o-z)
    subcrs7=9.d0/2.d0*(3.d0-t*u-u/t**2-t/u**2)
    RETURN

Here is the caller graph for this function:

title()

subroutine title

(

)

Definition at line 5981 of file hijing1.383.f.

    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

Here is the caller graph for this function:

vegas()

subroutine vegas	(	external	FXN,
			AVGI,
			SD,
			CHI2A
	)

Definition at line 5568 of file hijing1.383.f.

      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

Here is the call graph for this function:

Here is the caller graph for this function:

wdsax()

function wdsax

(

X

)

Definition at line 4619 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

wdsax1()

function wdsax1

(

X

)

Definition at line 4639 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

wdsax2()

function wdsax2

(

X

)

Definition at line 4655 of file hijing1.383.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: