C. mymat.f
C. Time-stamp: <2002-07-18 15:38:06 takeuchi>
C. special matrix element for t-W-b decay

C. I1=I2=I3=I4=0 => IR: mode ( 0: original 1: custom)
C. when you select custom: WT=> pm   WTMAX=> pz
      Subroutine mymat(I1, I2, I3, I4, IR, WT, WTMAX)

      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
C
C...PYTHIA Commonblocks.

      integer mode
      data mode/-1/

      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      SAVE /PYJETS/, mode, pm, pz, ymax

      real*8 wt, wtmax
      real meb2, mw2, mt2, costh
C.
      real pm, pz, ymax, y
      real pp(4)

      y(x)=pm*(1.0-x)**2+2.0*pz*(1.0-x**2)

      if (i1.eq.0 .and. i2.eq.0 .and. i3.eq.0 .and. i4.eq.0 ) then
         if (ir.eq.2) then
            mode=2
            print *, 'mymat: Flat ME for t-W-b selected'
         else if (ir.eq.1) then
            mode=1
            pm=wt
            pz=wtmax
            if (pm.le.pz) then
               ymax=4.0*pz**2/(2.0*pz-pm)
            else
               ymax=4.0*pm
            endif
            print *, 'mymat: Custom ME for t-W-b selected'
            print *, '(P-, P0)=', pm, pz
            print *, 'ymax=', ymax
         else
            mode=0
            print *, 'mymat: Original ME for t-W-b selected'
         endif
      endif

      if (mode.eq.-1) then
         mode=0
         print *, 'mymat: Original ME for t-W-b selected as default'
      endif

*** main ***
      if (mode.eq.1) then

C.    meb
         do i=1,4
            pp(i)=P(I4,i)+P(I2,i)
         enddo
         meb2=pp(4)**2-pp(1)**2-pp(2)**2-pp(3)**2

C.    mw
         do i=1,4
            pp(i)=P(I2,i)+P(I3,i)
         enddo
         mw2=pp(4)**2-pp(1)**2-pp(2)**2-pp(3)**2

C.    mbenu
         do i=1,4
            pp(i)=P(I2,i)+P(I3,i)+P(I4,i)
         enddo
         mt2=pp(4)**2-pp(1)**2-pp(2)**2-pp(3)**2
         
         costh=2.0*meb2/(mt2-mw2)-1.0
         
         wt=y(costh)
         wtmax=ymax
***         print *, 'custom'
      else if (mode.eq.0) then
***         print *, 'orig'
         WT=(P(I1,4)*P(I2,4)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
     $        P(I1,3)*P(I2,3))*(P(I3,4)*P(I4,4)-P(I3,1)*P(I4,1)-
     $        P(I3,2)*P(I4,2)-P(I3,3)*P(I4,3))
         WTMAX=(P(I1,5)**4-P(IR,5)**4)/8D0
      else if (mode.eq.2) then
         WT=1D0
         WTMAX=1D0
      endif

      END