CDECK ID>, HWCBCT. *CMZ :- -20/07/99 10:56:12 by Peter Richardson *-- Author : Peter Richardson C----------------------------------------------------------------------- SUBROUTINE HWCBCT(JHEP,KHEP,THEP,PCL,SPLIT) C----------------------------------------------------------------------- C Subroutine to split a baryonic cluster containing two heavy quarks C Based on HWCCUT C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,QM3,QM4, & PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5), & VCLUS(4),DQM,EMX,EMY,SKAPPA,RKAPPA,VTMP(4), & DELTM,PDIQUK(5),AY(5) INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,THEP,ID1,ID2,ID3,ID4,NTRY, & NTRYMX,J,IB LOGICAL SPLIT EXTERNAL HWUPCM,HWR,HWVDOT PARAMETER(SKAPPA=1.,NTRYMX=100) IF(IERROR.NE.0) RETURN EMC=PCL(5) ID1=IDHW(JHEP) ID2=IDHW(KHEP) ID3=IDHW(THEP) QM1=RMASS(ID1) QM2=RMASS(ID2) QM3=RMASS(ID3) 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.OR. & ID3.EQ.5.OR.ID3.EQ.11) THEN IB=2 ELSE IB=1 ENDIF C-- Set the positon of the cluster to be that of the heavy quark CALL HWVEQU(4,VHEP(1,THEP),VCLUS) C--SPLIT THE BARYONIC CLUSTER INTO A HEAVY FLAVOUR MESON AND A HEAVY C--FLAVOUR BARYON PXY=EMC-QM1-QM2-QM3 20 NTRY=NTRY+1 IF(NTRY.GT.NTRYMX) RETURN 30 EMX=QM1+QM2+PXY*HWR()**PSPLT(IB) EMY= QM3+PXY*HWR()**PSPLT(IB) IF(EMX+EMY.GE.EMC) GOTO 30 C--PULL A LIGHT QUARK PAIR OUT OF THE VACUUM 40 ID4=HWRINT(1,3) IF(QWT(ID4).LT.HWR()) GOTO 40 QM4=RMASS(ID4) C--Now combine particles 3 & 4 into a diquark C--If three also heavy this diquark doesn't exist in HERWIG C--just assume mass is sum of quark masses,as for other diquarks DQM=QM3+QM4 C--Now obtain the masses for the cluster splitting PCX=HWUPCM(EMX,QM1,DQM) IF(PCX.LT.ZERO) GOTO 20 PCY=HWUPCM(EMY,QM2,QM4) IF(PCY.LT.ZERO) GOTO 20 SPLIT=.TRUE. C--Now we've decided which light quark to pull out of the vacuum C--Find the direction of the second heavy quark CALL HWULOF(PCL,PHEP(1,THEP),AX) RCM=1./SQRT(HWVDOT(3,AX,AX)) CALL HWVSCA(3,RCM,AX,AX) C--Construct the new CoM momenta(collinear) PXY=HWUPCM(EMC,EMX,EMY) CALL HWVSCA(3,PXY,AX,PC) C--pc is momenta of Y cluster along 2nd quark dirn in cluster frame PC(4)=SQRT(PXY**2+EMY**2) PC(5)=EMY C--pa is momenta of 2nd quark in Y frame CALL HWVSCA(3,PCY,AX,PA) PA(4)=SQRT(PCY**2+QM3**2) PA(5)=QM3 C--pb is momenta of 2nd quark in cluster frame,pa now momenta of antiquark CALL HWULOB(PC,PA,PB) CALL HWVDIF(4,PC,PB,PA) PA(5)=QM4 LHEP=NHEP+1 MHEP=NHEP+2 C--boost these momenta back to lab frame CALL HWULOB(PCL,PB,PHEP(1,THEP)) CALL HWULOB(PCL,PA,PHEP(1,MHEP)) C--pc now becomes momenta of X cluster in cluster frame CALL HWVSCA(3,-ONE,PC,PC) PC(4)=EMC-PC(4) PC(5)=EMX C--find the dirn of the 1st heavy quark in the X frame C--transform to cluster frame CALL HWULOF(PCL,PHEP(1,JHEP),AY) C--transform to X-frame CALL HWULOF(PC,AY,AY) RCM=1./SQRT(HWVDOT(3,AY,AY)) CALL HWVSCA(3,RCM,AY,AY) C--pa now momenta of 1st havy quark along this dirn CALL HWVSCA(3,PCX,AY,PA) PA(4)=SQRT(PCX**2+QM1**2) PA(5)=QM1 C--pb now momenta of 1st heavy quark in cluster frame then to lab CALL HWULOB(PC,PA,PB) CALL HWULOB(PCL,PB,PHEP(1,JHEP)) C--now find the diquark momenta by momentum conservation DO 50 J=1,4 50 PDIQUK(J)=PCL(J)-PHEP(J,THEP)-PHEP(J,MHEP)-PHEP(J,JHEP) PDIQUK(5)=DQM C--Now obtain the quark momenta from the diquark DO 60 J=1,3 60 PA(J) = 0 PA(4) = QM2 PA(5) = QM2 CALL HWULOB(PDIQUK,PA,PHEP(1,KHEP)) CALL HWVDIF(4,PDIQUK,PHEP(1,KHEP),PHEP(1,LHEP)) 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,VHEP(1,LHEP)) CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP)) C--Relabel the colours of the quarks IDHEP(LHEP) = IDPDG(ID4) IDHEP(MHEP) = IDPDG(ID4) IF(IDHEP(JHEP).GT.0) THEN IDHW(LHEP) = ID4+6 IDHEP(LHEP) = -IDHEP(LHEP) IDHW(MHEP) = ID4 JDAHEP(2,LHEP) = JHEP JMOHEP(2,LHEP) = MHEP JMOHEP(2,MHEP) = JMOHEP(2,JHEP) JDAHEP(2,MHEP) = LHEP JMOHEP(2,JHEP) = LHEP ELSE IDHW(LHEP) = ID4 IDHW(MHEP) = ID4+6 IDHEP(MHEP) = -IDHEP(MHEP) JMOHEP(2,LHEP) = JHEP JDAHEP(2,MHEP) = JDAHEP(2,JHEP) JDAHEP(2,LHEP) = MHEP JMOHEP(2,MHEP) = LHEP JDAHEP(2,JHEP) = LHEP ENDIF ISTHEP(LHEP) = 151 ISTHEP(MHEP) = 151 JMOHEP(1,LHEP) = JMOHEP(1,KHEP) JDAHEP(1,LHEP) = 0 JMOHEP(1,MHEP) = JMOHEP(1,JHEP) JDAHEP(1,MHEP) = 0 NHEP = NHEP+2 999 END CDECK ID>, HWCBVI. *CMZ :- *-- Author : Mark Gibbs modified by Peter Richardson C----------------------------------------------------------------------- SUBROUTINE HWCBVI C----------------------------------------------------------------------- C FINDS UNPAIRED PARTONS AFTER BARYON-NUMBER VIOLATION C MODIFIED FOR RPARITY VIOLATING SUSY C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' COMMON/HWBVIC/NBV,IBV(18) DOUBLE PRECISION HWR,PDQ(5) INTEGER NBV,IBV,JBV,KBV,LBV,IHEP,IP1,IP2,IP3,JP1,JP2,JP3, & HWCBVT,NBR,MBV,IQ1,IQ2,IQ3,ID1,ID2,IDQ,IDIQK(3,3) LOGICAL SPLIT,DUNBV(18) DATA IDIQK/111,110,113,110,109,112,113,112,114/ C---Check for errors IF (IERROR.NE.0) RETURN C---Correct colour connections are gluon splitting CALL HWCCCC C---Reset bvi clustering flag HVFCEN = .FALSE. C---LIST PARTONS WITH WRONG COLOUR PARTNERS-QUARKS ONLY 5 NBV=0 DO 10 IHEP=1,NHEP IF (ISTHEP(IHEP).GT.149.AND.ISTHEP(IHEP).LT.155) THEN IF (QORQQB(IDHW(IHEP))) THEN IF (.NOT.QORQQB(IDHW(JMOHEP(2,IHEP))). & AND.JMOHEP(2,IHEP).GT.6) GOTO 10 ELSE C---Extra check for Gamma's IF (IDHW(IHEP).EQ.59) GO TO 10 C---End of bug fix. IF (QORQQB(IDHW(JDAHEP(2,IHEP)))) GO TO 10 GO TO 10 ENDIF IF(JMOHEP(2,IHEP).LT.6.AND. & .NOT.QBORQQ(IDHW(JMOHEP(2,IHEP)))) GOTO 10 C--new for hard process NBV=NBV+1 IF (NBV.GT.18) CALL HWWARN('HWCBVI',100,*999) IBV(NBV)=IHEP DUNBV(NBV)=.FALSE. ENDIF 10 CONTINUE C--NOW FIND THE ANTIQUARKS WITH WRONG COLOUR CONNECTIONS DO 11 IHEP=1,NHEP IF(ISTHEP(IHEP).GT.149.AND.ISTHEP(IHEP).LT.155) THEN IF(QBORQQ(IDHW(IHEP))) THEN IF(.NOT.QBORQQ(IDHW(JDAHEP(2,IHEP))).AND. & JDAHEP(2,IHEP).GT.6) GO TO 11 ELSE C--Extra check for gamma's IF(IDHW(IHEP).EQ.59) GO TO 11 IF(QBORQQ(IDHW(JMOHEP(2,IHEP)))) GO TO 11 GO TO 11 ENDIF IF(JDAHEP(2,IHEP).LT.6.AND. & .NOT.QORQQB(IDHW(JDAHEP(2,IHEP)))) GOTO 11 NBV=NBV+1 IF(NBV.GT.18) CALL HWWARN('HWCBVI',100,*999) IBV(NBV)=IHEP DUNBV(NBV)=.FALSE. ENDIF 11 CONTINUE IF (NBV.EQ.0) RETURN IF(MOD(NBV,3).NE.0) CALL HWWARN('HWCBVI',101,*999) C---PROCESS FOUND PARTONS, STARTING AT RANDOM POINT IN LIST NBR=NBV*HWR() DO 100 MBV=1,NBV JBV=MBV+NBR IF (JBV.GT.NBV) JBV=JBV-NBV IF (.NOT.DUNBV(JBV)) THEN DUNBV(JBV)=.TRUE. IP1=IBV(JBV) JP1=HWCBVT(IP1) C---FIND ASSOCIATED PARTONS DO 20 KBV=1,NBV IF (.NOT.DUNBV(KBV)) THEN IP2=IBV(KBV) JP2=HWCBVT(IP2) IF (JP2.EQ.JP1) THEN DUNBV(KBV)=.TRUE. DO 15 LBV=1,NBV IF (.NOT.DUNBV(LBV)) THEN IP3=IBV(LBV) JP3=HWCBVT(IP3) IF (JP3.EQ.JP2) THEN DUNBV(LBV)=.TRUE. GO TO 25 ENDIF ENDIF 15 CONTINUE ENDIF ENDIF 20 CONTINUE CALL HWWARN('HWCBVI',102,*999) 25 IQ1=0 C---LOOK FOR DIQUARK IF (ABS(IDHEP(IP1)).GT.100) THEN IQ1=IP1 IQ2=IP2 IQ3=IP3 ELSEIF (ABS(IDHEP(IP2)).GT.100) THEN IQ1=IP2 IQ2=IP3 IQ3=IP1 ELSEIF (ABS(IDHEP(IP3)).GT.100) THEN IQ1=IP3 IQ2=IP1 IQ3=IP2 ENDIF IF (IQ1.EQ.0) THEN C---NO DIQUARKS: COMBINE TWO (ANTI)QUARKS IF (ABS(IDHEP(IP1)).GT.3) THEN IQ1=IP2 IQ2=IP3 IQ3=IP1 ELSEIF (ABS(IDHEP(IP2)).GT.3) THEN IQ1=IP3 IQ2=IP1 IQ3=IP2 ELSE IQ1=IP1 IQ2=IP2 IQ3=IP3 ENDIF ID1=IDHEP(IQ1) ID2=IDHEP(IQ2) C---CHECK FLAVOURS IF (ID1.GT.0.AND.ID1.LT.4.AND. & ID2.GT.0.AND.ID2.LT.4) THEN IDQ=IDIQK(ID1,ID2) ELSEIF (ID1.LT.0.AND.ID1.GT.-4.AND. & ID1.LT.0.AND.ID2.GT.-4) THEN IDQ=IDIQK(-ID1,-ID2)+6 ELSE C---CANT MAKE DIQUARKS WITH HEAVY QUARKS: TRY CLUSTER SPLITTING CALL HWVSUM(4,PHEP(1,IQ1),PHEP(1,IQ2),PDQ) CALL HWUMAS(PDQ) C--Use the original splitting procedure CALL HWCCUT(IQ1,IQ2,PDQ,.FALSE.,SPLIT) IF(SPLIT) GOTO 5 C--If it fails try the new procedure CALL HWVSUM(4,PDQ,PHEP(1,IQ3),PDQ) CALL HWUMAS(PDQ) IF(ABS(ID1).GT.3) THEN CALL HWCBCT(IQ3,IQ2,IQ1,PDQ,SPLIT) ELSEIF(ABS(ID2).GT.3) THEN CALL HWCBCT(IQ3,IQ1,IQ2,PDQ,SPLIT) ELSE CALL HWWARN('HWCBVI',100,*999) ENDIF IF (SPLIT) GO TO 5 C---Unable to form cluster; dispose of event CALL HWWARN('HWCBVI',-3,*999) ENDIF C---OVERWRITE FIRST AND CANCEL SECOND IDHW(IQ1)=IDQ IDHEP(IQ1)=IDPDG(IDQ) CALL HWVSUM(4,PHEP(1,IQ1),PHEP(1,IQ2),PHEP(1,IQ1)) CALL HWUMAS(PHEP(1,IQ1)) ISTHEP(IQ2)=0 C---REMAKE COLOUR CONNECTIONS IF (QORQQB(IDQ)) THEN JMOHEP(2,IQ1)=IQ3 JDAHEP(2,IQ3)=IQ1 ELSE JDAHEP(2,IQ1)=IQ3 JMOHEP(2,IQ3)=IQ1 ENDIF ELSE C---SPLIT A DIQUARK NHEP=NHEP+1 CALL HWVSCA(5,HALF,PHEP(1,IQ1),PHEP(1,IQ1)) CALL HWVEQU(5,PHEP(1,IQ1),PHEP(1,NHEP)) ISTHEP(NHEP)=150 JMOHEP(1,NHEP)=JMOHEP(1,IQ1) JDAHEP(1,NHEP)=0 C---FIND FLAVOURS IDQ=IDHW(IQ1) DO 30 ID2=1,3 DO 30 ID1=1,3 IF (IDIQK(ID1,ID2).EQ.IDQ) THEN IDHW(IQ1)=ID1 IDHW(NHEP)=ID2 C---REMAKE COLOUR CONNECTIONS (DIQUARK) JMOHEP(2,IQ1)=IQ2 JMOHEP(2,IQ2)=NHEP JMOHEP(2,IQ3)=IQ1 JMOHEP(2,NHEP)=IQ3 JDAHEP(2,IQ1)=IQ3 JDAHEP(2,IQ2)=IQ1 JDAHEP(2,IQ3)=NHEP JDAHEP(2,NHEP)=IQ2 GO TO 35 ELSEIF (IDIQK(ID1,ID2).EQ.IDQ-6) THEN IDHW(IQ1)=ID1+6 IDHW(NHEP)=ID2+6 C---REMAKE COLOUR CONNECTIONS (ANTIDIQUARK) JMOHEP(2,IQ1)=IQ3 JMOHEP(2,IQ2)=IQ1 JMOHEP(2,IQ3)=NHEP JMOHEP(2,NHEP)=IQ2 JDAHEP(2,IQ1)=IQ2 JDAHEP(2,IQ2)=NHEP JDAHEP(2,IQ3)=IQ1 JDAHEP(2,NHEP)=IQ3 GO TO 35 ENDIF 30 CONTINUE CALL HWWARN('HWCBVI',104,*999) 35 IDHEP(IQ1)=IDPDG(IDHW(IQ1)) IDHEP(NHEP)=IDPDG(IDHW(NHEP)) ENDIF ENDIF 100 CONTINUE RETURN 999 END CDECK ID>, HWCBVT. *CMZ :- *-- Author : Peter Richardson C----------------------------------------------------------------------- FUNCTION HWCBVT(IP) C----------------------------------------------------------------------- C Function to find the baryon number violating vertex a parton came from C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER HWCBVT,IP,JP(2),KP,I,J,ID,TYPE,IDM,IDM2,IDM3,IDM4 JP(1) = IP ID = IDHW(IP) IF(ID.LE.6.OR.(ID.GE.115.AND.ID.LE.120)) THEN JP(2) = JMOHEP(2,IP) ELSE JP(2) = JDAHEP(2,IP) ENDIF DO I=1,2 IDM = JMOHEP(1,JMOHEP(1,JMOHEP(1,JMOHEP(1,JP(I))))) IF(IDHW(IDM).EQ.6.OR.IDHW(IDM).EQ.12) THEN JP(I)=IDM ENDIF ENDDO DO J=1,7 DO I=1,2 KP = JMOHEP(1,JP(I)) IDM = IDHW(KP) IDM2 = IDHW(JDAHEP(1,KP)) IDM3 = IDHW(JDAHEP(2,KP)) IDM4 = IDHW(JDAHEP(1,KP)+1) IF((ISTHEP(KP).EQ.155.AND. & ((IDM.GE.449.AND.IDM.LE.457.AND.IDM2.LE.12.AND. & IDM3.LE.12.AND.IDM4.LE.12).OR. & (((IDM.GE.411.AND.IDM.LE.424).OR.IDM.EQ.405.OR.IDM.EQ.406) & .AND.IDM2.LE.12.AND.IDM3.LE.12))) & .OR.(IDM.EQ.15.AND.IDM2.LE.12.AND. & IDHW(JMOHEP(1,KP)).LE.12.AND. & IDHW(JMOHEP(2,KP)).LE.12.AND.IDM3.GE.449.AND. & IDM3.LE.457).OR. & (IDM.EQ.15.AND.IDM2.GE.198.AND.IDM2.LE.200. & AND.ABS(IDPDG(IDM3)).GT.1000000)) THEN IF(IDHW(KP).EQ.449.AND.JDAHEP(1,KP).EQ.JP(I)) THEN KP = JMOHEP(1,KP) ELSEIF(IDHW(KP).EQ.15) THEN TYPE=IDHW(JDAHEP(1,KP)) IF(TYPE.GE.7.AND.TYPE.LE.12.AND. & JMOHEP(2,JDAHEP(2,KP)).EQ.JP(I)) THEN KP=IP ELSEIF(TYPE.LE.6.AND. & JDAHEP(2,JDAHEP(2,KP)).EQ.JP(I)) THEN KP=IP ELSE HWCBVT = KP RETURN ENDIF ELSE HWCBVT = KP RETURN ENDIF ENDIF JP(I) =KP ENDDO ENDDO HWCBVT = 0 999 END CDECK ID>, HWCCCC. *CMZ :- *-- Author : Peter Richardson C----------------------------------------------------------------------- SUBROUTINE HWCCCC C----------------------------------------------------------------------- C Subroutine to correct colour connections after the gluon splitting C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' INTEGER IHEP,STFSPT,LHEP,MHEP,RHEP IF(IERROR.NE.0) RETURN C--Find the first particle in the event record with status 150 DO IHEP=1,NHEP IF(ISTHEP(IHEP).GE.150.AND.ISTHEP(IHEP).LE.154) THEN STFSPT = IHEP GOTO 10 ENDIF ENDDO 10 CONTINUE C--Now find any that are colour connected to earlier particles C--in the event record DO IHEP=STFSPT,NHEP C--First the quarks and antidiquarks IF(IDHW(IHEP).LT.6.OR. & (IDHW(IHEP).GE.115.AND.IDHW(IHEP).LE.120)) THEN IF(JMOHEP(2,IHEP).LT.STFSPT) THEN LHEP = IHEP MHEP = JMOHEP(2,IHEP) RHEP = MHEP IF(MHEP.GT.6) RHEP = JDAHEP(1,MHEP) C--As from Rparity connect to particle not to antiparticle IF(IDHW(MHEP).NE.13) THEN JMOHEP(2,LHEP) = RHEP ELSE RHEP = RHEP+1 JMOHEP(2,LHEP) = RHEP ENDIF ENDIF ENDIF C--Now the antiquarks IF((IDHW(IHEP).GT.6.AND.IDHW(IHEP).LE.12).OR. & (IDHW(IHEP).GE.109.AND.IDHW(IHEP).LE.114)) THEN IF(JDAHEP(2,IHEP).LT.STFSPT) THEN LHEP = IHEP MHEP = JDAHEP(2,IHEP) RHEP = MHEP IF(MHEP.GT.6) RHEP = JDAHEP(1,MHEP) C--As from Rparity connect to antiparticle not particle IF(IDHW(MHEP).NE.13) THEN JDAHEP(2,LHEP) = RHEP ELSE JDAHEP(2,LHEP) = RHEP ENDIF ENDIF ENDIF ENDDO END