#include "isajet/pilot.h"
      SUBROUTINE JETGEN(J)
C
C          FRAGMENT JET J IN /JETSET/ INTO PRIMARY HADRONS USING THE
C          ALGORITHM OF FIELD AND FEYNMAN WITH
C          F(X)=1-XGEN(1)+XGEN(1)*(XGEN(2)+1)*(1-X)**XGEN(2)
C          FOR LIGHT QUARKS AND THE PETERSON F(X) WITH
C          EPSILON=XGEN(I)*AMASS(I)**2
C          FOR HEAVY QUARKS.
C          INCLUDE BARYONS USING DIQUARKS WITH PROBABILITY PBARY.
C          PROBABILITY PSPIN1 FOR SPIN 1 DEPENDS ON HEAVIEST FLAVOR.
C          FRAGMENT A GLUON LIKE A RANDOM QUARK.
C
C          Ver 7.30: Use delta function fragmentation for top quark.
C
C
#include "isajet/itapes.inc"
#include "isajet/jetset.inc"
#include "isajet/partcl.inc"
#include "isajet/frgpar.inc"
#include "isajet/const.inc"
#include "isajet/mbpar.inc"
C
      LOGICAL HEAVY
      NBEGIN=NPTCL+1
      PSUM=0.
      IFLBEG=JTYPE(J)
      HEAVY=.FALSE.
      IF(IABS(IFLBEG).GT.3.AND.IFLBEG.NE.9) HEAVY=.TRUE.
      PBEG=SQRT(PJSET(1,J)**2+PJSET(2,J)**2+PJSET(3,J)**2)
C          TOP QUARK...
      IF(IABS(IFLBEG).GE.6.AND.IABS(IFLBEG).LE.8) THEN
        NPTCL=NPTCL+1
        PPTCL(1,NPTCL)=0
        PPTCL(2,NPTCL)=0
        PPTCL(3,NPTCL)=PBEG
        PPTCL(4,NPTCL)=PJSET(4,J)
        PPTCL(5,NPTCL)=PJSET(5,J)
        IORIG(NPTCL)=-J
        IDCAY(NPTCL)=0
        IDENT(NPTCL)=JTYPE(J)
        RETURN
      ENDIF
C          EQUIVALENT QUARK FOR GLUON
      IF(IFLBEG.NE.9) GO TO 200
      IFLBEG=INT(RANF()/PUD)+1
      IF(RANF().GT..5) IFLBEG=-IFLBEG
C          CONSTRUCT FIRST QUARK
200   LOOP=0
      IFL1=IFLBEG
      CALL GETPT(PT1,SIGQT)
      PHI1=2.*PI*RANF()
      PX1=PT1*COS(PHI1)
      PY1=PT1*SIN(PHI1)
      PPLUS=PBEG+PJSET(4,J)
      PTRUE=PPLUS
935   CONTINUE
C          CONSTRUCT NEXT QUARK
300   LOOP=LOOP+1
      IF(PPLUS.LT.PEND.OR.LOOP.GT.10000) RETURN
C          IFL2 CAN BE DIQUARK ONLY IF IFL1 IS NOT
      IF(MOD(IFL1,100).EQ.0) GO TO 305
      IF(RANF().LT.PBARY) GO TO 310
      IFL2=ISIGN(INT(RANF()/PUD)+1,-IFL1)
      GO TO 320
305   IFL2=ISIGN(INT(RANF()/PUD)+1,+IFL1)
      GO TO 320
310   IQ1=INT(RANF()/PUD)+1
      IQ2=INT(RANF()/PUD)+1
      IF(IQ1.LE.IQ2) GO TO 315
      ISWAP=IQ1
      IQ1=IQ2
      IQ2=ISWAP
315   IFL2=ISIGN(1000*IQ1+100*IQ2,IFL1)
320   CONTINUE
      CALL GETPT(PT2,SIGQT)
      PHI2=2.*PI*RANF()
      PX2=PT2*COS(PHI2)
      PY2=PT2*SIN(PHI2)
C          CONSTRUCT MESON WITH FLAVOR MIXING
C          SPECIAL CASE - SUPERSYM
      IFLABS=IABS(IFL1)
      IF(IFLABS.GT.20.AND.IFLABS.LT.30) THEN
        IDHAD=IFL1
        GOTO 470
      ENDIF
      IF(MOD(IFL1,100).EQ.0) GO TO 420
      IF(MOD(IFL2,100).EQ.0) GO TO 425
      IHIGH=MAX0(IABS(IFL1),IABS(IFL2))
      JSPIN=INT(RANF()+PSPIN1(IHIGH))
      ID1=IFL1
      ID2=IFL2
      IF(ID1+ID2.NE.0) GO TO 400
      RND=RANF()
      ID1=IABS(ID1)
      ID1=INT(PMIX1(ID1,JSPIN+1)+RND)+INT(PMIX2(ID1,JSPIN+1)+RND)+1
      ID2=-ID1
400   IF(IABS(ID1).LE.IABS(ID2)) GO TO 410
      ISAVE=ID1
      ID1=ID2
      ID2=ISAVE
410   IDHAD=ISIGN(100*IABS(ID1)+10*IABS(ID2)+JSPIN,ID1)
      GO TO 470
C          CONSTRUCT BARYON IDENT.
420   ID3=MOD(IFL1/100,10)
      ID2=IFL1/1000
      ID1=IFL2
      GO TO 430
425   ID3=MOD(IFL2/100,10)
      ID2=IFL2/1000
      ID1=IFL1
430   IF(IABS(ID1).LE.IABS(ID2)) GO TO 431
      ISWAP=ID1
      ID1=ID2
      ID2=ISWAP
431   IF(IABS(ID2).LE.IABS(ID3)) GO TO 432
      ISWAP=ID2
      ID2=ID3
      ID3=ISWAP
432   IF(IABS(ID1).LE.IABS(ID2)) GO TO 440
      ISWAP=ID1
      ID1=ID2
      ID2=ISWAP
440   JSPIN=1
      IF(ID1.EQ.ID2.AND.ID2.EQ.ID3) GO TO 450
      IHIGH=IABS(ID3)
      JSPIN=INT(RANF()+PSPIN1(IHIGH))
450   IF(JSPIN.EQ.1.OR.ID1.EQ.ID2.OR.ID2.EQ.ID3) GO TO 460
      IF(RANF().GT.PISPN) GO TO 460
      ISWAP=ID1
      ID1=ID2
      ID2=ISWAP
460   IDHAD=1000*IABS(ID1)+100*IABS(ID2)+10*IABS(ID3)+JSPIN
      IDHAD=ISIGN(IDHAD,IFL1)
470   CONTINUE
      AM=AMASS(IDHAD)
      PX=PX1+PX2
      PY=PY1+PY2
      AMT2=PX**2+PY**2+AM**2
C          IF LEADING PARTICLE, FIND MINIMUM X
      XMIN=0.
      IF(LOOP.EQ.1) XMIN=AMIN1(SQRT(AMT2)/PPLUS,1.)
C          SELECT X
C          USE FIELD-FEYNMAN FUNCTION FOR LIGHT QUARKS.
C          USE PETERSON FRAGMENTATION FOR HEAVY QUARKS.
C          USE DISTRIBUTION FOR HEAVIER QUARK FOR DIQUARKS.
      II1=IABS(IFL1)
      IF(MOD(II1,100).EQ.0) II1=MOD(II1/100,10)
      IF(II1.LE.3) THEN
        X=RANF()
        IF(RANF().LT.XGEN(1)) X=1.-X**(1./(XGEN(2)+1.))
      ELSEIF(II1.LE.9) THEN
        CALL HEAVYX(X,XGEN(II1)/AM**2)
      ELSEIF(II1.GT.20.AND.II1.LT.30) THEN
        CALL HEAVYX(X,XGENSS(II1-20)/AM**2)
      ENDIF
      X=XMIN+(1.-XMIN)*X
      QPLUS=X*PPLUS
      QPLUS=AMAX1(QPLUS,1.E-6)
      QMINUS=AMT2/QPLUS
      P0=.5*(QPLUS+QMINUS)
      PZ=.5*(QPLUS-QMINUS)
C          DISCARD PARTICLE IF PZ<0
      IF(PZ.LT.0..AND..NOT.(HEAVY.AND.LOOP.EQ.1)) GO TO 500
C          ADD PARTICLE TO /PARTCL/
      IF(NPTCL.GE.MXPTCL) GO TO 9999
      NPTCL=NPTCL+1
      PPTCL(1,NPTCL)=PX
      PPTCL(2,NPTCL)=PY
      PPTCL(3,NPTCL)=PZ
      PPTCL(4,NPTCL)=P0
      PPTCL(5,NPTCL)=AM
      IORIG(NPTCL)=-J
      IDCAY(NPTCL)=0
      IDENT(NPTCL)=IDHAD
      PSUM=PSUM+QPLUS
C          SWAP QUARKS AND CONTINUE IF SUFFICIENT PPLUS
500   CONTINUE
      PX1=-PX2
      PY1=-PY2
      IFL1=-IFL2
      PPLUS=(1.-X)*PPLUS
      GO TO 300
C
9999  CALL PRTEVT(0)
      WRITE(ITLIS,10) NPTCL
10    FORMAT(//5X,'ERROR IN JETGEN...NPTCL >',I5)
      RETURN
      END