Coding rule violations corrected.

[u/mrichter/AliRoot.git] / HERWIG / src / hwcbct.f

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