Coding rule violations corrected.

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

CDECK  ID>, HWBTIM.

*CMZ :-        -26/04/91  14.27.17  by  Federico Carminati

*-- Author :    Ian Knowles

C-----------------------------------------------------------------------

      SUBROUTINE HWBTIM(INITBR,INTERF)

C-----------------------------------------------------------------------

C     Constructs full 4-momentum & production vertices in time-like jet

C     initiated by INITBR, interference partner INTERF and spin density

C     RHOPAR(INITBR). DECPAR(INITBR) returns jet's spin density matrix.

C     Includes azimuthal angular correlations between branching planes

C     due to spin (if AZSPIN) using the algorithm of Knowles & Collins.

C     Ses Nucl. Phys. B304 (1988) 794 & Comp. Phys. Comm. 58 (1990) 271.

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      DOUBLE PRECISION HWR,DMIN,PT,EIKON,EINUM,EIDEN1,EIDEN2,EISCR,

     & WT,SPIN,Z1,Z2,PRMAX,CAZ,CX,SX,ROHEP(3),RMAT(3,3),ZERO2(2)

      INTEGER INITBR,INTERF,IPAR,JPAR,KPAR,LPAR,MPAR,NTRY,JOLD

      LOGICAL EICOR,SWAP

      EXTERNAL HWR

      DATA ZERO2,DMIN/ZERO,ZERO,1.D-15/

      IF (IERROR.NE.0) RETURN

      JPAR=INITBR

      KPAR=INTERF

      IF ((JDAPAR(1,JPAR).NE.0).OR.(IDPAR(JPAR).EQ.13)) GOTO 30

C No branching, assign decay matrix

      CALL HWVZRO(2,DECPAR(1,JPAR))

      RETURN

C Advance up the leader

C     Find the parent and partner of J

  10  IPAR=JMOPAR(1,JPAR)

      KPAR=JPAR+1

C Generate new Rho

      IF (JMOPAR(1,KPAR).EQ.IPAR) THEN

C        Generate Rho'

         CALL HWBAZF(IPAR,JPAR,PHIPAR(1,IPAR),RHOPAR(1,IPAR),

     &                                   ZERO2,RHOPAR(1,JPAR))

      ELSE

         KPAR=JPAR-1

         IF (JMOPAR(1,KPAR).NE.IPAR)

     &   CALL HWWARN('HWBTIM',100,*999)

C        Generate Rho''

         CALL HWBAZF(IPAR,KPAR,RHOPAR(1,IPAR),PHIPAR(1,IPAR),

     &                         DECPAR(1,KPAR),RHOPAR(1,JPAR))

      ENDIF

C Generate azimuthal angle of J's branching

  30  IF (JDAPAR(1,JPAR).EQ.0) THEN

C        Final state gluon

         CALL HWVZRO(2,DECPAR(1,JPAR))

         IF (JPAR.EQ.INITBR) RETURN

         GOTO 70

      ELSE

C Assign an angle to a branching using an M-function

C        Find the daughters of J

         LPAR=JDAPAR(1,JPAR)

         MPAR=JDAPAR(2,JPAR)

C Soft correlations

         CALL HWUROT(PPAR(1,JPAR), ONE,ZERO,RMAT)

         CALL HWUROF(RMAT,PPAR(1,KPAR),ROHEP)

         PT=MAX(SQRT(ROHEP(1)*ROHEP(1)+ROHEP(2)*ROHEP(2)),DMIN)

         EIKON=1.

         SWAP=.FALSE.

         EICOR=AZSOFT.AND.((IDPAR(LPAR).EQ.13).OR.(IDPAR(MPAR).EQ.13))

         IF (EICOR) THEN

C           Rearrange s.t. LPAR is the (softest) gluon

            IF (IDPAR(MPAR).EQ.13) THEN

               IF (IDPAR(LPAR).NE.13.OR.

     &             PPAR(4,MPAR).LT.PPAR(4,LPAR)) THEN

                  SWAP=.TRUE.

                  LPAR=MPAR

                  MPAR=LPAR-1

               ENDIF

            ENDIF

            EINUM=(PPAR(4,KPAR)*PPAR(4,LPAR))

     &        *ABS(PPAR(2,LPAR)-PPAR(2,MPAR))

            EIDEN1=(PPAR(4,KPAR)*PPAR(4,LPAR))-ROHEP(3)*PPAR(3,LPAR)

            EIDEN2=PT*ABS(PPAR(1,LPAR))

            EISCR=1.-(PPAR(5,MPAR)/PPAR(4,MPAR))**2

     &           /MIN(PPAR(2,LPAR),PPAR(2,MPAR))

            EIKON=EISCR+EINUM/MAX(EIDEN1-EIDEN2,DMIN)

         ENDIF

C Spin correlations

         WT=0.

         SPIN=1.

         IF (AZSPIN) THEN

            Z1=PPAR(4,LPAR)/PPAR(4,JPAR)

            Z2=1.-Z1

            IF (IDPAR(JPAR).EQ.13.AND.IDPAR(LPAR).EQ.13) THEN

               WT=Z1*Z2/(Z1/Z2+Z2/Z1+Z1*Z2)

            ELSEIF (IDPAR(JPAR).EQ.13.AND.IDPAR(LPAR).LT.13) THEN

               WT=-2.*Z1*Z2/(Z1*Z1+Z2*Z2)

            ENDIF

         ENDIF

C Assign the azimuthal angle

         PRMAX=(1.+ABS(WT))*EIKON

         NTRY=0

   50    NTRY=NTRY+1

         IF (NTRY.GT.NBTRY) CALL HWWARN('HWBTIM',101,*999)

         CALL HWRAZM( ONE,CX,SX)

         CALL HWUROT(PPAR(1,JPAR),CX,SX,RMAT)

C Determine the angle between the branching planes

         CALL HWUROF(RMAT,PPAR(1,KPAR),ROHEP)

         CAZ=ROHEP(1)/PT

         PHIPAR(1,JPAR)=2.*CAZ*CAZ-1.

         PHIPAR(2,JPAR)=2.*CAZ*ROHEP(2)/PT

         IF (EICOR) EIKON=EISCR+EINUM/MAX(EIDEN1-EIDEN2*CAZ,DMIN)

         IF (AZSPIN) SPIN=1.+WT*(RHOPAR(1,JPAR)*PHIPAR(1,JPAR)

     &                          +RHOPAR(2,JPAR)*PHIPAR(2,JPAR))

         IF (SPIN*EIKON.LT.HWR()*PRMAX) GOTO 50

C Construct full 4-momentum of L and M

         JOLD=JPAR

         IF (SWAP) THEN

           PPAR(1,LPAR)=-PPAR(1,LPAR)

           PPAR(1,MPAR)=-PPAR(1,MPAR)

           JPAR=MPAR

         ELSE

           JPAR=LPAR

         ENDIF

         PPAR(2,LPAR)=0.

         CALL HWUROB(RMAT,PPAR(1,LPAR),PPAR(1,LPAR))

         PPAR(2,MPAR)=0.

         CALL HWUROB(RMAT,PPAR(1,MPAR),PPAR(1,MPAR))

C Assign production vertex to L and M

         CALL HWUDKL(IDPAR(JOLD),PPAR(1,JOLD),VPAR(1,LPAR))

         CALL HWVSUM(4,VPAR(1,JOLD),VPAR(1,LPAR),VPAR(1,LPAR))

         CALL HWVEQU(4,VPAR(1,LPAR),VPAR(1,MPAR))

      ENDIF

  60  IF (JDAPAR(1,JPAR).NE.0) GOTO 10

C Assign decay matrix

      CALL HWVZRO(2,DECPAR(1,JPAR))

C Backtrack down the leader

  70  IPAR=JMOPAR(1,JPAR)

      KPAR=JDAPAR(1,IPAR)

      IF (KPAR.EQ.JPAR) THEN

C        Develop the side branch

         JPAR=JDAPAR(2,IPAR)

         GOTO 60

      ELSE

C        Construct decay matrix

         CALL HWBAZF(IPAR,KPAR,DECPAR(1,JPAR),DECPAR(1,KPAR),

     &                         PHIPAR(1,IPAR),DECPAR(1,IPAR))

      ENDIF

      IF (IPAR.EQ.INITBR) RETURN

      JPAR=IPAR

      GOTO 70

  999 END

CDECK  ID>, HWBTOP.

*CMZ :-        -14/10/99  18.04.56  by  Mike Seymour

*-- Author :    Gennaro Corcella

C-----------------------------------------------------------------------

      SUBROUTINE HWBTOP

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      DOUBLE PRECISION HWBVMC,HWR,HWUALF,HWUSQR,X(3),W,

     & X3MIN,X3MAX,X1MIN,X1MAX,QSCALE,GLUFAC,R(3,3),M(3),

     & E(3),AW,PTSQ,EM,EPS,MASDEP,A,B,C,GAMDEP,LAMBDA,

     & PW(5),PT(5),PW1(5),CS,SN,EPG,QQ,RR,CC

      INTEGER ID,ID3,IHEP,KHEP,WHEP,ICMF,K

      EXTERNAL HWBVMC,HWUALF,HWUSQR,HWR

      LAMBDA(A,B,C)=(A**2+B**2+C**2-2*A*B-2*B*C-2*C*A)/(4*A)

C---FIND AN UNTREATED CMF

      ICMF=0

      DO 10 IHEP=1,NHEP

C----FIND A DECAYING TOP QUARK

 10     IF (ISTHEP(IHEP).EQ.155.AND.ISTHEP(JDAHEP(1,IHEP)).EQ.113

     &       .AND.(IDHW(IHEP).EQ.6.OR.IDHW(IHEP).EQ.12))

     &       ICMF=IHEP

      IF (ICMF.EQ.0) RETURN

      EM=PHEP(5,ICMF)

      X3MIN=2*GCUTME/EM

C---GENERATE X(1),X(3) ACCORDING TO 1/((1-X(1))*X(3)**2)

 100  CONTINUE

C-----AW=(MW/MT)**2

      AW=(PHEP(5,JDAHEP(1,ICMF))/EM)**2

C---CHOOSE X3

      X3MAX=1-AW

      X(3)=X3MIN*X3MAX/(X3MIN+(X3MAX-X3MIN)*HWR())

C--CC, QQ AND RR ARE THE VARIABLE DEFINED IN OUR PAPER

C--IN ORDER TO SOLVE THE CUBIC EQUATION

      CC=(1-AW)**2/4

      QQ=(AW**2-4*(1-X(3))*(2-CC-X(3))-2*AW*(3+2*X(3)))/3

     &     -((3+2*AW-4*X(3))**2)/9

      RR=((3+2*AW-4*X(3))*(AW**2-4*(1-X(3))*(2-CC-X(3))

     &     -2*AW*(3+2*X(3)))-3*(AW*(4-AW)*(2-CC)+(1-CC)

     &     *(2*(1-X(3))-AW)**2))/6-(ONE/27)*(3+2*AW-4*X(3))**3

C---CHOOSE X1

      X1MAX=2*(-QQ**3)**(ONE/6)*COS(ACOS(RR/SQRT(-QQ**3))/3)

     &     -(3+2*AW-4*X(3))/3

      X1MIN=1-X(3)+(AW*X(3))/(1-X(3))

      IF (X1MAX.GE.1.OR.X1MIN.GE.1.OR.X1MAX.LE.X1MIN) GOTO 100

      X(1)=1-(1-X1MAX)*((1-X1MIN)/(1-X1MAX))**HWR()

C---CALCULATE WEIGHT

      W=((1+1/AW-2*AW)*((1-AW)*X(3)-(1-X(1))*(1-X(3))-X(3)**2)

     &     +(1+1/(2*AW))*X(3)*(X(1)+X(3)-1)**2+2*X(3)**2*(1-X(1)))

     &     *(1/X3MIN-1/X3MAX)*LOG((1-X1MIN)/(1-X1MAX))

C---QSCALE=DURHAM-LIKE TRANSVERSE MOMENTUM OF THE GLUON

      QSCALE=EM*HWUSQR(X(3)*(1-X(1))/(2-X(1)-X(3)-AW))

C---FACTOR FOR GLUON EMISSION

      ID=IDHW(JDAHEP(2,ICMF))

      GLUFAC=0

      IF (QSCALE.GT.HWBVMC(13)) GLUFAC=CFFAC*HWUALF(1,QSCALE)

     &     /(PIFAC*(1-AW)*(1-2*AW+1/AW))

C---IN FRACTION GLUFAC*W OF EVENTS ADD A GLUON

      IF (GLUFAC*W.GT.HWR()) THEN

        ID3=13

      ELSE

        RETURN

      ENDIF

C---CHECK INFRA-RED CUT-OFF FOR GLUON

      M(1)=PHEP(5,JDAHEP(1,ICMF))

      M(2)=HWBVMC(ID)

      M(3)=HWBVMC(ID3)

      E(1)=HALF*EM*(X(1)+AW+(-M(2)**2-M(3)**2)/EM**2)

      E(3)=HALF*EM*X(3)

      E(2)=EM-E(1)-E(3)

      PTSQ=-LAMBDA(E(1)**2-M(1)**2,E(3)**2-M(3)**2,

     &     E(2)**2-M(2)**2)

      IF (PTSQ.LE.0.OR.E(1).LE.M(1).OR.E(2).LE.M(2).OR.E(3).LE.M(3))

     $     RETURN

C---CALCULATE MASS-DEPENDENT SUPPRESSION

      EPS=(RMASS(ID)/EM)**2

      EPG=(RMASS(ID3)/EM)**2

      GAMDEP=(1-AW)*(1+1/AW-2*AW)/(SQRT(1+AW**2+EPS**2

     &     -2*AW-2*EPS-2*AW*EPS)*(1+EPS+(1-EPS)**2/AW-2*AW))

      MASDEP=GAMDEP/(1-X(1))*((1+EPS+(1-EPS)**2/AW-2*AW)

     &     *((1-AW+EPS)*X(3)*(1-X(1))-(1-X(1))**2*(1-X(3))

     &     -X(3)**2*(1-X(1)+EPS))+(1+(1+EPS)/(2*AW))*X(3)

     &     *(1-X(1))*(X(1)+X(3)-1)**2+2*X(3)**2*(1-X(1))**2)

      IF (MASDEP.LT.HWR()*((1+1/AW-2*AW)*((1-AW)*X(3)

     &     -(1-X(1))*(1-X(3))-X(3)**2)+(1+1/(2*AW))*X(3)

     &     *(X(1)+X(3)-1)**2+2*X(3)**2*(1-X(1)))) RETURN

C---STORE OLD MOMENTA

c---PT = TOP MOMENTUM, PW= W MOMENTUM

      CALL HWVEQU(5,PHEP(1,ICMF),PT)

      CALL HWVEQU(5,PHEP(1,JDAHEP(1,ICMF)),PW)

C--------GET THE NON-EMITTING PARTON CMF DIRECTION

      CALL HWULOF(PHEP(1,ICMF),PW,PW)

      CALL HWRAZM(ONE,CS,SN)

      CALL HWUROT(PW,CS,SN,R)

      CALL HWUROF(R,PW,PW)

      CALL HWUMAS(PW)

C---REORDER ENTRIES: IHEP=EMITTER,  KHEP=EMITTED

      NHEP=NHEP+1

      IHEP=JDAHEP(2,ICMF)

      WHEP=JDAHEP(1,ICMF)

      KHEP=NHEP

C---SET UP MOMENTA IN TOP REST FRAME

      PHEP(1,ICMF)=0

      PHEP(2,ICMF)=0

      PHEP(3,ICMF)=0

      PHEP(4,ICMF)=EM

      PHEP(5,ICMF)=EM

      PHEP(4,IHEP)=HALF*EM*(2-X(1)-X(3)+EPS-AW+EPG)

      PHEP(4,KHEP)=HALF*EM*X(3)

      PHEP(5,IHEP)=RMASS(ID)

      PHEP(5,KHEP)=RMASS(ID3)

      PHEP(3,KHEP)=HALF*EM*((X(1)+AW-EPS-EPG)*X(3)-2*(1+EPS-AW

     $     -EPG-(2+EPS+EPG-AW-X(1)-X(3))))/HWUSQR((X(1)+AW

     $     -EPS-EPG)**2-4*AW)

      PHEP(3,IHEP)=-PHEP(3,KHEP)-HALF*EM

     $     *HWUSQR((X(1)+AW-EPS-EPG)**2-4*AW)

      PHEP(2,IHEP)=0

      PHEP(1,KHEP)=HWUSQR(PHEP(4,KHEP)**2-PHEP(5,KHEP)**2

     $     -PHEP(3,KHEP)**2)

      PHEP(1,IHEP)=-PHEP(1,KHEP)

      PHEP(2,KHEP)=0

      CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,KHEP),PW1)

      CALL HWVDIF(4,PHEP(1,ICMF),PW1,PW1)

      CALL HWUMAS(PW1)

      DO K=1,5

        PHEP(K,WHEP)=PW1(K)

      ENDDO

C---ORIENT IN CMF, THEN BOOST TO LAB

      CALL HWUROB(R,PHEP(1,ICMF),PHEP(1,ICMF))

      CALL HWUROB(R,PHEP(1,IHEP),PHEP(1,IHEP))

      CALL HWUROB(R,PHEP(1,WHEP),PHEP(1,WHEP))

      CALL HWUROB(R,PHEP(1,KHEP),PHEP(1,KHEP))

      CALL HWULOB(PT,PHEP(1,IHEP),PHEP(1,IHEP))

      CALL HWULOB(PT,PHEP(1,KHEP),PHEP(1,KHEP))

      CALL HWULOB(PT,PHEP(1,ICMF),PHEP(1,ICMF))

      CALL HWULOB(PT,PHEP(1,WHEP),PHEP(1,WHEP))

C---STATUS AND COLOUR CONNECTION

      ISTHEP(KHEP)=114

      IDHW(KHEP)=ID3

      IDHEP(KHEP)=IDPDG(ID3)

      JDAHEP(2,ICMF)=KHEP

      JMOHEP(1,KHEP)=ICMF

      JMOHEP(1,IHEP)=ICMF

      JDAHEP(1,KHEP)=0

      JMOHEP(2,IHEP)=ICMF

      JDAHEP(2,IHEP)=KHEP

      JMOHEP(2,KHEP)=IHEP

      JDAHEP(2,KHEP)=ICMF

 999  END

CDECK  ID>, HWBVMC.

*CMZ :-        -26/04/91  11.11.54  by  Bryan Webber

*-- Author :    Bryan Webber

C-----------------------------------------------------------------------

      FUNCTION HWBVMC(ID)

C-----------------------------------------------------------------------

C     VIRTUAL MASS CUTOFF FOR PARTON TYPE ID

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      DOUBLE PRECISION HWBVMC

      INTEGER ID

      IF (ID.EQ.13) THEN

        HWBVMC=RMASS(ID)+VGCUT

      ELSEIF (ID.LT.13) THEN

        HWBVMC=RMASS(ID)+VQCUT

      ELSEIF (ID.EQ.59) THEN

        HWBVMC=RMASS(ID)+VPCUT

      ELSE

        HWBVMC=RMASS(ID)

      ENDIF

      END