pyhiresd.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhiresd.f

    
C*********************************************************************  
    
      SUBROUTINE pyhiresd 
    
C...Allows resonances to decay (including parton showers for hadronic   
C...channels).  
      IMPLICIT DOUBLE PRECISION(d)  
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2) 
      SAVE /pyhiint2/ 
      common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3) 
      SAVE /pyhiint4/ 
      dimension iref(10,6),kdcy(2),kfl1(2),kfl2(2),nsd(2),ilin(6),  
     &coup(6,4),pk(6,4),pkk(6,6),cthe(2),phi(2),wdtp(0:40), 
     &wdte(0:40,0:5)    
      COMPLEX fgk,ha(6,6),hc(6,6)   
    
C...The F, Xi and Xj functions of Gunion and Kunszt 
C...(Phys. Rev. D33, 665, plus errata from the authors).    
      fgk(i1,i2,i3,i4,i5,i6)=4.*ha(i1,i3)*hc(i2,i6)*(ha(i1,i5)* 
     &hc(i1,i4)+ha(i3,i5)*hc(i3,i4))    
      digk(dt,du)=-4.*d34*d56+dt*(3.*dt+4.*du)+dt**2*(dt*du/(d34*d56)-  
     &2.*(1./d34+1./d56)*(dt+du)+2.*(d34/d56+d56/d34))  
      djgk(dt,du)=8.*(d34+d56)**2-8.*(d34+d56)*(dt+du)-6.*dt*du-    
     &2.*dt*du*(dt*du/(d34*d56)-2.*(1./d34+1./d56)*(dt+du)+ 
     &2.*(d34/d56+d56/d34)) 
    
C...Define initial two objects, initialize loop.    
      isub=mint(1)  
      sh=vint(44)   
      iref(1,5)=0   
      iref(1,6)=0   
      IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN   
        iref(1,1)=mint(84)+2+iset(isub) 
        iref(1,2)=0 
        iref(1,3)=mint(83)+6+iset(isub) 
        iref(1,4)=0 
      ELSEIF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN   
        iref(1,1)=mint(84)+1+iset(isub) 
        iref(1,2)=mint(84)+2+iset(isub) 
        iref(1,3)=mint(83)+5+iset(isub) 
        iref(1,4)=mint(83)+6+iset(isub) 
      ENDIF 
      np=1  
      ip=0  
  100 ip=ip+1   
      ninh=0    
    
C...Loop over one/two resonances; reset decay rates.    
      jtmax=2   
      IF(ip.EQ.1.AND.(iset(isub).EQ.1.OR.iset(isub).EQ.3)) jtmax=1  
      DO 140 jt=1,jtmax 
      kdcy(jt)=0    
      kfl1(jt)=0    
      kfl2(jt)=0    
      nsd(jt)=iref(ip,jt)   
      id=iref(ip,jt)    
      IF(id.EQ.0) goto 140  
      kfa=iabs(k(id,2)) 
      IF(kfa.LT.23.OR.kfa.GT.40) goto 140   
      IF(mdcy(kfa,1).NE.0) THEN 
        IF(isub.EQ.1.OR.isub.EQ.141) mint(61)=1 
        CALL pyhiwidt(kfa,p(id,5),wdtp,wdte)  
        IF(kchg(kfa,3).EQ.0) THEN   
          ipm=2 
        ELSE    
          ipm=(5+isign(1,k(id,2)))/2    
        ENDIF   
        IF(jtmax.EQ.1.OR.iabs(k(iref(ip,1),2)).NE.iabs(k(iref(ip,2),2)))    
     &  THEN    
          i12=4 
        ELSE    
          IF(jt.EQ.1) i12=int(4.5+rlu(0))   
          i12=9-i12 
        ENDIF   
        rkfl=(wdte(0,1)+wdte(0,ipm)+wdte(0,i12))*rlu(0) 
        DO 120 i=1,mdcy(kfa,3)  
        idc=i+mdcy(kfa,2)-1 
        kfl1(jt)=kfdp(idc,1)*isign(1,k(id,2))   
        kfl2(jt)=kfdp(idc,2)*isign(1,k(id,2))   
        rkfl=rkfl-(wdte(i,1)+wdte(i,ipm)+wdte(i,i12))   
        IF(rkfl.LE.0.) goto 130 
  120   CONTINUE    
  130   CONTINUE    
      ENDIF 
    
C...Summarize result on decay channel chosen.   
      IF((kfa.EQ.23.OR.kfa.EQ.24).AND.kfl1(jt).EQ.0) ninh=ninh+1    
      IF(kfl1(jt).EQ.0) goto 140    
      kdcy(jt)=2    
      IF(iabs(kfl1(jt)).LE.10.OR.kfl1(jt).EQ.21) kdcy(jt)=1 
      IF((iabs(kfl1(jt)).GE.23.AND.iabs(kfl1(jt)).LE.25).OR.    
     &(iabs(kfl1(jt)).EQ.37)) kdcy(jt)=3    
      nsd(jt)=n 
    
C...Fill decay products, prepared for parton showers for quarks.    
      IF(kdcy(jt).EQ.1) THEN    
        CALL lu2ent(-(n+1),kfl1(jt),kfl2(jt),p(id,5))   
      ELSE  
        CALL lu2ent(n+1,kfl1(jt),kfl2(jt),p(id,5))  
      ENDIF 
      IF(jtmax.EQ.1) THEN   
        cthe(jt)=vint(13)+(vint(33)-vint(13)+vint(34)-vint(14))*rlu(0)  
        IF(cthe(jt).GT.vint(33)) cthe(jt)=cthe(jt)+vint(14)-vint(33)    
        phi(jt)=vint(24)    
      ELSE  
        cthe(jt)=2.*rlu(0)-1.   
        phi(jt)=paru(2)*rlu(0)  
      ENDIF 
  140 CONTINUE  
      IF(mint(3).EQ.1.AND.ip.EQ.1) THEN 
        mint(25)=kfl1(1)    
        mint(26)=kfl2(1)    
      ENDIF 
      IF(jtmax.EQ.1.AND.kdcy(1).EQ.0) goto 530  
      IF(jtmax.EQ.2.AND.kdcy(1).EQ.0.AND.kdcy(2).EQ.0) goto 530 
      IF(mstp(45).LE.0.OR.iref(ip,2).EQ.0.OR.ninh.GE.1) goto 500    
      IF(k(iref(1,1),2).EQ.25.AND.ip.EQ.1) goto 500 
      IF(k(iref(1,1),2).EQ.25.AND.kdcy(1)*kdcy(2).EQ.0) goto 500    
    
C...Order incoming partons and outgoing resonances. 
      ilin(1)=mint(84)+1    
      IF(k(mint(84)+1,2).GT.0) ilin(1)=mint(84)+2   
      IF(k(ilin(1),2).EQ.21) ilin(1)=2*mint(84)+3-ilin(1)   
      ilin(2)=2*mint(84)+3-ilin(1)  
      imin=1    
      IF(iref(ip,5).EQ.25) imin=3   
      imax=2    
      iord=1    
      IF(k(iref(ip,1),2).EQ.23) iord=2  
      IF(k(iref(ip,1),2).EQ.24.AND.k(iref(ip,2),2).EQ.-24) iord=2   
      IF(iabs(k(iref(ip,iord),2)).EQ.25) iord=3-iord    
      IF(kdcy(iord).EQ.0) iord=3-iord   
    
C...Order decay products of resonances. 
      DO 390 jt=iord,3-iord,3-2*iord    
      IF(kdcy(jt).EQ.0) THEN    
        ilin(imax+1)=nsd(jt)    
        imax=imax+1 
      ELSEIF(k(nsd(jt)+1,2).GT.0) THEN  
        ilin(imax+1)=n+2*jt-1   
        ilin(imax+2)=n+2*jt 
        imax=imax+2 
        k(n+2*jt-1,2)=k(nsd(jt)+1,2)    
        k(n+2*jt,2)=k(nsd(jt)+2,2)  
      ELSE  
        ilin(imax+1)=n+2*jt 
        ilin(imax+2)=n+2*jt-1   
        imax=imax+2 
        k(n+2*jt-1,2)=k(nsd(jt)+1,2)    
        k(n+2*jt,2)=k(nsd(jt)+2,2)  
      ENDIF 
  390 CONTINUE  
    
C...Find charge, isospin, left- and righthanded couplings.  
      xw=paru(102)  
      DO 410 i=imin,imax    
      DO 400 j=1,4  
  400 coup(i,j)=0.  
      kfa=iabs(k(ilin(i),2))    
      IF(kfa.GT.20) goto 410    
      coup(i,1)=luchge(kfa)/3.  
      coup(i,2)=(-1)**mod(kfa,2)    
      coup(i,4)=-2.*coup(i,1)*xw    
      coup(i,3)=coup(i,2)+coup(i,4) 
  410 CONTINUE  
      sqmz=pmas(23,1)**2    
      gzmz=pmas(23,1)*pmas(23,2)    
      sqmw=pmas(24,1)**2    
      gzmw=pmas(24,1)*pmas(24,2)    
      sqmzp=pmas(32,1)**2   
      gzmzp=pmas(32,1)*pmas(32,2)   
    
C...Select random angles; construct massless four-vectors.  
  420 DO 430 i=n+1,n+4  
      k(i,1)=1  
      DO 430 j=1,5  
  430 p(i,j)=0. 
      DO 440 jt=1,jtmax 
      IF(kdcy(jt).EQ.0) goto 440    
      id=iref(ip,jt)    
      p(n+2*jt-1,3)=0.5*p(id,5) 
      p(n+2*jt-1,4)=0.5*p(id,5) 
      p(n+2*jt,3)=-0.5*p(id,5)  
      p(n+2*jt,4)=0.5*p(id,5)   
      cthe(jt)=2.*rlu(0)-1. 
      phi(jt)=paru(2)*rlu(0)    
      CALL ludbrb(n+2*jt-1,n+2*jt,acos(cthe(jt)),phi(jt),   
     &dble(p(id,1)/p(id,4)),dble(p(id,2)/p(id,4)),dble(p(id,3)/p(id,4)))    
  440 CONTINUE  
    
C...Store incoming and outgoing momenta, with random rotation to    
C...avoid accidental zeroes in HA expressions.  
      DO 450 i=1,imax   
      k(n+4+i,1)=1  
      p(n+4+i,4)=sqrt(p(ilin(i),1)**2+p(ilin(i),2)**2+p(ilin(i),3)**2+  
     &p(ilin(i),5)**2)  
      p(n+4+i,5)=p(ilin(i),5)   
      DO 450 j=1,3  
  450 p(n+4+i,j)=p(ilin(i),j)   
      therr=acos(2.*rlu(0)-1.)  
      phirr=paru(2)*rlu(0)  
      CALL ludbrb(n+5,n+4+imax,therr,phirr,0d0,0d0,0d0) 
      DO 460 i=1,imax   
      DO 460 j=1,4  
  460 pk(i,j)=p(n+4+i,j)    
    
C...Calculate internal products.    
      IF(isub.EQ.22.OR.isub.EQ.23.OR.isub.EQ.25) THEN   
        DO 470 i1=imin,imax-1   
        DO 470 i2=i1+1,imax 
        ha(i1,i2)=sqrt((pk(i1,4)-pk(i1,3))*(pk(i2,4)+pk(i2,3))/ 
     &  (1e-20+pk(i1,1)**2+pk(i1,2)**2))*cmplx(pk(i1,1),pk(i1,2))-  
     &  sqrt((pk(i1,4)+pk(i1,3))*(pk(i2,4)-pk(i2,3))/   
     &  (1e-20+pk(i2,1)**2+pk(i2,2)**2))*cmplx(pk(i2,1),pk(i2,2))   
        hc(i1,i2)=conjg(ha(i1,i2))  
        IF(i1.LE.2) ha(i1,i2)=cmplx(0.,1.)*ha(i1,i2)    
        IF(i1.LE.2) hc(i1,i2)=cmplx(0.,1.)*hc(i1,i2)    
        ha(i2,i1)=-ha(i1,i2)    
  470   hc(i2,i1)=-hc(i1,i2)    
      ENDIF 
      DO 480 i=1,2  
      DO 480 j=1,4  
  480 pk(i,j)=-pk(i,j)  
      DO 490 i1=imin,imax-1 
      DO 490 i2=i1+1,imax   
      pkk(i1,i2)=2.*(pk(i1,4)*pk(i2,4)-pk(i1,1)*pk(i2,1)-   
     &pk(i1,2)*pk(i2,2)-pk(i1,3)*pk(i2,3))  
  490 pkk(i2,i1)=pkk(i1,i2) 
    
      IF(iref(ip,5).EQ.25) THEN 
C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons 
        wt=16.*pkk(3,5)*pkk(4,6)    
        IF(ip.EQ.1) wtmax=sh**2 
        IF(ip.GE.2) wtmax=p(iref(ip,6),5)**4    
    
      ELSEIF(isub.EQ.1) THEN    
        IF(kfa.NE.37) THEN  
C...Angular weight for gamma*/Z0 -> 2 quarks/leptons    
          ei=kchg(iabs(mint(15)),1)/3.  
          ai=sign(1.,ei+0.1)    
          vi=ai-4.*ei*xw    
          ef=kchg(kfa,1)/3. 
          af=sign(1.,ef+0.1)    
          vf=af-4.*ef*xw    
          gg=1. 
          gz=1./(8.*xw*(1.-xw))*sh*(sh-sqmz)/((sh-sqmz)**2+gzmz**2) 
          zz=1./(16.*xw*(1.-xw))**2*sh**2/((sh-sqmz)**2+gzmz**2)    
          IF(mstp(43).EQ.1) THEN    
C...Only gamma* production included 
            gz=0.   
            zz=0.   
          ELSEIF(mstp(43).EQ.2) THEN    
C...Only Z0 production included 
            gg=0.   
            gz=0.   
          ENDIF 
          asym=2.*(ei*ai*gz*ef*af+4.*vi*ai*zz*vf*af)/(ei**2*gg*ef**2+   
     &    ei*vi*gz*ef*vf+(vi**2+ai**2)*zz*(vf**2+af**2))    
          wt=1.+asym*cthe(jt)+cthe(jt)**2   
          wtmax=2.+abs(asym)    
        ELSE    
C...Angular weight for gamma*/Z0 -> H+ + H- 
          wt=1.-cthe(jt)**2 
          wtmax=1.  
        ENDIF   
    
      ELSEIF(isub.EQ.2) THEN    
C...Angular weight for W+/- -> 2 quarks/leptons 
        wt=(1.+cthe(jt))**2 
        wtmax=4.    
    
      ELSEIF(isub.EQ.15.OR.isub.EQ.19) THEN 
C...Angular weight for f + fb -> gluon/gamma + Z0 ->    
C...-> gluon/gamma + 2 quarks/leptons   
        wt=((coup(1,3)*coup(3,3))**2+(coup(1,4)*coup(3,4))**2)* 
     &  (pkk(1,3)**2+pkk(2,4)**2)+((coup(1,3)*coup(3,4))**2+    
     &  (coup(1,4)*coup(3,3))**2)*(pkk(1,4)**2+pkk(2,3)**2) 
        wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*  
     &  ((pkk(1,3)+pkk(1,4))**2+(pkk(2,3)+pkk(2,4))**2) 
    
      ELSEIF(isub.EQ.16.OR.isub.EQ.20) THEN 
C...Angular weight for f + fb' -> gluon/gamma + W+/- -> 
C...-> gluon/gamma + 2 quarks/leptons   
        wt=pkk(1,3)**2+pkk(2,4)**2  
        wtmax=(pkk(1,3)+pkk(1,4))**2+(pkk(2,3)+pkk(2,4))**2 
    
      ELSEIF(isub.EQ.22) THEN   
C...Angular weight for f + fb -> Z0 + Z0 -> 4 quarks/leptons    
        s34=p(iref(ip,iord),5)**2   
        s56=p(iref(ip,3-iord),5)**2 
        ti=pkk(1,3)+pkk(1,4)+s34    
        ui=pkk(1,5)+pkk(1,6)+s56    
        wt=coup(1,3)**4*((coup(3,3)*coup(5,3)*abs(fgk(1,2,3,4,5,6)/ 
     &  ti+fgk(1,2,5,6,3,4)/ui))**2+(coup(3,4)*coup(5,3)*abs(   
     &  fgk(1,2,4,3,5,6)/ti+fgk(1,2,5,6,4,3)/ui))**2+(coup(3,3)*    
     &  coup(5,4)*abs(fgk(1,2,3,4,6,5)/ti+fgk(1,2,6,5,3,4)/ui))**2+ 
     &  (coup(3,4)*coup(5,4)*abs(fgk(1,2,4,3,6,5)/ti+fgk(1,2,6,5,4,3)/  
     &  ui))**2)+coup(1,4)**4*((coup(3,3)*coup(5,3)*abs(    
     &  fgk(2,1,5,6,3,4)/ti+fgk(2,1,3,4,5,6)/ui))**2+(coup(3,4)*    
     &  coup(5,3)*abs(fgk(2,1,6,5,3,4)/ti+fgk(2,1,3,4,6,5)/ui))**2+ 
     &  (coup(3,3)*coup(5,4)*abs(fgk(2,1,5,6,4,3)/ti+fgk(2,1,4,3,5,6)/  
     &  ui))**2+(coup(3,4)*coup(5,4)*abs(fgk(2,1,6,5,4,3)/ti+   
     &  fgk(2,1,4,3,6,5)/ui))**2)   
        wtmax=4.*s34*s56*(coup(1,3)**4+coup(1,4)**4)*(coup(3,3)**2+ 
     &  coup(3,4)**2)*(coup(5,3)**2+coup(5,4)**2)*4.*(ti/ui+ui/ti+  
     &  2.*sh*(s34+s56)/(ti*ui)-s34*s56*(1./ti**2+1./ui**2))    
    
      ELSEIF(isub.EQ.23) THEN   
C...Angular weight for f + fb' -> Z0 + W +/- -> 4 quarks/leptons    
        d34=p(iref(ip,iord),5)**2   
        d56=p(iref(ip,3-iord),5)**2 
        dt=pkk(1,3)+pkk(1,4)+d34    
        du=pkk(1,5)+pkk(1,6)+d56    
        cawz=coup(2,3)/sngl(dt)-2.*(1.-xw)*coup(1,2)/(sh-sqmw)  
        cbwz=coup(1,3)/sngl(du)+2.*(1.-xw)*coup(1,2)/(sh-sqmw)  
        wt=coup(5,3)**2*abs(cawz*fgk(1,2,3,4,5,6)+cbwz* 
     &  fgk(1,2,5,6,3,4))**2+coup(5,4)**2*abs(cawz* 
     &  fgk(1,2,3,4,6,5)+cbwz*fgk(1,2,6,5,3,4))**2  
        wtmax=4.*d34*d56*(coup(5,3)**2+coup(5,4)**2)*(cawz**2*  
     &  digk(dt,du)+cbwz**2*digk(du,dt)+cawz*cbwz*djgk(dt,du))  
    
      ELSEIF(isub.EQ.24) THEN   
C...Angular weight for f + fb -> Z0 + H0 -> 2 quarks/leptons + H0   
        wt=((coup(1,3)*coup(3,3))**2+(coup(1,4)*coup(3,4))**2)* 
     &  pkk(1,3)*pkk(2,4)+((coup(1,3)*coup(3,4))**2+(coup(1,4)* 
     &  coup(3,3))**2)*pkk(1,4)*pkk(2,3)    
        wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*  
     &  (pkk(1,3)+pkk(1,4))*(pkk(2,3)+pkk(2,4)) 
    
      ELSEIF(isub.EQ.25) THEN   
C...Angular weight for f + fb -> W+ + W- -> 4 quarks/leptons    
        d34=p(iref(ip,iord),5)**2   
        d56=p(iref(ip,3-iord),5)**2 
        dt=pkk(1,3)+pkk(1,4)+d34    
        du=pkk(1,5)+pkk(1,6)+d56    
        cdww=(coup(1,3)*sqmz/(sh-sqmz)+coup(1,2))/sh    
        caww=cdww+0.5*(coup(1,2)+1.)/sngl(dt)   
        cbww=cdww+0.5*(coup(1,2)-1.)/sngl(du)   
        ccww=coup(1,4)*sqmz/(sh-sqmz)/sh    
        wt=abs(caww*fgk(1,2,3,4,5,6)-cbww*fgk(1,2,5,6,3,4))**2+ 
     &  ccww**2*abs(fgk(2,1,5,6,3,4)-fgk(2,1,3,4,5,6))**2   
        wtmax=4.*d34*d56*(caww**2*digk(dt,du)+cbww**2*digk(du,dt)-caww* 
     &  cbww*djgk(dt,du)+ccww**2*(digk(dt,du)+digk(du,dt)-djgk(dt,du))) 
    
      ELSEIF(isub.EQ.26) THEN   
C...Angular weight for f + fb' -> W+/- + H0 -> 2 quarks/leptons + H0    
        wt=pkk(1,3)*pkk(2,4)    
        wtmax=(pkk(1,3)+pkk(1,4))*(pkk(2,3)+pkk(2,4))   
    
      ELSEIF(isub.EQ.30) THEN   
C...Angular weight for f + g -> f + Z0 -> f + 2 quarks/leptons  
        IF(k(ilin(1),2).GT.0) wt=((coup(1,3)*coup(3,3))**2+ 
     &  (coup(1,4)*coup(3,4))**2)*(pkk(1,4)**2+pkk(3,5)**2)+    
     &  ((coup(1,3)*coup(3,4))**2+(coup(1,4)*coup(3,3))**2)*    
     &  (pkk(1,3)**2+pkk(4,5)**2)   
        IF(k(ilin(1),2).LT.0) wt=((coup(1,3)*coup(3,3))**2+ 
     &  (coup(1,4)*coup(3,4))**2)*(pkk(1,3)**2+pkk(4,5)**2)+    
     &  ((coup(1,3)*coup(3,4))**2+(coup(1,4)*coup(3,3))**2)*    
     &  (pkk(1,4)**2+pkk(3,5)**2)   
        wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*  
     &  ((pkk(1,3)+pkk(1,4))**2+(pkk(3,5)+pkk(4,5))**2) 
    
      ELSEIF(isub.EQ.31) THEN   
C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons  
        IF(k(ilin(1),2).GT.0) wt=pkk(1,4)**2+pkk(3,5)**2    
        IF(k(ilin(1),2).LT.0) wt=pkk(1,3)**2+pkk(4,5)**2    
        wtmax=(pkk(1,3)+pkk(1,4))**2+(pkk(3,5)+pkk(4,5))**2 
    
      ELSEIF(isub.EQ.141) THEN  
C...Angular weight for gamma*/Z0/Z'0 -> 2 quarks/leptons    
        ei=kchg(iabs(mint(15)),1)/3.    
        ai=sign(1.,ei+0.1)  
        vi=ai-4.*ei*xw  
        api=sign(1.,ei+0.1) 
        vpi=api-4.*ei*xw    
        ef=kchg(kfa,1)/3.   
        af=sign(1.,ef+0.1)  
        vf=af-4.*ef*xw  
        apf=sign(1.,ef+0.1) 
        vpf=apf-4.*ef*xw    
        gg=1.   
        gz=1./(8.*xw*(1.-xw))*sh*(sh-sqmz)/((sh-sqmz)**2+gzmz**2)   
        gzp=1./(8.*xw*(1.-xw))*sh*(sh-sqmzp)/((sh-sqmzp)**2+gzmzp**2)   
        zz=1./(16.*xw*(1.-xw))**2*sh**2/((sh-sqmz)**2+gzmz**2)  
        zzp=2./(16.*xw*(1.-xw))**2* 
     &  sh**2*((sh-sqmz)*(sh-sqmzp)+gzmz*gzmzp)/    
     &  (((sh-sqmz)**2+gzmz**2)*((sh-sqmzp)**2+gzmzp**2))   
        zpzp=1./(16.*xw*(1.-xw))**2*sh**2/((sh-sqmzp)**2+gzmzp**2)  
        IF(mstp(44).EQ.1) THEN  
C...Only gamma* production included 
          gz=0. 
          gzp=0.    
          zz=0. 
          zzp=0.    
          zpzp=0.   
        ELSEIF(mstp(44).EQ.2) THEN  
C...Only Z0 production included 
          gg=0. 
          gz=0. 
          gzp=0.    
          zzp=0.    
          zpzp=0.   
        ELSEIF(mstp(44).EQ.3) THEN  
C...Only Z'0 production included    
          gg=0. 
          gz=0. 
          gzp=0.    
          zz=0. 
          zzp=0.    
        ELSEIF(mstp(44).EQ.4) THEN  
C...Only gamma*/Z0 production included  
          gzp=0.    
          zzp=0.    
          zpzp=0.   
        ELSEIF(mstp(44).EQ.5) THEN  
C...Only gamma*/Z'0 production included 
          gz=0. 
          zz=0. 
          zzp=0.    
        ELSEIF(mstp(44).EQ.6) THEN  
C...Only Z0/Z'0 production included 
          gg=0. 
          gz=0. 
          gzp=0.    
        ENDIF   
        asym=2.*(ei*ai*gz*ef*af+ei*api*gzp*ef*apf+4.*vi*ai*zz*vf*af+    
     &  (vi*api+vpi*ai)*zzp*(vf*apf+vpf*af)+4.*vpi*api*zpzp*vpf*apf)/   
     &  (ei**2*gg*ef**2+ei*vi*gz*ef*vf+ei*vpi*gzp*ef*vpf+   
     &  (vi**2+ai**2)*zz*(vf**2+af**2)+(vi*vpi+ai*api)*zzp* 
     &  (vf*vpf+af*apf)+(vpi**2+api**2)*zpzp*(vpf**2+apf**2))   
        wt=1.+asym*cthe(jt)+cthe(jt)**2 
        wtmax=2.+abs(asym)  
    
      ELSE  
        wt=1.   
        wtmax=1.    
      ENDIF 
C...Obtain correct angular distribution by rejection techniques.    
      IF(wt.LT.rlu(0)*wtmax) goto 420   
    
C...Construct massive four-vectors using angles chosen. Mark decayed    
C...resonances, add documentation lines. Shower evolution.  
  500 DO 520 jt=1,jtmax 
      IF(kdcy(jt).EQ.0) goto 520    
      id=iref(ip,jt)    
      CALL ludbrb(nsd(jt)+1,nsd(jt)+2,acos(cthe(jt)),phi(jt),   
     &dble(p(id,1)/p(id,4)),dble(p(id,2)/p(id,4)),dble(p(id,3)/p(id,4)))    
      k(id,1)=k(id,1)+10    
      k(id,4)=nsd(jt)+1 
      k(id,5)=nsd(jt)+2 
      idoc=mint(83)+mint(4) 
      DO 510 i=nsd(jt)+1,nsd(jt)+2  
      mint(4)=mint(4)+1 
      i1=mint(83)+mint(4)   
      k(i,3)=i1 
      k(i1,1)=21    
      k(i1,2)=k(i,2)    
      k(i1,3)=iref(ip,jt+2) 
      DO 510 j=1,5  
  510 p(i1,j)=p(i,j)    
      IF(jtmax.EQ.1) THEN   
        mint(7)=mint(83)+6+2*iset(isub) 
        mint(8)=mint(83)+7+2*iset(isub) 
      ENDIF 
      IF(mstp(71).GE.1.AND.kdcy(jt).EQ.1) CALL lushow(nsd(jt)+1,    
     &nsd(jt)+2,p(id,5))    
    
C...Check if new resonances were produced, loop back if needed. 
      IF(kdcy(jt).NE.3) goto 520    
      np=np+1   
      iref(np,1)=nsd(jt)+1  
      iref(np,2)=nsd(jt)+2  
      iref(np,3)=idoc+1 
      iref(np,4)=idoc+2 
      iref(np,5)=k(iref(ip,jt),2)   
      iref(np,6)=iref(ip,jt)    
  520 CONTINUE  
  530 IF(ip.LT.np) goto 100 
    
      RETURN    
      END