4 *CMZ :- -20/07/99 10:56:12 by Peter Richardson
6 *-- Author : Peter Richardson
8 C-----------------------------------------------------------------------
10 SUBROUTINE HWCBCT(JHEP,KHEP,THEP,PCL,SPLIT)
12 C-----------------------------------------------------------------------
14 C Subroutine to split a baryonic cluster containing two heavy quarks
18 C-----------------------------------------------------------------------
20 INCLUDE 'HERWIG61.INC'
22 DOUBLE PRECISION HWUPCM,HWR,HWVDOT,EMC,QM1,QM2,QM3,QM4,
24 & PXY,PCX,PCY,RCM,PCL(5),AX(5),PA(5),PB(5),PC(5),
26 & VCLUS(4),DQM,EMX,EMY,SKAPPA,RKAPPA,VTMP(4),
28 & DELTM,PDIQUK(5),AY(5)
30 INTEGER HWRINT,JHEP,KHEP,LHEP,MHEP,THEP,ID1,ID2,ID3,ID4,NTRY,
36 EXTERNAL HWUPCM,HWR,HWVDOT
38 PARAMETER(SKAPPA=1.,NTRYMX=100)
40 IF(IERROR.NE.0) RETURN
60 C Decide if cluster contains a b-(anti)quark
62 IF (ID1.EQ.5.OR.ID1.EQ.11.OR.ID2.EQ.5.OR.ID2.EQ.11.OR.
64 & ID3.EQ.5.OR.ID3.EQ.11) THEN
74 C-- Set the positon of the cluster to be that of the heavy quark
76 CALL HWVEQU(4,VHEP(1,THEP),VCLUS)
78 C--SPLIT THE BARYONIC CLUSTER INTO A HEAVY FLAVOUR MESON AND A HEAVY
86 IF(NTRY.GT.NTRYMX) RETURN
88 30 EMX=QM1+QM2+PXY*HWR()**PSPLT(IB)
90 EMY= QM3+PXY*HWR()**PSPLT(IB)
92 IF(EMX+EMY.GE.EMC) GOTO 30
94 C--PULL A LIGHT QUARK PAIR OUT OF THE VACUUM
98 IF(QWT(ID4).LT.HWR()) GOTO 40
102 C--Now combine particles 3 & 4 into a diquark
104 C--If three also heavy this diquark doesn't exist in HERWIG
106 C--just assume mass is sum of quark masses,as for other diquarks
110 C--Now obtain the masses for the cluster splitting
112 PCX=HWUPCM(EMX,QM1,DQM)
114 IF(PCX.LT.ZERO) GOTO 20
116 PCY=HWUPCM(EMY,QM2,QM4)
118 IF(PCY.LT.ZERO) GOTO 20
122 C--Now we've decided which light quark to pull out of the vacuum
124 C--Find the direction of the second heavy quark
126 CALL HWULOF(PCL,PHEP(1,THEP),AX)
128 RCM=1./SQRT(HWVDOT(3,AX,AX))
130 CALL HWVSCA(3,RCM,AX,AX)
132 C--Construct the new CoM momenta(collinear)
134 PXY=HWUPCM(EMC,EMX,EMY)
136 CALL HWVSCA(3,PXY,AX,PC)
138 C--pc is momenta of Y cluster along 2nd quark dirn in cluster frame
140 PC(4)=SQRT(PXY**2+EMY**2)
144 C--pa is momenta of 2nd quark in Y frame
146 CALL HWVSCA(3,PCY,AX,PA)
148 PA(4)=SQRT(PCY**2+QM3**2)
152 C--pb is momenta of 2nd quark in cluster frame,pa now momenta of antiquark
154 CALL HWULOB(PC,PA,PB)
156 CALL HWVDIF(4,PC,PB,PA)
164 C--boost these momenta back to lab frame
166 CALL HWULOB(PCL,PB,PHEP(1,THEP))
168 CALL HWULOB(PCL,PA,PHEP(1,MHEP))
170 C--pc now becomes momenta of X cluster in cluster frame
172 CALL HWVSCA(3,-ONE,PC,PC)
178 C--find the dirn of the 1st heavy quark in the X frame
180 C--transform to cluster frame
182 CALL HWULOF(PCL,PHEP(1,JHEP),AY)
184 C--transform to X-frame
186 CALL HWULOF(PC,AY,AY)
188 RCM=1./SQRT(HWVDOT(3,AY,AY))
190 CALL HWVSCA(3,RCM,AY,AY)
192 C--pa now momenta of 1st havy quark along this dirn
194 CALL HWVSCA(3,PCX,AY,PA)
196 PA(4)=SQRT(PCX**2+QM1**2)
200 C--pb now momenta of 1st heavy quark in cluster frame then to lab
202 CALL HWULOB(PC,PA,PB)
204 CALL HWULOB(PCL,PB,PHEP(1,JHEP))
206 C--now find the diquark momenta by momentum conservation
210 50 PDIQUK(J)=PCL(J)-PHEP(J,THEP)-PHEP(J,MHEP)-PHEP(J,JHEP)
214 C--Now obtain the quark momenta from the diquark
224 CALL HWULOB(PDIQUK,PA,PHEP(1,KHEP))
226 CALL HWVDIF(4,PDIQUK,PHEP(1,KHEP),PHEP(1,LHEP))
228 C--Construct new vertex positions
232 CALL HWVSCA(3,RKAPPA,AX,AX)
234 DELTM=(EMX-EMY)*(EMX+EMY)/(TWO*EMC)
236 CALL HWVSCA(3,DELTM,AX,VTMP)
238 VTMP(4)=(HALF*EMC-PXY)*RKAPPA
240 CALL HWULB4(PCL,VTMP,VTMP)
242 CALL HWVSUM(4,VTMP,VCLUS,VHEP(1,LHEP))
244 CALL HWVEQU(4,VHEP(1,LHEP),VHEP(1,MHEP))
246 C--Relabel the colours of the quarks
248 IDHEP(LHEP) = IDPDG(ID4)
250 IDHEP(MHEP) = IDPDG(ID4)
252 IF(IDHEP(JHEP).GT.0) THEN
256 IDHEP(LHEP) = -IDHEP(LHEP)
260 JDAHEP(2,LHEP) = JHEP
262 JMOHEP(2,LHEP) = MHEP
264 JMOHEP(2,MHEP) = JMOHEP(2,JHEP)
266 JDAHEP(2,MHEP) = LHEP
268 JMOHEP(2,JHEP) = LHEP
276 IDHEP(MHEP) = -IDHEP(MHEP)
278 JMOHEP(2,LHEP) = JHEP
280 JDAHEP(2,MHEP) = JDAHEP(2,JHEP)
282 JDAHEP(2,LHEP) = MHEP
284 JMOHEP(2,MHEP) = LHEP
286 JDAHEP(2,JHEP) = LHEP
294 JMOHEP(1,LHEP) = JMOHEP(1,KHEP)
298 JMOHEP(1,MHEP) = JMOHEP(1,JHEP)
310 *-- Author : Mark Gibbs modified by Peter Richardson
312 C-----------------------------------------------------------------------
316 C-----------------------------------------------------------------------
318 C FINDS UNPAIRED PARTONS AFTER BARYON-NUMBER VIOLATION
320 C MODIFIED FOR RPARITY VIOLATING SUSY
322 C-----------------------------------------------------------------------
324 INCLUDE 'HERWIG61.INC'
326 COMMON/HWBVIC/NBV,IBV(18)
328 DOUBLE PRECISION HWR,PDQ(5)
330 INTEGER NBV,IBV,JBV,KBV,LBV,IHEP,IP1,IP2,IP3,JP1,JP2,JP3,
332 & HWCBVT,NBR,MBV,IQ1,IQ2,IQ3,ID1,ID2,IDQ,IDIQK(3,3)
334 LOGICAL SPLIT,DUNBV(18)
336 DATA IDIQK/111,110,113,110,109,112,113,112,114/
340 IF (IERROR.NE.0) RETURN
342 C---Correct colour connections are gluon splitting
346 C---Reset bvi clustering flag
350 C---LIST PARTONS WITH WRONG COLOUR PARTNERS-QUARKS ONLY
356 IF (ISTHEP(IHEP).GT.149.AND.ISTHEP(IHEP).LT.155) THEN
358 IF (QORQQB(IDHW(IHEP))) THEN
360 IF (.NOT.QORQQB(IDHW(JMOHEP(2,IHEP))).
362 & AND.JMOHEP(2,IHEP).GT.6) GOTO 10
366 C---Extra check for Gamma's
368 IF (IDHW(IHEP).EQ.59) GO TO 10
372 IF (QORQQB(IDHW(JDAHEP(2,IHEP)))) GO TO 10
378 IF(JMOHEP(2,IHEP).LT.6.AND.
380 & .NOT.QBORQQ(IDHW(JMOHEP(2,IHEP)))) GOTO 10
382 C--new for hard process
386 IF (NBV.GT.18) CALL HWWARN('HWCBVI',100,*999)
396 C--NOW FIND THE ANTIQUARKS WITH WRONG COLOUR CONNECTIONS
400 IF(ISTHEP(IHEP).GT.149.AND.ISTHEP(IHEP).LT.155) THEN
402 IF(QBORQQ(IDHW(IHEP))) THEN
404 IF(.NOT.QBORQQ(IDHW(JDAHEP(2,IHEP))).AND.
406 & JDAHEP(2,IHEP).GT.6) GO TO 11
410 C--Extra check for gamma's
412 IF(IDHW(IHEP).EQ.59) GO TO 11
414 IF(QBORQQ(IDHW(JMOHEP(2,IHEP)))) GO TO 11
420 IF(JDAHEP(2,IHEP).LT.6.AND.
422 & .NOT.QORQQB(IDHW(JDAHEP(2,IHEP)))) GOTO 11
426 IF(NBV.GT.18) CALL HWWARN('HWCBVI',100,*999)
438 IF(MOD(NBV,3).NE.0) CALL HWWARN('HWCBVI',101,*999)
440 C---PROCESS FOUND PARTONS, STARTING AT RANDOM POINT IN LIST
448 IF (JBV.GT.NBV) JBV=JBV-NBV
450 IF (.NOT.DUNBV(JBV)) THEN
458 C---FIND ASSOCIATED PARTONS
462 IF (.NOT.DUNBV(KBV)) THEN
474 IF (.NOT.DUNBV(LBV)) THEN
498 CALL HWWARN('HWCBVI',102,*999)
504 IF (ABS(IDHEP(IP1)).GT.100) THEN
512 ELSEIF (ABS(IDHEP(IP2)).GT.100) THEN
520 ELSEIF (ABS(IDHEP(IP3)).GT.100) THEN
532 C---NO DIQUARKS: COMBINE TWO (ANTI)QUARKS
534 IF (ABS(IDHEP(IP1)).GT.3) THEN
542 ELSEIF (ABS(IDHEP(IP2)).GT.3) THEN
566 IF (ID1.GT.0.AND.ID1.LT.4.AND.
568 & ID2.GT.0.AND.ID2.LT.4) THEN
572 ELSEIF (ID1.LT.0.AND.ID1.GT.-4.AND.
574 & ID1.LT.0.AND.ID2.GT.-4) THEN
576 IDQ=IDIQK(-ID1,-ID2)+6
580 C---CANT MAKE DIQUARKS WITH HEAVY QUARKS: TRY CLUSTER SPLITTING
582 CALL HWVSUM(4,PHEP(1,IQ1),PHEP(1,IQ2),PDQ)
586 C--Use the original splitting procedure
588 CALL HWCCUT(IQ1,IQ2,PDQ,.FALSE.,SPLIT)
592 C--If it fails try the new procedure
594 CALL HWVSUM(4,PDQ,PHEP(1,IQ3),PDQ)
598 IF(ABS(ID1).GT.3) THEN
600 CALL HWCBCT(IQ3,IQ2,IQ1,PDQ,SPLIT)
602 ELSEIF(ABS(ID2).GT.3) THEN
604 CALL HWCBCT(IQ3,IQ1,IQ2,PDQ,SPLIT)
608 CALL HWWARN('HWCBVI',100,*999)
614 C---Unable to form cluster; dispose of event
616 CALL HWWARN('HWCBVI',-3,*999)
620 C---OVERWRITE FIRST AND CANCEL SECOND
624 IDHEP(IQ1)=IDPDG(IDQ)
626 CALL HWVSUM(4,PHEP(1,IQ1),PHEP(1,IQ2),PHEP(1,IQ1))
628 CALL HWUMAS(PHEP(1,IQ1))
632 C---REMAKE COLOUR CONNECTIONS
634 IF (QORQQB(IDQ)) THEN
654 CALL HWVSCA(5,HALF,PHEP(1,IQ1),PHEP(1,IQ1))
656 CALL HWVEQU(5,PHEP(1,IQ1),PHEP(1,NHEP))
660 JMOHEP(1,NHEP)=JMOHEP(1,IQ1)
672 IF (IDIQK(ID1,ID2).EQ.IDQ) THEN
678 C---REMAKE COLOUR CONNECTIONS (DIQUARK)
698 ELSEIF (IDIQK(ID1,ID2).EQ.IDQ-6) THEN
704 C---REMAKE COLOUR CONNECTIONS (ANTIDIQUARK)
728 CALL HWWARN('HWCBVI',104,*999)
730 35 IDHEP(IQ1)=IDPDG(IDHW(IQ1))
732 IDHEP(NHEP)=IDPDG(IDHW(NHEP))
748 *-- Author : Peter Richardson
750 C-----------------------------------------------------------------------
754 C-----------------------------------------------------------------------
756 C Function to find the baryon number violating vertex a parton came from
758 C-----------------------------------------------------------------------
760 INCLUDE 'HERWIG61.INC'
762 INTEGER HWCBVT,IP,JP(2),KP,I,J,ID,TYPE,IDM,IDM2,IDM3,IDM4
768 IF(ID.LE.6.OR.(ID.GE.115.AND.ID.LE.120)) THEN
780 IDM = JMOHEP(1,JMOHEP(1,JMOHEP(1,JMOHEP(1,JP(I)))))
782 IF(IDHW(IDM).EQ.6.OR.IDHW(IDM).EQ.12) THEN
798 IDM2 = IDHW(JDAHEP(1,KP))
800 IDM3 = IDHW(JDAHEP(2,KP))
802 IDM4 = IDHW(JDAHEP(1,KP)+1)
804 IF((ISTHEP(KP).EQ.155.AND.
806 & ((IDM.GE.449.AND.IDM.LE.457.AND.IDM2.LE.12.AND.
808 & IDM3.LE.12.AND.IDM4.LE.12).OR.
810 & (((IDM.GE.411.AND.IDM.LE.424).OR.IDM.EQ.405.OR.IDM.EQ.406)
812 & .AND.IDM2.LE.12.AND.IDM3.LE.12)))
814 & .OR.(IDM.EQ.15.AND.IDM2.LE.12.AND.
816 & IDHW(JMOHEP(1,KP)).LE.12.AND.
818 & IDHW(JMOHEP(2,KP)).LE.12.AND.IDM3.GE.449.AND.
822 & (IDM.EQ.15.AND.IDM2.GE.198.AND.IDM2.LE.200.
824 & AND.ABS(IDPDG(IDM3)).GT.1000000)) THEN
826 IF(IDHW(KP).EQ.449.AND.JDAHEP(1,KP).EQ.JP(I)) THEN
830 ELSEIF(IDHW(KP).EQ.15) THEN
832 TYPE=IDHW(JDAHEP(1,KP))
834 IF(TYPE.GE.7.AND.TYPE.LE.12.AND.
836 & JMOHEP(2,JDAHEP(2,KP)).EQ.JP(I)) THEN
840 ELSEIF(TYPE.LE.6.AND.
842 & JDAHEP(2,JDAHEP(2,KP)).EQ.JP(I)) THEN
878 *-- Author : Peter Richardson
880 C-----------------------------------------------------------------------
884 C-----------------------------------------------------------------------
886 C Subroutine to correct colour connections after the gluon splitting
888 C-----------------------------------------------------------------------
890 INCLUDE 'HERWIG61.INC'
892 INTEGER IHEP,STFSPT,LHEP,MHEP,RHEP
894 IF(IERROR.NE.0) RETURN
896 C--Find the first particle in the event record with status 150
900 IF(ISTHEP(IHEP).GE.150.AND.ISTHEP(IHEP).LE.154) THEN
912 C--Now find any that are colour connected to earlier particles
914 C--in the event record
918 C--First the quarks and antidiquarks
920 IF(IDHW(IHEP).LT.6.OR.
922 & (IDHW(IHEP).GE.115.AND.IDHW(IHEP).LE.120)) THEN
924 IF(JMOHEP(2,IHEP).LT.STFSPT) THEN
928 MHEP = JMOHEP(2,IHEP)
932 IF(MHEP.GT.6) RHEP = JDAHEP(1,MHEP)
934 C--As from Rparity connect to particle not to antiparticle
936 IF(IDHW(MHEP).NE.13) THEN
938 JMOHEP(2,LHEP) = RHEP
944 JMOHEP(2,LHEP) = RHEP
952 C--Now the antiquarks
954 IF((IDHW(IHEP).GT.6.AND.IDHW(IHEP).LE.12).OR.
956 & (IDHW(IHEP).GE.109.AND.IDHW(IHEP).LE.114)) THEN
958 IF(JDAHEP(2,IHEP).LT.STFSPT) THEN
962 MHEP = JDAHEP(2,IHEP)
966 IF(MHEP.GT.6) RHEP = JDAHEP(1,MHEP)
968 C--As from Rparity connect to antiparticle not particle
970 IF(IDHW(MHEP).NE.13) THEN
972 JDAHEP(2,LHEP) = RHEP
976 JDAHEP(2,LHEP) = RHEP