C####################################################################
C
      SUBROUTINE MBSET 
C         
C          SET PARAMETERS FOR GENERATING MINBIAS EVENTS OR BEAM JETS, 
C          ALLOWING DIFFERENT PARAMETERS FOR TWO CASES.     
C         
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS       
      COMMON/MBPAR/PUD0,PJSPN,PISPN,SIGQT0,XGEN0(2),PMIX01(3,2)       
     1,PMIX02(3,2),PBARY0     
      PARAMETER (LIMPOM=20)   
      COMMON/MBGEN/POMWT(LIMPOM),POMGEN(LIMPOM),MNPOM,MXPOM,PDIFFR,   
     $NPOM,XBARY(2),DXBARY(2),XPOM(LIMPOM,2)      
      COMMON/PRIMAR/NJET,SCM,HALFE,ECM,IDIN(2),NEVENT,NTRIES,NSIGMA   
      COMMON/TOTALS/NKINPT,NWGEN,NKEEP,SUMWT,WT   
      COMMON/KEYS/IKEYS,KEYON,KEYS(10)  
      COMMON/XKEYS/REAC       
      LOGICAL KEYS  
      LOGICAL KEYON 
      CHARACTER*8 REAC        
C         
C         
C          DN/DY INCREASES WITH LOG(S). INCLUDED IN SPLITTING FUNCTION
C          BECAUSE AVERAGE MULTIPLICITY COMES FROM SINGLE CHAIN GRAPH.
      XGEN0(1)=.9   
      XGEN0(2)=1.+0.35*ALOG(ECM/60.)    
C         
C          POMWT ARE (RELATIVE) PROBABILITIES FOR N CUT POMERONS.     
C          PDIFFR IS DIFFRACTIVE PROBABILITY.     
C          SIGQT0 IS MEAN PT. 
      IF(KEYS(4)) THEN        
          PDIFFR=.33
          SIGQT0=.35  
          PSUM=0.     
          DO 100 I=1,LIMPOM     
          POMWT(I)=(1.+4.*I**2)*EXP(-1.8*I)         
          PSUM=PSUM+POMWT(I)    
  100     CONTINUE
      ELSE
          PDIFFR=0.   
          SIGQT0=.45  
          PSUM=0.     
          DO 110 I=1,LIMPOM     
          POMWT(I)=(1.+4.*I**2)*EXP(-1.8*I)         
          PSUM=PSUM+POMWT(I)    
  110     CONTINUE    
          POMWT(1)=.1*POMWT(1)  
          POMWT(2)=.2*POMWT(2)  
          POMWT(3)=.5*POMWT(3)  
      ENDIF         
C         
C          RENORMALIZE POMWT. 
      PSUM=1./PSUM  
      DO 200 I=1,LIMPOM       
      POMWT(I)=PSUM*POMWT(I)  
  200 CONTINUE      
      PSUM=0.       
      DO 210 I=MNPOM,MXPOM    
      PSUM=PSUM+POMWT(I)      
  210 CONTINUE      
C         
C          POMGEN IS USED TO SELECT NUMBER OF POMERONS.     
      PGEN=0.       
      PSUM=1./PSUM  
      DO 300 I=1,LIMPOM       
      POMGEN(I)=0.  
  300 CONTINUE      
      DO 310 I=MNPOM,MXPOM    
      PGEN=PGEN+PSUM*POMWT(I) 
      POMGEN(I)=PGEN
  310 CONTINUE      
      POMGEN(MXPOM)=1.        
C         
C          SET /TOTALS/ FOR MINBIAS EVENTS USING LOG**2(S) FIT TO     
C          TOTAL CROSS SECTION.         
      IF(KEYS(4)) THEN        
           SIGTOT=25.65*(1.+.0102*ALOG(SCM/1.76)**2) 
           SIGTOT=PSUM*SIGTOT    
           NKINPT=NEVENT         
           SUMWT=SIGTOT*NKINPT   
      ENDIF         
C         
      RETURN        
      END
        SUBROUTINE MOVLEV(A,B,N)
        DIMENSION A(N),B(N)
        DO 10 I=1,N
        B(I)=A(I)
10      CONTINUE
        RETURN
        END
      SUBROUTINE DECAY(IP)  
C          THIS SUBROUTINE DECAYS PARTICLE IP FROM /PARTCL/ USING THE   
C          BRANCHING RATIOS FROM /DKYTAB/ AND ADDS THE DECAY PRODUCTS   
C          TO /PARTCL/ WITH IORIG=IP.   
C          QUARK-BASED IDENT CODE   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/WCON/SIN2W,WMASS(4),WGAM(4),AQ(12,4),BQ(12,4),COUT(4), 
     1MATCH(25,4),WCBR(25,4),CUTOFF,CUTPOW,TBRWW(4,2),RBRWW(12,4,2),EZ, 
     2AQDP(12,4),BQDP(12,4),EZDP    
      DOUBLE PRECISION AQDP,BQDP,EZDP   
      PARAMETER (MXPTCL=4000,IPACK=10000)   
      COMMON/PARTCL/NPTCL,PPTCL(5,MXPTCL),IORIG(MXPTCL),IDENT(MXPTCL)   
     1,IDCAY(MXPTCL)    
C          LOOK MUST BE DIMENSIONED TO THE MAXIMUM VALUE OF INDEX.  
      PARAMETER (MXLOOK=400)    
      PARAMETER (MXDKY=1000)    
      COMMON/DKYTAB/LOOK(MXLOOK),CBR(MXDKY),MODE(5,MXDKY)   
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      COMMON/JWORK/ZZC(MXJSET),JMATCH(MXJSET),TNEW,P1CM(4), 
     1J1,J2,J3,J4,J5,E1CM,E2CM,E3CM,E4CM,E5CM   
      LOGICAL TNEW  
      DIMENSION JJ(5),EE(5) 
      EQUIVALENCE (J1,JJ(1)),(E1CM,EE(1))   
      COMMON/CONST/PI,SQRT2,ALFA,GF,UNITS   
      DIMENSION PGEN(5,5),RND(5),U(3),BETA(3),IDQK(3),ROT(3,3),PSAVE(3) 
     1,REDUCE(5)    
      DIMENSION PSUM(5),POLD(4),PNEW(4) 
      DIMENSION PREST(4,5)  
      DATA REDUCE/1.,1.,2.,5.,15./  
      DATA PSUM/5*0./   
      DATA TWOME/1.022006E-3/   
C          FUNCTION DEFINITIONS 
C          USE DOUBLE PRECISION FOR PCM ON 32-BIT MACHINES  
      PCM(A,B,C)=DBLPCM(A,B,C)  
      DOT(I1,I2)=PREST(4,I1)*PREST(4,I2)-PREST(1,I1)*PREST(1,I2)    
     $-PREST(2,I1)*PREST(2,I2)-PREST(3,I1)*PREST(3,I2)  
C          CHARGED W PROPAGATOR.    
      WPROP(Z)=(Z-WMASS(2)**2)**2+(WMASS(2)*WGAM(2))**2 
C          SELECT DECAY MODE    
      IDLV1=IDENT(IP)   
      CALL FLAVOR(IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX) 
1     CONTINUE  
      IPOINT=LOOK(INDEX)-1  
      IF(IPOINT.LT.0) RETURN    
      TRY=RANF()    
100   IPOINT=IPOINT+1   
      IF(TRY.GT.CBR(IPOINT)) GO TO 100  
      NADD=0    
      SUM=0.    
      NSTART=NPTCL+1    
      DO 110 I=1,5  
      IF(MODE(I,IPOINT).EQ.0) GO TO 110 
      IF(NPTCL+NADD+1.GT.MXPTCL) GO TO 9999 
      NADD=NADD+1   
      NEW=NPTCL+NADD    
      IDENT(NEW)=MODE(I,IPOINT) 
      IDLV1=IDENT(NEW)  
      PPTCL(5,NEW)=AMASS(IDLV1) 
      SUM=SUM+PPTCL(5,NEW)  
110   CONTINUE  
      NADD1=NADD-1  
      DO 120 J=1,5  
      PGEN(J,1)=PPTCL(J,IP) 
120   CONTINUE  
      PGEN(5,NADD)=PPTCL(5,NPTCL+NADD)  
      IF(NADD.EQ.1) GO TO 700   
      IF(NADD.EQ.2) GO TO 400   
C          USE KROLL-WADA DISTRIBUTION FOR DALITZ DECAYS.   
      IF(.NOT.((IDENT(IP).EQ.110.OR.IDENT(IP).EQ.220).AND.  
     1IABS(IDENT(NPTCL+2)).EQ.12)) GO TO 130    
125   AMEE=TWOME*(PPTCL(5,IP)/TWOME)**RANF()    
      REE=(TWOME/AMEE)**2   
      WTEE=(1.-(AMEE/PPTCL(5,IP))**2)**3*SQRT(1.-REE)*(1.+.5*REE)   
      IF(WTEE.LT.RANF()) GO TO 125  
      PGEN(5,2)=AMEE    
      GO TO 400 
130   CONTINUE  
C          CALCULATE MAXIMUM PHASE-SPACE WEIGHT 
      WTMAX=1./REDUCE(NADD) 
      SUM1=PGEN(5,1)    
      SUM2=SUM-PPTCL(5,NPTCL+1) 
      DO 200 I=1,NADD1  
      WTMAX=WTMAX*PCM(SUM1,SUM2,PPTCL(5,NPTCL+I))   
      SUM1=SUM1-PPTCL(5,NPTCL+I)    
      SUM2=SUM2-PPTCL(5,NPTCL+I+1)  
200   CONTINUE  
C          GENERATE UNIFORM NADD-BODY PHASE SPACE   
300   CONTINUE  
      RND(1)=1. 
      DO 310 I=2,NADD1  
      RNEW=RANF()   
      I1=I-1    
      DO 320 JJ1=1,I1   
      J=I-JJ1   
      JSAVE=J+1 
      IF(RNEW.LE.RND(J)) GO TO 310  
      RND(JSAVE)=RND(J) 
320   CONTINUE  
310   RND(JSAVE)=RNEW   
      RND(NADD)=0.  
      WT=1. 
      SUM1=SUM  
      DO 330 I=2,NADD   
      SUM1=SUM1-PPTCL(5,NPTCL+I-1)  
      PGEN(5,I)=SUM1+RND(I)*(PGEN(5,1)-SUM) 
      WT=WT*PCM(PGEN(5,I-1),PGEN(5,I),PPTCL(5,NPTCL+I-1))   
330   CONTINUE  
      IF(WT.LT.RANF()*WTMAX) GO TO 300  
C          CARRY OUT TWO-BODY DECAYS IN PGEN FRAMES 
400   CONTINUE  
      DO 410 I=1,NADD1  
      QCM=PCM(PGEN(5,I),PGEN(5,I+1),PPTCL(5,NPTCL+I))   
      U(3)=2.*RANF()-1. 
      PHI=2.*PI*RANF()  
      U(1)=SQRT(1.-U(3)**2)*COS(PHI)    
      U(2)=SQRT(1.-U(3)**2)*SIN(PHI)    
      DO 420 J=1,3  
      PPTCL(J,NPTCL+I)=QCM*U(J) 
      PGEN(J,I+1)=-PPTCL(J,NPTCL+I) 
420   CONTINUE  
      PPTCL(4,NPTCL+I)=SQRT(QCM**2+PPTCL(5,NPTCL+I)**2) 
      PGEN(4,I+1)=SQRT(QCM**2+PGEN(5,I+1)**2)   
410   CONTINUE  
      DO 430 J=1,4  
      PPTCL(J,NPTCL+NADD)=PGEN(J,NADD)  
430   CONTINUE  
C          BOOST PGEN FRAMES TO LAB FRAME   
C          ALSO SAVE MOMENTA IN REST FRAME (LAST FRAME) 
      DO 500 II=1,NADD1 
      I=NADD-II 
      DO 510 J=1,3  
      BETA(J)=PGEN(J,I)/PGEN(4,I)   
510   CONTINUE  
      GAMMA=PGEN(4,I)/PGEN(5,I) 
      DO 520 K=I,NADD   
      K1=NPTCL+K    
      BP=BETA(1)*PPTCL(1,K1)+BETA(2)*PPTCL(2,K1)+BETA(3)*PPTCL(3,K1)    
      DO 530 J=1,3  
      PREST(J,K)=PPTCL(J,K1)    
      PPTCL(J,K1)=PPTCL(J,K1)+GAMMA*BETA(J)*(PPTCL(4,K1)    
     1+BP*GAMMA/(GAMMA+1.)) 
530   CONTINUE  
      PREST(4,K)=PPTCL(4,K1)    
      PPTCL(4,K1)=GAMMA*(PPTCL(4,K1)+BP)    
520   CONTINUE  
500   CONTINUE  
C          MATRIX ELEMENTS  
      IF(NADD.EQ.3.AND.(IDENT(IP).EQ.221.OR.IDENT(IP).EQ.331)) GO TO 610    
      IF(NADD.EQ.3.AND.IABS(IDENT(NPTCL+1)).LT.20.AND.  
     1IDENT(NPTCL+1).NE.10) GO TO 620   
      GO TO 800 
C          OMEG AND PHI DECAY   
C          USE VECTORS IN REST FRAME    
610   WT=(PPTCL(5,NPTCL+1)*PPTCL(5,NPTCL+2)*PPTCL(5,NPTCL+3))**2    
     1-(PPTCL(5,NPTCL+1)*DOT(2,3))**2   
     2-(PPTCL(5,NPTCL+2)*DOT(1,3))**2   
     3-(PPTCL(5,NPTCL+3)*DOT(1,2))**2   
     4+2.*DOT(1,2)*DOT(2,3)*DOT(1,3)    
      IF(WT.LT.RANF()*PPTCL(5,IP)**6/108.) GO TO 300    
      GO TO 800 
C          SEMILEPTONIC AND QUARK DECAYS    
C          USE VECTORS IN REST FRAME, WHERE IP HAS (M,0,0,0)    
C          INCLUDE W PROPAGATOR 
620   WT=(PPTCL(5,IP)*PREST(4,2))*DOT(1,3)  
      S12=PPTCL(5,NPTCL+1)**2+PPTCL(5,NPTCL+2)**2+2.*DOT(1,2)   
      S12MAX=PPTCL(5,IP)**2 
      WT=WT*WPROP(S12MAX)/WPROP(S12)    
      IF(WT.LT.RANF()*PPTCL(5,IP)**4/16.) GO TO 300 
      GO TO 800 
C          ONE-PARTICLE DECAYS  
700   CONTINUE  
      DO 710 J=1,5  
      PPTCL(J,NPTCL+1)=PPTCL(J,IP)  
710   CONTINUE  
C          SWAP PARTICLES AND ANTIPARTICLES IF IDENT(IP)<0  
800   CONTINUE  
      IF(IDENT(IP).GE.0.OR.IABS(IDENT(IP)).EQ.20) GO TO 900 
      DO 810 I=1,NADD   
      IDABS=IABS(IDENT(NPTCL+I))    
      IFL1=IDABS/1000   
      IFL2=MOD(IDABS/100,10)    
      IFL3=MOD(IDABS/10,10) 
      IF(IFL1.EQ.0.AND.IFL2.NE.0.AND.IFL2.EQ.IFL3) GO TO 810    
      IF(IDABS.EQ.9.OR.IDABS.EQ.10.OR.IDABS.EQ.20) GO TO 810    
      IF(IDABS.EQ.29.OR.IDABS.EQ.30.OR.IDABS.EQ.40) GO TO 810   
      IDENT(NPTCL+I)=-IDENT(NPTCL+I)    
810   CONTINUE  
C          REMOVE QUARKS FROM /PARTCL/ AND TRANSFORM BACK TO REST FRAME 
900   CONTINUE  
      NPTCL=NPTCL+NADD  
      NQK=0 
      IF(IABS(IDENT(NPTCL)).GE.10.AND.MOD(IDENT(NPTCL),100).NE.0)   
     1GO TO 1000    
      NOFF=NPTCL-NSTART+1   
      DO 910 II=1,NOFF  
      I=NPTCL+1-II  
      IF(IABS(IDENT(NPTCL)).LT.10.OR.MOD(IDENT(NPTCL),100).NE.0)    
     1NQK=NQK+1 
      IREST=NADD+1-II   
      DO 911 J=1,4  
911   PPTCL(J,I)=PREST(J,IREST) 
910   CONTINUE  
C          COPY DECAY PRODUCTS INTO /JETSET/    
      IF(NJSET+NADD.GT.MXJSET) GO TO 9998   
      NJSAVE=NJSET  
      NPTCL=NPTCL-NADD  
      DO 920 I=1,NADD   
      NJSET=NJSET+1 
      DO 921 K=1,5  
921   PJSET(K,NJSET)=PPTCL(K,NPTCL+I)   
      JORIG(NJSET)=0    
      JTYPE(NJSET)=IDENT(NPTCL+I)   
      JDCAY(NJSET)=0    
      JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD    
C          QCD EVOLUTION STARTS FROM PARENT MASS    
C          BUT USE NADD*ENERGY TO PRESERVE TP --> W+ BT 
      IF(IABS(JTYPE(NJSET)).GE.10.AND.MOD(JTYPE(NJSET),100).NE.0)   
     1GO TO 920 
      JDCAY(NJSET)=-1   
      PJSET(5,NJSET)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,NJSET)) 
920   CONTINUE  
C          PERFORM QCD JET EVOLUTION    
      CALL QCDJET(NJSAVE+1) 
C          DECAY QUARKS AND ROTATE TO PROPER ANGLES 
C          HADRONIZE JETS   
      NJ1=NJSAVE+1  
      DO 931 I=NJ1,NJSET    
      IF(JDCAY(I).NE.0) GO TO 931   
      IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0)   
     1GO TO 935 
      NEXT=NPTCL+1  
      PJET=SQRT(PJSET(1,I)**2+PJSET(2,I)**2+PJSET(3,I)**2)  
      CTHQK=PJSET(3,I)/PJET 
      STHQK=SQRT(1.-CTHQK**2)   
      CPHIQK=PJSET(1,I)/(PJET*STHQK)    
      SPHIQK=PJSET(2,I)/(PJET*STHQK)    
      CALL JETGEN(I)    
      IF(NEXT.GT.NPTCL) GO TO 931   
      ROT(1,1)=CPHIQK*CTHQK 
      ROT(2,1)=SPHIQK*CTHQK 
      ROT(3,1)=-STHQK   
      ROT(1,2)=-SPHIQK  
      ROT(2,2)=CPHIQK   
      ROT(3,2)=0.   
      ROT(1,3)=CPHIQK*STHQK 
      ROT(2,3)=SPHIQK*STHQK 
      ROT(3,3)=CTHQK    
      DO 932 K=NEXT,NPTCL   
      DO 933 J=1,3  
      PSAVE(J)=PPTCL(J,K)   
      PPTCL(J,K)=0. 
933   CONTINUE  
      DO 932 J=1,3  
      DO 932 JJ1=1,3    
      PPTCL(J,K)=PPTCL(J,K)+ROT(J,JJ1)*PSAVE(JJ1)   
932   CONTINUE  
      GOTO 931  
C          ADD LEPTON TO /PARTCL/   
  935 NPTCL=NPTCL+1 
      DO 936 K=1,5  
  936 PPTCL(K,NPTCL)=PJSET(K,I) 
      IDENT(NPTCL)=JTYPE(I) 
931   CONTINUE  
C          RESET NJSET SO DECAY JETS DO NOT APPEAR IN /JETSET/  
      NJSET=NJSAVE  
C          CHECK FOR AT LEAST TWO PARTICLES 
      IF(NPTCL.GT.NSTART) GO TO 939 
      NPTCL=NSTART-1    
      GO TO 1   
939   CONTINUE  
C          CONSERVE CHARGE  
      SUMQ=0.   
      DO 960 I=NSTART,NPTCL 
      IDLV1=IDENT(I)    
      SUMQ=SUMQ+CHARGE(IDLV1)   
960   CONTINUE  
      IDLV1=IDENT(IP)   
      SUMQ=SUMQ-CHARGE(IDLV1)   
      IF(SUMQ.EQ.0.) GO TO 970  
      DO 961 I=NSTART,NPTCL 
      ID1=IDENT(I)  
      IF(IABS(ID1).GT.1000) GO TO 961   
      I1=MOD(IABS(ID1)/100,10)  
      I2=MOD(IABS(ID1)/10,10)   
      I3=MOD(IABS(ID1),10)  
      IF(I1.EQ.1.AND.I2.GT.2.AND.SUMQ*ID1.GT.0.) GO TO 962  
      IF(I1.EQ.2.AND.I2.GT.2.AND.SUMQ*ID1.LT.0.) GO TO 963  
      IF(I1.EQ.1.AND.I2.EQ.2.AND.SUMQ*ID1.GT.0.) GO TO 964  
      IF(I1.EQ.1.AND.I2.EQ.1) GO TO 965 
      GO TO 961 
962   IDENT(I)=ISIGN(200+10*I2+I3,ID1)  
      GO TO 969 
963   IDENT(I)=ISIGN(100+10*I2+I3,ID1)  
      GO TO 969 
964   IDENT(I)=110+I3   
      GO TO 969 
965   IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ))    
969   SUMQ=SIGN(ABS(SUMQ)-1.,SUMQ)  
      IDLV1=IDENT(I)    
      PPTCL(5,I)=AMASS(IDLV1)   
      PPTCL(4,I)=SQRT(PPTCL(1,I)**2+PPTCL(2,I)**2+PPTCL(3,I)**2 
     1+PPTCL(5,I)**2)   
C          SUM CANNOT VANISH FOR FRACTIONALLY CHARGED INITIAL PARTICLE. 
      IF(ABS(SUMQ).LT.1.) GO TO 970 
961   CONTINUE  
      NPTCL=NSTART-1    
      GO TO 1   
C          RESCALE MOMENTA FOR CORRECT MASS 
970   CONTINUE  
      PSUM(4)=PPTCL(5,IP)   
      PSUM(5)=PSUM(4)   
      NPTLV1=NPTCL  
      CALL RESCAL(NSTART,NPTLV1,PSUM,IFAIL) 
      IF(IFAIL.EQ.0) GO TO 940  
      NPTCL=NSTART-1    
      GO TO 1   
940   CONTINUE  
C          BOOST BACK TO LAB FRAME  
      DO 950 I=NSTART,NPTCL 
      BP=BETA(1)*PPTCL(1,I)+BETA(2)*PPTCL(2,I)+BETA(3)*PPTCL(3,I)   
      DO 951 J=1,3  
      PPTCL(J,I)=PPTCL(J,I)+GAMMA*BETA(J)*(PPTCL(4,I)   
     1+BP*GAMMA/(GAMMA+1.)) 
951   CONTINUE  
      PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP)  
950   CONTINUE  
C          SET IORIG AND IDCAY  
1000  CONTINUE  
      IDCAY(IP)=IPACK*NSTART+NPTCL  
      DO 1010 I=NSTART,NPTCL    
      IORIG(I)=IP   
      IDCAY(I)=0    
1010  CONTINUE  
      RETURN    
9999  CALL PRTEVT(0)    
      WRITE(ITLIS,10) NPTCL 
10    FORMAT(//5X,'ERROR IN DECAY...NPTCL > ',I5)   
      RETURN    
9998  CALL PRTEVT(0)    
      WRITE(ITLIS,20) NJSET 
20    FORMAT(//5X,'ERROR IN DECAY...NJSET > ',I5)   
      RETURN    
      END   
      SUBROUTINE FLAVOR(ID,IFL1,IFL2,IFL3,JSPIN,INDEX)  
C   
C          THIS SUBROUTINE UNPACKS THE IDENT CODE ID=+/-IJKL    
C   
C          MESONS-- 
C          I=0, J<=K, +/- IS SIGN FOR J 
C          ID=110 FOR PI0, ID=220 FOR ETA, ETC. 
C   
C          BARYONS--    
C          I<=J<=K IN GENERAL   
C          J<I<K FOR SECOND STATE ANTISYMMETRIC IN (I,J), EG. L = 2130  
C   
C          OTHER--  
C          ID=1,...,6 FOR QUARKS    
C          ID=9 FOR GLUON   
C          ID=10 FOR PHOTON 
C          ID=11,...,16 FOR LEPTONS 
C          ID=20 FOR KS, ID=-20 FOR KL  
C   
C          I=21...26 FOR SCALAR QUARKS  
C          I=29 FOR GLUINO  
C          I=30 FOR PHOTINO 
C          I=31...36 FOR SCALAR LEPTONS 
C          I=39 FOR WINO    
C          I=40 FOR ZINO    
C   
C          ID=80 FOR W+ 
C          ID=81,...,89 FOR HIGGS MESONS    
C          ID=90 FOR Z0 
C   
C          DIQUARKS--   
C          ID=+/-IJ00, I<J FOR DIQUARK COMPOSED OF I,J. 
C   
C          INDEX IS A SEQUENCE NUMBER USED INTERNALLY   
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/QLMASS/AMLEP(52),NQLEP,NMES,NBARY  
      IDABS=IABS(ID)    
      I=IDABS/1000  
      J=MOD(IDABS/100,10)   
      K=MOD(IDABS/10,10)    
      JSPIN=MOD(IDABS,10)   
      IF(ID.NE.0.AND.MOD(ID,100).EQ.0) GO TO 300    
      IF(J.EQ.0) GO TO 200  
      IF(I.EQ.0) GO TO 100  
C          BARYONS  
C          ONLY X,Y BARYONS ARE QQX, QQY, Q=U,D,S.  
      IFL1=ISIGN(I,ID)  
      IFL2=ISIGN(J,ID)  
      IFL3=ISIGN(K,ID)  
      IF(K.LE.6) THEN   
        INDEX=MAX0(I-1,J-1)**2+I+MAX0(I-J,0)+(K-1)*K*(2*K-1)/6  
     1  +109*JSPIN+36*NMES+NQLEP+11 
      ELSE  
        INDEX=MAX0(I-1,J-1)**2+I+MAX0(I-J,0)+9*(K-7)+91 
     1  +109*JSPIN+36*NMES+NQLEP+11 
      ENDIF 
      RETURN    
C          MESONS   
100   CONTINUE  
      IFL1=0    
      IFL2=ISIGN(J,ID)  
      IFL3=ISIGN(K,-ID) 
      INDEX=J+K*(K-1)/2+36*JSPIN+NQLEP  
      INDEX=INDEX+11    
      RETURN    
C          QUARKS, LEPTONS, ETC 
200   CONTINUE  
      IFL1=0    
      IFL2=0    
      IFL3=0    
      JSPIN=0   
      INDEX=IDABS   
      IF(IDABS.LT.20) RETURN    
C          DEFINE INDEX=20 FOR KS, INDEX=21 FOR KL  
      INDEX=IDABS+1 
      IF(ID.EQ.20) INDEX=20 
C          INDEX=NQLEP+1,...,NQLEP+11 FOR W+, HIGGS, Z0 
      IF(IDABS.LT.80) RETURN    
      INDEX=NQLEP+IDABS-79  
      RETURN    
300   IFL1=ISIGN(I,ID)  
      IFL2=ISIGN(J,ID)  
      IFL3=0    
      JSPIN=0   
      INDEX=0   
      RETURN    
      END   
      FUNCTION LABEL(ID)    
C   
C          RETURN THE CHARACTER*8 LABEL FOR THE PARTICLE ID.    
C          QUARK-BASED IDENT CODE.  
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/QLMASS/AMLEP(52),NQLEP,NMES,NBARY  
C   
      CHARACTER*8 LABEL 
      CHARACTER*8 LLEP,LMES0,LMES1,LBAR0,LABAR0,LBAR1,LABAR1    
      CHARACTER*8 LQQ,LAQQ  
      DIMENSION LLEP(104)   
      DIMENSION LMES0(64),LMES1(64) 
      DIMENSION LBAR0(109),LABAR0(109),LBAR1(109),LABAR1(109)   
      DIMENSION LQQ(21),LAQQ(21)    
C          DIQUARK LABELS   
      DATA LQQ/ 
     1'UU0. ','UD0. ','DD0. ','US0. ','DS0. ','SS0. ','UC0. ','DC0. ',  
     2'SC0. ','CC0. ','UB0. ','DB0. ','SB0. ','CB0. ','BB0. ','UT0. ',  
     3'DT0. ','ST0. ','CT0. ','BT0. ','TT0. '/  
      DATA LAQQ/    
     1'AUU0.','AUD0.','ADD0.','AUS0.','ADS0.','ASS0.','AUC0.','ADC0.',  
     2'ASC0.','ACC0.','AUB0.','ADB0.','ASB0.','ACB0.','ABB0.','AUT0.',  
     3'ADT0.','AST0.','ACT0.','ABT0.','ATT0.'/  
C          QUARK AND LEPTON LABELS  
      DATA LLEP/    
     $'     ','UP   ','UB   ','DN   ','DB   ','ST   ','SB   ','CH   ',  
     $'CB   ','BT   ','BB   ','TP   ','TB   ','Y    ','YB   ','X    ',  
     $'XB   ','GL   ','ERR  ','GM   ','ERR  ','NUE  ','ANUE ','E-   ',  
     $'E+   ','NUM  ','ANUM ','MU-  ','MU+  ','NUT  ','ANUT ','TAU- ',  
     $'TAU+ ','ERR  ','ERR  ','ERR  ','ERR  ','ERR  ','ERR  ','KS   ',  
     $'ERR  ','ERR  ','KL   ',  
     $'UPSS ','UBSS ','DNSS ','DBSS ','STSS ','SBSS ','CHSS ','CBSS ',  
     $'BTSS ','BBSS ','TPSS ','TBSS ','ERR  ','ERR  ','ERR  ','ERR  ',  
     $'GLSS ','ERR  ','GMSS ','ERR  ','NESS ','ANESS','E-SS ','E+SS ',  
     $'NMSS ','ANMSS','MU-SS','MU+SS','NTSS ','ANTSS','T-SS ','T+SS ',  
     $'ERR  ','ERR  ','ERR  ','ERR  ','W+SS ','W-SS ','Z0SS ','ERR  ',  
     $'W+   ','W-   ','H10  ','AH10 ','H20  ','AH20 ','H30  ','AH30 ',  
     $'H4+  ','H4-  ','H5+  ','H5-  ','H6+  ','H6-  ','H7++ ','H7-- ',  
     $'H8++ ','H8-- ','H9++ ','H9-- ','Z0   '/  
C          0- MESON LABELS  
      DATA LMES0/   
     1'PI0  ','PI+  ','ETA  ','PI-  ','K+   ','K0   ','ETAP ','AK0  ',  
     2'K-   ','AD0  ','D-   ','F-   ','ETAC ','F+   ','D+   ','D0   ',  
     2'UB.  ','DB.  ','SB.  ','CB.  ','BB.  ','BC.  ','BS.  ','BD.  ',  
     3'BU.  ','UT.  ','DT.  ','ST.  ','CT.  ','BT.  ','TT.  ','TB.  ',  
     4'TC.  ','TS.  ','TD.  ','TU.  ','UY.  ','DY.  ','SY.  ','CY.  ',  
     5'BY.  ','TY.  ','YY.  ','YT.  ','YB.  ','YC.  ','YS.  ','YD.  ',  
     6'YU.  ','UX.  ','DX.  ','SX.  ','CX.  ','BX.  ','TX.  ','YX.  ',  
     7'XX.  ','XY.  ','XT.  ','XB.  ','XC.  ','XS.  ','XD.  ','XU.  '/  
C          1- MESON LABELS  
      DATA LMES1/   
     1'RHO0 ','RHO+ ','OMEG ','RHO- ','K*+  ','K*0  ','PHI  ','AK*0 ',  
     2'K*-  ','AD*0 ','D*-  ','F*-  ','JPSI ','F*+  ','D*+  ','D*0  ',  
     3'UB*  ','DB*  ','SB*  ','CB*  ','UPSL ','BC*  ','BS*  ','BD*  ',  
     4'BU*  ','UT*  ','DT*  ','ST*  ','CT*  ','BT*  ','TT*  ','TB*  ',  
     5'TC*  ','TS*  ','TD*  ','TU*  ','UY*  ','DY*  ','SY*  ','CY*  ',  
     6'BY*  ','TY*  ','YY*  ','YT*  ','YB*  ','YC*  ','YS*  ','YD*  ',  
     7'YU*  ','UX*  ','DX*  ','SX*  ','CX*  ','BX*  ','TX*  ','YX*  ',  
     8'XX*  ','XY*  ','XT*  ','XB*  ','XC*  ','XS*  ','XD*  ','XU*  '/  
C          1/2+ BARYON LABELS   
      DATA LBAR0/   
     1'ERR  ','P    ','N    ','ERR  ','ERR  ','S+   ','S0   ','S-   ',  
     2'L    ','XI0  ','XI-  ','ERR  ','ERR  ','ERR  ','SC++ ','SC+  ',  
     3'SC0  ','LC+  ','USC. ','DSC. ','SSC. ','SDC. ','SUC. ','UCC. ',  
     4'DCC. ','SCC. ','ERR  ','ERR  ','ERR  ','ERR  ','UUB. ','UDB. ',  
     5'DDB. ','DUB. ','USB. ','DSB. ','SSB. ','SDB. ','SUB. ','UCB. ',  
     6'DCB. ','SCB. ','CCB. ','CSB. ','CDB. ','CUB. ','UBB. ','DBB. ',  
     7'SBB. ','CBB. ','ERR  ','ERR  ','ERR  ','ERR  ','ERR  ','UTT. ',  
     8'UDT. ','DDT. ','DUT. ','UST. ','DST. ','SST. ','SDT. ','SUT. ',  
     9'UCT. ','DCT. ','SCT. ','CCT. ','CST. ','CDT. ','CUT. ','UBT. ',  
     1'DBT. ','SBT. ','CBT. ','BBT. ','BCT. ','BST. ','BDT. ','BUT. ',  
     2'UTT. ','DTT. ','STT. ','CTT. ','BTT. ','ERR  ','ERR  ','ERR  ',  
     3'ERR  ','ERR  ','ERR  ','UUY. ','UDY. ','DDY. ','DUY. ','USY. ',  
     4'DSY. ','SSY. ','SDY. ','SUY. ','UUX. ','UDX. ','DDX. ','DUX. ',  
     5'USX. ','DSX. ','SSX. ','SDX. ','SUX. '/  
      DATA LABAR0/  
     1'ERR  ','AP   ','AN   ','ERR  ','ERR  ','AS-  ','AS0  ','AS+  ',  
     2'AL   ','AXI0 ','AXI+ ','ERR  ','ERR  ','ERR  ','ASC--','ASC- ',  
     3'ASC0 ','ALC- ','AUSC.','ADSC.','ASSC.','ASDC.','ASUC.','AUCC.',  
     4'ADCC.','ASCC.','ERR  ','ERR  ','ERR  ','ERR  ','AUUB.','AUDB.',  
     5'ADDB.','ADUB.','AUSB.','ADSB.','ASSB.','ASDB.','ASUB.','AUCB.',  
     6'ADCB.','ASCB.','ACCB.','ACSB.','ACDB.','ACUB.','AUBB.','ADBB.',  
     7'ASBB.','ACBB.','ERR  ','ERR  ','ERR  ','ERR  ','ERR  ','AUTT.',  
     8'AUDT.','ADDT.','ADUT.','AUST.','ADST.','ASST.','ASDT.','ASUT.',  
     9'AUCT.','ADCT.','ASCT.','ACCT.','ACST.','ACDT.','ACUT.','AUBT.',  
     1'ADBT.','ASBT.','ACBT.','ABBT.','ABCT.','ABST.','ABDT.','ABUT.',  
     2'AUTT.','ADTT.','ASTT.','ACTT.','ABTT.','ERR  ','ERR  ','ERR  ',  
     3'ERR  ','ERR  ','ERR  ','AUUY.','AUDY.','ADDY.','ADUY.','AUSY.',  
     4'ADSY.','ASSY.','ASDY.','ASUY.','AUUX.','AUDX.','ADDX.','ADUX.',  
     5'AUSX.','ADSX.','ASSX.','ASDX.','ASUX.'/  
C          3/2+ BARYON LABELS   
      DATA LBAR1/   
     1'DL++ ','DL+  ','DL0  ','DL-  ','ERR  ','S*+  ','S*0  ','S*-  ',  
     2'ERR  ','XI*0 ','XI*- ','OM-  ','ERR  ','ERR  ','UUC* ','UDC* ',  
     3'DDC* ','ERR  ','USC* ','DSC* ','SSC* ','ERR  ','ERR  ','UCC* ',  
     4'DCC* ','SCC* ','CCC* ','ERR  ','ERR  ','ERR  ','UUB* ','UDB* ',  
     5'DDB* ','ERR  ','USB* ','DSB* ','SSB* ','ERR  ','ERR  ','UCB* ',  
     6'DCB* ','SCB* ','CCB* ','ERR  ','ERR  ','ERR  ','UBB* ','DBB* ',  
     7'SBB* ','CBB* ','BBB* ','ERR  ','ERR  ','ERR  ','ERR  ','UTT* ',  
     8'UDT* ','DDT* ','ERR  ','UST* ','DST* ','SST* ','ERR  ','ERR  ',  
     9'UCT* ','DCT* ','SCT* ','CCT* ','ERR  ','ERR  ','ERR  ','UBT* ',  
     1'DBT* ','SBT* ','CBT* ','BBT* ','ERR  ','ERR  ','ERR  ','ERR  ',  
     2'UTT* ','DTT* ','STT* ','CTT* ','BTT* ','TTT* ','ERR  ','ERR  ',  
     3'ERR  ','ERR  ','ERR  ','UUY* ','UDY* ','DDY* ','ERR  ','USY* ',  
     4'DSY* ','SSY* ','ERR  ','ERR  ','UUX* ','UDX* ','DDX* ','ERR  ',  
     5'USX* ','DSX* ','SSX* ','ERR  ','ERR  '/  
      DATA LABAR1/  
     1'ADL--','ADL- ','ADL0 ','ADL+ ','ERR  ','AS*- ','AS*0 ','AS*+ ',  
     2'ERR  ','AXI*0','AXI*+','AOM+ ','ERR  ','ERR  ','AUUC*','AUDC*',  
     3'ADDC*','ERR  ','AUSC*','ADSC*','ASSC*','ERR  ','ERR  ','AUCC*',  
     4'ADCC*','ASCC*','ACCC*','ERR  ','ERR  ','ERR  ','AUUB*','AUDB*',  
     5'ADDB*','ERR  ','AUSB*','ADSB*','ASSB*','ERR  ','ERR  ','AUCB*',  
     6'ADCB*','ASCB*','ACCB*','ERR  ','ERR  ','ERR  ','AUBB*','ADBB*',  
     7'ASBB*','ACBB*','ABBB*','ERR  ','ERR  ','ERR  ','ERR  ','AUTT*',  
     8'AUDT*','ADDT*','ERR  ','AUST*','ADST*','ASST*','ERR  ','ERR  ',  
     9'AUCT*','ADCT*','ASCT*','ACCT*','ERR  ','ERR  ','ERR  ','AUBT*',  
     1'ADBT*','ASBT*','ACBT*','ABBT*','ERR  ','ERR  ','ERR  ','ERR  ',  
     2'AUTT*','ADTT*','ASTT*','ACTT*','ABTT*','ATTT*','ERR  ','ERR  ',  
     3'ERR  ','ERR  ','ERR  ','AUUY*','AUDY*','ADDY*','ERR  ','AUSY*',  
     4'ADSY*','ASSY*','ERR  ','ERR  ','AUUX*','AUDX*','ADDX*','ERR  ',  
     5'AUSX*','ADSX*','ASSX*','ERR  ','ERR  '/  
C          ENTRY
      CALL FLAVOR(ID,IFL1,IFL2,IFL3,JSPIN,INDEX)    
      IF(IABS(ID).LT.100) GO TO 200 
      IF(IABS(ID).LT.1000) GO TO 100    
      IF(ID.NE.0.AND.MOD(ID,100).EQ.0) GO TO 300    
C          BARYONS  
      INDEX=INDEX-109*JSPIN-36*NMES-NQLEP   
      INDEX=INDEX-11    
      IF(JSPIN.EQ.0.AND.ID.GT.0) LABEL=LBAR0(INDEX) 
      IF(JSPIN.EQ.0.AND.ID.LT.0) LABEL=LABAR0(INDEX)    
      IF(JSPIN.EQ.1.AND.ID.GT.0) LABEL=LBAR1(INDEX) 
      IF(JSPIN.EQ.1.AND.ID.LT.0) LABEL=LABAR1(INDEX)    
      RETURN    
C          MESONS   
100   CONTINUE  
      I=MAX0(IFL2,IFL3) 
      J=-MIN0(IFL2,IFL3)    
      INDEX=MAX0(I-1,J-1)**2+I+MAX0(I-J,0)  
      IF(JSPIN.EQ.0) LABEL=LMES0(INDEX) 
      IF(JSPIN.EQ.1) LABEL=LMES1(INDEX) 
      RETURN    
C          QUARKS, LEPTONS, ETC.    
200   CONTINUE  
      INDEX=2*INDEX 
      IF(ID.LE.0) INDEX=INDEX+1 
      LABEL=LLEP(INDEX) 
      RETURN    
300   I=IABS(IFL1)  
      J=IABS(IFL2)  
      INDEX=I+J*(J-1)/2 
      IF(ID.GT.0) LABEL=LQQ(INDEX)  
      IF(ID.LT.0) LABEL=LAQQ(INDEX) 
      RETURN    
      END   

      BLOCK DATA ALDATA 
C          INITIALIZE ALL COMMON BLOCKS 
C.......................................................................    
C          WARNINGS: MANY VARIABLES SET IN ALDATA ARE ALSO SET BY      .    
C          SUBROUTINE RESET.                                           .    
C                                                                      .    
C          ALDATA SHOULD ALWAYS BE LOADED WHEN USING ISAJET OR WHEN    .    
C          READING AN ISAJET TAPE.                                     .    
C.......................................................................    
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/PJETS/PJETS(5,3),IDJETS(3),QWJET(5),IDENTW 
     $,PPAIR(5,4),IDPAIR(4),JPAIR(4),NPAIR  
      COMMON/PINITS/PINITS(5,2),IDINIT(2)   
      COMMON/LSTPRT/LSTPRT  
C          LOOK MUST BE DIMENSIONED TO THE MAXIMUM VALUE OF INDEX.  
      PARAMETER (MXLOOK=400)    
      PARAMETER (MXDKY=1000)    
      COMMON/DKYTAB/LOOK(MXLOOK),CBR(MXDKY),MODE(5,MXDKY)   
      COMMON/DYLIM/QMIN,QMAX,QTMIN,QTMAX,YWMIN,YWMAX,XWMIN,XWMAX,THWMIN,    
     2  THWMAX,PHWMIN,PHWMAX    
     3  ,SETLMQ(12) 
      LOGICAL SETLMQ    
      DIMENSION BLIM1(12)   
      EQUIVALENCE(BLIM1(1),QMIN)    
      COMMON/FRGPAR/PUD,PBARY,SIGQT,PEND,XGEN(8),PSPIN1(8), 
     $PMIX1(3,2),PMIX2(3,2),XGENSS(9)   
      DIMENSION PMIXX1(6),PMIXX2(6) 
      EQUIVALENCE (PMIX1(1,1),PMIXX1(1))    
      EQUIVALENCE (PMIX2(1,1),PMIXX2(1))    
      DIMENSION FRPAR(32)   
      EQUIVALENCE(FRPAR(1),PUD) 
      COMMON/IDRUN/IDVER,IDG(2),IEVT    
      COMMON/JETLIM/PMIN(3),PMAX(3),PTMIN(3),PTMAX(3),YJMIN(3),YJMAX(3) 
     1 ,PHIMIN(3),PHIMAX(3),XJMIN(3),XJMAX(3),THMIN(3),THMAX(3) 
     2 ,SETLMJ(36)  
      DIMENSION BLIMS(36)   
      EQUIVALENCE(BLIMS(1),PMIN(1)) 
      LOGICAL SETLMJ    
      COMMON/FIXPAR/FIXP(3),FIXPT(3),FIXYJ(3),FIXPHI(3),FIXXJ(3)    
     2   ,FIXQM,FIXQT,FIXYW,FIXXW,FIXPHW    
      LOGICAL FIXQM,FIXQT,FIXYW,FIXXW,FIXPHW    
      LOGICAL FIXP,FIXPT,FIXYJ,FIXPHI,FIXXJ 
      COMMON/SGNPAR/CTHS(2,3),THS(2,3),YJS(2,3),XJS(2,3)    
      COMMON/JETPAR/P(3),PT(3),YJ(3),PHI(3),XJ(3),TH(3),CTH(3),STH(3)   
     1 ,JETTYP(3),SHAT,THAT,UHAT,QSQ,X1,X2,PBEAM(2) 
     2 ,QMW,QW,QTW,YW,XW,THW,QTMW,PHIW,SHAT1,THAT1,UHAT1,JWTYP  
     3 ,ALFQSQ,CTHW,STHW,Q0W    
     4 ,INITYP(2),ISIGS,PBEAMS(5)   
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      PARAMETER (MXSIGS=800,IOPAK=100)  
      COMMON/JETSIG/SIGMA,SIGS(MXSIGS),NSIGS,INOUT(MXSIGS),SIGEVT   
      COMMON/MBPAR/PUD0,PJSPN,PISPN,SIGQT0,XGEN0(2),PMIX01(3,2) 
     1,PMIX02(3,2),PBARY0   
      COMMON/NODCAY/NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR 
      LOGICAL NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR   
      PARAMETER (MXPTCL=4000,IPACK=10000)   
      COMMON/PARTCL/NPTCL,PPTCL(5,MXPTCL),IORIG(MXPTCL),IDENT(MXPTCL)   
     1,IDCAY(MXPTCL)    
      COMMON/PRIMAR/NJET,SCM,HALFE,ECM,IDIN(2),NEVENT,NTRIES,NSIGMA 
      COMMON/PRTOUT/NEVPRT,NJUMP    
      COMMON/QCDPAR/ALAM,ALAM2,CUTJET,ISTRUC    
      COMMON/QLMASS/AMLEP(52),NQLEP,NMES,NBARY  
      COMMON/Q1Q2/GOQ(25,3),GOALL(3),GODY(4),STDDY,GOWW(25,2),ALLWW(2)  
      LOGICAL GOQ,GOALL,GODY,STDDY,GOWW,ALLWW   
      COMMON/SEED/XSEED 
      CHARACTER*24 XSEED    
      COMMON/TOTALS/NKINPT,NWGEN,NKEEP,SUMWT,WT 
      COMMON/TYPES/LOC(50),NTYP,NJTTYP(3),NWWTYP(2) 
      COMMON/XTYPES/PARTYP(40),TITLE(10),JETYP(30,3),WWTYP(30,2)    
      CHARACTER*8 JETYP,WWTYP,TITLE,PARTYP  
      COMMON/PASCHAR/JETYP11,JETYP12,JETYP13,WWTYP11,WWTYP12  
      CHARACTER*8 JETYP11,JETYP12,JETYP13
      CHARACTER*8 WWTYP11,WWTYP12  
      COMMON/WCON/SIN2W,WMASS(4),WGAM(4),AQ(12,4),BQ(12,4),COUT(4), 
     1MATCH(25,4),WCBR(25,4),CUTOFF,CUTPOW,TBRWW(4,2),RBRWW(12,4,2),EZ, 
     2AQDP(12,4),BQDP(12,4),EZDP    
      DOUBLE PRECISION AQDP,BQDP,EZDP   
      PARAMETER (LIMPOM=20) 
      COMMON/MBGEN/POMWT(LIMPOM),POMGEN(LIMPOM),MNPOM,MXPOM,PDIFFR, 
     $NPOM,XBARY(2),DXBARY(2),XPOM(LIMPOM,2)    
      COMMON/FORCE/NFORCE,IFORCE(20),MFORCE(5,20)   
      PARAMETER (MAXLEN=1024)   
      COMMON/ZEVEL/IZEVEL(MAXLEN)   
      DIMENSION ZEVEL(MAXLEN)   
      EQUIVALENCE(ZEVEL(1),IZEVEL(1))   
      EQUIVALENCE (IZVL1,IZEVEL(1)) 
      EQUIVALENCE (IZVL2,IZEVEL(2)) 
      COMMON/FINAL/NKINF,SIGF,ALUM,ACCEPT,NRECS 
      COMMON/KEYS/IKEYS,KEYON,KEYS(10)  
      COMMON/XKEYS/REAC 
      LOGICAL KEYS  
      LOGICAL KEYON 
      CHARACTER*8 REAC  
C   
C   
C          DATA FOR IDRUN   
C          IDVER=100*VERSION+CYCLE  
      DATA IDVER/524/   
C   
C          DATA FOR ITAPES  
      DATA ITDKY,ITEVT,ITCOM,ITLIS/1,2,5,6/ 
C   
C          DATA FOR QLMASS  
C          AMLEP LABELED BY INDEX...SEE FLAVOR  
C          SETW RESETS W+- AND Z0 MASSES    
      DATA AMLEP/.3,.3,.5,1.6,4.9,30.,-1.,-1.,0.,0.,    
     $0.,.511003E-3,0.,.105661,0.,1.807,3*-1.,.49767,.49767,    
     $100.3,100.3,100.5,101.6,104.9,130.,2*-1.,100.,0., 
     $100.,100.005,100.,100.1,100.,101.8,2*-1.,100.,100.,   
     $11*0./    
      DATA NQLEP,NMES,NBARY/41,2,2/ 
C   
C          DATA FOR PJETS   
      DATA IDJETS/3*0/,IDENTW/0/    
C   
C          DATA FOR PINITS  
      DATA IDINIT/2*0/  
C   
C          DATA FOR LSTPRT  
      DATA LSTPRT/0/    
C   
C          DATA FOR MBPAR   
      DATA PUD0/.45/,PJSPN,PISPN/2*.5/,SIGQT0/.35/,XGEN0/.9,1./,PMIX01/ 
     $.25,.25,.5,0.,.5,1./,PMIX02/.5,.5,1.,0.,0.,1./    
      DATA PBARY0/.075/ 
C   
C          DATA FOR MBGEN   
      DATA MNPOM,MXPOM/1,LIMPOM/    
C   
C          DATA FOR SEED    
      DATA XSEED/'0'/   
C   
C          DATA FOR FRGPAR  
C          F(X)=1-XGEN(1)+XGEN(1)*(XGEN(2)+1)*(1-X)**XGEN(2) FOR U,D,S  
C          PETERSON FRAGMENTATION, EPSILON=XGEN(I)*M(I)**2 FOR C,B,T    
      DATA PUD,PBARY/.43,.10/   
      DATA SIGQT,PEND/.35,.14/  
      DATA XGEN/.96,3.,0.,.8,.5,.5,.5,.5/   
      DATA PSPIN1/.5,.5,.5,.75,.75,.75,.75,.75/ 
      DATA PMIX1/.25,.25,.5,0.,.5,1./,PMIX2/.5,.5,1.,0.,0.,1./  
      DATA XGENSS/9*.5/ 
C   
C          DATA FOR JETLIM  
      DATA BLIMS/36*-1.E9/  
C   
C          DATA FOR NODCAY  
      DATA NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR/6*.FALSE./   
C   
C          DATA FOR TYPES   
      DATA LOC/50*0/,NTYP/50/   
      DATA NJTTYP/3*0/  
      DATA JETYP/90*'        '/ 
      DATA JETYP11,JETYP12,JETYP13/3*'ALL     '/   
      DATA NWWTYP/2*0/  
      DATA WWTYP/60*'        '/ 
      DATA WWTYP11,WWTYP12/2*'ALL     '/  
      DATA JWTYP/4/ 
C   
C          DATA FOR PRIMAR  
      DATA IDIN/1120,1120/  
      DATA NTRIES/200/  
      DATA NSIGMA/20/   
C   
C          DATA FOR DKYTAB  
      DATA LOOK/MXLOOK*0/   
      DATA CBR/MXDKY*0./    
      DATA MODE/MXDKY*0,MXDKY*0,MXDKY*0,MXDKY*0,MXDKY*0/    
C   
C          DATA FOR Q1Q2    
      DATA GOQ/75*.TRUE./   
      DATA GOALL/3*.TRUE./  
      DATA GODY/.TRUE.,.FALSE.,.FALSE.,.TRUE./  
      DATA GOWW/50*.TRUE./,ALLWW/2*.TRUE./  
C   
C          DATA FOR WCON    
      DATA MATCH/   
     $0,3,2,5,4,7,6,9,8,11,10,13,12,0,0,17,16,0,0,21,20,0,0,25,24,  
     $0,5,0,0,2,0,8,7,0,0,12,11,0,17,0,0,14,21,0,0,18,25,0,0,22,    
     $0,0,4,3,0,9,0,0,6,13,0,0,10,0,16,15,0,0,20,19,0,0,24,23,0,    
     $0,3,2,5,4,7,6,9,8,11,10,13,12,15,14,17,16,19,18,21,20,23,22,25,24/    
      DATA CUTOFF,CUTPOW/.200,1.0/  
C   
C          DATA FOR TOTALS  
      DATA NKINPT,NWGEN,NKEEP/3*0/,SUMWT/0./    
C   
C          DATA FOR DYLIM   
      DATA BLIM1/12*-1.E9/  
C   
C          DATA FOR PARTCL  
      DATA NPTCL/0/ 
C   
C            DATA FOR PRTOUT    
      DATA NEVPRT,NJUMP/1,1/    
C   
C          DATA FOR JETSET  
      DATA NJSET/0/ 
C   
C          DATA FOR QCDPAR  
      DATA ALAM,ALAM2/.2,.04/,CUTJET/6./,ISTRUC/3/  
C   
C          DATA FOR FORCE   
      DATA NFORCE/0/    
C   
C          DATA FOR NRECS   
      DATA NRECS/0/ 
C   
C          DATA FOR KEYS    
      DATA KEYS/10*.FALSE./ 
C   
      END   
      FUNCTION AMASS(ID)    
C          THIS FUNCTION RETURNS THE MASS OF THE PARTICLE WITH  
C          IDENT CODE ID.   
C          QUARK-BASED IDENT CODE   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/QLMASS/AMLEP(52),NQLEP,NMES,NBARY  
      DIMENSION AMMES0(10),AMMES1(10),AMBAR0(30),AMBAR1(30) 
C          0- MESON MASS TABLE  
      DATA AMMES0/.13496,.13957,.5488,.49367,.49767,.9576,1.8633    
     1,1.8683,2.030,2.976/  
C          1- MESON MASS TABLE  
      DATA AMMES1/.776,.776,.7826,.8881,.8922,1.0196,2.006,2.0086   
     1,2.140,3.097/ 
C          1/2+ BARYON MASS TABLE   
      DATA AMBAR0/-1.,.93828,.93957,2*-1.,1.1894,1.1925,1.1974  
     1,1.1156,1.3149,1.3213,3*-1.,2.43,2.43,2.43,2.26   
     2,2.50,2.50,2.60,2.40,2.40,3.55,3.55,3.70,4*-1./   
C          3/2+ BARYON MASS TABLE   
      DATA AMBAR1/1.232,1.232,1.232,1.232,-1.,1.3823,1.3820 
     1,1.3875,-1.,1.5318,1.5350,1.6722,2*-1.    
     2,2.63,2.63,2.63,-1.,2.70,2.70,2.80,2*-1.,3.75,3.75    
     3,3.90,4.80,3*-1./ 
C          ENTRY    
      CALL FLAVOR(ID,IFL1,IFL2,IFL3,JSPIN,INDEX)    
      IF(ID.NE.0.AND.MOD(ID,100).EQ.0) GO TO 400    
      IF(IABS(IFL1).GT.4.OR.IABS(IFL2).GT.4.OR.IABS(IFL3).GT.4) 
     1GO TO 300 
      IF(IFL2.EQ.0) GO TO 200   
      IF(IFL1.EQ.0) GO TO 100   
C          BARYONS  
      INDEX=INDEX-109*JSPIN-36*NMES-NQLEP   
      INDEX=INDEX-11    
      AMASS=(1-JSPIN)*AMBAR0(INDEX)+JSPIN*AMBAR1(INDEX) 
      RETURN    
C          MESONS   
100   CONTINUE  
      INDEX=INDEX-36*JSPIN-NQLEP    
      INDEX=INDEX-11    
      AMASS=(1-JSPIN)*AMMES0(INDEX)+JSPIN*AMMES1(INDEX) 
      RETURN    
C          QUARKS AND LEPTONS   
200   CONTINUE  
      AMASS=AMLEP(INDEX)    
      RETURN    
C          B AND T PARTICLES    
300   CONTINUE  
      AMASS=AMLEP(IABS(IFL2))+AMLEP(IABS(IFL3))-.03+.04*JSPIN   
      IF(IFL1.NE.0) AMASS=AMASS+AMLEP(IABS(IFL1))   
      RETURN    
C          DIQUARKS 
400   AMASS=AMLEP(IABS(IFL1))+AMLEP(IABS(IFL2)) 
      RETURN    
      END   
      FUNCTION CHARGE(ID)   
C   
C          COMPUTE CHARGE OF PARTICLE WITH IDENT CODE ID    
C          ICHRG MUST BE DIMENSIONED NQLEP+12   
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      DIMENSION ICHRG(53),IFL(3)    
      DATA ICHRG/0,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,0,-3,3*0  
     $,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,0,-3,3,0  
     $,3,0,0,0,3,3,3,6,6,6,0/   
      IDABS=IABS(ID)    
      CALL FLAVOR(ID,IFL(1),IFL(2),IFL(3),JSPIN,INDEX)  
      IF(IDABS.LT.100) GO TO 200    
      ISUM=0    
      DO 100 I=1,3  
      ISUM=ISUM+ICHRG(IABS(IFL(I))+1)*ISIGN(1,IFL(I))   
  100 CONTINUE  
      CHARGE=ISUM/3.    
      RETURN    
200   CHARGE=ICHRG(INDEX+1)*ISIGN(1,ID) 
      CHARGE=CHARGE/3.  
      RETURN    
      END   
      FUNCTION DBLPCM(A,B,C)    
C          CALCULATE COM MOMENTUM FOR A-->B+C WITH DOUBLE PRECISION.    
C          NEEDED TO FIX BUG ON 32-BIT MACHINES AT HIGH ENERGY. 
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      DOUBLE PRECISION DA,DB,DC,DVAL    
C          CONVERT TO DOUBLE PRECISION  
      DA=A  
      DB=B  
      DC=C  
      DVAL=(DA**2-DB**2-DC**2)**2-(2.*DB*DC)**2 
C          CONVERT BACK TO SINGLE PRECISION 
      VAL=DVAL  
      DBLPCM=SQRT(VAL)/(2.*A)   
      RETURN    
      END   
      SUBROUTINE GETPT(PT0,PTMEAN)  
C          GENERATE PT WITH 1/(1+B*PT**2)**4 DISTRIBUTION   
C          (APPROXIMATELY AN EXPONENTIAL FOR PT < 2 GEV.)   
C          CON1=16/(3*PI)   
C          CON2=-1/3    
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      DATA CON1/1.697652726/,CON2/-.3333333333/ 
      R=RANF()  
      ARG=AMAX1(R**CON2-1.,0.)  
      PT0=PTMEAN*CON1*SQRT(ARG) 
      RETURN    
      END   
      SUBROUTINE IDGEN  
C   
C          CALL SYSTEM DATE AND TIME ROUTINES (NON-STANDARD) TO SET UP  
C          RUN IDENTIFICATION:  
C          IDVER=100*VERSN     (INTEGER ISAJET VERSION NUMBER)  
C          IDG(1)=YYMMDD       (INTEGER YEAR-MONTH-DAY) 
C          IDG(2)=HHMMSS       (INTEGER HOUR-MINUTE-SECOND) 
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/IDRUN/IDVER,IDG(2),IEVT    
C          DEFAULT RUN ID IS ZERO.  
      IYMD=0.   
      IHMS=0.   
C          CALL DATIME FOR DATE AND TIME. (IN CERN LIBRARY) 
CCCCC      CALL DATIME(IYMD,IHMS)    
      IDG(1)=IYMD   
      IDG(2)=IHMS   
      RETURN    
      END   
      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   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      PARAMETER (MXPTCL=4000,IPACK=10000)   
      COMMON/PARTCL/NPTCL,PPTCL(5,MXPTCL),IORIG(MXPTCL),IDENT(MXPTCL)   
     1,IDCAY(MXPTCL)    
      COMMON/FRGPAR/PUD,PBARY,SIGQT,PEND,XGEN(8),PSPIN1(8), 
     $PMIX1(3,2),PMIX2(3,2),XGENSS(9)   
      DIMENSION PMIXX1(6),PMIXX2(6) 
      EQUIVALENCE (PMIX1(1,1),PMIXX1(1))    
      EQUIVALENCE (PMIX2(1,1),PMIXX2(1))    
      DIMENSION FRPAR(32)   
      EQUIVALENCE(FRPAR(1),PUD) 
      COMMON/CONST/PI,SQRT2,ALFA,GF,UNITS   
      COMMON/MBPAR/PUD0,PJSPN,PISPN,SIGQT0,XGEN0(2),PMIX01(3,2) 
     1,PMIX02(3,2),PBARY0   
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)  
      IF(IFLBEG.NE.9) GO TO 200 
C          EQUIVALENT QUARK FOR GLUON   
      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   
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=ID1   
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    
      IF(PPLUS.GT.PEND.AND.LOOP.LT.10000) GO TO 300 
      RETURN    
9999  CALL PRTEVT(0)    
      WRITE(ITLIS,10) NPTCL 
10    FORMAT(//5X,'ERROR IN JETGEN...NPTCL >',I5)   
      RETURN    
      END   
      SUBROUTINE PRTEVT(IPRT)   
C   
C          PRINT THE EVENT STORED IN /PJETS/, /JETSET/, AND /PARTCL/    
C          IF IPRT IS SELECTED BY NEVPRT AND NJUMP. 
C          IPRT=0 ALWAYS PRINTS EVENT   
C          IPRT<0 PRINTS ONLY JET PARAMETERS    
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (LIMPOM=20) 
      COMMON/MBGEN/POMWT(LIMPOM),POMGEN(LIMPOM),MNPOM,MXPOM,PDIFFR, 
     $NPOM,XBARY(2),DXBARY(2),XPOM(LIMPOM,2)    
      COMMON/PJETS/PJETS(5,3),IDJETS(3),QWJET(5),IDENTW 
     $,PPAIR(5,4),IDPAIR(4),JPAIR(4),NPAIR  
      COMMON/PINITS/PINITS(5,2),IDINIT(2)   
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      COMMON/IDRUN/IDVER,IDG(2),IEVT    
      PARAMETER (MXSIGS=800,IOPAK=100)  
      COMMON/JETSIG/SIGMA,SIGS(MXSIGS),NSIGS,INOUT(MXSIGS),SIGEVT   
      COMMON/KEYS/IKEYS,KEYON,KEYS(10)  
      COMMON/XKEYS/REAC 
      LOGICAL KEYS  
      LOGICAL KEYON 
      CHARACTER*8 REAC  
      COMMON/JETPAR/P(3),PT(3),YJ(3),PHI(3),XJ(3),TH(3),CTH(3),STH(3)   
     1 ,JETTYP(3),SHAT,THAT,UHAT,QSQ,X1,X2,PBEAM(2) 
     2 ,QMW,QW,QTW,YW,XW,THW,QTMW,PHIW,SHAT1,THAT1,UHAT1,JWTYP  
     3 ,ALFQSQ,CTHW,STHW,Q0W    
     4 ,INITYP(2),ISIGS,PBEAMS(5)   
      COMMON/LSTPRT/LSTPRT  
      PARAMETER (MXPTCL=4000,IPACK=10000)   
      COMMON/PARTCL/NPTCL,PPTCL(5,MXPTCL),IORIG(MXPTCL),IDENT(MXPTCL)   
     1,IDCAY(MXPTCL)    
      COMMON/PRIMAR/NJET,SCM,HALFE,ECM,IDIN(2),NEVENT,NTRIES,NSIGMA 
      COMMON/PRTOUT/NEVPRT,NJUMP    
      COMMON/WSIG/SIGLLQ    
      COMMON/SEED/XSEED 
      CHARACTER*24 XSEED    
C   
C          LABELS ARE CHARACTER*8   
      CHARACTER*8 LABEL,LW,LJET,LPTCL   
C   
C          DECIDE WHETHER TO PRINT  
      IF(IPRT.GT.NJUMP*NEVPRT) THEN 
        IF(NJUMP.GT.0) THEN 
          IF(MOD(IPRT,NJUMP).EQ.0) WRITE(ITLIS,5) IDG,IEVT,XSEED    
5         FORMAT(/6X,'RUN ID',2I10,5X,'EVENT NO',I8,5X,'SEED',2X,A24)   
        ENDIF   
        RETURN  
      ENDIF 
      IF(IPRT.GT.1.AND.MOD(IPRT,NJUMP).NE.0) RETURN 
      IF(IEVT.EQ.LSTPRT) RETURN 
      PI=4.*ATAN(1.)    
      LSTPRT=IEVT   
C          PRINT EVENT NUMBER   
      WRITE(ITLIS,10) IDG,IEVT,XSEED    
10    FORMAT('1',5X,'RUN ID',2I10,5X,'EVENT NO',I8,5X,'SEED',2X,A24)    
C   
C          PRINT JET PARAMETERS 
      IF(NJET.EQ.0) GO TO 300   
      WRITE(ITLIS,20)   
20    FORMAT(//20X,'JET PARAMETERS'//3X,'JET',4X,'TYPE ',   
     18X,'PX',8X,'PY',8X,'PZ',8X,'P0',8X,'PT',  
     25X,'THETA',7X,'PHI',9X,'X',9X,'Y')    
      IF(KEYS(3)) THEN  
        LW=LABEL(IDENTW)    
        WRITE(ITLIS,31) LW,(QWJET(K),K=1,4),QTW,THW,PHIW,XW,YW  
31      FORMAT(5X,'W',4X,A5,5F10.2,4F10.4)  
      ENDIF 
      DO 100 I=1,NJET   
      LJET=LABEL(IDJETS(I)) 
      WRITE(ITLIS,30) I,LJET,(PJETS(K,I),K=1,4),PT(I),TH(I),PHI(I), 
     $XJ(I),YJ(I)   
30    FORMAT(1X,I5,4X,A5,5F10.2,4F10.4) 
100   CONTINUE  
C          PRINT WPAIR DECAY PARAMETERS 
      IF(KEYS(6)) THEN  
        WRITE(ITLIS,101)    
101     FORMAT(//20X,'WPAIR DECAY PARAMETERS'//3X,'JET',4X,'TYPE ', 
     1  8X,'PX',8X,'PY',8X,'PZ',8X,'P0',8X,'PT',    
     2  5X,'THETA',7X,'PHI',9X,'X',9X,'Y')  
        DO 102 I=1,NPAIR    
        JET=JPAIR(I)    
        LJET=LABEL(IDPAIR(I))   
        PTPRT=SQRT(PPAIR(1,I)**2+PPAIR(2,I)**2) 
        THPRT=ACOS(PPAIR(3,I)/SQRT(PTPRT**2+PPAIR(3,I)**2)) 
        PHIPRT=ATAN2(PPAIR(2,I),PPAIR(1,I)) 
        XPRT=2*PPAIR(3,I)/ECM   
        YPRT=-ALOG(TAN(THPRT/2.))   
        WRITE(ITLIS,30) JET,LJET,(PPAIR(K,I),K=1,4),PTPRT,THPRT,PHIPRT, 
     $  XPRT,YPRT   
102     CONTINUE    
      ENDIF 
      IF(IPRT.LT.0) RETURN  
C   
C          PRINT JET CROSS SECTIONS 
201   CONTINUE  
      IF(KEYS(1)) THEN  
        WRITE(ITLIS,39) SIGEVT,SIGMA    
39      FORMAT(/'  D(SIGMA)/D(PT**2)D(Y1)D(Y2)=',E12.4,5X,  
     $  'SUMMED OVER QUARK TYPES=',E12.4)   
      ENDIF 
      IF(KEYS(3)) THEN  
        IF(NJET.EQ.3) THEN  
          WRITE(ITLIS,38) LW,QMW,SIGEVT,SIGLLQ  
38        FORMAT(//5X,'MASS OF ',A2,' = ',F10.3//   
     C    5X,'D(SIGMA)/D(Q**2)D(QT**2)D(YW)D(YJ) = ',E12.4/ 
     C    5X,'D(SIGMA)/D(Q**2)D(QT**2)D(YW)D(YJ)D(OMEGA) = ',E12.4) 
        ELSE    
          WRITE(ITLIS,37) LW,QMW,SIGEVT,SIGLLQ  
37        FORMAT(//5X,'MASS OF ',A2,' = ',F10.3//   
     1    5X,'D(SIGMA)/D(Q**2)D(YW) =',E12.4/   
     2    5X,'D(SIGMA)/D(Q**2)D(YW)D(OMEGA) =',E12.4)   
        ENDIF   
      ENDIF 
      IF(KEYS(2)) THEN  
        WRITE(ITLIS,32)  SIGEVT,SIGMA   
32      FORMAT(/'  D(SIGMA)/D(COS TH) = ',E12.4,5X, 
     C  'SUMMED OVER TYPES = ',E12.4)   
      ENDIF 
C   
C          PRINT /JETSET/ PARAMETERS    
      IF(KEYS(4)) GO TO 300 
      IF(NJSET.EQ.0) GO TO 300  
      WRITE(ITLIS,70)   
70    FORMAT(//20X,'PARTON CASCADE PARAMETERS'//    
     C  6X,'I',3X,'JET',4X,'ORIG',4X,'TYPE',9X,'DECAY', 
     C  8X,'PX',8X,'PY',8X,'PZ',8X,'P0',6X,'MASS')  
      DO 310 I=1,NJSET  
      JET=JORIG(I)/JPACK    
      I1=MOD(JORIG(I),JPACK)    
      JTLV1=JTYPE(I)    
      LJET=LABEL(JTLV1) 
      J1=JDCAY(I)/JPACK 
      J2=MOD(JDCAY(I),JPACK)    
      IF(JDCAY(I).EQ.0) THEN    
        WRITE(ITLIS,71) I,JET,I1,LJET,(PJSET(K,I),K=1,5)    
71      FORMAT(1X,I6,I6,I8,4X,A5,8X,'FINAL',5F10.2) 
      ELSEIF(J1.NE.J2) THEN 
        WRITE(ITLIS,72) I,JET,I1,LJET,J1,J2,(PJSET(K,I),K=1,5)  
72      FORMAT(1X,I6,I6,I8,4X,A5,4X,I4,'-',I4,5F10.2)   
      ELSE  
        WRITE(ITLIS,73) I,JET,I1,LJET,(PJSET(K,I),K=1,5)    
73      FORMAT(1X,I6,I6,I8,4X,A5,6X,'INITIAL',5F10.2)   
      ENDIF 
310   CONTINUE  
C   
C          PRINT HADRON PARAMETERS  
300   IF(NPTCL.EQ.0) RETURN 
      IF(.NOT.KEYS(2)) WRITE(ITLIS,45) NPOM 
45    FORMAT(//' NUMBER OF POMERONS =',I5)  
      WRITE(ITLIS,40)   
40    FORMAT(//20X,'HADRON PARAMETERS'//7X,'I',3X,'JET',5X,'ORIG'   
     C  ,4X,'TYPE',11X,'DECAY',8X,'PX',8X,'PY',8X,'PZ',8X,'P0'  
     C  ,8X,'PT',5X,'THETA',7X,'PHI')   
      DO 200 I=1,NPTCL  
      I1=IABS(IORIG(I)) 
      JET=I1/IPACK  
      I1=I1-IPACK*JET   
      I1=ISIGN(I1,IORIG(I)) 
      IDLV1=IDENT(I)    
      LPTCL=LABEL(IDLV1)    
      J1=IDCAY(I)/IPACK 
      J2=MOD(IDCAY(I),IPACK)    
      PTHAD=SQRT(PPTCL(1,I)**2+PPTCL(2,I)**2)   
      PHAD=SQRT(PPTCL(3,I)**2+PTHAD**2) 
      PHIHAD=ATAN2(PPTCL(2,I),PPTCL(1,I))   
      IF(PHIHAD.LT.0.) PHIHAD=2.*PI+PHIHAD  
      THHAD=ACOS(PPTCL(3,I)/PHAD)   
      IF(IDCAY(I).EQ.0) THEN    
        WRITE(ITLIS,50) I,JET,I1,LPTCL,(PPTCL(K,I),K=1,4),PTHAD,THHAD,  
     $  PHIHAD  
50      FORMAT(1X,I7,I6,I9,4X,A5,9X,'STABLE',5F10.2,2F10.4) 
      ELSE  
        WRITE(ITLIS,60) I,JET,I1,LPTCL,J1,J2,(PPTCL(K,I),K=1,4),    
     $  PTHAD,THHAD,PHIHAD  
60      FORMAT(1X,I7,I6,I9,4X,A5,4X,I5,'-',I5,5F10.2,2F10.4)    
      ENDIF 
200   CONTINUE  
      RETURN    
      END   
      SUBROUTINE QCDJET(NJMIN)  
C   
C          CARRY OUT FINAL STATE QCD JET EVOLUTION USING THE ALGORITHM  
C          OF FOX AND WOLFRAM.  EVOLVE EACH PARTON IN T WITH FIXED ZC   
C          AND ITERATE AS FOLLOWS-- 
C   
C          (0) EVOLVE INITIAL PARTONS.  
C          (1) PICK I AND FIND MATCHING J>I.    
C          (2) SOLVE KINEMATICS.    
C          (3) FOR K=I,J, GENERATE Z(K) AND EVOLVE T(K1), T(K2).  IF NO 
C              GOOD, EVOLVE T(K). OTHERWISE, ADD K1 AND K2 TO /JETSET/. 
C          (4) IF I OR J NO GOOD, THEN (2). 
C          (5) THEN (1).    
C   
C          USE Z=(E+P)/(E0+P0) AND A LARGE TCUT.    
C          JMATCH(J1)=J2 IF J1 AND J2 MATCH.    
C          JMATCH(J)=JPACK*J1+J2 IF J1,...,J2 MATCH. USED FOR MULTIPLE  
C          INITIAL PARTONS. 
C          JMATCH(J)=0 FOR INITIAL JET PARTONS  
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (MXPTCL=4000,IPACK=10000)   
      COMMON/PARTCL/NPTCL,PPTCL(5,MXPTCL),IORIG(MXPTCL),IDENT(MXPTCL)   
     1,IDCAY(MXPTCL)    
      COMMON/QCDPAR/ALAM,ALAM2,CUTJET,ISTRUC    
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      COMMON/JWORK/ZZC(MXJSET),JMATCH(MXJSET),TNEW,P1CM(4), 
     1J1,J2,J3,J4,J5,E1CM,E2CM,E3CM,E4CM,E5CM   
      LOGICAL TNEW  
      DIMENSION JJ(5),EE(5) 
      EQUIVALENCE (J1,JJ(1)),(E1CM,EE(1))   
      COMMON/CONST/PI,SQRT2,ALFA,GF,UNITS   
C   
      DIMENSION PSUM(5) 
      DATA PSUM/5*0./   
C   
C          (0) EVOLVE INITIAL PARTON MASSES.    
C   
      DO 100 J=NJMIN,NJSET  
      J1=J  
      J2=JMATCH(J)  
      IF(J2.GT.JPACK) GO TO 150 
      IF(J2.LE.J1) GO TO 100    
C          TWO PARTONS  
      IF(JDCAY(J1).EQ.-1) CALL QCDT(J1) 
      IF(JDCAY(J2).EQ.-1) CALL QCDT(J2) 
      JPRNT=MOD(JORIG(J),JPACK) 
      IF(JPRNT.EQ.0) THEN   
        AM0=PJSET(4,J1)+PJSET(4,J2) 
      ELSE  
        AM0=PJSET(5,JPRNT)  
      ENDIF 
110   AM1=PJSET(5,J1)   
      AM2=PJSET(5,J2)   
      IF(AM0.LE.AM1+AM2) GO TO 120  
      GO TO 100 
120   J3=J1 
      IF(RANF().GT..5) J3=J2    
      IF(JDCAY(J3).EQ.-1) CALL QCDT(J3) 
      GO TO 110 
C          MORE THAN TWO PARTONS    
150   JI1=JMATCH(J)/JPACK   
      IF(J.NE.JI1) GO TO 100    
      JI2=JMATCH(J)-JPACK*JI1   
      NJI=JI2-JI1+1 
      AM0=0.    
      AMSUM=0.  
      DO 160 JI=JI1,JI2 
      IF(JDCAY(JI).EQ.-1) CALL QCDT(JI) 
      AM0=AM0+PJSET(4,JI)   
160   AMSUM=AMSUM+PJSET(5,JI)   
170   IF(AM0.LT.AMSUM) THEN 
        J3=NJI*RANF()+JI1   
        AMSUM=AMSUM-PJSET(5,J3) 
        IF(JDCAY(J3).EQ.-1) CALL QCDT(J3)   
        AMSUM=AMSUM+PJSET(5,J3) 
        GO TO 170   
      ENDIF 
100   CONTINUE  
C   
C          (1) LOOP OVER ACTIVE PARTONS 
C   
      NJ1=NJMIN 
1     NJ2=NJSET 
      DO 200 J=NJ1,NJ2  
      J1=J  
      J2=JMATCH(J1) 
      NJI=2 
      IF(J2.LE.J1) GO TO 200    
C   
C          (2) SOLVE KINEMATICS.    
C   
C          INITIAL PARTONS--KEEP DIRECTIONS FIXED.  
210   IF(MOD(JORIG(J),JPACK).NE.0) GO TO 230    
      IF(JMATCH(J).GT.JPACK) GO TO 400  
      AM0=PJSET(4,J1)+PJSET(4,J2)   
      AM1=PJSET(5,J1)   
      AM2=PJSET(5,J2)   
      PJSET(4,J1)=(AM0**2+AM1**2-AM2**2)/(2*AM0)    
      PJSET(4,J2)=(AM0**2+AM2**2-AM1**2)/(2*AM0)    
      PCM2=((AM0**2-AM1**2-AM2**2)**2-(2*AM1*AM2)**2)/(4*AM0**2)    
      DO 220 L=1,2  
      POLD2=PJSET(1,JJ(L))**2+PJSET(2,JJ(L))**2+PJSET(3,JJ(L))**2   
      RATIO=SQRT(PCM2/POLD2)    
      DO 225 K=1,3  
225   PJSET(K,JJ(L))=RATIO*PJSET(K,JJ(L))   
220   CONTINUE  
      GO TO 300 
C   
C          NJI.LE.5 INITIAL PARTONS 
400   CONTINUE  
      JI1=JMATCH(J)/JPACK   
      IF(J.NE.JI1) GO TO 200    
      JI2=JMATCH(J)-JPACK*JI1   
      NJI=JI2-JI1+1 
      AM0=0.    
      DO 410 JI=JI1,JI2 
      AM0=AM0+PJSET(4,JI)   
      JJ(JI-JI1+1)=JI   
      PJSET(4,JI)=SQRT(PJSET(1,JI)**2+PJSET(2,JI)**2+PJSET(3,JI)**2 
     1+PJSET(5,JI)**2)  
      DO 420 K=1,5  
420   PPTCL(K,NPTCL+JI-JI1+1)=PJSET(K,JI)   
410   CONTINUE  
      PSUM(4)=AM0   
      PSUM(5)=PSUM(4)   
      NPTLV1=NPTCL  
      CALL RESCAL(NPTLV1+1,NPTLV1+NJI,PSUM,IFAIL)   
      DO 430 JI=JI1,JI2 
      DO 430 K=1,5  
      PJSET(K,JI)=PPTCL(K,NPTCL+JI-JI1+1)   
430   CONTINUE  
      GO TO 300 
C   
C          SOLVE KINEMATICS FOR GENERAL PARTONS.    
230   J0=MOD(JORIG(J),JPACK)    
      AM0=PJSET(5,J0)   
      AM1=PJSET(5,J1)   
      AM2=PJSET(5,J2)   
      E1CM=(AM0**2+AM1**2-AM2**2)/(2*AM0)   
      E2CM=(AM0**2+AM2**2-AM1**2)/(2*AM0)   
      P12CM=SQRT((AM0**2-AM1**2-AM2**2)**2-(2*AM1*AM2)**2)/(2*AM0)  
      NJI=2 
C          DETERMINE E1, P1, AND COS(THCM) FROM Z(J0).  
C          OCCASIONALLY COS(TH)>1.  IF SO THEN RESET Z. 
      E0=PJSET(4,J0)    
      P0=SQRT(PJSET(1,J0)**2+PJSET(2,J0)**2+PJSET(3,J0)**2) 
      Z1=ZZC(J0)    
      IF(Z1.GT.0.5) THEN    
        E1MAX=(E0*E1CM+P0*P12CM)/AM0    
        P1MAX=(P0*E1CM+E0*P12CM)/AM0    
        ZMAX=(E1MAX+P1MAX)/(E0+P0)  
        E1MIN=(E0*E1CM-P0*P12CM)/AM0    
        P1MIN=(P0*E1CM-E0*P12CM)/AM0    
        P1MIN=ABS(P1MIN)    
        ZMIN=(E1MIN+P1MIN)/(E0+P0)  
        IF(Z1.LT.ZMIN.OR.Z1.GT.ZMAX) Z1=ZMIN+Z1*(ZMAX-ZMIN) 
        ZZC(J0)=Z1  
        ZEP=Z1*(E0+P0)  
        P1=(ZEP**2-AM1**2)/(2.*ZEP) 
        E1=(ZEP**2+AM1**2)/(2.*ZEP) 
        CTHCM=(E1*AM0-E0*E1CM)/(P0*P12CM)   
      ELSE  
        Z2=1.-Z1    
        E2MAX=(E0*E2CM+P0*P12CM)/AM0    
        P2MAX=(P0*E2CM+E0*P12CM)/AM0    
        ZMAX=(E2MAX+P2MAX)/(E0+P0)  
        E2MIN=(E0*E2CM-P0*P12CM)/AM0    
        P2MIN=(P0*E2CM-E0*P12CM)/AM0    
        P2MIN=ABS(P2MIN)    
        ZMIN=(E2MIN+P2MIN)/(E0+P0)  
        IF(Z2.LT.ZMIN.OR.Z2.GT.ZMAX) Z2=ZMIN+Z2*(ZMAX-ZMIN) 
        ZZC(J0)=Z2  
        ZEP=Z2*(E0+P0)  
        P2=(ZEP**2-AM2**2)/(2.*ZEP) 
        E2=(ZEP**2+AM2**2)/(2.*ZEP) 
        CTHCM=-(E2*AM0-E0*E2CM)/(P0*P12CM)  
      ENDIF 
C          AVOID DISASTER   
      IF(ABS(CTHCM).GT.1.) CTHCM=SIGN(RANF(),CTHCM) 
      STHCM=SQRT(1.-CTHCM**2)   
      PHICM=2*PI*RANF() 
      CPHICM=COS(PHICM) 
      SPHICM=SIN(PHICM) 
C   
C          CONSTRUCT CM MOMENTA.    
      PT0=SQRT(PJSET(1,J0)**2+PJSET(2,J0)**2)   
      CTH0=PJSET(3,J0)/P0   
      STH0=PT0/P0   
      CPHI0=PJSET(1,J0)/PT0 
      SPHI0=PJSET(2,J0)/PT0 
      P1CM(1)=P12CM*(CPHI0*(CTH0*CPHICM*STHCM+STH0*CTHCM)   
     1-SPHI0*SPHICM*STHCM)  
      P1CM(2)=P12CM*(SPHI0*(CTH0*CPHICM*STHCM+STH0*CTHCM)   
     1+CPHI0*SPHICM*STHCM)  
      P1CM(3)=P12CM*(-STH0*CPHICM*STHCM+CTH0*CTHCM) 
C          BOOST WITH P0 TO GET LAB MOMENTA 
      DO 240 L=1,2  
      SGN=3-2*L 
      BP=0  
      DO 241 K=1,3  
241   BP=BP+PJSET(K,J0)*SGN*P1CM(K) 
      BP=BP/AM0 
      PJSET(4,JJ(L))=PJSET(4,J0)*EE(L)/PJSET(5,J0)+BP   
      DO 242 K=1,3  
242   PJSET(K,JJ(L))=SGN*P1CM(K)+PJSET(K,J0)*EE(L)/PJSET(5,J0)  
     1+PJSET(K,J0)*BP/(PJSET(4,J0)+PJSET(5,J0)) 
240   CONTINUE  
C   
C          (3) PICK Z AND DECAY PARTONS. CHECK. 
C   
300   CONTINUE  
      TNEW=.FALSE.  
      DO 310 L=1,NJI    
      IF(JDCAY(JJ(L)).GE.0) GO TO 310   
      IF(NJSET+2.GT.MXJSET) GO TO 9999  
      CALL QCDZ(JJ(L))  
      CALL QCDT(NJSET+1)    
      CALL QCDT(NJSET+2)    
C   
C          CHECK WHETHER MASSES ALLOWED.    
      AM0=PJSET(5,JJ(L))    
      AM1=PJSET(5,NJSET+1)  
      AM2=PJSET(5,NJSET+2)  
      IF(AM1+AM2.GE.AM0) GO TO 320  
C   
C          CHECK WHETHER Z ALLOWED. 
      E1CM=(AM0**2+AM1**2-AM2**2)/(2*AM0)   
      E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)  
      P12CM=SQRT((AM0**2-AM1**2-AM2**2)**2-(2*AM1*AM2)**2)/(2*AM0)  
      E0=PJSET(4,JJ(L)) 
      P0=SQRT(PJSET(1,JJ(L))**2+PJSET(2,JJ(L))**2+PJSET(3,JJ(L))**2)    
      IF(ZZC(JJ(L)).GT.0.5) THEN    
        ZEP=ZZC(JJ(L))*(E0+P0)  
        P1=(ZEP**2-AM1**2)/(2.*ZEP) 
        E1=(ZEP**2+AM1**2)/(2.*ZEP) 
        CTHCM=(E1*AM0-E0*E1CM)/(P0*P12CM)   
        IF(ABS(CTHCM).GE.1..OR.P1.LE.0.) GO TO 320  
      ELSE  
        ZEP=(1.-ZZC(JJ(L)))*(E0+P0) 
        P2=(ZEP**2-AM2**2)/(2.*ZEP) 
        E2=(ZEP**2+AM2**2)/(2.*ZEP) 
        CTHCM=-(E2*AM0-E0*E2CM)/(P0*P12CM)  
        IF(ABS(CTHCM).GE.1..OR.P2.LE.0.) GO TO 320  
      ENDIF 
C   
C          REQUIRE Z AND 1-Z WITHIN KINEMATIC LIMITS.   
      ZLIM=(AM0/(E0+P0))**2 
      ZLIM1=CUTJET/(E0+P0)  
      ZLIM=AMAX1(ZLIM,ZLIM1)    
      IF(ZZC(JJ(L)).LE.ZLIM.OR.ZZC(JJ(L)).GE.(1.-ZLIM)) GO TO 320   
C   
C          ADD NEW PARTONS TO /JETSET/. 
      JDCAY(JJ(L))=JPACK*(NJSET+1)+(NJSET+2)    
      NJSET=NJSET+2 
      GO TO 310 
C   
C          DISCARD PARTONS AND EVOLVE JJ(L) AGAIN.  
320   TNEW=.TRUE.   
      CALL QCDT(JJ(L))  
310   CONTINUE  
C   
C          (4) RESOLVE KINEMATICS IF ANY PARTON MASS IS CHANGED.    
C   
      IF(TNEW) GO TO 210    
200   CONTINUE  
C   
C          (5) ITERATE ENTIRE PROCEEDURE.   
C   
      NJ1=NJ2+1 
      IF(NJ1.LE.NJSET) GO TO 1  
      RETURN    
C   
C          ERROR MESSAGE    
C   
9999  CALL PRTEVT(0)    
      WRITE(ITLIS,10) NJSET 
10    FORMAT(//' ERROR IN QCDJET...NJSET > ',I4)    
      RETURN    
      END   
      SUBROUTINE READIN(IFL)    
C   
C       READ IN USER DATA   
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (LIMPOM=20) 
      COMMON/MBGEN/POMWT(LIMPOM),POMGEN(LIMPOM),MNPOM,MXPOM,PDIFFR, 
     $NPOM,XBARY(2),DXBARY(2),XPOM(LIMPOM,2)    
      COMMON/FORCE/NFORCE,IFORCE(20),MFORCE(5,20)   
      COMMON/QCDPAR/ALAM,ALAM2,CUTJET,ISTRUC    
      COMMON/IDRUN/IDVER,IDG(2),IEVT    
      COMMON/FRGPAR/PUD,PBARY,SIGQT,PEND,XGEN(8),PSPIN1(8), 
     $PMIX1(3,2),PMIX2(3,2),XGENSS(9)   
      DIMENSION PMIXX1(6),PMIXX2(6) 
      EQUIVALENCE (PMIX1(1,1),PMIXX1(1))    
      EQUIVALENCE (PMIX2(1,1),PMIXX2(1))    
      DIMENSION FRPAR(32)   
      EQUIVALENCE(FRPAR(1),PUD) 
      COMMON/KEYS/IKEYS,KEYON,KEYS(10)  
      COMMON/XKEYS/REAC 
      LOGICAL KEYS  
      LOGICAL KEYON 
      CHARACTER*8 REAC  
      COMMON/PRTOUT/NEVPRT,NJUMP    
      COMMON/SEED/XSEED 
      CHARACTER*24 XSEED    
      COMMON/TYPES/LOC(50),NTYP,NJTTYP(3),NWWTYP(2) 
      COMMON/XTYPES/PARTYP(40),TITLE(10),JETYP(30,3),WWTYP(30,2)    
      CHARACTER*8 JETYP,WWTYP,TITLE,PARTYP  
      COMMON/PASCHAR/JETYP11,JETYP12,JETYP13,WWTYP11,WWTYP12  
      CHARACTER*8 JETYP11,JETYP12,JETYP13
      CHARACTER*8 WWTYP11,WWTYP12  
      COMMON/PRIMAR/NJET,SCM,HALFE,ECM,IDIN(2),NEVENT,NTRIES,NSIGMA 
      COMMON/JETLIM/PMIN(3),PMAX(3),PTMIN(3),PTMAX(3),YJMIN(3),YJMAX(3) 
     1 ,PHIMIN(3),PHIMAX(3),XJMIN(3),XJMAX(3),THMIN(3),THMAX(3) 
     2 ,SETLMJ(36)  
      DIMENSION BLIMS(36)   
      EQUIVALENCE(BLIMS(1),PMIN(1)) 
      LOGICAL SETLMJ    
      COMMON/FIXPAR/FIXP(3),FIXPT(3),FIXYJ(3),FIXPHI(3),FIXXJ(3)    
     2   ,FIXQM,FIXQT,FIXYW,FIXXW,FIXPHW    
      LOGICAL FIXQM,FIXQT,FIXYW,FIXXW,FIXPHW    
      LOGICAL FIXP,FIXPT,FIXYJ,FIXPHI,FIXXJ 
      COMMON/SGNPAR/CTHS(2,3),THS(2,3),YJS(2,3),XJS(2,3)    
      COMMON/NODCAY/NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR 
      LOGICAL NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR   
      COMMON/WCON/SIN2W,WMASS(4),WGAM(4),AQ(12,4),BQ(12,4),COUT(4), 
     1MATCH(25,4),WCBR(25,4),CUTOFF,CUTPOW,TBRWW(4,2),RBRWW(12,4,2),EZ, 
     2AQDP(12,4),BQDP(12,4),EZDP    
      DOUBLE PRECISION AQDP,BQDP,EZDP   
      COMMON/DYLIM/QMIN,QMAX,QTMIN,QTMAX,YWMIN,YWMAX,XWMIN,XWMAX,THWMIN,    
     2  THWMAX,PHWMIN,PHWMAX    
     3  ,SETLMQ(12) 
      LOGICAL SETLMQ    
      DIMENSION BLIM1(12)   
      EQUIVALENCE(BLIM1(1),QMIN)    
      COMMON/QLMASS/AMLEP(52),NQLEP,NMES,NBARY  
      COMMON/Q1Q2/GOQ(25,3),GOALL(3),GODY(4),STDDY,GOWW(25,2),ALLWW(2)  
      LOGICAL GOQ,GOALL,GODY,STDDY,GOWW,ALLWW   
      COMMON/JETPAR/P(3),PT(3),YJ(3),PHI(3),XJ(3),TH(3),CTH(3),STH(3)   
     1 ,JETTYP(3),SHAT,THAT,UHAT,QSQ,X1,X2,PBEAM(2) 
     2 ,QMW,QW,QTW,YW,XW,THW,QTMW,PHIW,SHAT1,THAT1,UHAT1,JWTYP  
     3 ,ALFQSQ,CTHW,STHW,Q0W    
     4 ,INITYP(2),ISIGS,PBEAMS(5)   
C   
      LOGICAL SETTYP,DUMY   
      CHARACTER*8 TTL(10),WORD,LSTRUC,BLANK 
      CHARACTER*8 WTYP(4),RDID(2)   
      DIMENSION NLAP(3,17)  
      DOUBLE PRECISION SEED 
C   
C          OVERLAPPING VARIABLE FLAGS.  
      DATA NPROC/6/ 
      DATA NLAP/1,2,3, 1,2,7 ,1,2,8, 1,3,5, 1,3,6, 1,3,7, 1,3,8, 1,5,7, 
     X          1,5,8, 1,6,7, 1,6,8, 2,3,7, 2,3,8, 3,5,7, 3,6,7, 3,5,8, 
     X          3,6,8/  
      DATA BLANK/'        '/    
C   
C          ENTRY    
      JETYP(1,1)=JETYP11
      JETYP(1,2)=JETYP12
      JETYP(1,3)=JETYP13
      WWTYP(1,1)=WWTYP11
      WWTYP(1,2)=WWTYP12  
      IFL=0 
      VERSN=IDVER/100.  
      WRITE(ITLIS,10) VERSN 
10    FORMAT('1',//5X,'*****  ISAJET VERSION ',F6.2,'  *****')  
      WRITE(ITLIS,11)   
   11 FORMAT(////30X,' COMMANDS READ BY READIN')    
C   
C          READ TITLE   
      READ(ITCOM,1) TTL 
    1 FORMAT(10A8)  
      WRITE(ITLIS,2) TTL    
    2 FORMAT(' ',10A8)  
      IF(TTL(1).EQ.'STOP    ') IFL=1001 
      IF(IFL.EQ.1001) GOTO 999  
C   
C          READ ENERGY AND NO. OF EVENTS    
      READ(ITCOM,*) ECM,NEVENT,NEVPRT,NJUMP 
      WRITE(ITLIS,*) ECM,NEVENT,NEVPRT,NJUMP    
C          RESET ALL VARIABLES IF TITLE IS NOT 'SAME'   
      IF(TTL(1).EQ.'SAME    ') GO TO 21 
      DO 20 I=1,10  
   20 TITLE(I)=TTL(I)   
      CALL RESET    
      KEYON=.FALSE. 
C   
C          READ REACTION AND SET KEYS AND NJET  
C   
      READ(ITCOM,3) REAC    
    3 FORMAT(A8)    
      WRITE(ITLIS,4) REAC   
   4  FORMAT(1X,A8) 
      DO 18 I=1,NPROC   
18    KEYS(I)=.FALSE.   
      KEYON=.FALSE. 
C   
C   
      IF(REAC.EQ.'TWOJET  ') THEN   
        KEYS(1)=.TRUE.  
        IKEY=1  
        NJET=2  
      ELSEIF(REAC.EQ.'E+E-    ') THEN   
        KEYS(2)=.TRUE.  
        IKEY=2  
        NJET=2  
      ELSEIF(REAC.EQ.'DRELLYAN') THEN   
        KEYS(3)=.TRUE.  
        IKEY=3  
        NJET=3  
      ELSEIF(REAC.EQ.'MINBIAS ') THEN   
        KEYS(4)=.TRUE.  
        IKEY=4  
        NJET=0  
      ELSEIF(REAC.EQ.'SUPERSYM'.OR.REAC.EQ.'SUSY    ') THEN 
        KEYS(5)=.TRUE.  
        IKEY=5  
        NJET=2  
      ELSEIF(REAC.EQ.'WPAIR   ') THEN   
        KEYS(6)=.TRUE.  
        IKEY=6  
        NJET=2  
      ELSE  
        KEYON=.FALSE.   
        GO TO 890   
      ENDIF 
C   
C   
   21 CONTINUE  
      SCM=ECM**2    
      HALFE=ECM/2   
      NSEL=0    
C   
C          READ KEYWORD. FOR EACH RECOGNIZED KEYWORD READ CORRESPONDING 
C          VARIABLES AND SET LOC FLAG.  
C   
      NSEL=0    
100   CONTINUE  
      READ(ITCOM,3) WORD    
      WRITE(ITLIS,4) WORD   
      NSEL=NSEL+1   
C          KEYWORD END  
      IF(WORD.EQ.'END     ') THEN   
        GO TO 101   
C          KEYWORD P    
      ELSEIF(WORD.EQ.'P       ') THEN   
        READ(ITCOM,*)  (PMIN(K),PMAX(K),K=1,NJET)   
        WRITE(ITLIS,*) (PMIN(K),PMAX(K),K=1,NJET)   
        LOC(1)=NSEL 
C          KEYWORD Y    
      ELSEIF(WORD.EQ.'Y       ') THEN   
        READ(ITCOM,*)  (YJMIN(K),YJMAX(K),K=1,NJET) 
        WRITE(ITLIS,*) (YJMIN(K),YJMAX(K),K=1,NJET) 
        LOC(2)=NSEL 
C          KEYWORD X    
      ELSEIF(WORD.EQ.'X       ') THEN   
        READ(ITCOM,*)  (XJMIN(K),XJMAX(K),K=1,NJET) 
        WRITE(ITLIS,*) (XJMIN(K),XJMAX(K),K=1,NJET) 
        LOC(3)=NSEL 
C          KEYWORD PHI  
      ELSEIF(WORD.EQ.'PHI     ') THEN   
        READ(ITCOM,*)  (PHIMIN(K),PHIMAX(K),K=1,NJET)   
        WRITE(ITLIS,*) (PHIMIN(K),PHIMAX(K),K=1,NJET)   
        LOC(4)=NSEL 
C          KEYWORD TH   
        ELSEIF(WORD.EQ.'TH      '.OR.WORD.EQ.'THETA   ') THEN   
        READ(ITCOM,*)  (THMIN(K),THMAX(K),K=1,NJET) 
        WRITE(ITLIS,*) (THMIN(K),THMAX(K),K=1,NJET) 
        LOC(5)=NSEL 
        LOC(6)=NSEL 
C          KEYWORD PT   
      ELSEIF(WORD.EQ.'PT      '.OR.WORD.EQ.'PPERP   ') THEN 
        READ(ITCOM,*)  (PTMIN(K),PTMAX(K),K=1,NJET) 
        WRITE(ITLIS,*) (PTMIN(K),PTMAX(K),K=1,NJET) 
        LOC(7)=NSEL 
        LOC(8)=NSEL 
C          KEYWORD NODECAY  
      ELSEIF(WORD.EQ.'NODECAY ') THEN   
        READ(ITCOM,571) NODCAY  
571     FORMAT(L1)  
        WRITE(ITLIS,572) NODCAY 
572     FORMAT(' ',L1)  
        LOC(9)=NSEL 
C          KEYWORD NOETA    
      ELSEIF(WORD.EQ.'NOETA   ') THEN   
        READ(ITCOM,571) NOETA   
        WRITE(ITLIS,572) NOETA  
        LOC(10)=NSEL    
C          KEYWORD NOPI0    
      ELSEIF(WORD.EQ.'NOPI0   ') THEN   
        READ(ITCOM,571) NOPI0   
        WRITE(ITLIS,572) NOPI0  
        LOC(11)=NSEL    
C          KEYWORD BEAMS    
      ELSEIF(WORD.EQ.'BEAMS   ') THEN   
        READ(ITCOM,*) RDID(1),RDID(2)   
        WRITE(ITLIS,*) RDID(1),RDID(2)  
        IDIN(1)=0   
        IDIN(2)=0   
        DO 583 K=1,2    
        IF(RDID(K).EQ.'P       ') IDIN(K)=+1120 
        IF(RDID(K).EQ.'AP      ') IDIN(K)=-1120 
        IF(RDID(K).EQ.'N       ') IDIN(K)=+1220 
        IF(RDID(K).EQ.'AN      ') IDIN(K)=-1220 
583     CONTINUE    
        IF(IDIN(1)*IDIN(2).EQ.0) GO TO 920  
        LOC(12)=NSEL    
C          KEYWORD FRAGMENT 
      ELSEIF(WORD.EQ.'FRAGMENT') THEN   
        READ(ITCOM,*)  FRPAR    
        WRITE(ITLIS,*)  FRPAR   
        LOC(13)=NSEL    
C          KEYWORD SEED 
      ELSEIF(WORD.EQ.'SEED    ') THEN   
        READ(ITCOM,*) SEED  
        WRITE(ITLIS,*) SEED 
CCCCC        CALL RANSET(SEED)   
        WRITE(XSEED,'(E24.15)') SEED    
        LOC(14)=NSEL    
C          KEYWORDS JETTYPE1, JETTYPE2, JETTYPE3    
      ELSEIF(WORD.EQ.'JETTYPE1'.OR.WORD.EQ.'JETTYPE2'.OR.WORD.EQ.   
     $'JETTYPE3') THEN  
        IF(WORD.EQ.'JETTYPE1') IJ=1 
        IF(WORD.EQ.'JETTYPE2') IJ=2 
        IF(WORD.EQ.'JETTYPE3') IJ=3 
        READ(ITCOM,*) (JETYP(K,IJ),K=1,25)  
        DO 601 K=1,25   
601     IF(JETYP(K,IJ).NE.BLANK) NJTTYP(IJ)=NJTTYP(IJ)+1    
        WRITE(ITLIS,*) (JETYP(K,IJ),K=1,NJTTYP(IJ)) 
        LOC(15)=NSEL    
C          KEYWORD SIN2W    
      ELSEIF(WORD.EQ.'SIN2W   ') THEN   
        READ(ITCOM,*)  SIN2W    
        WRITE(ITLIS,*) SIN2W    
        LOC(17)=NSEL    
C          KEYWORD TMASS    
      ELSEIF(WORD.EQ.'TMASS   ') THEN   
        READ(ITCOM,*)  AMLEP(6),AMLEP(7),AMLEP(8)   
        WRITE(ITLIS,*) AMLEP(6),AMLEP(7),AMLEP(8)   
        LOC(18)=NSEL    
C          KEYWORD QMW  
      ELSEIF(WORD.EQ.'QMW     ') THEN   
        READ(ITCOM,*) QMIN,QMAX 
        WRITE(ITLIS,*) QMIN,QMAX    
        LOC(19)=NSEL    
C          KEYWORD QTW  
      ELSEIF(WORD.EQ.'QTW     ') THEN   
        READ(ITCOM,*)  QTMIN,QTMAX  
        WRITE(ITLIS,*) QTMIN,QTMAX  
        LOC(20)=NSEL    
C          KEYWORD YW   
      ELSEIF(WORD.EQ.'YW      ') THEN   
        READ(ITCOM,*)  YWMIN,YWMAX  
        WRITE(ITLIS,*) YWMIN,YWMAX  
        LOC(21)=NSEL    
C          KEYWORD XW   
      ELSEIF(WORD.EQ.'XW      ') THEN   
        READ(ITCOM,*)  XWMIN,XWMAX  
        WRITE(ITLIS,*) XWMIN,XWMAX  
        LOC(22)=NSEL    
C          KEYWORD THW  
      ELSEIF(WORD.EQ.'THW     ') THEN   
        READ(ITCOM,*)  THWMIN,THWMAX    
        WRITE(ITLIS,*) THWMIN,THWMAX    
        LOC(23)=NSEL    
C          KEYWORD PHIW 
      ELSEIF(WORD.EQ.'PHIW    ') THEN   
        READ(ITCOM,*)  PHWMIN,PHWMAX    
        WRITE(ITLIS,*) PHWMIN,PHWMAX    
        LOC(24)=NSEL    
C          KEYWORD NONUNU   
      ELSEIF(WORD.EQ.'NONUNU  ') THEN   
        READ(ITCOM,571) NONUNU  
        WRITE(ITLIS,572) NONUNU 
        LOC(25)=NSEL    
C          KEYWORD WTYPE    
      ELSEIF(WORD.EQ.'WTYPE   ') THEN   
        DO 191 J=1,4    
        WTYP(J)=BLANK   
191     GODY(J)=.FALSE. 
        READ(ITCOM,*) WTYP  
        WRITE(ITLIS,*) WTYP 
        DO 192 K=1,4    
        IF(WTYP(K).EQ.'GM      ') GODY(1)=.TRUE.    
        IF(WTYP(K).EQ.'W+      ') GODY(2)=.TRUE.    
        IF(WTYP(K).EQ.'W-      ') GODY(3)=.TRUE.    
        IF(WTYP(K).EQ.'Z0      ') GODY(4)=.TRUE.    
  192   CONTINUE    
        IF(GODY(1)) JWTYP=1 
        IF(GODY(2).OR.GODY(3)) JWTYP=3  
        IF(GODY(4)) JWTYP=4 
        IF(GODY(3).AND.(GODY(1).OR.GODY(4))) GOTO 921   
        IF(GODY(2).AND.(GODY(1).OR.GODY(4))) GOTO 921   
        LOC(26)=NSEL    
C          KEYWORD LAMBDA   
      ELSEIF(WORD.EQ.'LAMBDA  ') THEN   
    
       READ(ITCOM,*)  ALAM  
        WRITE(ITLIS,*)  ALAM    
        ALAM2=ALAM**2   
        LOC(27)=NSEL    
C          KEYWORD NTRIES   
      ELSEIF(WORD.EQ.'NTRIES  ') THEN   
        READ(ITCOM,*) NTRIES    
        WRITE(ITLIS,*) NTRIES   
        LOC(28)=NSEL    
C          KEYWORD CUTOFF   
      ELSEIF(WORD.EQ.'CUTOFF  ') THEN   
        READ(ITCOM,*) CUTOFF,CUTPOW 
        WRITE(ITLIS,*) CUTOFF,CUTPOW    
        LOC(29)=NSEL    
C          KEYWORD XGEN 
      ELSEIF(WORD.EQ.'XGEN    ') THEN   
        READ(ITCOM,*) XGEN  
        WRITE(ITLIS,*) XGEN 
        LOC(30)=NSEL    
C          KEYWORD SIGQT    
      ELSEIF(WORD.EQ.'SIGQT   ') THEN   
        READ(ITCOM,*) SIGQT 
        WRITE(ITLIS,*) SIGQT    
        LOC(31)=NSEL    
C          KEYWORD CUTJET   
      ELSEIF(WORD.EQ.'CUTJET  ') THEN   
        READ(ITCOM,*)  CUTJET   
        WRITE(ITLIS,*)  CUTJET  
        LOC(32)=NSEL    
C          KEYWORD STRUC    
      ELSEIF(WORD.EQ.'STRUC   ') THEN   
        ISTRUC=0    
        READ(ITCOM,*) LSTRUC    
        WRITE(ITLIS,4) LSTRUC   
        IF(LSTRUC.EQ.'OWENS   ') ISTRUC=1   
        IF(LSTRUC.EQ.'BAIER   ') ISTRUC=2   
        IF(LSTRUC.EQ.'EICHTEN '.OR.LSTRUC.EQ.'EHLQ    ') ISTRUC=3   
        IF(LSTRUC.EQ.'DUKE    '.OR.LSTRUC.EQ.'DO      ') ISTRUC=4   
        IF(ISTRUC.EQ.0) GO TO 920   
        LOC(33)=NSEL    
C          KEYWORD NPOMERON 
      ELSEIF(WORD.EQ.'NPOMERON') THEN   
        READ(ITCOM,*) MNPOM,MXPOM   
        WRITE(ITLIS,*) MNPOM,MXPOM  
        IF(MNPOM.LT.1.OR.MNPOM.GT.MXPOM.OR.MXPOM.GT.LIMPOM) 
     1  GO TO 930   
        LOC(34)=NSEL    
C          KEYWORD FORCE    
      ELSEIF(WORD.EQ.'FORCE   ') THEN   
        NFORCE=NFORCE+1 
        IF(NFORCE.GT.20) GO TO 930  
        DO 721 K=1,5    
721     MFORCE(K,NFORCE)=0  
        READ(ITCOM,*) IFORCE(NFORCE),(MFORCE(K,NFORCE),K=1,5)   
        WRITE(ITLIS,*) IFORCE(NFORCE),(MFORCE(K,NFORCE),K=1,5)  
        LOC(35)=NSEL    
C          KEYWORD HMASS    
      ELSEIF(WORD.EQ.'HMASS   ') THEN   
        CALL FLAVOR(80,I1,I2,I3,J1,INDEX)   
        READ(ITCOM,*) (AMLEP(INDEX+K),K=1,9)    
        WRITE(ITLIS,*) (AMLEP(INDEX+K),K=1,9)   
        LOC(36)=NSEL    
C          KEYWORDS WMODE1,WMODE2   
      ELSEIF(WORD.EQ.'WMODE1  '.OR.WORD.EQ.'WMODE2  ') THEN 
        IF(WORD.EQ.'WMODE1  ') IJ=1 
        IF(WORD.EQ.'WMODE2  ') IJ=2 
        READ(ITCOM,*) (WWTYP(K,IJ),K=1,25)  
        DO 602 K=1,25   
602     IF(WWTYP(K,IJ).NE.BLANK) NWWTYP(IJ)=NWWTYP(IJ)+1    
        WRITE(ITLIS,*) (WWTYP(K,IJ),K=1,NWWTYP(IJ)) 
        LOC(37)=NSEL    
      ELSEIF(WORD.EQ.'NOEVOLVE') THEN   
        READ (ITCOM,571) NOEVOL 
        WRITE(ITLIS,572) NOEVOL 
        LOC(38)=NSEL    
      ELSEIF(WORD.EQ.'NOHADRON') THEN   
        READ (ITCOM,571) NOHADR 
        WRITE(ITLIS,572) NOHADR 
        LOC(39)=NSEL    
C          KEYWORD GAUGINO  
      ELSEIF(WORD.EQ.'GAUGINO ') THEN   
        CALL FLAVOR(29,I1,I2,I3,J1,IDG1)    
        CALL FLAVOR(30,I1,I2,I3,J1,IDG2)    
        CALL FLAVOR(39,I1,I2,I3,J1,IDG3)    
        CALL FLAVOR(40,I1,I2,I3,J1,IDG4)    
        READ(ITCOM,*) AMLEP(IDG1),AMLEP(IDG2),AMLEP(IDG3),AMLEP(IDG4)   
        WRITE(ITLIS,*) AMLEP(IDG1),AMLEP(IDG2),AMLEP(IDG3),AMLEP(IDG4)  
        LOC(40)=NSEL    
C          KEYWORD SQUARK   
      ELSEIF(WORD.EQ.'SQUARK  ') THEN   
        CALL FLAVOR(21,I1,I2,I3,J1,IDXQK)   
        READ(ITCOM,*) (AMLEP(IDXQK+K-1),K=1,6)  
        WRITE(ITLIS,*) (AMLEP(IDXQK+K-1),K=1,6) 
        LOC(41)=NSEL    
C          KEYWORD SLEPTON  
      ELSEIF(WORD.EQ.'SLEPTON ') THEN   
        CALL FLAVOR(31,I1,I2,I3,J1,IDXLEP)  
        READ(ITCOM,*) (AMLEP(IDXLEP+K-1),K=1,6) 
        WRITE(ITLIS,*) (AMLEP(IDXLEP+K-1),K=1,6)    
        LOC(42)=NSEL    
C          KEYWORD NSIGMA   
      ELSEIF(WORD.EQ.'NSIGMA  ') THEN   
        READ(ITCOM,*) NSIGMA    
        WRITE(ITLIS,*) NSIGMA   
        LOC(43)=NSEL    
      ELSEIF(WORD.EQ.'XGENSS  ') THEN   
        READ(ITCOM,*) XGENSS(9),(XGENSS(KK),KK=1,8) 
        WRITE(ITLIS,*) XGENSS(9),(XGENSS(KK),KK=1,8)    
        LOC(44)=NSEL    
C          NONE OF THE ABOVE    
      ELSE  
        GO TO 900   
      ENDIF 
C   
      GO TO 100 
C   
C          CHECK FOR ERRORS AND OVERLAPPING VARIABLES.  
C          PRINT ERROR MESSAGE OR WARNINGS  
  101 CONTINUE  
      IF(LOC(2)*LOC(5).NE.0) GOTO 910   
      IF(LOC(2)*LOC(6).NE.0) GOTO 910   
      IF(LOC(15).NE.0.OR.LOC(37).NE.0) DUMY=SETTYP(0)   
      DO 120 I=1,17 
      I1=NLAP(1,I)  
      I2=NLAP(2,I)  
      I3=NLAP(3,I)  
      IF(LOC(I1)*LOC(I2)*LOC(I3).NE.0) WRITE(ITLIS,1001)    
  120 CONTINUE  
 1001 FORMAT(//2X,'YOU HAVE GIVEN LIMITS FOR AN OVERLAPPING SET OF',    
     C  ' VARIABLES. SET MINIMIZING PPERP INTERVAL WILL BE USED')   
      GOTO 999  
  890 WRITE(ITLIS,1999) 
 1999 FORMAT(//' YOU FORGOT TO SELECT A PROCESS FOR GENERATION.'    
     X  /' AVAILABLE AT PRESENT ARE ',  
     X  /' TWOJET  , E+E-    , DRELLYAN, MINBIAS , WPAIR')  
      IFL=9 
      GOTO 999  
  900 WRITE(ITLIS,2000) WORD    
 2000 FORMAT(//2X,A8,' IS NOT A RECOGNIZABLE PARAMETER. JOB TERMINATED')    
      IFL=10    
      GOTO 999  
  910 WRITE(ITLIS,2001) 
 2001 FORMAT(//2X,' YOU CANNOT GIVE LIMITS FOR BOTH THETA AND Y. MAKE', 
     C  ' UP YOUR MIND. JOB TERMINATED')    
      IFL=11    
      GOTO 999  
  920 WRITE(ITLIS,2002) 
 2002 FORMAT(/'  WHAT IS THAT SUPPOSED TO BE')  
      IFL=12    
      GOTO 999  
  921 WRITE(ITLIS,2003) 
 2003 FORMAT(/'  YOU CANNOT RUN WS AND Z0 OR GAMMAS AT THE SAME TIME')  
      IFL=13    
      GOTO 999  
930   WRITE(ITLIS,2004) 
2004  FORMAT(//'  PARAMETER OUT OF RANGE. JOB TERMINATED.') 
      IFL=14    
 999  RETURN    
      END   
      SUBROUTINE RESCAL(N1,N2,PSUM,IFAIL)   
C          RESCALE MOMENTA OF PARTICLES N1...N2 TO GIVE TOTAL   
C          FOUR-MOMENTUM PSUM.  
C          RETURN IFAIL=0 IF OK, IFAIL=1 IF NO GOOD.    
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (MXPTCL=4000,IPACK=10000)   
      COMMON/PARTCL/NPTCL,PPTCL(5,MXPTCL),IORIG(MXPTCL),IDENT(MXPTCL)   
     1,IDCAY(MXPTCL)    
      DIMENSION PSUM(5),PADD(5),BETA(3) 
      DATA ERRLIM/.0001/    
C          ORIGIONAL MOMENTUM IS PADD.  
      IFAIL=1   
      IF(N1.GE.N2) RETURN   
      DO 100 K=1,5  
100   PADD(K)=0.    
      DO 110 IP=N1,N2   
      DO 110 K=1,5  
      PADD(K)=PADD(K)+PPTCL(K,IP)   
110   CONTINUE  
      IF(PADD(5).GE.PSUM(5)) RETURN 
      PADD(5)=PADD(4)**2-PADD(1)**2-PADD(2)**2-PADD(3)**2   
      IF(PADD(5).LE.0) RETURN   
      PADD(5)=SQRT(PADD(5)) 
      DO 120 K=1,3  
120   BETA(K)=-PADD(K)/PADD(5)  
      GAMMA=PADD(4)/PADD(5) 
C          BOOST PARTICLES TO REST. 
200   CONTINUE  
      DO 210 IP=N1,N2   
      BP=0. 
      DO 220 K=1,3  
220   BP=BP+PPTCL(K,IP)*BETA(K) 
      DO 230 K=1,3  
230   PPTCL(K,IP)=PPTCL(K,IP)+BETA(K)*PPTCL(4,IP)   
     $+BETA(K)*BP/(GAMMA+1.)    
      PPTCL(4,IP)=GAMMA*PPTCL(4,IP)+BP  
210   CONTINUE  
      IF(IFAIL.EQ.0) RETURN 
C          RESCALE MOMENTA IN REST FRAME.   
      SCAL=1.   
      DO 301 IPASS=1,200    
      SUM=0.    
      DO 310 IP=N1,N2   
      DO 320 K=1,3  
320   PPTCL(K,IP)=SCAL*PPTCL(K,IP)  
      PPTCL(4,IP)=SQRT(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(3,IP)**2 
     $+PPTCL(5,IP)**2)  
      SUM=SUM+PPTCL(4,IP)   
310   CONTINUE  
      SCAL=PSUM(5)/SUM  
301   IF(ABS(SCAL-1.).LE.ERRLIM) GO TO 300  
300   CONTINUE  
C          BOOST BACK WITH PSUM.    
      BMAG=0.   
      DO 400 K=1,3  
      BETA(K)=PSUM(K)/PSUM(5)   
      BMAG=BMAG+ABS(BETA(K))    
400   CONTINUE  
      GAMMA=PSUM(4)/PSUM(5) 
      IFAIL=0   
      IF(BMAG.EQ.0.) RETURN 
      GO TO 200 
      END   
      SUBROUTINE RESET  
C          RESET ALL USER DEFINED VARIABLES 
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/QCDPAR/ALAM,ALAM2,CUTJET,ISTRUC    
      COMMON/DYLIM/QMIN,QMAX,QTMIN,QTMAX,YWMIN,YWMAX,XWMIN,XWMAX,THWMIN,    
     2  THWMAX,PHWMIN,PHWMAX    
     3  ,SETLMQ(12) 
      LOGICAL SETLMQ    
      DIMENSION BLIM1(12)   
      EQUIVALENCE(BLIM1(1),QMIN)    
      COMMON/FRGPAR/PUD,PBARY,SIGQT,PEND,XGEN(8),PSPIN1(8), 
     $PMIX1(3,2),PMIX2(3,2),XGENSS(9)   
      DIMENSION PMIXX1(6),PMIXX2(6) 
      EQUIVALENCE (PMIX1(1,1),PMIXX1(1))    
      EQUIVALENCE (PMIX2(1,1),PMIXX2(1))    
      DIMENSION FRPAR(32)   
      EQUIVALENCE(FRPAR(1),PUD) 
      COMMON/JETLIM/PMIN(3),PMAX(3),PTMIN(3),PTMAX(3),YJMIN(3),YJMAX(3) 
     1 ,PHIMIN(3),PHIMAX(3),XJMIN(3),XJMAX(3),THMIN(3),THMAX(3) 
     2 ,SETLMJ(36)  
      DIMENSION BLIMS(36)   
      EQUIVALENCE(BLIMS(1),PMIN(1)) 
      LOGICAL SETLMJ    
      COMMON/FIXPAR/FIXP(3),FIXPT(3),FIXYJ(3),FIXPHI(3),FIXXJ(3)    
     2   ,FIXQM,FIXQT,FIXYW,FIXXW,FIXPHW    
      LOGICAL FIXQM,FIXQT,FIXYW,FIXXW,FIXPHW    
      LOGICAL FIXP,FIXPT,FIXYJ,FIXPHI,FIXXJ 
      COMMON/SGNPAR/CTHS(2,3),THS(2,3),YJS(2,3),XJS(2,3)    
      COMMON/JETPAR/P(3),PT(3),YJ(3),PHI(3),XJ(3),TH(3),CTH(3),STH(3)   
     1 ,JETTYP(3),SHAT,THAT,UHAT,QSQ,X1,X2,PBEAM(2) 
     2 ,QMW,QW,QTW,YW,XW,THW,QTMW,PHIW,SHAT1,THAT1,UHAT1,JWTYP  
     3 ,ALFQSQ,CTHW,STHW,Q0W    
     4 ,INITYP(2),ISIGS,PBEAMS(5)   
      COMMON/NODCAY/NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR 
      LOGICAL NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR   
      COMMON/PRIMAR/NJET,SCM,HALFE,ECM,IDIN(2),NEVENT,NTRIES,NSIGMA 
      COMMON/QLMASS/AMLEP(52),NQLEP,NMES,NBARY  
      COMMON/Q1Q2/GOQ(25,3),GOALL(3),GODY(4),STDDY,GOWW(25,2),ALLWW(2)  
      LOGICAL GOQ,GOALL,GODY,STDDY,GOWW,ALLWW   
      COMMON/SEED/XSEED 
      CHARACTER*24 XSEED    
      COMMON/TYPES/LOC(50),NTYP,NJTTYP(3),NWWTYP(2) 
      COMMON/XTYPES/PARTYP(40),TITLE(10),JETYP(30,3),WWTYP(30,2)    
      CHARACTER*8 JETYP,WWTYP,TITLE,PARTYP  
      COMMON/WCON/SIN2W,WMASS(4),WGAM(4),AQ(12,4),BQ(12,4),COUT(4), 
     1MATCH(25,4),WCBR(25,4),CUTOFF,CUTPOW,TBRWW(4,2),RBRWW(12,4,2),EZ, 
     2AQDP(12,4),BQDP(12,4),EZDP    
      DOUBLE PRECISION AQDP,BQDP,EZDP   
      COMMON/FORCE/NFORCE,IFORCE(20),MFORCE(5,20)   
      PARAMETER (LIMPOM=20) 
      COMMON/MBGEN/POMWT(LIMPOM),POMGEN(LIMPOM),MNPOM,MXPOM,PDIFFR, 
     $NPOM,XBARY(2),DXBARY(2),XPOM(LIMPOM,2)    
      CHARACTER*8 BLANK 
      DATA BLANK/'        '/    
      DATA UNDEF/-1.E9/ 
C          RESET DYLIM  
      DO 110 I=1,12 
      BLIM1(I)=UNDEF    
      SETLMQ(I)=.TRUE.  
110   CONTINUE  
C          RESET FRGPAR 
      PUD=.43   
      PBARY=.10 
      SIGQT=.35 
      PEND=.14  
      XGEN(1)=.96   
      XGEN(2)=3.    
      XGEN(3)=0.    
      XGEN(4)=.8    
      XGEN(5)=.5    
      XGEN(6)=.5    
      XGEN(7)=.5    
      XGEN(8)=.5    
      DO 111 K=1,9  
111   XGENSS(K)=.5  
      PSPIN1(1)=.5  
      PSPIN1(2)=.5  
      PSPIN1(3)=.5  
      PSPIN1(4)=.75 
      PSPIN1(5)=.75 
      PSPIN1(6)=.75 
      PSPIN1(7)=.75 
      PSPIN1(8)=.75 
      PMIXX1(1)=.25 
      PMIXX1(2)=.25 
      PMIXX1(3)=.5  
      PMIXX1(4)=0.  
      PMIXX1(5)=.5  
      PMIXX1(6)=1.  
      PMIXX2(1)=.5  
      PMIXX2(2)=.5  
      PMIXX2(3)=1.  
      PMIXX2(4)=0.  
      PMIXX2(5)=0.  
      PMIXX2(6)=1.  
C          RESET JETLIM 
      DO 120 I=1,36 
      BLIMS(I)=UNDEF    
      SETLMJ(I)=.TRUE.  
120   CONTINUE  
C          RESET NODCAY 
      NODCAY=.FALSE.    
      NOETA=.FALSE. 
      NOPI0=.FALSE. 
      NONUNU=.FALSE.    
      NOEVOL=.FALSE.    
      NOHADR=.FALSE.    
C          RESET PRIMAR 
      IDIN(1)=1120  
      IDIN(2)=1120  
      NTRIES=200    
      NSIGMA=20 
C          RESET QCDPAR 
      ALAM=.2   
      ALAM2=ALAM**2 
      CUTJET=6. 
      ISTRUC=3  
C          RESET QLMASS 
      AMLEP(6)=30.  
      AMLEP(7)=-1.  
      AMLEP(8)=-1.  
      DO 125 I=1,9  
      CALL FLAVOR(80+I,I1,I2,I3,J1,INDEX)   
125   AMLEP(INDEX)=0.   
      CALL FLAVOR(29,I1,I2,I3,J1,INDEX) 
      AMLEP(INDEX)=100. 
      CALL FLAVOR(30,I1,I2,I3,J1,INDEX) 
      AMLEP(INDEX)=0.   
      CALL FLAVOR(39,I1,I2,I3,J1,INDEX) 
      AMLEP(INDEX)=100. 
      CALL FLAVOR(40,I1,I2,I3,J1,INDEX) 
      AMLEP(INDEX)=100. 
      DO 126 I=1,6  
      CALL FLAVOR(20+I,I1,I2,I3,J1,INDEX)   
      AMLEP(INDEX)=100.+AMASS(I)    
      CALL FLAVOR(30+I,I1,I2,I3,J1,INDEX)   
      AMLEP(INDEX)=100.+AMASS(I+10) 
126   CONTINUE  
C          RESET Q1Q2   
      DO 130 I=1,25 
      DO 130 J=1,3  
130   GOQ(I,J)=.TRUE.   
      DO 131 I=1,3  
131   GOALL(I)=.TRUE.   
      GODY(1)=.TRUE.    
      GODY(2)=.FALSE.   
      GODY(3)=.FALSE.   
      GODY(4)=.TRUE.    
      DO 132 I=1,2  
      ALLWW(I)=.TRUE.   
      DO 132 J=1,25 
132   GOWW(J,I)=.TRUE.  
C          RESET SEED   
      ISEED=0   
C          RESET TYPES  
      DO 140 I=1,NTYP   
140   LOC(I)=0  
      DO 141 I=1,3  
      NJTTYP(I)=0   
      JETYP(1,I)='ALL     ' 
      DO 141 K=2,25 
141   JETYP(K,I)=BLANK  
      JWTYP=4   
      DO 142 I=1,2  
      NWWTYP(I)=0   
      WWTYP(1,I)='ALL     ' 
      DO 142 K=2,4  
142   WWTYP(K,I)=BLANK  
C          RESET WCON   
      SIN2W=.215    
      CUTOFF=.200   
      CUTPOW=1.0    
C          RESET MBGEN  
      MNPOM=1   
      MXPOM=LIMPOM  
C          RESET FORCE  
      NFORCE=0  
      RETURN    
      END   
      SUBROUTINE SETCON 
C          THIS SUBROUTINE SETS THE CONSTANTS IN /CONST/.   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/CONST/PI,SQRT2,ALFA,GF,UNITS   
      PI=4.*ATAN(1.)    
      SQRT2=SQRT(2.)    
      ALFA=1./137.036   
      GF=1.16570E-5 
      UNITS=1./2.56815  
      RETURN    
      END   
      SUBROUTINE SETDKY(LPRINT) 
C          THIS SUBROUTINE READS IN THE DECAY TABLE FROM TAPE ITDKY 
C          AND SETS UP /DKYTAB/.    
C          QUARK-BASED IDENT CODE   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/FORCE/NFORCE,IFORCE(20),MFORCE(5,20)   
C          LOOK MUST BE DIMENSIONED TO THE MAXIMUM VALUE OF INDEX.  
      PARAMETER (MXLOOK=400)    
      PARAMETER (MXDKY=1000)    
      COMMON/DKYTAB/LOOK(MXLOOK),CBR(MXDKY),MODE(5,MXDKY)   
      COMMON/NODCAY/NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR 
      LOGICAL NODCAY,NOETA,NOPI0,NONUNU,NOEVOL,NOHADR   
      DIMENSION IMODE(6)    
      CHARACTER*8 LABEL,LMODE(6),LRES   
      CHARACTER*8 IQUIT,IBLANK,LREAD(10)    
      LOGICAL LPRINT    
      DATA IQUIT/'////'/,IBLANK/' '/    
      IF(LPRINT) WRITE(ITLIS,10)    
10    FORMAT('1',30('*')/' *',28X,'*'/  
     1' *',5X,'ISAJET DECAY TABLE',5X,'*'/  
     2' *',28X,'*'/' ',30('*')//    
     36X,'PART',18X,'DECAY MODE',19X,'CUM BR',15X,'IDENT',17X,  
     4'DECAY IDENT'/)   
      LOOP=0    
      IOLD=0    
      DO 100 I=1,400    
      LOOK(I)=0 
100   CONTINUE  
      IF(NODCAY) RETURN 
      REWIND ITDKY  
200   LOOP=LOOP+1   
      IF(LOOP.GT.MXDKY) GO TO 9999  
220   DO 210 I=1,5  
      IMODE(I)=0    
      LMODE(I)=IBLANK   
210   CONTINUE  
      READ(ITDKY,*) IRES,ITYPE,BR,IMODE 
      IF(IRES.EQ.0) GO TO 300   
      IF(NOPI0.AND.IRES.EQ.110) GO TO 220   
      IF(NOETA.AND.IRES.EQ.220) GO TO 220   
      IF(IRES.EQ.IOLD) GO TO 230    
      CALL FLAVOR(IRES,IFL1,IFL2,IFL3,JSPIN,INDEX)  
      LOOK(INDEX)=LOOP  
230   IOLD=IRES 
      CBR(LOOP)=BR  
      DO 240 I=1,5  
      MODE(I,LOOP)=IMODE(I) 
      IF(IMODE(I).NE.0) LMODE(I)=LABEL(IMODE(I))    
240   CONTINUE  
      LRES=LABEL(IRES)  
      IF(LPRINT) WRITE(ITLIS,20) LRES,(LMODE(K),K=1,5), 
     1BR,IRES,(IMODE(K),K=1,5)  
20    FORMAT(6X,A5,6X,5(A5,2X),3X,F8.5,15X,I5,4X,5(I5,2X))  
      GO TO 200 
C          SET FORCED DECAY MODES   
300   IF(NFORCE.EQ.0) GO TO 400 
      DO 310 I=1,NFORCE 
      LOOP=LOOP+1   
      IF(LOOP.GT.MXDKY) GO TO 9999  
      CALL FLAVOR(IFORCE(I),IFL1,IFL2,IFL3,JSPIN,INDEX) 
      LOOK(INDEX)=LOOP  
      DO 320 K=1,5  
320   MODE(K,LOOP)=MFORCE(K,I)  
      CBR(LOOP)=1.  
310   CONTINUE  
C          READ AND PRINT NOTES FROM DECAYTABLE FILE    
400   IF(.NOT.LPRINT) RETURN    
      WRITE(ITLIS,1001) 
1001  FORMAT('1',24('*')/' *',22X,'*'/  
     1' *',5X,'ISAJET NOTES',5X,'*'/    
     2' *',22X,'*'/' ',24('*')//)   
410   READ(ITDKY,1002) LREAD    
1002  FORMAT(10A8)  
      IF(LREAD(1).EQ.IQUIT) RETURN  
      WRITE(ITLIS,1003) LREAD   
1003  FORMAT(1X,10A8)   
      GO TO 410 
9999  WRITE(ITLIS,3001) LOOP    
3001  FORMAT(//' ***** ERROR IN SETDKY ... LOOP = ',I6,' *****')    
      CALL EXIT 
      END   
      LOGICAL FUNCTION SETTYP(LPRT) 
C   
C          SET JETTYPE FLAGS AND WMODE FLAGS FOR WPAIR  
C          NOTE THAT FLAGS ARE NOW CHARACTER*8  
C   
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      COMMON/KEYS/IKEYS,KEYON,KEYS(10)  
      COMMON/XKEYS/REAC 
      LOGICAL KEYS  
      LOGICAL KEYON 
      CHARACTER*8 REAC  
      COMMON/TYPES/LOC(50),NTYP,NJTTYP(3),NWWTYP(2) 
      COMMON/XTYPES/PARTYP(40),TITLE(10),JETYP(30,3),WWTYP(30,2)    
      CHARACTER*8 JETYP,WWTYP,TITLE,PARTYP  
      COMMON/Q1Q2/GOQ(25,3),GOALL(3),GODY(4),STDDY,GOWW(25,2),ALLWW(2)  
      LOGICAL GOQ,GOALL,GODY,STDDY,GOWW,ALLWW   
C   
      CHARACTER*8 WORD,BLANK,LIST(29),LISTW(4),LISTXY(4)    
      DATA BLANK/'        '/    
      DATA LIST/'GL','UP','UB','DN','DB','ST','SB','CH','CB','BT','BB', 
     $'TP','TB','NUE','ANUE','E-','E+','NUM','ANUM','MU-','MU+',    
     $'NUT','ANUT','TAU-','TAU+','ALL','QUARKS','LEPTONS','NUS'/    
      DATA LISTW/'GM','W+','W-','Z0'/   
      DATA LISTXY/'Y','YB','X','XB'/    
C   
C          JETTYPE FLAGS FOR TWOJET, E+E-, DRELLYAN 
C          NOTE THAT NJTTYP IS SET IN READIN TO THE NUMBER OF NON-BLANK 
C          VALUES READ. 
C   
      SETTYP=.FALSE.    
C   
C   
C          CHECK FOR LEGAL JET TYPE NAMES AND SET APPROPRIATE FLAGS 
C   
      IF(KEYS(6)) GO TO 6   
C   
      DO 1000 JET=1,3   
      IF(NJTTYP(JET).EQ.0) GO TO 1000   
C           INITIALIZE EVERYTHING TO .FALSE.    
      GOALL(JET)=.FALSE.    
      DO 1100 K=1,25    
      GOQ(K,JET)=.FALSE.    
1100  CONTINUE  
C          LOOP OVER NON-BLANK JETTYPE ENTRIES  
      DO 1200 I=1,NJTTYP(JET)   
      WORD=JETYP(I,JET) 
C         BLANK 
      IF(WORD.EQ.BLANK) THEN    
        GO TO 1200  
      ENDIF 
C          ALL  
      IF(WORD.EQ.'ALL     ') THEN   
        GOALL(JET)=.TRUE.   
        DO 1210 K=1,25  
1210    GOQ(K,JET)=.TRUE.   
        GO TO 1000  
      ENDIF 
C          QUARKS   
      IF(WORD.EQ.'QUARKS  ') THEN   
        DO 1220 K=2,13  
1220    GOQ(K,JET)=.TRUE.   
        GO TO 1200  
      ENDIF 
C          CHARGED LEPTONS  
      IF(WORD.EQ.'LEPTONS ') THEN   
        DO 1230 K=16,24,4   
        GOQ(K,JET)=.TRUE.   
1230    GOQ(K+1,JET)=.TRUE. 
        GO TO 1200  
      ENDIF 
C          NEUTRINOS    
      IF(WORD.EQ.'NUS     ') THEN   
        DO 1240 K=14,22,4   
        GOQ(K,JET)=.TRUE.   
1240    GOQ(K+1,JET)=.TRUE. 
        GO TO 1200  
      ENDIF 
C          EXPLICIT TYPES   
      DO 1250 K=1,25    
      IF(WORD.EQ.LIST(K)) THEN  
        GOQ(K,JET)=.TRUE.   
        GO TO 1200  
      ENDIF 
1250  CONTINUE  
      DO 1260 K=1,4 
      IF(KEYS(5).AND.WORD.EQ.LISTW(K)) THEN 
        GOQ(K+13,JET)=.TRUE.    
        GO TO 1200  
      ENDIF 
1260  CONTINUE  
      DO 1270 K=1,4 
      IF(KEYS(1).AND.WORD.EQ.LISTXY(K)) THEN    
        GOQ(13+K,JET)=.TRUE.    
        GO TO 1200  
      ENDIF 
1270  CONTINUE  
C          ERROR    
      WRITE(ITLIS,1300) WORD,JET    
1300  FORMAT(1X,A8,' IS NOT RECOGNIZABLE FOR JETTYPE',I1)   
      SETTYP=.TRUE. 
1200  CONTINUE  
1000  CONTINUE  
      RETURN    
C   
C          JETTYPE AND WMODE FLAGS FOR WPAIR    
C          NOTE THAT NJTTYP AND NWWTYP ARE THE NUMBER OF NON-BLANK  
C          VALUES.  
C   
6     DO 2000 JET=1,2   
      IF(NJTTYP(JET).EQ.0) GO TO 2300   
C          INITIALIZE TO FALSE  
      GOALL(JET)=.FALSE.    
      DO 2100 K=1,4 
2100  GOQ(K,JET)=.FALSE.    
C   
C          LOOP OVER NON-BLANK JETTYPE FLAGS    
C   
      DO 2200 I=1,NJTTYP(JET)   
      WORD=JETYP(I,JET) 
C          BLANK    
      IF(WORD.EQ.BLANK) THEN    
        GO TO 2200  
      ENDIF 
C          ALL  
      IF(WORD.EQ.'ALL     ') THEN   
        GOALL(JET)=.TRUE.   
        DO 2210 K=1,4   
2210    GOQ(K,JET)=.TRUE.   
        GO TO 2300  
      ENDIF 
C          EXPLICIT TYPES   
      DO 2220 K=1,4 
      IF(WORD.EQ.LISTW(K)) THEN 
        GOQ(K,JET)=.TRUE.   
        GO TO 2200  
      ENDIF 
2220  CONTINUE  
C          ERROR    
      WRITE(ITLIS,1300) WORD,JET    
      SETTYP=.TRUE. 
2200  CONTINUE  
C   
C          LOOP OVER NONBLANK WMODE FLAGS   
C   
2300  IF(NWWTYP(JET).EQ.0) GO TO 2000   
      ALLWW(JET)=.FALSE.    
C         INITIALIZE EVERYTHING TO FALSE    
      DO 2350 K=1,25    
2350  GOWW(K,JET)=.FALSE.   
C   
      DO 2400 I=1,NWWTYP(JET)   
      WORD=WWTYP(I,JET) 
      IF(WORD.NE.BLANK) NWWTYP(JET)=I   
C          BLANK    
      IF(WORD.EQ.BLANK) THEN    
      GO TO 2400    
      ENDIF 
C          ALL  
      IF(WORD.EQ.'ALL     ') THEN   
      ALLWW(JET)=.TRUE. 
      DO 2410 K=1,25    
2410  GOWW(K,JET)=.TRUE.    
      GO TO 2000    
      ENDIF 
C          QUARKS   
      IF(WORD.EQ.'QUARKS  ') THEN   
      DO 2420 K=2,13    
2420  GOWW(K,JET)=.TRUE.    
      GO TO 2400    
      ENDIF 
C          CHARGED LEPTONS  
      IF(WORD.EQ.'LEPTONS ') THEN   
      DO 2430 K=16,24,4 
      GOWW(K,JET)=.TRUE.    
2430  GOWW(K+1,JET)=.TRUE.  
      GO TO 2400    
      ENDIF 
C          NEUTRINOS    
      IF(WORD.EQ.'NUS     ') THEN   
      DO 2440 K=14,22,4 
      GOWW(K,JET)=.TRUE.    
2440  GOWW(K+1,JET)=.TRUE.  
      GO TO 2400    
      ENDIF 
C          EXPLICIT TYPES   
      DO 2450 K=1,25    
      IF(WORD.EQ.LIST(K)) THEN  
      GOWW(K,JET)=.TRUE.    
      GO TO 2400    
      ENDIF 
2450  CONTINUE  
C          ERROR    
      WRITE(ITLIS,2500) WORD,JET    
2500  FORMAT(1X,A8,' IS NOT A VALID CODE FOR WMODE',I1) 
      SETTYP=.TRUE. 
2400  CONTINUE  
2000  CONTINUE  
C   
      RETURN    
      END   
      SUBROUTINE HEAVYX(X,EPS)  
C   
C          GENERATE X FOR HEAVY PARTICLE FRAGMENTATION ACCORDING TO 
C          THE PETERSON FORM    
C          D(X)=1/(X*(1-1/X-EPS/(1-X))**2)  
C              =D0(X)*D1(X)*D2(X)   
C          D0(X)=(1-X)**2/((1-X)**2+EPS)**2 
C          D1(X)=X  
C          D2(X)=(((1-X)**2+EPS)/((1-X)**2+EPS*X))**2   
C          USING X=1-Y**POW 
C   
      DATA ALN4/1.3863/ 
C   
C          CHOOSE POW FOR X=1-Y**POW.   
C          GENERATE FLAT IN X IF EPS>1. 
      IF(EPS.LT.1.) THEN    
        POW=ALOG((3.+EPS)/EPS)/ALN4 
        YMX=(EPS*(3.*POW-1.)/(POW+1.))**(.5/POW)    
        ZMX=1-YMX**POW  
        D0MX=(1-ZMX)**2/((1.-ZMX)**2+EPS)**2*POW*YMX**(POW-1.)  
        D2MX=2./(2.-SQRT(EPS))  
      ELSE  
        POW=1.  
        ZMX=0.  
        D0MX=(1.-ZMX)**2/((1.-ZMX)**2+EPS)**2   
        D2MX=1.+EPS 
      ENDIF 
C   
C          GENERATE Z ACCORDING TO (1-Z)**2/((1-Z)**2+EPS*Z)**2 
1     CONTINUE  
      Y=RANF()  
      Z=1.-Y**POW   
C   
      D0Z=(1.-Z)**2/((1.-Z)**2+EPS)**2*POW*Y**(POW-1.)  
      IF(D0Z.LT.RANF()*D0MX) GO TO 1    
C   
C          CHECK REMAINING FACTORS  
      D1=Z  
      D2=(((1.-Z)**2+EPS)/((1.-Z)**2+EPS*Z))**2 
      IF(D1*D2.LT.RANF()*D2MX) GO TO 1  
C   
C          GOOD X   
      X=Z   
      RETURN    
      END   
      SUBROUTINE QCDT(J)    
C          AUXILIARY ROUTINE FOR QCDJET.  CALCULATE ZC AND STORE    
C          IN ZZC(J).  GENERATE NEW MASS WITH THIS ZC AND STORE 
C          IN PJSET(5,J).   
C          MUST INCLUDE 1/2 SYMMETRY FACTOR IN GAMGG. NO FIX IS NEEDED  
C          IN QCDZ SINCE GAMGG+2*GAMQQ IS USED AS THE NORMALIZATION.    
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      COMMON/JWORK/ZZC(MXJSET),JMATCH(MXJSET),TNEW,P1CM(4), 
     1J1,J2,J3,J4,J5,E1CM,E2CM,E3CM,E4CM,E5CM   
      LOGICAL TNEW  
      DIMENSION JJ(5),EE(5) 
      EQUIVALENCE (J1,JJ(1)),(E1CM,EE(1))   
      COMMON/QCDPAR/ALAM,ALAM2,CUTJET,ISTRUC    
C          CALCULATE ZC 
      AM0=PJSET(5,J)    
      JTLV1=JTYPE(J)    
      AM1=AMASS(JTLV1)+CUTJET   
      AM2=CUTJET    
      IF(AM1+AM2.GE.AM0) GO TO 300  
      T0=AM0**2 
      T1=AM1**2 
      T2=AM2**2 
      ZC=(T0-T1+T2-SQRT((T0-T1-T2)**2-4*T1*T2))/(2*T0)  
      ZZC(J)=ZC 
C          COUNT LIGHT FERMIONS 
      NF=3  
      DO 110 IQ=4,6 
      IF(AM0.LT.2*AMASS(IQ)) GO TO 120  
      NF=NF+1   
110   CONTINUE  
120   B0=11.-2.*NF/3.   
C          CALCULATE GAMMA(ZC)  
C   
      JTABS=IABS(JTYPE(J))  
C          INITIAL GLUON    
      IF(JTABS.EQ.9) THEN   
        GAMQQ=(1.-2.*ZC)*(1.-ZC*(1.-ZC))/3. 
        GAMGG=12.*ALOG((1.-ZC)/ZC)-9.*(1.-2.*ZC)-6.*GAMQQ   
        GAMGG=0.5*GAMGG 
        GAM=GAMGG+NF*GAMQQ  
C          INITIAL QUARK OR DIQUARK 
      ELSEIF(JTABS.LT.9.OR.MOD(JTABS,100).EQ.0) THEN    
        GAM=8./3.*ALOG((1-ZC)/ZC)-2.*(1.-2.*ZC) 
C          INITIAL GLUINO   
      ELSEIF(JTABS.EQ.29) THEN  
        GAM=6.*ALOG((1.-ZC)/ZC)-9./2.*(1.-2.*ZC)    
C          INITIAL SQUARK   
      ELSEIF(JTABS.GT.20.AND.JTABS.LT.29) THEN  
        GAM = 8./3.*(ALOG((1.-ZC)/ZC)-(1.-2.*ZC))   
      ENDIF 
C          GENERATE NEW MASS    
200   GB=2*GAM/B0   
      PROB=(ALOG(AM1/ALAM)/ALOG(AM0/ALAM))**GB  
      IF(PROB.GT.RANF()) GO TO 300  
      POW=(1-(1-PROB)*RANF())**(1./GB)  
      AMNEW=ALAM*(AM0/ALAM)**POW    
      IF(AMNEW.LE.AM1) GO TO 300    
      PJSET(5,J)=AMNEW  
      RETURN    
C          FINAL PARTON 
300   PJSET(5,J)=AM1-CUTJET 
      JDCAY(J)=0    
      RETURN    
      END   
      SUBROUTINE QCDZ(J)    
C          AUXILIARY ROUTINE FOR QCDJET.  GENERATE Z FOR PARTON J   
C          AND STORE IN ZZC(J).  ADD POSSIBLE NEW PARTONS TO /JETSET/.  
      COMMON/ITAPES/ITDKY,ITEVT,ITCOM,ITLIS 
      PARAMETER (MXJSET=400,JPACK=1000) 
      COMMON/JETSET/NJSET,PJSET(5,MXJSET),JORIG(MXJSET),JTYPE(MXJSET),  
     $JDCAY(MXJSET) 
      COMMON/JWORK/ZZC(MXJSET),JMATCH(MXJSET),TNEW,P1CM(4), 
     1J1,J2,J3,J4,J5,E1CM,E2CM,E3CM,E4CM,E5CM   
      LOGICAL TNEW  
      DIMENSION JJ(5),EE(5) 
      EQUIVALENCE (J1,JJ(1)),(E1CM,EE(1))   
      COMMON/QCDPAR/ALAM,ALAM2,CUTJET,ISTRUC    
C          ALTARELLI-PARISI FUNCTIONS.  
      PQQ(Z)=4*(1+Z**2)/(3*(1-Z))   
      PGQ(Z)=.5*(Z**2+(1-Z)**2) 
      PGG(Z)=6*(1-Z*(1-Z))**2/(Z*(1-Z)) 
      PGSGS(Z)=3.*(1.+Z**2)/(1.-Z)  
      PQSQS(Z)=8./3.*Z/(1.-Z)   
C          INITIALIZE.  
      AM0=PJSET(5,J)    
      ZC=ZZC(J) 
C   
      JTABS=IABS(JTYPE(J))  
C          INITIAL GLUON    
      IF (JTABS.EQ.9) THEN  
200     CONTINUE    
        GAMQQ=(1-2*ZC)*(1-ZC*(1-ZC))    
        GAMGG=12*ALOG((1-ZC)/ZC)-9*(1-2*ZC)-6*GAMQQ 
        NF=3    
        DO 210 IQ=4,6   
        IF(AM0.LT.2*AMASS(IQ))GO TO 220 
210     NF=NF+1 
220     PROBG=GAMGG/(GAMGG+2*NF*GAMQQ)  
        PROBQ=GAMQQ/(GAMGG+2*NF*GAMQQ)  
        RND=RANF()  
        IF(PROBG.LT.RND) GO TO 240  
C          GL--->GL+GL  
230     ZGEN=(ZC/(1-ZC))**(1-2*RANF())  
        Z=ZGEN/(1.+ZGEN)    
        GZ=6./(Z*(1.-Z))    
        IF(PGG(Z).LT.GZ*RANF()) GO TO 230   
        JTYPE(NJSET+1)=9    
        JTYPE(NJSET+2)=9    
        ZZC(J)=Z    
        GO TO 300   
C          GL--->QK+QB  
240     Z=RANF()    
        IF(PGQ(Z).LT.0.5*RANF()) GO TO 240  
        IFL=(RND-PROBG)/PROBQ+1.    
        IF(IFL.GT.NF) IFL=NF-IFL    
        JTYPE(NJSET+1)=IFL  
        JTYPE(NJSET+2)=-IFL 
        ZZC(J)=Z    
C          INITIAL QUARK    
      ELSEIF (JTABS.LT.9.OR.MOD(JTABS,100).EQ.0) THEN   
100     CONTINUE    
        Z=1-(ZC/(1-ZC))**RANF()*(1-ZC)  
        GZ=8./(3.*(1-Z))    
        IF(PQQ(Z).LT.GZ*RANF()) GO TO 100   
        JTYPE(NJSET+1)=JTYPE(J) 
        JTYPE(NJSET+2)=9    
        ZZC(J)=Z    
C          INITIAL GLUINO   
       ELSEIF (JTABS.EQ.29) THEN    
400      Z=1.-(ZC/(1.-ZC))**RANF()*(1.-ZC)  
         GZ=6./(1.-Z)   
         IF(PGSGS(Z) .LT. GZ*RANF()) GOTO 400   
         JTYPE(NJSET+1)=JTYPE(J)    
         JTYPE(NJSET+2)=9   
         ZZC(J)=Z   
C          INITIAL SQUARK   
      ELSEIF(JTABS.GT.20.AND.JTABS.LT.29) THEN  
500     CONTINUE    
        Z=1-(ZC/(1-ZC))**RANF()*(1-ZC)  
        GZ=8./(3.*(1-Z))    
        IF(PQSQS(Z).LT.GZ*RANF()) GO TO 500 
        JTYPE(NJSET+1)=JTYPE(J) 
        JTYPE(NJSET+2)=9    
        ZZC(J)=Z    
      ENDIF 
C   
C          SET MASSES AND FLAGS.    
300   CONTINUE  
      PJSET(5,NJSET+1)=Z*AM0    
      PJSET(5,NJSET+2)=(1.-Z)*AM0   
      JET=IABS(JORIG(J))/JPACK  
      JORIG(NJSET+1)=JPACK*JET+J    
      JORIG(NJSET+2)=JPACK*JET+J    
      JDCAY(NJSET+1)=-1 
      JDCAY(NJSET+2)=-1 
      JMATCH(NJSET+1)=NJSET+2   
      JMATCH(NJSET+2)=NJSET+1   
      RETURN    
      END   

