* * $Id: hwbjco.F,v 1.1.1.1 1996/03/08 17:02:10 mclareni Exp $ * * $Log: hwbjco.F,v $ * Revision 1.1.1.1 1996/03/08 17:02:10 mclareni * Herwig58 * * *CMZ : 29/08/94 11.51.46 by Unknown *-- Author : CDECK ID>, HWBJCO. *CMZ :- -26/04/91 14.25.31 by Federico Carminati *-- Author : Bryan Webber C------------------------------------------------------------------------ SUBROUTINE HWBJCO C COMBINES JETS WITH REQUIRED KINEMATICS C----------------------------------------------------------------------- C REVISED 15/10/90 TO PRESERVE LEPTON MOMENTA IN DIS C Modified to preserve photon beam momentum in point-like C interactions: IGK 1/10/92 C REVISED 9/12/92 TO INCLUDE NMXJET C 16/10/93 TO USE BREIT FRAME FOR DIS IF BREIT=.TRUE. AND C Q**2>10**-4 C----------------------------------------------------------------------- #include "herwig58/herwig58.inc" INTEGER LJET,IJ1,IST,IP,ICM,IP1,IP2,NP,IHEP,MHEP,JP,KP,LP,KHEP, & JHEP,NE,IJT,IEND(2) INTEGER IJET(NMXJET),IPAR(NMXJET) 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) LOGICAL AZCOR,JETRAD,DISPRO,DISLOW,BREITF PARAMETER (EPS=1.D-4) IF (IERROR.NE.0) RETURN AZCOR=AZSOFT.OR.AZSPIN BREITF=.FALSE. C---FIRST LOOK FOR SPACELIKE JETS LJET=131 1 IJET(1)=1 5 CONTINUE IJ1=IJET(1) DO 20 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 10 JHEP=IP1,IP2 NP=NP+1 IPAR(NP)=JHEP 10 IJET(NP)=JDAHEP(1,JHEP) GO TO 30 ENDIF 20 CONTINUE C---NO MORE JETS? IF (LJET.EQ.131) THEN LJET=133 GO TO 1 ENDIF RETURN 30 CONTINUE 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---BUT ONLY IF Q**2>10**-4 IF (PB(5).LT.-1.E-2) THEN BREITF=.TRUE. 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) ELSE BREITF=.FALSE. ENDIF ENDIF PTX=0. PTY=0. PF=1.D0 DO 50 IP=1,2 MHEP=IJET(IP) IF (JDAHEP(1,MHEP).EQ.0) THEN C---SPECIAL FOR NON-PARTON JETS IHEP=MHEP GO TO 40 ELSE IST=134+IP DO 35 IHEP=MHEP,NHEP IF (ISTHEP(IHEP).EQ.IST) GO TO 40 35 CONTINUE C---COULDN'T FIND SPACELIKE PARTON CALL HWWARN('HWBJCO',101,*999) ENDIF 40 CALL HWVSCA(3,PF,PHEP(1,IHEP),PS) IF (PTINT(3,IP).GT.0.) 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.BREITF) 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.BREITF) 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.0..AND.PTINF.NE.0.) 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 C (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 46 JHEP=IHEP,KHEP CALL HWUROF(RS,PHEP(1,JHEP),PHEP(1,JHEP)) 46 CONTINUE 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) 50 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 (BREITF) 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.0.) 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 55 IHEP=IJET(2),IEND(2) CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP)) C---BOOST FROM BREIT FRAME IF NECESSARY IF (BREITF) THEN CALL HWUROB(RBR,PHEP(1,IHEP),PHEP(1,IHEP)) CALL HWULOB(PBR,PHEP(1,IHEP),PHEP(1,IHEP)) ENDIF 55 ISTHEP(IHEP)=ISTHEP(IHEP)+10 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.0.) 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 57 IHEP=IJT,IEND(2) CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP)) 57 ISTHEP(IHEP)=ISTHEP(IHEP)+10 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.0.) 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 62 IP=1,2 PB(3)=AL(IP)-(1./AL(IP)) PB(4)=AL(IP)+(1./AL(IP)) IJT=IJET(IP) DO 60 IHEP=IJT,IEND(IP) CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP)) 60 ISTHEP(IHEP)=ISTHEP(IHEP)+10 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 62 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 66 KP=3,NP CALL HWVSUM(4,PHEP(1,IPAR(KP)),PQ,PQ) 66 CONTINUE ENDIF LP=1 ENDIF IF (.NOT.DISLOW) THEN C---FIND JET CM MOMENTA ECM=PQ(5) EMS=0. JETRAD=.FALSE. DO 67 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.0.) CALL HWWARN('HWBJCO',104,*999) 67 CONTINUE PF=1. IF (JETRAD) THEN C---JETS DID RADIATE IF (EMS.GE.ECM) THEN FROST=.TRUE. RETURN ENDIF DO 70 NE=1,NETRY EMS=-ECM DMS=0. DO 68 KP=LP,NP ES=SQRT(PF*PJ(KP)+PM(KP)) EMS=EMS+ES 68 DMS=DMS+PJ(KP)/ES DPF=2.*EMS/DMS IF (DPF.GT.PF) DPF=0.9*PF PF=PF-DPF IF (ABS(DPF).LT.EPS) GO TO 71 70 CONTINUE CALL HWWARN('HWBJCO',105,*999) ENDIF 71 CONTINUE ENDIF C---BOOST PC AND PQ TO BREIT FRAME IF NECESSARY IF (DISPRO.AND.BREITF) 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 80 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.BREITF) 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.BREITF) 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.BREITF) 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.BREITF) THEN CALL HWULOF(PBR,PB,PB) CALL HWUROF(RBR,PB,PB) ENDIF PTCON=PB(1)**2+PB(2)**2 IF (PTCON.NE.0..AND.PTINF.NE.0.) 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 75 JHEP=IHEP,KHEP CALL HWULOF(PB,PHEP(1,JHEP),PHEP(1,JHEP)) CALL HWUROB(RS,PHEP(1,JHEP),PHEP(1,JHEP)) C---BOOST FROM BREIT FRAME IF NECESSARY IF (DISPRO.AND.BREITF) THEN CALL HWUROB(RBR,PHEP(1,JHEP),PHEP(1,JHEP)) CALL HWULOB(PBR,PHEP(1,JHEP),PHEP(1,JHEP)) ENDIF 75 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 80 CONTINUE IF (ISTHEP(ICM).EQ.110) ISTHEP(ICM)=120 ENDIF GO TO 5 999 END