[u/mrichter/AliRoot.git] / ISAJET / code / decjet.F

#include "isajet/pilot.h"
      LOGICAL FUNCTION DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA)
C
C          Auxiliary routine for DECAY. Evolve and hadronize partons.
C          Check conservation laws. Return TRUE if OK, FALSE otherwise.
C
C          IP = particle to be decayed.
C          NADD = number of products (NPTCL+1, ..., NPTCL+NADD).
C          IDABS = absolute values of decay IDENT's.
C          PREST = 4-momenta in rest frame.
C          WDECAY = logical flag for real W decay.
C          BETA,GAMMA = boost parameters.
C
#if defined(CERNLIB_IMPNONE)
      IMPLICIT NONE
#endif
#include "isajet/itapes.inc"
#include "isajet/wcon.inc"
#include "isajet/partcl.inc"
#include "isajet/dkytab.inc"
#include "isajet/jetset.inc"
#include "isajet/jwork.inc"
#include "isajet/const.inc"
C   
      REAL PGEN(5,5),RND(5),U(3),BETA(3),IDQK(3),ROT(3,3),PSAVE(3)  
     1,REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME,CHARGE    
      REAL PSUM(5),POLD(4),PNEW(4),SUM,WTMAX,SUM1,SUM2  
      REAL PREST(4,6),PWREST(5),BETAW(3),DOT,PCM    
      REAL AMEE,REE,WTEE,SWAP,RNEW,WT,QCM,PHI,S12,S12MAX,GAMMAW,BP  
      REAL PJET,CTHQK,STHQK,CPHIQK,SPHIQK,SUMQ,A,B,C,GAMMA  
      REAL CHARGW   
      LOGICAL WDECAY    
      INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2,I3  
      INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5),NEXT    
      INTEGER JJ1,II,K1,K,NJSAVE,NJSAV1,NJSAV2,NJ1,NPRTN,KK,NHDRN1  
      INTEGER IFAIL,JSAVE,JETIP,JET,NJADD,NPTLV1,IDANTI,NPJET(5)    
      INTEGER NHDRN,NPJET3,NPTCLW,NPBEG(5)  
C   
C          Copy decay products into /JETSET/ and do QCD evolution.  
C   
      IF(NJSET+NADD.GT.MXJSET) GO TO 9998   
      NJSAVE=NJSET  
      NSTART=NPTCL-NADD+1   
      NPTCL=NSTART-1    
      DO 100 I=1,NADD   
        NJSET=NJSET+1   
        DO 110 K=1,4    
110     PJSET(K,NJSET)=PREST(K,I)   
        PJSET(5,NJSET)=PPTCL(5,NPTCL+I) 
        JORIG(NJSET)=JPACK*I    
        JTYPE(NJSET)=IDENT(NPTCL+I) 
        JDCAY(NJSET)=0  
        JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD  
100   CONTINUE  
C   
C          For heavy quarks match 1+2 and 3+(1+2). Boost 1+2 to rest.   
C   
      IF(WDECAY) THEN   
        JMATCH(NJSAVE+1)=NJSAVE+2   
        JMATCH(NJSAVE+2)=NJSAVE+1   
        NJSET=NJSET+1   
        DO 120 K=1,4    
          PWREST(K)=PJSET(K,NJSAVE+1)+PJSET(K,NJSAVE+2) 
          PJSET(K,NJSET)=PWREST(K)  
120     CONTINUE    
        PWREST(5)=SQRT(PWREST(4)**2-PWREST(1)**2-PWREST(2)**2   
     $  -PWREST(3)**2)  
        PJSET(5,NJSET)=PWREST(5)    
        JMATCH(NJSAVE+3)=NJSAVE+4   
        JMATCH(NJSAVE+4)=NJSAVE+3   
        JORIG(NJSAVE+4)=-1  
        IDLV1=JTYPE(NJSAVE+1)   
        CHARGW=CHARGE(IDLV1)    
        IDLV1=JTYPE(NJSAVE+2)   
        CHARGW=CHARGW+CHARGE(IDLV1) 
        JTYPE(NJSAVE+4)=80*SIGN(1.,CHARGW)  
        JDCAY(NJSAVE+4)=0   
C          Boost W vectors to rest. 
        DO 130 K=1,3    
130     BETAW(K)=PWREST(K)/PWREST(4)    
        GAMMAW=PWREST(4)/PWREST(5)  
        NJSAV1=NJSAVE+1 
        NJSAV2=NJSAVE+2 
        DO 140 J=NJSAV1,NJSAV2  
          BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)+BETAW(3)*PJSET(3,J)    
          DO 141 K=1,3  
141       PJSET(K,J)=PJSET(K,J)-GAMMAW*BETAW(K)*(PJSET(4,J) 
     $    -BP*GAMMAW/(GAMMAW+1.))   
          PJSET(4,J)=GAMMAW*(PJSET(4,J)-BP) 
140     CONTINUE    
      ENDIF 
C   
C          Do evolution and save new W momentum. Start from parent  
C          mass or NADD*energy. 
      NJSAV1=NJSAVE+1   
      DO 150 J=NJSAV1,NJSET 
        IF(IABS(JTYPE(J)).LT.10.OR.MOD(JTYPE(J),100).EQ.0) THEN 
          JDCAY(J)=-1   
          PJSET(5,J)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,J)) 
        ENDIF   
150   CONTINUE  
C   
      CALL QCDJET(NJSAVE+1) 
C   
      IF(WDECAY) THEN   
        PWREST(4)=PJSET(4,NJSAVE+4) 
        GAMMAW=PWREST(4)/PWREST(5)  
        DO 200 K=1,3    
          PWREST(K)=PJSET(K,NJSAVE+4)   
          BETAW(K)=PWREST(K)/PWREST(4)  
200     CONTINUE    
      ENDIF 
C   
C          Put final partons in particle table - temporary IORIG.   
C          Also include virtual or real W momentum for quark decays.    
C   
      NJ1=NJSAVE+1  
      IF(WDECAY) THEN   
C          Real or virtual W.   
        NPTCL=NPTCL+1   
        NPTCLW=NPTCL    
        DO 210 K=1,5    
210     PPTCL(K,NPTCL)=PJSET(K,NJSAVE+4)    
        IORIG(NPTCL)=IP 
        IDENT(NPTCL)=JTYPE(NJSAVE+4)    
        IDCAY(NPTCL)=0  
C          Jet 3    
        NPBEG(3)=NPTCL+1    
        DO 220 J=NJ1,NJSET  
          IF(JDCAY(J).NE.0) GO TO 220   
          IF(JORIG(J)/JPACK.NE.3) GO TO 220 
          NPTCL=NPTCL+1 
          DO 221 K=1,5  
221       PPTCL(K,NPTCL)=PJSET(K,J) 
          IORIG(NPTCL)=3*IPACK+IP   
          IDENT(NPTCL)=JTYPE(J) 
          IDCAY(NPTCL)=0    
220     CONTINUE    
C          Jets 1 and 2 
        NPJET3=NPTCL    
        DO 230 JET=1,2  
          NPBEG(JET)=NPTCL+1    
          DO 240 J=NJ1,NJSET    
            IF(JDCAY(J).NE.0) GO TO 240 
            IF(JORIG(J)/JPACK.NE.JET) GO TO 240 
            NPTCL=NPTCL+1   
            BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)  
     $      +BETAW(3)*PJSET(3,J)    
            DO 241 K=1,3    
241         PPTCL(K,NPTCL)=PJSET(K,J)+GAMMAW*BETAW(K)*(PJSET(4,J)   
     $      +BP*GAMMAW/(GAMMAW+1.)) 
            PPTCL(4,NPTCL)=GAMMAW*(PJSET(4,J)+BP)   
            PPTCL(5,NPTCL)=PJSET(5,J)   
            IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+NPTCLW  
            IDENT(NPTCL)=JTYPE(J)   
            IDCAY(NPTCL)=0  
240       CONTINUE  
230     CONTINUE    
C          Quark decays to W plus jet 3; then W decays. 
        IDCAY(IP)=IPACK*NPTCLW+NPJET3   
        IDCAY(NPTCLW)=IPACK*(NPJET3+1)+NPTCL    
      ELSE  
C          Not quark decay, so just copy partons.   
        DO 250 JET=1,NADD   
          NPBEG(JET)=NPTCL+1    
          DO 260 J=NJ1,NJSET    
            IF(JDCAY(J).NE.0) GO TO 260 
            IF(JORIG(J)/JPACK.NE.JET) GO TO 260 
            NPTCL=NPTCL+1   
            DO 261 K=1,5    
261         PPTCL(K,NPTCL)=PJSET(K,J)   
            IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+IP  
            IDENT(NPTCL)=JTYPE(J)   
            IDCAY(NPTCL)=0  
260       CONTINUE  
250     CONTINUE    
        IDCAY(IP)=NSTART*IPACK+NPTCL    
      ENDIF 
      NHDRN=NPTCL   
C   
C          Hadronize quarks and rotate to proper angles.    
C   
      DO 300 JET=1,NADD 
        NPRTN=NPBEG(JET)-1  
        DO 310 I=NJ1,NJSET  
          IF(JDCAY(I).NE.0) GO TO 310   
          IF(JORIG(I)/JPACK.NE.JET) GO TO 310   
          NPRTN=NPRTN+1 
          IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0)   
     $    GO TO 330 
C   
C          Fragment parton: 
          NEXT=NPTCL+1  
          PJET=SQRT(PJSET(1,I)**2+PJSET(2,I)**2+PJSET(3,I)**2)  
          CTHQK=PJSET(3,I)/PJET 
          STHQK=1.-CTHQK**2
          IF(STHQK.LT.1) THEN
            STHQK=SQRT(STHQK)
            CPHIQK=PJSET(1,I)/(PJET*STHQK)
            SPHIQK=PJSET(2,I)/(PJET*STHQK)
          ELSE
            STHQK=0
            CPHIQK=1
            SPHIQK=0
          ENDIF
          CALL JETGEN(I)    
          IF(NEXT.GT.NPTCL) GO TO 310   
          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    
C   
          DO 320 II=NEXT,NPTCL  
            DO 321 K=1,3    
              PSAVE(K)=PPTCL(K,II)  
              PPTCL(K,II)=0.    
321         CONTINUE    
            DO 322 K=1,3    
            DO 322 KK=1,3   
322         PPTCL(K,II)=PPTCL(K,II)+ROT(K,KK)*PSAVE(KK) 
            IORIG(II)=IPACK*JET+NPRTN   
            IDCAY(II)=0 
320       CONTINUE  
          IDCAY(NPRTN)=NEXT*IPACK+NPTCL 
          GO TO 310 
C   
C          or add lepton:   
330       NPTCL=NPTCL+1 
          DO 331 K=1,5  
331       PPTCL(K,NPTCL)=PJSET(K,I) 
          IORIG(NPTCL)=IPACK*JET+NPRTN  
          IDENT(NPTCL)=JTYPE(I) 
          IDCAY(NPTCL)=0    
          IDCAY(NPRTN)=NPTCL*IPACK+NPTCL    
310     CONTINUE    
        NPJET(JET)=NPTCL    
300   CONTINUE  
C   
C          Reset NJSET so decay jets do not appear in /JETSET/  
      NJADD=NJSET   
      NJSET=NJSAVE  
C   
C          Check for at least two particles 
      IF(NPTCL.LT.NHDRN+2) THEN 
        NPTCL=NSTART-1  
        DECJET=.FALSE.  
        RETURN  
      ENDIF 
C   
C          Conserve charge  
C   
      SUMQ=0.   
      NHDRN1=NHDRN+1    
      DO 400 I=NHDRN1,NPTCL 
        IDLV1=IDENT(I)  
        SUMQ=SUMQ+CHARGE(IDLV1) 
400   CONTINUE  
      IDLV1=IDENT(IP)   
      SUMQ=SUMQ-CHARGE(IDLV1)   
C   
      IF(SUMQ.EQ.0.) GO TO 500  
C   
C          Charge wrong--fix it by swapping UP and DN quarks.   
      DO 410 I=NHDRN1,NPTCL 
        ID1=IDENT(I)    
        IF(IABS(ID1).GT.1000) GO TO 410 
        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.) THEN 
          IDENT(I)=ISIGN(200+10*I2+I3,ID1)  
        ELSEIF(I1.EQ.2.AND.I2.GT.2.AND.SUMQ*ID1.LT.0.) THEN 
          IDENT(I)=ISIGN(100+10*I2+I3,ID1)  
        ELSEIF(I1.EQ.1.AND.I2.EQ.2.AND.SUMQ*ID1.GT.0.) THEN 
          IDENT(I)=110+I3   
        ELSEIF(I1.EQ.1.AND.I2.EQ.1) THEN    
          IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ))    
        ELSE    
          GO TO 410 
        ENDIF   
        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   
     $  +PPTCL(5,I)**2) 
C          Sum cannot vanish for fractionally charged initial particle. 
        IF(ABS(SUMQ).LT.1.) GO TO 500   
410   CONTINUE  
C          Failed to conserve charge.   
      NPTCL=NSTART-1    
      DECJET=.FALSE.    
      RETURN    
C   
C          Rescale momenta for correct mass 
C   
500   CONTINUE  
      IF(WDECAY) THEN   
C          Quark decay. First rescale jet3 + W  
        DO 510 K=1,5    
510     PPTCL(K,NPTCL+1)=PPTCL(K,NPTCLW)    
        NPTLV1=NPTCL+1  
        DO 520 K=1,3    
520     PSUM(K)=0.  
        PSUM(4)=PPTCL(5,IP) 
        PSUM(5)=PSUM(4) 
        CALL RESCAL(NPJET(2)+1,NPTLV1,PSUM,IFAIL) 
        IF(IFAIL.NE.0) THEN 
          NPTCL=NSTART-1    
          DECJET=.FALSE.    
          RETURN    
        ENDIF   
        DO 530 K=1,3    
530     BETAW(K)=PPTCL(K,NPTCL+1)/PPTCL(4,NPTCL+1)  
        GAMMAW=PPTCL(4,NPTCL+1)/PPTCL(5,NPTCL+1)    
C          Then rescale W   
        PSUM(4)=PPTCL(5,NPTCLW) 
        PSUM(5)=PSUM(4) 
        CALL RESCAL(NHDRN1,NPJET(2),PSUM,IFAIL) 
        IF(IFAIL.NE.0) THEN 
          NPTCL=NSTART-1    
          DECJET=.FALSE.    
          RETURN    
        ENDIF   
      ELSE  
C          General decay with no W. 
        DO 550 K=1,3    
550     PSUM(K)=0.  
        PSUM(4)=PPTCL(5,IP) 
        PSUM(5)=PSUM(4) 
        NPTLV1=NPTCL    
        CALL RESCAL(NHDRN1,NPTLV1,PSUM,IFAIL)   
        IF(IFAIL.NE.0) THEN 
          NPTCL=NSTART-1    
          DECJET=.FALSE.    
          RETURN    
        ENDIF   
      ENDIF 
C   
C          Boost back to lab frame. Reset IORIG.    
C   
      IF(WDECAY) THEN   
        DO 600 I=NHDRN1,NPTCL  
          JET=IORIG(I)/IPACK    
          IF(JET.NE.1.AND.JET.NE.2) GO TO 600   
          BP=BETAW(1)*PPTCL(1,I)+BETAW(2)*PPTCL(2,I)+BETAW(3)*PPTCL(3,I)    
          DO 610 J=1,3  
610       PPTCL(J,I)=PPTCL(J,I)+GAMMAW*BETAW(J)*(PPTCL(4,I) 
     $    +BP*GAMMAW/(GAMMAW+1.))   
          PPTCL(4,I)=GAMMAW*(PPTCL(4,I)+BP) 
600     CONTINUE    
      ENDIF 
C   
      DO 620 I=NSTART,NPTCL 
        IORIG(I)=MOD(IORIG(I),IPACK)    
        BP=BETA(1)*PPTCL(1,I)+BETA(2)*PPTCL(2,I)+BETA(3)*PPTCL(3,I) 
        DO 621 J=1,3    
          PPTCL(J,I)=PPTCL(J,I)+GAMMA*BETA(J)*(PPTCL(4,I)   
     $    +BP*GAMMA/(GAMMA+1.)) 
621     CONTINUE    
        PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP)    
620   CONTINUE  
C   
C          Normal exit  
C   
      DECJET=.TRUE. 
      RETURN    
C   
C          Error messages.  
C   
9998  DECJET=.FALSE.
      CALL PRTEVT(0)    
      WRITE(ITLIS,99980) NJSET  
99980 FORMAT(//5X,'ERROR IN DECJET...NJSET > ',I5)  
      RETURN    
      END
Commit	Line	Data
0795afa3	1	#include "isajet/pilot.h"
	2	LOGICAL FUNCTION DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA)
	3	C
	4	C Auxiliary routine for DECAY. Evolve and hadronize partons.
	5	C Check conservation laws. Return TRUE if OK, FALSE otherwise.
	6	C
	7	C IP = particle to be decayed.
	8	C NADD = number of products (NPTCL+1, ..., NPTCL+NADD).
	9	C IDABS = absolute values of decay IDENT's.
	10	C PREST = 4-momenta in rest frame.
	11	C WDECAY = logical flag for real W decay.
	12	C BETA,GAMMA = boost parameters.
	13	C
	14	#if defined(CERNLIB_IMPNONE)
	15	IMPLICIT NONE
	16	#endif
	17	#include "isajet/itapes.inc"
	18	#include "isajet/wcon.inc"
	19	#include "isajet/partcl.inc"
	20	#include "isajet/dkytab.inc"
	21	#include "isajet/jetset.inc"
	22	#include "isajet/jwork.inc"
	23	#include "isajet/const.inc"
	24	C
	25	REAL PGEN(5,5),RND(5),U(3),BETA(3),IDQK(3),ROT(3,3),PSAVE(3)
	26	1,REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME,CHARGE
	27	REAL PSUM(5),POLD(4),PNEW(4),SUM,WTMAX,SUM1,SUM2
	28	REAL PREST(4,6),PWREST(5),BETAW(3),DOT,PCM
	29	REAL AMEE,REE,WTEE,SWAP,RNEW,WT,QCM,PHI,S12,S12MAX,GAMMAW,BP
	30	REAL PJET,CTHQK,STHQK,CPHIQK,SPHIQK,SUMQ,A,B,C,GAMMA
	31	REAL CHARGW
	32	LOGICAL WDECAY
	33	INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2,I3
	34	INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5),NEXT
	35	INTEGER JJ1,II,K1,K,NJSAVE,NJSAV1,NJSAV2,NJ1,NPRTN,KK,NHDRN1
	36	INTEGER IFAIL,JSAVE,JETIP,JET,NJADD,NPTLV1,IDANTI,NPJET(5)
	37	INTEGER NHDRN,NPJET3,NPTCLW,NPBEG(5)
	38	C
	39	C Copy decay products into /JETSET/ and do QCD evolution.
	40	C
	41	IF(NJSET+NADD.GT.MXJSET) GO TO 9998
	42	NJSAVE=NJSET
	43	NSTART=NPTCL-NADD+1
	44	NPTCL=NSTART-1
	45	DO 100 I=1,NADD
	46	NJSET=NJSET+1
	47	DO 110 K=1,4
	48	110 PJSET(K,NJSET)=PREST(K,I)
	49	PJSET(5,NJSET)=PPTCL(5,NPTCL+I)
	50	JORIG(NJSET)=JPACK*I
	51	JTYPE(NJSET)=IDENT(NPTCL+I)
	52	JDCAY(NJSET)=0
	53	JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD
	54	100 CONTINUE
	55	C
	56	C For heavy quarks match 1+2 and 3+(1+2). Boost 1+2 to rest.
	57	C
	58	IF(WDECAY) THEN
	59	JMATCH(NJSAVE+1)=NJSAVE+2
	60	JMATCH(NJSAVE+2)=NJSAVE+1
	61	NJSET=NJSET+1
	62	DO 120 K=1,4
	63	PWREST(K)=PJSET(K,NJSAVE+1)+PJSET(K,NJSAVE+2)
	64	PJSET(K,NJSET)=PWREST(K)
65	120 CONTINUE
66	PWREST(5)=SQRT(PWREST(4)2-PWREST(1)2-PWREST(2)**2
67	$ -PWREST(3)**2)
68	PJSET(5,NJSET)=PWREST(5)
69	JMATCH(NJSAVE+3)=NJSAVE+4
70	JMATCH(NJSAVE+4)=NJSAVE+3
71	JORIG(NJSAVE+4)=-1
72	IDLV1=JTYPE(NJSAVE+1)
73	CHARGW=CHARGE(IDLV1)
74	IDLV1=JTYPE(NJSAVE+2)
75	CHARGW=CHARGW+CHARGE(IDLV1)
76	JTYPE(NJSAVE+4)=80*SIGN(1.,CHARGW)
77	JDCAY(NJSAVE+4)=0
78	C Boost W vectors to rest.
79	DO 130 K=1,3
80	130 BETAW(K)=PWREST(K)/PWREST(4)
81	GAMMAW=PWREST(4)/PWREST(5)
82	NJSAV1=NJSAVE+1
83	NJSAV2=NJSAVE+2
84	DO 140 J=NJSAV1,NJSAV2
85	BP=BETAW(1)PJSET(1,J)+BETAW(2)PJSET(2,J)+BETAW(3)*PJSET(3,J)
86	DO 141 K=1,3
87	141 PJSET(K,J)=PJSET(K,J)-GAMMAWBETAW(K)(PJSET(4,J)
88	$ -BP*GAMMAW/(GAMMAW+1.))
89	PJSET(4,J)=GAMMAW*(PJSET(4,J)-BP)
90	140 CONTINUE
91	ENDIF
92	C
93	C Do evolution and save new W momentum. Start from parent
94	C mass or NADD*energy.
95	NJSAV1=NJSAVE+1
96	DO 150 J=NJSAV1,NJSET
97	IF(IABS(JTYPE(J)).LT.10.OR.MOD(JTYPE(J),100).EQ.0) THEN
98	JDCAY(J)=-1
99	PJSET(5,J)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,J))
100	ENDIF
101	150 CONTINUE
102	C
103	CALL QCDJET(NJSAVE+1)
104	C
105	IF(WDECAY) THEN
106	PWREST(4)=PJSET(4,NJSAVE+4)
107	GAMMAW=PWREST(4)/PWREST(5)
108	DO 200 K=1,3
109	PWREST(K)=PJSET(K,NJSAVE+4)
110	BETAW(K)=PWREST(K)/PWREST(4)
111	200 CONTINUE
112	ENDIF
113	C
114	C Put final partons in particle table - temporary IORIG.
115	C Also include virtual or real W momentum for quark decays.
116	C
117	NJ1=NJSAVE+1
118	IF(WDECAY) THEN
119	C Real or virtual W.
120	NPTCL=NPTCL+1
121	NPTCLW=NPTCL
122	DO 210 K=1,5
123	210 PPTCL(K,NPTCL)=PJSET(K,NJSAVE+4)
124	IORIG(NPTCL)=IP
125	IDENT(NPTCL)=JTYPE(NJSAVE+4)
126	IDCAY(NPTCL)=0
127	C Jet 3
128	NPBEG(3)=NPTCL+1
129	DO 220 J=NJ1,NJSET
130	IF(JDCAY(J).NE.0) GO TO 220
131	IF(JORIG(J)/JPACK.NE.3) GO TO 220
132	NPTCL=NPTCL+1
133	DO 221 K=1,5
134	221 PPTCL(K,NPTCL)=PJSET(K,J)
135	IORIG(NPTCL)=3*IPACK+IP
136	IDENT(NPTCL)=JTYPE(J)
137	IDCAY(NPTCL)=0
138	220 CONTINUE
139	C Jets 1 and 2
140	NPJET3=NPTCL
141	DO 230 JET=1,2
142	NPBEG(JET)=NPTCL+1
143	DO 240 J=NJ1,NJSET
144	IF(JDCAY(J).NE.0) GO TO 240
145	IF(JORIG(J)/JPACK.NE.JET) GO TO 240
146	NPTCL=NPTCL+1
147	BP=BETAW(1)PJSET(1,J)+BETAW(2)PJSET(2,J)
148	$ +BETAW(3)*PJSET(3,J)
149	DO 241 K=1,3
150	241 PPTCL(K,NPTCL)=PJSET(K,J)+GAMMAWBETAW(K)(PJSET(4,J)
151	$ +BP*GAMMAW/(GAMMAW+1.))
152	PPTCL(4,NPTCL)=GAMMAW*(PJSET(4,J)+BP)
153	PPTCL(5,NPTCL)=PJSET(5,J)
154	IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+NPTCLW
155	IDENT(NPTCL)=JTYPE(J)
156	IDCAY(NPTCL)=0
157	240 CONTINUE
158	230 CONTINUE
159	C Quark decays to W plus jet 3; then W decays.
160	IDCAY(IP)=IPACK*NPTCLW+NPJET3
161	IDCAY(NPTCLW)=IPACK*(NPJET3+1)+NPTCL
162	ELSE
163	C Not quark decay, so just copy partons.
164	DO 250 JET=1,NADD
165	NPBEG(JET)=NPTCL+1
166	DO 260 J=NJ1,NJSET
167	IF(JDCAY(J).NE.0) GO TO 260
168	IF(JORIG(J)/JPACK.NE.JET) GO TO 260
169	NPTCL=NPTCL+1
170	DO 261 K=1,5
171	261 PPTCL(K,NPTCL)=PJSET(K,J)
172	IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+IP
173	IDENT(NPTCL)=JTYPE(J)
174	IDCAY(NPTCL)=0
175	260 CONTINUE
176	250 CONTINUE
177	IDCAY(IP)=NSTART*IPACK+NPTCL
178	ENDIF
179	NHDRN=NPTCL
180	C
181	C Hadronize quarks and rotate to proper angles.
182	C
183	DO 300 JET=1,NADD
184	NPRTN=NPBEG(JET)-1
185	DO 310 I=NJ1,NJSET
186	IF(JDCAY(I).NE.0) GO TO 310
187	IF(JORIG(I)/JPACK.NE.JET) GO TO 310
188	NPRTN=NPRTN+1
189	IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0)
190	$ GO TO 330
191	C
192	C Fragment parton:
193	NEXT=NPTCL+1
194	PJET=SQRT(PJSET(1,I)2+PJSET(2,I)2+PJSET(3,I)**2)
195	CTHQK=PJSET(3,I)/PJET
196	STHQK=1.-CTHQK**2
197	IF(STHQK.LT.1) THEN
198	STHQK=SQRT(STHQK)
199	CPHIQK=PJSET(1,I)/(PJET*STHQK)
200	SPHIQK=PJSET(2,I)/(PJET*STHQK)
201	ELSE
202	STHQK=0
203	CPHIQK=1
204	SPHIQK=0
205	ENDIF
206	CALL JETGEN(I)
207	IF(NEXT.GT.NPTCL) GO TO 310
208	ROT(1,1)=CPHIQK*CTHQK
209	ROT(2,1)=SPHIQK*CTHQK
210	ROT(3,1)=-STHQK
211	ROT(1,2)=-SPHIQK
212	ROT(2,2)=CPHIQK
213	ROT(3,2)=0.
214	ROT(1,3)=CPHIQK*STHQK
215	ROT(2,3)=SPHIQK*STHQK
216	ROT(3,3)=CTHQK
217	C
218	DO 320 II=NEXT,NPTCL
219	DO 321 K=1,3
220	PSAVE(K)=PPTCL(K,II)
221	PPTCL(K,II)=0.
222	321 CONTINUE
223	DO 322 K=1,3
224	DO 322 KK=1,3
225	322 PPTCL(K,II)=PPTCL(K,II)+ROT(K,KK)*PSAVE(KK)
226	IORIG(II)=IPACK*JET+NPRTN
227	IDCAY(II)=0
228	320 CONTINUE
229	IDCAY(NPRTN)=NEXT*IPACK+NPTCL
230	GO TO 310
231	C
232	C or add lepton:
233	330 NPTCL=NPTCL+1
234	DO 331 K=1,5
235	331 PPTCL(K,NPTCL)=PJSET(K,I)
236	IORIG(NPTCL)=IPACK*JET+NPRTN
237	IDENT(NPTCL)=JTYPE(I)
238	IDCAY(NPTCL)=0
239	IDCAY(NPRTN)=NPTCL*IPACK+NPTCL
240	310 CONTINUE
241	NPJET(JET)=NPTCL
242	300 CONTINUE
243	C
244	C Reset NJSET so decay jets do not appear in /JETSET/
245	NJADD=NJSET
246	NJSET=NJSAVE
247	C
248	C Check for at least two particles
249	IF(NPTCL.LT.NHDRN+2) THEN
250	NPTCL=NSTART-1
251	DECJET=.FALSE.
252	RETURN
253	ENDIF
254	C
255	C Conserve charge
256	C
257	SUMQ=0.
258	NHDRN1=NHDRN+1
259	DO 400 I=NHDRN1,NPTCL
260	IDLV1=IDENT(I)
261	SUMQ=SUMQ+CHARGE(IDLV1)
262	400 CONTINUE
263	IDLV1=IDENT(IP)
264	SUMQ=SUMQ-CHARGE(IDLV1)
265	C
266	IF(SUMQ.EQ.0.) GO TO 500
267	C
268	C Charge wrong--fix it by swapping UP and DN quarks.
269	DO 410 I=NHDRN1,NPTCL
270	ID1=IDENT(I)
271	IF(IABS(ID1).GT.1000) GO TO 410
272	I1=MOD(IABS(ID1)/100,10)
273	I2=MOD(IABS(ID1)/10,10)
274	I3=MOD(IABS(ID1),10)
275	IF(I1.EQ.1.AND.I2.GT.2.AND.SUMQ*ID1.GT.0.) THEN
276	IDENT(I)=ISIGN(200+10*I2+I3,ID1)
277	ELSEIF(I1.EQ.2.AND.I2.GT.2.AND.SUMQ*ID1.LT.0.) THEN
278	IDENT(I)=ISIGN(100+10*I2+I3,ID1)
279	ELSEIF(I1.EQ.1.AND.I2.EQ.2.AND.SUMQ*ID1.GT.0.) THEN
280	IDENT(I)=110+I3
281	ELSEIF(I1.EQ.1.AND.I2.EQ.1) THEN
282	IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ))
283	ELSE
284	GO TO 410
285	ENDIF
286	SUMQ=SIGN(ABS(SUMQ)-1.,SUMQ)
287	IDLV1=IDENT(I)
288	PPTCL(5,I)=AMASS(IDLV1)
289	PPTCL(4,I)=SQRT(PPTCL(1,I)2+PPTCL(2,I)2+PPTCL(3,I)**2
290	$ +PPTCL(5,I)**2)
291	C Sum cannot vanish for fractionally charged initial particle.
292	IF(ABS(SUMQ).LT.1.) GO TO 500
293	410 CONTINUE
294	C Failed to conserve charge.
295	NPTCL=NSTART-1
296	DECJET=.FALSE.
297	RETURN
298	C
299	C Rescale momenta for correct mass
300	C
301	500 CONTINUE
302	IF(WDECAY) THEN
303	C Quark decay. First rescale jet3 + W
304	DO 510 K=1,5
305	510 PPTCL(K,NPTCL+1)=PPTCL(K,NPTCLW)
306	NPTLV1=NPTCL+1
307	DO 520 K=1,3
308	520 PSUM(K)=0.
309	PSUM(4)=PPTCL(5,IP)
310	PSUM(5)=PSUM(4)
311	CALL RESCAL(NPJET(2)+1,NPTLV1,PSUM,IFAIL)
312	IF(IFAIL.NE.0) THEN
313	NPTCL=NSTART-1
314	DECJET=.FALSE.
315	RETURN
316	ENDIF
317	DO 530 K=1,3
318	530 BETAW(K)=PPTCL(K,NPTCL+1)/PPTCL(4,NPTCL+1)
319	GAMMAW=PPTCL(4,NPTCL+1)/PPTCL(5,NPTCL+1)
320	C Then rescale W
321	PSUM(4)=PPTCL(5,NPTCLW)
322	PSUM(5)=PSUM(4)
323	CALL RESCAL(NHDRN1,NPJET(2),PSUM,IFAIL)
324	IF(IFAIL.NE.0) THEN
325	NPTCL=NSTART-1
326	DECJET=.FALSE.
327	RETURN
328	ENDIF
329	ELSE
330	C General decay with no W.
331	DO 550 K=1,3
332	550 PSUM(K)=0.
333	PSUM(4)=PPTCL(5,IP)
334	PSUM(5)=PSUM(4)
335	NPTLV1=NPTCL
336	CALL RESCAL(NHDRN1,NPTLV1,PSUM,IFAIL)
337	IF(IFAIL.NE.0) THEN
338	NPTCL=NSTART-1
339	DECJET=.FALSE.
340	RETURN
341	ENDIF
342	ENDIF
343	C
344	C Boost back to lab frame. Reset IORIG.
345	C
346	IF(WDECAY) THEN
347	DO 600 I=NHDRN1,NPTCL
348	JET=IORIG(I)/IPACK
349	IF(JET.NE.1.AND.JET.NE.2) GO TO 600
350	BP=BETAW(1)PPTCL(1,I)+BETAW(2)PPTCL(2,I)+BETAW(3)*PPTCL(3,I)
351	DO 610 J=1,3
352	610 PPTCL(J,I)=PPTCL(J,I)+GAMMAWBETAW(J)(PPTCL(4,I)
353	$ +BP*GAMMAW/(GAMMAW+1.))
354	PPTCL(4,I)=GAMMAW*(PPTCL(4,I)+BP)
355	600 CONTINUE
356	ENDIF
357	C
358	DO 620 I=NSTART,NPTCL
359	IORIG(I)=MOD(IORIG(I),IPACK)
360	BP=BETA(1)PPTCL(1,I)+BETA(2)PPTCL(2,I)+BETA(3)*PPTCL(3,I)
361	DO 621 J=1,3
362	PPTCL(J,I)=PPTCL(J,I)+GAMMABETA(J)(PPTCL(4,I)
363	$ +BP*GAMMA/(GAMMA+1.))
364	621 CONTINUE
365	PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP)
366	620 CONTINUE
367	C
368	C Normal exit
369	C
370	DECJET=.TRUE.
371	RETURN
372	C
373	C Error messages.
374	C
375	9998 DECJET=.FALSE.
376	CALL PRTEVT(0)
377	WRITE(ITLIS,99980) NJSET
378	99980 FORMAT(//5X,'ERROR IN DECJET...NJSET > ',I5)
379	RETURN
380	END