CDECK ID>, HWCCUT. *CMZ :- -26/04/91 14.29.39 by Federico Carminati *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWCCUT(JHEP,KHEP,PCL,BTCLUS,SPLIT) C----------------------------------------------------------------------- C Cuts into 2 the cluster, momentum PCL, made of partons JHEP & KHEP C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWREXQ,HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,EMX,EMY, & QM3,PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5),SKAPPA,DELTM, & VSCA,VTMP(4),RKAPPA,VCLUS INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,ID1,ID2,ID3,NTRY,NTRYMX,J,IB LOGICAL BTCLUS,SPLIT EXTERNAL HWREXQ,HWUPCM,HWR,HWVDOT,HWRINT COMMON/HWCFRM/VCLUS(4,NMXHEP) PARAMETER (SKAPPA=1.,NTRYMX=100) IF (IERROR.NE.0) RETURN EMC=PCL(5) ID1=IDHW(JHEP) ID2=IDHW(KHEP) QM1=RMASS(ID1) QM2=RMASS(ID2) SPLIT=.FALSE. NTRY=0 C Decide if cluster contains a b-(anti)quark IF (ID1.EQ.5.OR.ID1.EQ.11.OR.ID2.EQ.5.OR.ID2.EQ.11) THEN IB=2 ELSE IB=1 ENDIF IF (BTCLUS) THEN C Split beam and target clusters as soft clusters C Both (remnant) children treated like soft clusters if IOPREM=0(1) 10 ID3=HWRINT(1,2) QM3=RMASS(ID3) IF (EMC.LE.QM1+QM2+2.*QM3) THEN ID3=3-ID3 QM3=RMASS(ID3) IF (EMC.LE.QM1+QM2+2.*QM3) RETURN ENDIF PXY=EMC-QM1-QM2-TWO*QM3 IF (ISTHEP(JHEP).EQ.153.OR.ISTHEP(JHEP).EQ.154.OR. & IOPREM.EQ.0) THEN EMX=QM1+QM3+HWREXQ(BTCLM,PXY) ELSE EMX=QM1+QM3+PXY*HWR()**PSPLT(IB) ENDIF IF (ISTHEP(KHEP).EQ.153.OR.ISTHEP(KHEP).EQ.154.OR. & IOPREM.EQ.0) THEN EMY=QM2+QM3+HWREXQ(BTCLM,PXY) ELSE EMY=QM2+QM3+PXY*HWR()**PSPLT(IB) ENDIF IF (EMX+EMY.GE.EMC) THEN NTRY=NTRY+1 IF (NTRY.GT.NTRYMX) RETURN GOTO 10 ENDIF PCX=HWUPCM(EMX,QM1,QM3) PCY=HWUPCM(EMY,QM2,QM3) ELSE C Choose fragment masses for ordinary cluster PXY=EMC-QM1-QM2 20 NTRY=NTRY+1 IF (NTRY.GT.NTRYMX) RETURN 30 EMX=QM1+PXY*HWR()**PSPLT(IB) EMY=QM2+PXY*HWR()**PSPLT(IB) IF (EMX+EMY.GE.EMC) GOTO 30 C u,d,s pair production with weights QWT 40 ID3=HWRINT(1,3) IF (QWT(ID3).LT.HWR()) GOTO 40 QM3=RMASS(ID3) PCX=HWUPCM(EMX,QM1,QM3) IF (PCX.LT.ZERO) GOTO 20 PCY=HWUPCM(EMY,QM2,QM3) IF (PCY.LT.ZERO) GOTO 20 SPLIT=.TRUE. ENDIF C Boost antiquark to CoM frame to find axis CALL HWULOF(PCL,PHEP(1,KHEP),AX) RCM=1./SQRT(HWVDOT(3,AX,AX)) CALL HWVSCA(3,RCM,AX,AX) C Construct new CoM momenta (collinear) PXY=HWUPCM(EMC,EMX,EMY) CALL HWVSCA(3,PXY,AX,PC) PC(4)=SQRT(PXY**2+EMY**2) PC(5)=EMY CALL HWVSCA(3,PCY,AX,PA) PA(4)=SQRT(PCY**2+QM2**2) PA(5)=QM2 CALL HWULOB(PC,PA,PB) CALL HWVDIF(4,PC,PB,PA) PA(5)=QM3 LHEP=NHEP+1 MHEP=NHEP+2 CALL HWULOB(PCL,PB,PHEP(1,KHEP)) CALL HWULOB(PCL,PA,PHEP(1,MHEP)) CALL HWVSCA(3,-ONE,PC,PC) PC(4)=EMC-PC(4) PC(5)=EMX CALL HWVSCA(3,PCX,AX,PA) PA(4)=SQRT(PCX**2+QM3**2) CALL HWULOB(PC,PA,PB) CALL HWULOB(PCL,PB,PHEP(1,LHEP)) DO 50 J=1,4 50 PHEP(J,JHEP)=PCL(J)-PHEP(J,KHEP)-PHEP(J,LHEP)-PHEP(J,MHEP) PHEP(5,JHEP)=QM1 CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP)) C Construct new vertex positions RKAPPA=GEV2MM/SKAPPA CALL HWVSCA(3,RKAPPA,AX,AX) DELTM=(EMX-EMY)*(EMX+EMY)/(TWO*EMC) CALL HWVSCA(3,DELTM,AX,VTMP) VTMP(4)=(HALF*EMC-PXY)*RKAPPA CALL HWULB4(PCL,VTMP,VTMP) CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VHEP(1,LHEP)) CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP)) VSCA=0.25*EMC+HALF*(PXY+DELTM) CALL HWVSCA(3,VSCA,AX,VTMP) VTMP(4)=(EMC-VSCA)*RKAPPA CALL HWULB4(PCL,VTMP,VTMP) CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,MHEP)) VSCA=-0.25*EMC+HALF*(DELTM-PXY) CALL HWVSCA(3,VSCA,AX,VTMP) VTMP(4)=(EMC+VSCA)*RKAPPA CALL HWULB4(PCL,VTMP,VTMP) CALL HWVSUM(4,VTMP,VCLUS(1,JHEP),VCLUS(1,JHEP)) C (Re-)label quarks IDHW(LHEP)=ID3+6 IDHW(MHEP)=ID3 IDHEP(MHEP)= IDPDG(ID3) IDHEP(LHEP)=-IDPDG(ID3) ISTHEP(LHEP)=151 ISTHEP(MHEP)=151 JMOHEP(2,JHEP)=LHEP JDAHEP(2,KHEP)=MHEP JMOHEP(1,LHEP)=JMOHEP(1,KHEP) JMOHEP(2,LHEP)=MHEP JDAHEP(1,LHEP)=0 JDAHEP(2,LHEP)=JHEP JMOHEP(1,MHEP)=JMOHEP(1,JHEP) JMOHEP(2,MHEP)=KHEP JDAHEP(1,MHEP)=0 JDAHEP(2,MHEP)=LHEP NHEP=NHEP+2 999 END CDECK ID>, HWCDEC. *CMZ :- -26/04/91 10.18.56 by Bryan Webber *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWCDEC C----------------------------------------------------------------------- C DECAYS CLUSTERS INTO PRIMARY HADRONS C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER JCL,KCL,IP,JP,KP,IST,ID1,ID2,ID3 IF (IERROR.NE.0) RETURN IF (IPROC/1000.EQ.9.OR.IPROC/1000.EQ.5) THEN C---RELABEL CLUSTER CONNECTED TO REMNANT IN DIS DO 10 JCL=2,NHEP IF (ISTHEP(JCL).EQ.164) GOTO 20 IF (ISTHEP(JCL).EQ.165) THEN IP=JMOHEP(1,JCL) JP=JMOHEP(2,JCL) KP=IP IF (ISTHEP(IP).EQ.162) THEN KP=JP JP=IP ENDIF IF (JMOHEP(2,KP).NE.JP) THEN IP=JMOHEP(2,KP) ELSE IP=JDAHEP(2,KP) ENDIF KCL=JDAHEP(1,IP) IF (ISTHEP(KCL)/10.NE.16) CALL HWWARN('HWCDEC',100,*999) ISTHEP(KCL)=164 GOTO 20 ENDIF 10 CONTINUE ENDIF 20 CONTINUE DO 30 JCL=1,NHEP IST=ISTHEP(JCL) IF (IST.GT.162.AND.IST.LT.166) THEN C---DON'T HADRONIZE BEAM/TARGET CLUSTERS IF (IST.EQ.163.OR..NOT.GENSOF) THEN C---SET UP FLAVOURS FOR CLUSTER DECAY CALL HWCFLA(IDHW(JMOHEP(1,JCL)),IDHW(JMOHEP(2,JCL)),ID1,ID3) CALL HWCHAD(JCL,ID1,ID3,ID2) ENDIF ENDIF 30 CONTINUE ISTAT=50 999 END CDECK ID>, HWCFLA. *CMZ :- -26/04/91 10.18.56 by Bryan Webber *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWCFLA(JD1,JD2,ID1,ID2) C----------------------------------------------------------------------- C SETS UP FLAVOURS FOR CLUSTER DECAY C----------------------------------------------------------------------- INTEGER JD1,JD2,ID1,ID2,JD,JDEC(12) DATA JDEC/1,2,3,10,11,12,4,5,6,7,8,9/ JD=JD1 IF (JD.GT.12) JD=JD-108 ID1=JDEC(JD) JD=JD2 IF (JD.GT.12) JD=JD-96 ID2=JDEC(JD-6) END