CDECK ID>, HWBFIN. *CMZ :- -26/04/91 10.18.56 by Bryan Webber *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWBFIN(IHEP) C----------------------------------------------------------------------- C DELETES INTERNAL LINES FROM SHOWER, MAKES COLOUR CONNECTION INDEX C AND COPIES INTO /HEPEVT/ IN COLOUR ORDER. C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER IHEP,ID,IJET,KHEP,IPAR,JPAR,NXPAR,IP,JP IF (IERROR.NE.0) RETURN C---SAVE VIRTUAL PARTON DATA NHEP=NHEP+1 IF(NHEP.GT.NMXHEP) CALL HWWARN('HWBFIN',100,*999) ID=IDPAR(2) IDHW(NHEP)=ID IDHEP(NHEP)=IDPDG(ID) ISTHEP(NHEP)=ISTHEP(IHEP)+20 JMOHEP(1,NHEP)=IHEP JMOHEP(2,NHEP)=JMOHEP(1,IHEP) JDAHEP(1,IHEP)=NHEP JDAHEP(1,NHEP)=0 JDAHEP(2,NHEP)=0 CALL HWVEQU(5,PPAR(1,2),PHEP(1,NHEP)) CALL HWVEQU(4,VPAR(1,2),VHEP(1,NHEP)) C---FINISHED FOR SPECTATOR OR NON-PARTON JETS IF (ISTHEP(NHEP).GT.136) RETURN IF (ID.GT.13.AND.ID.LT.209 .AND. ID.NE.59) RETURN IF (ID.GT.220.AND.ABS(IDPDG(ID)).LT.1000000) RETURN IF (ID.GT.424.AND.ID.NE.449) RETURN IF (.NOT.TMPAR(2).AND.ID.EQ.59) RETURN IDHEP(NHEP)=94 IJET=NHEP IF (NPAR.GT.2) THEN C---SAVE CONE DATA NHEP=NHEP+1 IF(NHEP.GT.NMXHEP) CALL HWWARN('HWBFIN',101,*999) IDHW(NHEP)=IDPAR(1) IDHEP(NHEP)=0 ISTHEP(NHEP)=100 JMOHEP(1,NHEP)=IHEP JMOHEP(2,NHEP)=JCOPAR(1,1) JDAHEP(1,NHEP)=0 JDAHEP(2,NHEP)=0 CALL HWVEQU(5,PPAR,PHEP(1,NHEP)) CALL HWVEQU(4,VPAR(1,2),VHEP(1,NHEP)) ENDIF KHEP=NHEP C---START WITH ANTICOLOUR DAUGHTER OF HARDEST PARTON IPAR=2 JPAR=JCOPAR(4,IPAR) NXPAR=NPAR/2 DO 20 IP=1,NXPAR DO 10 JP=1,NXPAR IF (JPAR.EQ.0) GOTO 15 IF (JCOPAR(2,JPAR).EQ.IPAR) THEN IPAR=JPAR JPAR=JCOPAR(4,IPAR) ELSE IPAR=JPAR JPAR=JCOPAR(1,IPAR) ENDIF 10 CONTINUE C---COULDN'T FIND COLOUR PARTNER CALL HWWARN('HWBFIN',1,*999) 15 JPAR=JCOPAR(1,IPAR) KHEP=KHEP+1 IF(KHEP.GT.NMXHEP) CALL HWWARN('HWBFIN',102,*999) ID=IDPAR(IPAR) IF (TMPAR(IPAR)) THEN IF (ID.LT.14) THEN ISTHEP(KHEP)=139 ELSEIF (ID.EQ.59) THEN ISTHEP(KHEP)=139 ELSEIF (ID.LT.109) THEN ISTHEP(KHEP)=130 ELSEIF (ID.LT.120) THEN ISTHEP(KHEP)=139 ELSEIF (ABS(IDPDG(ID)).LT.1000000) THEN ISTHEP(KHEP)=130 ELSEIF (ID.LT.425) THEN ISTHEP(KHEP)=139 ELSEIF (ID.EQ.449) THEN ISTHEP(KHEP)=139 ELSE ISTHEP(KHEP)=130 ENDIF ELSE ISTHEP(KHEP)=ISTHEP(IHEP)+24 ENDIF IDHW(KHEP)=ID IDHEP(KHEP)=IDPDG(ID) CALL HWVEQU(5,PPAR(1,IPAR),PHEP(1,KHEP)) CALL HWVEQU(4,VPAR(1,IPAR),VHEP(1,KHEP)) JMOHEP(1,KHEP)=IJET JMOHEP(2,KHEP)=KHEP+1 JDAHEP(1,KHEP)=0 JDAHEP(2,KHEP)=KHEP-1 20 CONTINUE JMOHEP(2,KHEP)=0 JDAHEP(2,NHEP+1)=0 JDAHEP(1,IJET)=NHEP+1 JDAHEP(2,IJET)=KHEP NHEP=KHEP 999 END CDECK ID>, HWBGEN. *CMZ :- -14/10/99 18.04.56 by Mike Seymour *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWBGEN C----------------------------------------------------------------------- C BRANCHING GENERATOR WITH INTERFERING GLUONS C HWBGEN EVOLVES QCD JETS ACCORDING TO THE METHOD OF C G.MARCHESINI & B.R.WEBBER, NUCL. PHYS. B238(1984)1 C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWULDO,HWRGAU,EINHEP,ERTXI,RTXI,XF INTEGER NTRY,LASHEP,IHEP,NRHEP,ID,IST,JHEP,KPAR,I,J,IRHEP(NMXJET), & IRST(NMXJET) LOGICAL HWRLOG EXTERNAL HWULDO,HWRGAU IF (IERROR.NE.0) RETURN IF (IPRO.EQ.80) RETURN C---CHECK THAT EMSCA IS SET IF (EMSCA.LE.ZERO) CALL HWWARN('HWBGEN',200,*999) IF (HARDME) THEN C---FORCE A BRANCH INTO THE `DEAD ZONE' IN E+E- IF (IPROC/10.EQ.10) CALL HWBDED(1) C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DIS IF (IPRO.EQ.90) CALL HWBDIS(1) C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DRELL-YAN PROCESSES IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(1) C---FORCE A BRANCH INTO THE `DEAD ZONE' IN TOP DECAYS CALL HWBTOP ENDIF C---GENERATE INTRINSIC PT ONCE AND FOR ALL DO 5 JNHAD=1,2 IF (PTRMS.NE.0.) THEN PTINT(1,JNHAD)=HWRGAU(1,ZERO,PXRMS) PTINT(2,JNHAD)=HWRGAU(2,ZERO,PXRMS) PTINT(3,JNHAD)=PTINT(1,JNHAD)**2+PTINT(2,JNHAD)**2 ELSE CALL HWVZRO(3,PTINT(1,JNHAD)) ENDIF 5 CONTINUE NTRY=0 LASHEP=NHEP 10 NTRY=NTRY+1 IF (NTRY.GT.NETRY) CALL HWWARN('HWBGEN',ISLENT*100,*999) NRHEP=0 NHEP=LASHEP FROST=.FALSE. DO 100 IHEP=1,LASHEP IST=ISTHEP(IHEP) IF (IST.GE.111.AND.IST.LE.115) THEN NRHEP=NRHEP+1 IRHEP(NRHEP)=IHEP IRST(NRHEP)=IST ID=IDHW(IHEP) IF (IST.NE.115) THEN C---FOUND A PARTON TO EVOLVE NEVPAR=IHEP NPAR=2 IDPAR(1)=17 IDPAR(2)=ID TMPAR(1)=.TRUE. PPAR(2,1)=0. PPAR(4,1)=1. DO 15 J=1,2 DO 15 I=1,2 JMOPAR(I,J)=0 15 JCOPAR(I,J)=0 C---SET UP EVOLUTION SCALE AND FRAME JHEP=JMOHEP(2,IHEP) IF (ID.EQ.13) THEN IF (HWRLOG(HALF)) JHEP=JDAHEP(2,IHEP) ELSEIF (IST.GT.112) THEN IF ((ID.GT.6.AND.ID.LT.13).OR. & (ID.GT.214.AND.ID.LT.221)) JHEP=JDAHEP(2,IHEP) ELSE IF (ID.LT.7.OR.(ID.GT.208.AND.ID.LT.215)) JHEP=JDAHEP(2,IHEP) ENDIF IF (JHEP.LE.0.OR.JHEP.GT.NHEP) THEN CALL HWWARN('HWBGEN',1,*999) JHEP=IHEP ENDIF JCOPAR(1,1)=JHEP EINHEP=PHEP(4,IHEP) ERTXI=HWULDO(PHEP(1,IHEP),PHEP(1,JHEP)) IF (ERTXI.LT.ZERO) ERTXI=0. IF (IST.LE.112.AND.IHEP.EQ.JHEP) ERTXI=0. IF (ISTHEP(JHEP).EQ.155) THEN ERTXI=ERTXI/PHEP(5,JHEP) RTXI=1. ELSE ERTXI=SQRT(ERTXI) RTXI=ERTXI/EINHEP ENDIF IF (RTXI.EQ.ZERO) THEN XF=1. PPAR(1,1)=0. PPAR(3,1)=1. PPAR(1,2)=EINHEP PPAR(2,2)=0. PPAR(4,2)=EINHEP ELSE XF=1./RTXI PPAR(1,1)=1. PPAR(3,1)=0. PPAR(1,2)=ERTXI PPAR(2,2)=1. PPAR(4,2)=ERTXI ENDIF IF (PPAR(4,2).LT.PHEP(5,IHEP)) PPAR(4,2)=PHEP(5,IHEP) C---STORE MASS PPAR(5,2)=PHEP(5,IHEP) CALL HWVZRO(4,VPAR(1,1)) CALL HWVZRO(4,VPAR(1,2)) IF (IST.GT.112) THEN TMPAR(2)=.TRUE. INHAD=0 JNHAD=0 XFACT=0. ELSE TMPAR(2)=.FALSE. JNHAD=IST-110 INHAD=JNHAD IF (JDAHEP(1,JNHAD).NE.0) INHAD=JDAHEP(1,JNHAD) XFACT=XF/PHEP(4,INHAD) ANOMSC(1,JNHAD)=ZERO ANOMSC(2,JNHAD)=ZERO ENDIF C---FOR QUARKS IN A COLOUR SINGLET, ALLOW SOFT MATRIX-ELEMENT CORRECTION HARDST=PPAR(4,2) IF (SOFTME.AND.IDHW(IHEP).LT.13.AND. $ ((JMOHEP(2,JHEP).EQ.IHEP.AND.JDAHEP(2,JHEP).EQ.IHEP).OR. $ ISTHEP(JHEP).EQ.155)) HARDST=0 C---CREATE BRANCHES AND COMPUTE ENERGIES DO 20 KPAR=2,NMXPAR IF (TMPAR(KPAR)) THEN CALL HWBRAN(KPAR) ELSE CALL HWSBRN(KPAR) ENDIF IF (IERROR.NE.0) RETURN IF (FROST) GOTO 100 IF (KPAR.EQ.NPAR) GOTO 30 20 CONTINUE C---COMPUTE MASSES AND 3-MOMENTA 30 CONTINUE CALL HWBMAS IF (AZSPIN) CALL HWBSPN IF (TMPAR(2)) THEN CALL HWBTIM(2,1) ELSE CALL HWBSPA ENDIF C---ENTER PARTON JET IN /HEPEVT/ CALL HWBFIN(IHEP) ELSE C---COPY SPECTATOR NHEP=NHEP+1 IF (ID.GT.120.AND.ID.LT.133 .OR. ID.GE.198.AND.ID.LE.201) THEN ISTHEP(NHEP)=190 ELSE ISTHEP(NHEP)=152 ENDIF IDHW(NHEP)=ID IDHEP(NHEP)=IDPDG(ID) JMOHEP(1,NHEP)=IHEP JMOHEP(2,NHEP)=0 JDAHEP(2,NHEP)=0 JDAHEP(1,IHEP)=NHEP CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP)) ENDIF ISTHEP(IHEP)=ISTHEP(IHEP)+10 ENDIF 100 CONTINUE IF (.NOT.FROST) THEN C---COMBINE JETS ISTAT=20 CALL HWBJCO ENDIF IF (.NOT.FROST) THEN C---ATTACH SPECTATORS ISTAT=30 CALL HWSSPC ENDIF IF (FROST) THEN C---BAD JET: RESTORE PARTONS AND RE-EVOLVE DO 120 I=1,NRHEP 120 ISTHEP(IRHEP(I))=IRST(I) GOTO 10 ENDIF C---CONNECT COLOURS CALL HWBCON ISTAT=40 LASHEP=NHEP IF (HARDME) THEN C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN E+E- IF (IPROC/10.EQ.10) CALL HWBDED(2) C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DIS IF (IPRO.EQ.90) CALL HWBDIS(2) C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DRELL-YAN PROC IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(2) ENDIF C---IF THE CLEAN-UP OPERATION ADDED ANY PARTONS TO THE EVENT RECORD C IT MIGHT NEED RESHOWERING IF (NHEP.GT.LASHEP) THEN LASHEP=NHEP GOTO 10 ENDIF 999 END CDECK ID>, HWBJCO. *CMZ :- -26/04/91 14.25.31 by Federico Carminati *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWBJCO C----------------------------------------------------------------------- C COMBINES JETS WITH REQUIRED KINEMATICS C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWULDO,EPS,PTX,PTY,PF,PTINF,PTCON,CN,CP,SP,PP0, & PM0,ET0,DET,ECM,EMJ,EMP,EMS,DMS,ES,DPF,ALF,AL(2),ET(2),PP(2), & PT(3),PB(5),PC(5),PQ(5),PR(5),PS(5),RR(3,3),RS(3,3),ETC, & PJ(NMXJET),PM(NMXJET),PBR(5),RBR(3,3),DISP(4) INTEGER LJET,IJ1,IST,IP,ICM,IP1,IP2,NP,IHEP,MHEP,JP,KP,LP,KHEP, & JHEP,NE,IJT,IEND(2),IJET(NMXJET),IPAR(NMXJET) LOGICAL AZCOR,JETRAD,DISPRO,DISLOW EXTERNAL HWULDO PARAMETER (EPS=1.D-4) IF (IERROR.NE.0) RETURN AZCOR=AZSOFT.OR.AZSPIN C---FIRST LOOK FOR SPACELIKE JETS LJET=131 10 IJET(1)=1 20 IJ1=IJET(1) DO 40 IHEP=IJ1,NHEP IST=ISTHEP(IHEP) IF (IST.EQ.137.OR.IST.EQ.138) IST=133 IF (IST.EQ.LJET) THEN C---FOUND AN UNBOOSTED JET - FIND PARTNERS IP=JMOHEP(1,IHEP) ICM=JMOHEP(1,IP) DISPRO=IPRO/10.EQ.9.AND.IDHW(ICM).EQ.15 DISLOW=DISPRO.AND.JDAHEP(1,ICM).EQ.JDAHEP(2,ICM)-1 IF (IST.EQ.131) THEN IP1=JMOHEP(1,ICM) IP2=JMOHEP(2,ICM) ELSE IP1=JDAHEP(1,ICM) IP2=JDAHEP(2,ICM) ENDIF IF (IP1.NE.IP) CALL HWWARN('HWBJCO',100,*999) NP=0 DO 30 JHEP=IP1,IP2 NP=NP+1 IPAR(NP)=JHEP 30 IJET(NP)=JDAHEP(1,JHEP) GOTO 50 ENDIF 40 CONTINUE C---NO MORE JETS? IF (LJET.EQ.131) THEN LJET=133 GOTO 10 ENDIF RETURN 50 IF (LJET.EQ.131) THEN C---SPACELIKE JETS: FIND SPACELIKE PARTONS IF (NP.NE.2) CALL HWWARN('HWBJCO',103,*999) C---special for DIS: FIND BOOST AND ROTATION FROM LAB TO BREIT FRAME IF (DISPRO.AND.BREIT) THEN IP=2 IF (JDAHEP(1,IP).NE.0) IP=JDAHEP(1,IP) CALL HWVDIF(4,PHEP(1,JMOHEP(1,ICM)),PHEP(1,JDAHEP(1,ICM)),PB) CALL HWUMAS(PB) C---IF Q**2<10**-2, SOMETHING MUST HAVE ALREADY GONE WRONG IF (PB(5)**2.LT.1.D-2) CALL HWWARN('HWBJCO',102,*999) CALL HWVSCA(4,PB(5)**2/HWULDO(PHEP(1,IP),PB),PHEP(1,IP),PBR) CALL HWVSUM(4,PB,PBR,PBR) CALL HWUMAS(PBR) CALL HWULOF(PBR,PB,PB) CALL HWUROT(PB,ONE,ZERO,RBR) ENDIF PTX=0. PTY=0. PF=1.D0 DO 90 IP=1,2 MHEP=IJET(IP) IF (JDAHEP(1,MHEP).EQ.0) THEN C---SPECIAL FOR NON-PARTON JETS IHEP=MHEP GOTO 70 ELSE IST=134+IP DO 60 IHEP=MHEP,NHEP 60 IF (ISTHEP(IHEP).EQ.IST) GOTO 70 C---COULDN'T FIND SPACELIKE PARTON CALL HWWARN('HWBJCO',101,*999) ENDIF 70 CALL HWVSCA(3,PF,PHEP(1,IHEP),PS) IF (PTINT(3,IP).GT.ZERO) THEN C---ADD INTRINSIC PT PT(1)=PTINT(1,IP) PT(2)=PTINT(2,IP) PT(3)=0. CALL HWUROT(PS, ONE,ZERO,RS) CALL HWUROB(RS,PT,PT) CALL HWVSUM(3,PS,PT,PS) ENDIF JP=IJET(IP)+1 IF (AZCOR.AND.JP.LE.NHEP.AND.IDHW(JP).EQ.17) THEN C---ALIGN CONE WITH INTERFERING PARTON CALL HWUROT(PS, ONE,ZERO,RS) CALL HWUROF(RS,PHEP(1,JP),PR) PTCON=PR(1)**2+PR(2)**2 KP=JMOHEP(2,JP) IF (KP.EQ.0) THEN CALL HWWARN('HWBJCO',1,*999) PTINF=0. ELSE CALL HWVEQU(4,PHEP(1,KP),PB) IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PB,PB) CALL HWUROF(RBR,PB,PB) ENDIF PTINF=PB(1)**2+PB(2)**2 IF (PTINF.LT.EPS) THEN C---COLLINEAR JETS: ALIGN CONES KP=JDAHEP(1,KP)+1 IF (ISTHEP(KP).EQ.100.AND.ISTHEP(KP-1)/10.EQ.14) THEN CALL HWVEQU(4,PHEP(1,KP),PB) IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PB,PB) CALL HWUROF(RBR,PB,PB) ENDIF PTINF=PB(1)**2+PB(2)**2 ELSE PTINF=0. ENDIF ENDIF ENDIF IF (PTCON.NE.ZERO.AND.PTINF.NE.ZERO) THEN CN=1./SQRT(PTINF*PTCON) CP=CN*(PR(1)*PB(1)+PR(2)*PB(2)) SP=CN*(PR(1)*PB(2)-PR(2)*PB(1)) ELSE CALL HWRAZM( ONE,CP,SP) ENDIF ELSE CALL HWRAZM( ONE,CP,SP) ENDIF C---ROTATE SO SPACELIKE IS ALONG AXIS (APART FROM INTRINSIC PT) CALL HWUROT(PS,CP,SP,RS) IHEP=IJET(IP) KHEP=JDAHEP(2,IHEP) IF (KHEP.LT.IHEP) KHEP=IHEP IEND(IP)=KHEP DO 80 JHEP=IHEP,KHEP CALL HWUROF(RS,PHEP(1,JHEP),PHEP(1,JHEP)) 80 CALL HWUROF(RS,VHEP(1,JHEP),VHEP(1,JHEP)) PP(IP)=PHEP(4,IHEP)+PF*PHEP(3,IHEP) ET(IP)=PHEP(1,IHEP)**2+PHEP(2,IHEP)**2-PHEP(5,IHEP)**2 C---REDEFINE HARD CM PTX=PTX+PHEP(1,IHEP) PTY=PTY+PHEP(2,IHEP) 90 PF=-PF PHEP(1,ICM)=PTX PHEP(2,ICM)=PTY C---special for DIS: keep lepton momenta fixed IF (DISPRO) THEN IP1=JMOHEP(1,ICM) IP2=JDAHEP(1,ICM) IJT=IJET(1) C---IJT will be used to store lepton momentum transfer CALL HWVDIF(4,PHEP(1,IP1),PHEP(1,IP2),PHEP(1,IJT)) CALL HWUMAS(PHEP(1,IJT)) IF (IDHEP(IP1).EQ.IDHEP(IP2)) THEN IDHW(IJT)=200 ELSEIF (IDHEP(IP1).LT.IDHEP(IP2)) THEN IDHW(IJT)=199 ELSE IDHW(IJT)=198 ENDIF IDHEP(IJT)=IDPDG(IDHW(IJT)) ISTHEP(IJT)=3 C---calculate boost for struck parton C PC is momentum of outgoing parton(s) IP2=JDAHEP(2,ICM) IF (.NOT.DISLOW) THEN C---FOR heavy QQbar PQ and PC are old and new QQbar momenta CALL HWVSUM(4,PHEP(1,IP2-1),PHEP(1,IP2),PQ) CALL HWUMAS(PQ) PC(5)=PQ(5) ELSE PC(5)=PHEP(5,JDAHEP(1,IP2)) ENDIF CALL HWVSUM(2,PHEP(1,IJT),PHEP(1,IJET(2)),PC) ET(1)=ET(2) C---USE BREIT FRAME BOSON MOMENTUM IF NECESSARY IF (BREIT) THEN ET(2)=ET(1)+PC(5)**2+PHEP(5,IJET(2))**2 PM0=PHEP(5,IJT) PP0=-PM0 ELSE ET(2)=PC(1)**2+PC(2)**2+PC(5)**2 PP0=PHEP(4,IJT)+PHEP(3,IJT) PM0=PHEP(4,IJT)-PHEP(3,IJT) ENDIF ET0=(PP0*PM0)+ET(1)-ET(2) DET=ET0**2-4.*(PP0*PM0)*ET(1) IF (DET.LT.ZERO) THEN FROST=.TRUE. RETURN ENDIF ALF=(SQRT(DET)-ET0)/(2.*PP0*PP(2)) PB(1)=0. PB(2)=0. PB(5)=2.D0 PB(3)=ALF-(1./ALF) PB(4)=ALF+(1./ALF) DO 100 IHEP=IJET(2),IEND(2) CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP)) CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP)) C---BOOST FROM BREIT FRAME IF NECESSARY IF (BREIT) THEN CALL HWUROB(RBR,PHEP(1,IHEP),PHEP(1,IHEP)) CALL HWULOB(PBR,PHEP(1,IHEP),PHEP(1,IHEP)) CALL HWUROB(RBR,VHEP(1,IHEP),VHEP(1,IHEP)) CALL HWULB4(PBR,VHEP(1,IHEP),VHEP(1,IHEP)) ENDIF 100 ISTHEP(IHEP)=ISTHEP(IHEP)+10 CALL HWVDIF(4,VHEP(1,IPAR(2)),VHEP(1,IJET(2)),DISP) DO 110 IHEP=IJET(2),IEND(2) 110 CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP)) IF (IEND(2).GT.IJET(2)+1) ISTHEP(IJET(2)+1)=100 CALL HWVSUM(4,PHEP(1,IJT),PHEP(1,IJET(2)),PC) CALL HWVSUM(4,PHEP(1,IP1),PHEP(1,IJET(2)),PHEP(1,ICM)) CALL HWUMAS(PHEP(1,ICM)) ELSEIF (IPRO/10.EQ.5) THEN C Special to preserve photon momentum ETC=PTX**2+PTY**2+PHEP(5,ICM)**2 ET0=ETC+ET(1)-ET(2) DET=ET0**2-4.*ETC*ET(1) IF (DET.LT.ZERO) THEN FROST=.TRUE. RETURN ENDIF ALF=(SQRT(DET)+ET0-2.*ET(1))/(2.*PP(1)*PP(2)) PB(1)=0. PB(2)=0. PB(3)=ALF-1./ALF PB(4)=ALF+1./ALF PB(5)=2. IJT=IJET(2) DO 120 IHEP=IJT,IEND(2) CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP)) CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP)) 120 ISTHEP(IHEP)=ISTHEP(IHEP)+10 CALL HWVDIF(4,VHEP(1,IPAR(2)),VHEP(1,IJT),DISP) DO 130 IHEP=IJT,IEND(2) 130 CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP)) IF (IEND(2).GT.IJT+1) ISTHEP(IJT+1)=100 ISTHEP(IJET(1))=ISTHEP(IJET(1))+10 CALL HWVSUM(2,PHEP(3,IPAR(1)),PHEP(3,IJT),PHEP(3,ICM)) ELSE PHEP(4,ICM)=SQRT(PTX**2+PTY**2+PHEP(3,ICM)**2+PHEP(5,ICM)**2) C---NOW BOOST TO REQUIRED Q**2 AND X-F PP0=PHEP(4,ICM)+PHEP(3,ICM) PM0=PHEP(4,ICM)-PHEP(3,ICM) ET0=(PP0*PM0)+ET(1)-ET(2) DET=ET0**2-4.*(PP0*PM0)*ET(1) IF (DET.LT.ZERO) THEN FROST=.TRUE. RETURN ENDIF DET=SQRT(DET)+ET0 AL(1)= 2.*PM0*PP(1)/DET AL(2)=(PM0/PP(2))*(1.-2.*ET(1)/DET) PB(1)=0. PB(2)=0. PB(5)=2. DO 160 IP=1,2 PB(3)=AL(IP)-(1./AL(IP)) PB(4)=AL(IP)+(1./AL(IP)) IJT=IJET(IP) DO 140 IHEP=IJT,IEND(IP) CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP)) CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP)) 140 ISTHEP(IHEP)=ISTHEP(IHEP)+10 CALL HWVDIF(4,VHEP(1,IPAR(IP)),VHEP(1,IJT),DISP) DO 150 IHEP=IJT,IEND(IP) 150 CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP)) IF (IEND(IP).GT.IJT+1) THEN ISTHEP(IJT+1)=100 ELSEIF (IEND(IP).EQ.IJT) THEN C---NON-PARTON JET ISTHEP(IJT)=3 ENDIF 160 CONTINUE ENDIF ISTHEP(ICM)=120 ELSE C---TIMELIKE JETS C special for DIS: preserve outgoing lepton momentum IF (DISPRO) THEN CALL HWVEQU(5,PHEP(1,IPAR(1)),PHEP(1,IJET(1))) ISTHEP(IJET(1))=1 LP=2 ELSE CALL HWVEQU(5,PHEP(1,ICM),PC) C--- PQ AND PC ARE OLD AND NEW PARTON CM CALL HWVSUM(4,PHEP(1,IPAR(1)),PHEP(1,IPAR(2)),PQ) PQ(5)=PHEP(5,ICM) IF (NP.GT.2) THEN DO 170 KP=3,NP 170 CALL HWVSUM(4,PHEP(1,IPAR(KP)),PQ,PQ) ENDIF LP=1 ENDIF IF (.NOT.DISLOW) THEN C---FIND JET CM MOMENTA ECM=PQ(5) EMS=0. JETRAD=.FALSE. DO 180 KP=LP,NP EMJ=PHEP(5,IJET(KP)) EMP=PHEP(5,IPAR(KP)) JETRAD=JETRAD.OR.EMJ.NE.EMP EMS=EMS+EMJ PM(KP)= EMJ**2 C---N.B. ROUNDING ERRORS HERE AT HIGH ENERGIES PJ(KP)=(HWULDO(PHEP(1,IPAR(KP)),PQ)/ECM)**2-EMP**2 IF (PJ(KP).LE.ZERO) CALL HWWARN('HWBJCO',104,*999) 180 CONTINUE PF=1. IF (JETRAD) THEN C---JETS DID RADIATE IF (EMS.GE.ECM) THEN FROST=.TRUE. RETURN ENDIF DO 200 NE=1,NETRY EMS=-ECM DMS=0. DO 190 KP=LP,NP ES=SQRT(PF*PJ(KP)+PM(KP)) EMS=EMS+ES 190 DMS=DMS+PJ(KP)/ES DPF=2.*EMS/DMS IF (DPF.GT.PF) DPF=0.9*PF PF=PF-DPF 200 IF (ABS(DPF).LT.EPS) GOTO 210 CALL HWWARN('HWBJCO',105,*999) ENDIF 210 CONTINUE ENDIF C---BOOST PC AND PQ TO BREIT FRAME IF NECESSARY IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PC,PC) CALL HWUROF(RBR,PC,PC) IF (.NOT.DISLOW) THEN CALL HWULOF(PBR,PQ,PQ) CALL HWUROF(RBR,PQ,PQ) ENDIF ENDIF DO 230 IP=LP,NP C---FIND CM ROTATION FOR JET IP IF (.NOT.DISLOW) THEN CALL HWVEQU(4,PHEP(1,IPAR(IP)),PR) IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PR,PR) CALL HWUROF(RBR,PR,PR) ENDIF CALL HWULOF(PQ,PR,PR) CALL HWUROT(PR, ONE,ZERO,RR) PR(1)=0. PR(2)=0. PR(3)=SQRT(PF*PJ(IP)) PR(4)=SQRT(PF*PJ(IP)+PM(IP)) PR(5)=PHEP(5,IJET(IP)) CALL HWUROB(RR,PR,PR) CALL HWULOB(PC,PR,PR) ELSE CALL HWVEQU(5,PC,PR) ENDIF C---NOW PR IS LAB/BREIT MOMENTUM OF JET IP KP=IJET(IP)+1 IF (AZCOR.AND.KP.LE.NHEP.AND.IDHW(KP).EQ.17) THEN C---ALIGN CONE WITH INTERFERING PARTON CALL HWUROT(PR, ONE,ZERO,RS) JP=JMOHEP(2,KP) IF (JP.EQ.0) THEN CALL HWWARN('HWBJCO',2,*999) PTINF=0. ELSE CALL HWVEQU(4,PHEP(1,JP),PS) IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PS,PS) CALL HWUROF(RBR,PS,PS) ENDIF CALL HWUROF(RS,PS,PS) PTINF=PS(1)**2+PS(2)**2 IF (PTINF.LT.EPS) THEN C---COLLINEAR JETS: ALIGN CONES JP=JDAHEP(1,JP)+1 IF (ISTHEP(JP).EQ.100.AND.ISTHEP(JP-1)/10.EQ.14) THEN CALL HWVEQU(4,PHEP(1,JP),PS) IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PS,PS) CALL HWUROF(RBR,PS,PS) ENDIF CALL HWUROF(RS,PS,PS) PTINF=PS(1)**2+PS(2)**2 ELSE PTINF=0. ENDIF ENDIF ENDIF CALL HWVEQU(4,PHEP(1,KP),PB) IF (DISPRO.AND.BREIT) THEN CALL HWULOF(PBR,PB,PB) CALL HWUROF(RBR,PB,PB) ENDIF PTCON=PB(1)**2+PB(2)**2 IF (PTCON.NE.ZERO.AND.PTINF.NE.ZERO) THEN CN=1./SQRT(PTINF*PTCON) CP=CN*(PS(1)*PB(1)+PS(2)*PB(2)) SP=CN*(PS(1)*PB(2)-PS(2)*PB(1)) ELSE CALL HWRAZM( ONE,CP,SP) ENDIF ELSE CALL HWRAZM( ONE,CP,SP) ENDIF CALL HWUROT(PR,CP,SP,RS) C---FIND BOOST FOR JET IP ALF=(PHEP(3,IJET(IP))+PHEP(4,IJET(IP)))/ & (PR(4)+SQRT((PR(4)+PR(5))*(PR(4)-PR(5)))) PB(1)=0. PB(2)=0. PB(3)=ALF-(1./ALF) PB(4)=ALF+(1./ALF) PB(5)=2. IHEP=IJET(IP) KHEP=JDAHEP(2,IHEP) IF (KHEP.LT.IHEP) KHEP=IHEP DO 220 JHEP=IHEP,KHEP CALL HWULOF(PB,PHEP(1,JHEP),PHEP(1,JHEP)) CALL HWUROB(RS,PHEP(1,JHEP),PHEP(1,JHEP)) CALL HWULF4(PB,VHEP(1,JHEP),VHEP(1,JHEP)) CALL HWUROB(RS,VHEP(1,JHEP),VHEP(1,JHEP)) C---BOOST FROM BREIT FRAME IF NECESSARY IF (DISPRO.AND.BREIT) THEN CALL HWUROB(RBR,PHEP(1,JHEP),PHEP(1,JHEP)) CALL HWULOB(PBR,PHEP(1,JHEP),PHEP(1,JHEP)) CALL HWUROB(RBR,VHEP(1,JHEP),VHEP(1,JHEP)) CALL HWULB4(PBR,VHEP(1,JHEP),VHEP(1,JHEP)) ENDIF CALL HWVSUM(4,VHEP(1,JHEP),VHEP(1,IPAR(IP)),VHEP(1,JHEP)) 220 ISTHEP(JHEP)=ISTHEP(JHEP)+10 IF (KHEP.GT.IHEP+1) THEN ISTHEP(IHEP+1)=100 ELSEIF (KHEP.EQ.IHEP) THEN C---NON-PARTON JET ISTHEP(IHEP)=190 ENDIF 230 CONTINUE IF (ISTHEP(ICM).EQ.110) ISTHEP(ICM)=120 ENDIF GOTO 20 999 END CDECK ID>, HWBMAS. *CMZ :- -26/04/91 11.11.54 by Bryan Webber *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWBMAS C----------------------------------------------------------------------- C Passes backwards through a jet cascade calculating the masses C and magnitudes of the longitudinal and transverse three momenta. C Components given relative to direction of parent for a time-like C vertex and with respect to z-axis for space-like vertices. C C On input PPAR(1-5,*) contains: C (E*sqrt(Xi),Xi,3-mom (if external),E,M-sq (if external)) C C On output PPAR(1-5,*) (if TMPAR(*)), containts: C (P-trans,Xi or Xilast,P-long,E,M) C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWUSQR,EXI,PISQ,PJPK,EJEK,PTSQ,Z,ZMIN,ZMAX, $ EMI,EMJ,EMK,C,NQ,HWBVMC,RHO,POLD,PNEW,EOLD,ENEW,A,B INTEGER IPAR,JPAR,KPAR,MPAR,I,J,K EXTERNAL HWUSQR IF (IERROR.NE.0) RETURN IF (NPAR.GT.2) THEN DO 30 MPAR=NPAR-1,3,-2 JPAR=MPAR C Find parent and partner of this branch IPAR=JMOPAR(1,JPAR) KPAR=JPAR+1 C Determine type of branching IF (TMPAR(IPAR)) THEN C Time-like branching C Compute mass of parent EXI=PPAR(1,JPAR)*PPAR(1,KPAR) PPAR(5,IPAR)=PPAR(5,JPAR)+PPAR(5,KPAR)+2.*EXI C Compute three momentum of parent PISQ=PPAR(4,IPAR)*PPAR(4,IPAR)-PPAR(5,IPAR) PPAR(3,IPAR)=HWUSQR(PISQ) C---SPECIAL FOR G-->QQBAR: READJUST ANGULAR DISTRIBUTION IF (IDPAR(IPAR).EQ.13 .AND. IDPAR(JPAR).LT.13) THEN Z=PPAR(4,JPAR)/PPAR(4,IPAR) ZMIN=HWBVMC(IDPAR(JPAR))/PPAR(1,JPAR)*Z RHO=(Z*(3-Z*(3-2*Z))-ZMIN*(3-ZMIN*(3-2*ZMIN))) $ /(2*(1-2*ZMIN)*(1-ZMIN*(1-ZMIN))) NQ=PPAR(3,IPAR)*(PPAR(3,IPAR)+PPAR(4,IPAR)) EMI=PPAR(5,IPAR) EMJ=PPAR(5,JPAR) EMK=PPAR(5,KPAR) ZMIN=MAX((EMI+EMJ-EMK)/(2*(EMI+NQ)), $ (EMI+EMJ-EMK-SQRT((EMI-EMJ-EMK)**2-4*EMJ*EMK))/(2*EMI)) ZMAX=1-MAX((EMI-EMJ+EMK)/(2*(EMI+NQ)), $ (EMI-EMJ+EMK-SQRT((EMI-EMJ-EMK)**2-4*EMJ*EMK))/(2*EMI)) C=2*RMASS(IDPAR(JPAR))**2/EMI Z=(4*ZMIN*(1.5*(1+C-ZMIN)+ZMIN**2)*(1-RHO) $ +4*ZMAX*(1.5*(1+C-ZMAX)+ZMAX**2)*RHO-2-3*C)/(1+2*C)**1.5 Z=SQRT(1+2*C)*SINH(LOG(Z+SQRT(Z**2+1))/3)+0.5 Z=(Z*NQ+(EMI+EMJ-EMK)/2)/(NQ+EMI) PPAR(4,JPAR)=Z*PPAR(4,IPAR) PPAR(4,KPAR)=PPAR(4,IPAR)-PPAR(4,JPAR) PPAR(3,JPAR)=HWUSQR(PPAR(4,JPAR)**2-EMJ) PPAR(3,KPAR)=HWUSQR(PPAR(4,KPAR)**2-EMK) PPAR(2,JPAR)=EXI/(PPAR(4,JPAR)*PPAR(4,KPAR)) IF(JDAPAR(2,JPAR).NE.0)PPAR(2,JDAPAR(2,JPAR))=PPAR(2,JPAR) IF(JDAPAR(2,KPAR).NE.0)PPAR(2,JDAPAR(2,KPAR))=PPAR(2,JPAR) C---FIND DESCENDENTS OF THIS SPLITTING AND READJUST THEIR MOMENTA TOO DO 20 J=JPAR+2,NPAR-1,2 I=J 10 I=JMOPAR(1,I) IF (I.GT.IPAR) GOTO 10 IF (I.EQ.IPAR) THEN I=JMOPAR(1,J) K=J+1 POLD=PPAR(3,J)+PPAR(3,K) EOLD=PPAR(4,J)+PPAR(4,K) PNEW=HWUSQR(PPAR(4,I)**2-PPAR(5,I)) ENEW=PPAR(4,I) A=(ENEW*EOLD-PNEW*POLD)/PPAR(5,I) B=(PNEW*EOLD-ENEW*POLD)/PPAR(5,I) PPAR(3,J)=A*PPAR(3,J)+B*PPAR(4,J) PPAR(4,J)=(PPAR(4,J)+B*PPAR(3,J))/A PPAR(3,K)=PNEW-PPAR(3,J) PPAR(4,K)=ENEW-PPAR(4,J) PPAR(2,J)=1-(PPAR(3,J)*PPAR(3,K)+PPAR(1,J)*PPAR(1,K)) $ /(PPAR(4,J)*PPAR(4,K)) IF (JDAPAR(2,J).NE.0) PPAR(2,JDAPAR(2,J))=PPAR(2,J) IF (JDAPAR(2,K).NE.0) PPAR(2,JDAPAR(2,K))=PPAR(2,J) ENDIF 20 CONTINUE ENDIF C Compute daughter' transverse and longitudinal momenta PJPK=PPAR(3,JPAR)*PPAR(3,KPAR) EJEK=PPAR(4,JPAR)*PPAR(4,KPAR)-EXI PTSQ=(PJPK+EJEK)*(PJPK-EJEK)/PISQ PPAR(1,JPAR)=HWUSQR(PTSQ) PPAR(3,JPAR)=HWUSQR(PPAR(3,JPAR)*PPAR(3,JPAR)-PTSQ) PPAR(1,KPAR)=-PPAR(1,JPAR) PPAR(3,KPAR)= PPAR(3,IPAR)-PPAR(3,JPAR) ELSE C Space-like branching C Re-arrange such that JPAR is time-like IF (TMPAR(KPAR)) THEN KPAR=JPAR JPAR=JPAR+1 ENDIF C Compute time-like branch PTSQ=(2.-PPAR(2,JPAR))*PPAR(1,JPAR)*PPAR(1,JPAR) & -PPAR(5,JPAR) PPAR(1,JPAR)=HWUSQR(PTSQ) PPAR(3,JPAR)=(1.-PPAR(2,JPAR))*PPAR(4,JPAR) PPAR(3,IPAR)=PPAR(3,KPAR)-PPAR(3,JPAR) PPAR(5,IPAR)=0. PPAR(1,KPAR)=0. ENDIF C Reset Xi to Xilast PPAR(2,KPAR)=PPAR(2,IPAR) 30 CONTINUE ENDIF DO 40 IPAR=2,NPAR 40 PPAR(5,IPAR)=HWUSQR(PPAR(5,IPAR)) PPAR(1,2)=0. PPAR(2,2)=0. END