ISAJET/code/dhelas.F

   1 #include "isajet/pilot.h"
   2 C  *********************************************************************
   3 C  ***                                                               ***
   4 C  ***               coded by H. Murayama & I. Watanabe              ***
   5 C  *** For the formalism and notations, see the following reference: ***
   6 C  ***           H. Murayama, I. Watanabe and K. Hagiwara            ***
   7 C  ***           "HELAS: HELicity Amplitude Subroutines              ***
   8 C  ***               for Feynman diagram evaluation"                 ***
   9 C  ***               KEK Report 91-11, December 1991                 ***
  10 C  ***                                                               ***
  11 C  *********************************************************************
  12 C
  13 C  Converted to double precision by W. Long and T. Seltzer for MadGraph.
  14 C
  15 C  Minor changes for portability by FEP, July 1999. The code is not ANSI
  16 C  standard, but that cannot be helped if MadGraph compatibility is to
  17 C  be maintained.
  18 C
  19 C ======================================================================
  20 C
  21       SUBROUTINE BOOSTX(P,Q , PBOOST)
  22 C
  23 C this subroutine performs the lorentz boost of a four-momentum.  the
  24 C momentum p is assumed to be given in the rest frame of q.  pboost is
  25 C the momentum p boosted to the frame in which q is given.  q must be a
  26 C timelike momentum.
  27 C
  28 C input:
  29 C       real    p(0:3)         : four-momentum p in the q rest  frame
  30 C       real    q(0:3)         : four-momentum q in the boosted frame
  31 C
  32 C output:
  33 C       real    pboost(0:3)    : four-momentum p in the boosted frame
  34 C
  35 #if defined(CERNLIB_IMPNONE)
  36       IMPLICIT NONE
  37 #endif
  38       REAL*8    P(0:3),Q(0:3),PBOOST(0:3),PQ,QQ,M,LF
  39       REAL*8 RXZERO
  40       PARAMETER( RXZERO=0.0D0 )
  41 C
  42       QQ=Q(1)**2+Q(2)**2+Q(3)**2
  43 C
  44       IF ( QQ .NE. RXZERO ) THEN
  45          PQ=P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3)
  46          M=SQRT(Q(0)**2-QQ)
  47          LF=((Q(0)-M)*PQ/QQ+P(0))/M
  48          PBOOST(0) = (P(0)*Q(0)+PQ)/M
  49          PBOOST(1) =  P(1)+Q(1)*LF
  50          PBOOST(2) =  P(2)+Q(2)*LF
  51          PBOOST(3) =  P(3)+Q(3)*LF
  52       ELSE
  53          PBOOST(0)=P(0)
  54          PBOOST(1)=P(1)
  55          PBOOST(2)=P(2)
  56          PBOOST(3)=P(3)
  57       ENDIF
  58 C
  59       RETURN
  60       END
  61 C
  62 C **********************************************************************
  63 C
  64       SUBROUTINE COUP1X(SW2 , GW,GWWA,GWWZ)
  65 C
  66 C this subroutine sets up the coupling constants of the gauge bosons in
  67 C the standard model.
  68 C
  69 C input:
  70 C       real    sw2            : square of sine of the weak angle
  71 C
  72 C output:
  73 C       real    gw             : weak coupling constant
  74 C       real    gwwa           : dimensionless coupling of w-,w+,a
  75 C       real    gwwz           : dimensionless coupling of w-,w+,z
  76 C
  77 #if defined(CERNLIB_IMPNONE)
  78       IMPLICIT NONE
  79 #endif
  80       REAL*8    SW2,GW,GWWA,GWWZ,ALPHA,FOURPI,EE,SW,CW
  81       REAL*8 RXONE, RXFOUR, RXOTE, RXPI, RIALPH
  82       PARAMETER( RXONE=1.0D0, RXFOUR=4.0D0, RXOTE=128.0D0 )
  83       PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
  84 C
  85       ALPHA = RXONE / RXOTE
  86 C      alpha = r_one / r_ialph
  87       FOURPI = RXFOUR * RXPI
  88       EE=SQRT( ALPHA * FOURPI )
  89       SW=SQRT( SW2 )
  90       CW=SQRT( RXONE - SW2 )
  91 C
  92       GW    =  EE/SW
  93       GWWA  =  EE
  94       GWWZ  =  EE*CW/SW
  95 C
  96       RETURN
  97       END
  98 C
  99 C ----------------------------------------------------------------------
 100 C
 101       SUBROUTINE COUP2X(SW2 , GAL,GAU,GAD,GWF,GZN,GZL,GZU,GZD,G1)
 102 C
 103 C this subroutine sets up the coupling constants for the fermion-
 104 C fermion-vector vertices in the standard model.  the array of the
 105 C couplings specifies the chirality of the flowing-in fermion.  g??(1)
 106 C denotes a left-handed coupling, and g??(2) a right-handed coupling.
 107 C
 108 C input:
 109 C       real    sw2            : square of sine of the weak angle
 110 C
 111 C output:
 112 C       real    gal(2)         : coupling with a of charged leptons
 113 C       real    gau(2)         : coupling with a of up-type quarks
 114 C       real    gad(2)         : coupling with a of down-type quarks
 115 C       real    gwf(2)         : coupling with w-,w+ of fermions
 116 C       real    gzn(2)         : coupling with z of neutrinos
 117 C       real    gzl(2)         : coupling with z of charged leptons
 118 C       real    gzu(2)         : coupling with z of up-type quarks
 119 C       real    gzd(2)         : coupling with z of down-type quarks
 120 C       real    g1(2)          : unit coupling of fermions
 121 C
 122 #if defined(CERNLIB_IMPNONE)
 123       IMPLICIT NONE
 124 #endif
 125       REAL*8 GAL(2),GAU(2),GAD(2),GWF(2),GZN(2),GZL(2),GZU(2),GZD(2),
 126      &     G1(2),SW2,ALPHA,FOURPI,EE,SW,CW,EZ,EY
 127 C
 128       REAL*8 RXZERO, RXHALF, RXONE, RXTWO, RTHREE, RXFOUR, RXOTE
 129       REAL*8 RXPI, RIALPH
 130       PARAMETER( RXZERO=0.0D0, RXHALF=0.5D0, RXONE=1.0D0, RXTWO=2.0D0,
 131      $     RTHREE=3.0D0 )
 132       PARAMETER( RXFOUR=4.0D0, RXOTE=128.0D0 )
 133       PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
 134 C
 135       ALPHA = RXONE / RXOTE
 136 C      alpha = r_one / r_ialph
 137       FOURPI = RXFOUR * RXPI
 138       EE=SQRT( ALPHA * FOURPI )
 139       SW=SQRT( SW2 )
 140       CW=SQRT( RXONE - SW2 )
 141       EZ=EE/(SW*CW)
 142       EY=EE*(SW/CW)
 143 C
 144       GAL(1) =  EE
 145       GAL(2) =  EE
 146       GAU(1) = -EE*RXTWO/RTHREE
 147       GAU(2) = -EE*RXTWO/RTHREE
 148       GAD(1) =  EE   /RTHREE
 149       GAD(2) =  EE   /RTHREE
 150       GWF(1) = -EE/SQRT(RXTWO*SW2)
 151       GWF(2) =  RXZERO
 152       GZN(1) = -EZ*  RXHALF
 153       GZN(2) =  RXZERO
 154       GZL(1) = -EZ*(-RXHALF+SW2)
 155       GZL(2) = -EY
 156       GZU(1) = -EZ*( RXHALF-SW2*RXTWO/RTHREE)
 157       GZU(2) =  EY*          RXTWO/RTHREE
 158       GZD(1) = -EZ*(-RXHALF+SW2   /RTHREE)
 159       GZD(2) = -EY             /RTHREE
 160       G1(1)  =  RXONE
 161       G1(2)  =  RXONE
 162 C
 163       RETURN
 164       END
 165 C
 166 C ----------------------------------------------------------------------
 167 C
 168       SUBROUTINE COUP3X(SW2,ZMASS,HMASS ,
 169      &                  GWWH,GZZH,GHHH,GWWHH,GZZHH,GHHHH)
 170 C
 171 C this subroutine sets up the coupling constants of the gauge bosons and
 172 C higgs boson in the standard model.
 173 C
 174 C input:
 175 C       real    sw2            : square of sine of the weak angle
 176 C       real    zmass          : mass of z
 177 C       real    hmass          : mass of higgs
 178 C
 179 C output:
 180 C       real    gwwh           : dimensionful  coupling of w-,w+,h
 181 C       real    gzzh           : dimensionful  coupling of z, z, h
 182 C       real    ghhh           : dimensionful  coupling of h, h, h
 183 C       real    gwwhh          : dimensionful  coupling of w-,w+,h, h
 184 C       real    gzzhh          : dimensionful  coupling of z, z, h, h
 185 C       real    ghhhh          : dimensionless coupling of h, h, h, h
 186 C
 187 #if defined(CERNLIB_IMPNONE)
 188       IMPLICIT NONE
 189 #endif
 190       REAL*8    SW2,ZMASS,HMASS,GWWH,GZZH,GHHH,GWWHH,GZZHH,GHHHH,
 191      &        ALPHA,FOURPI,EE2,SC2,V
 192 C
 193       REAL*8 RXHALF, RXONE, RXTWO, RTHREE, RXFOUR, RXOTE
 194       REAL*8 RXPI, RIALPH
 195       PARAMETER( RXHALF=0.5D0, RXONE=1.0D0, RXTWO=2.0D0, RTHREE=3.0D0 )
 196       PARAMETER( RXFOUR=4.0D0, RXOTE=128.0D0 )
 197       PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
 198 C
 199       ALPHA = RXONE / RXOTE
 200 C      alpha = r_one / r_ialph
 201       FOURPI = RXFOUR * RXPI
 202       EE2=ALPHA*FOURPI
 203       SC2=SW2*( RXONE - SW2 )
 204       V = RXTWO * ZMASS*SQRT(SC2)/SQRT(EE2)
 205 C
 206       GWWH  =   EE2/SW2*RXHALF*V
 207       GZZH  =   EE2/SC2*RXHALF*V
 208       GHHH  =  -HMASS**2/V*RTHREE
 209       GWWHH =   EE2/SW2*RXHALF
 210       GZZHH =   EE2/SC2*RXHALF
 211       GHHHH = -(HMASS/V)**2*RTHREE
 212 C
 213       RETURN
 214       END
 215 C
 216 C ----------------------------------------------------------------------
 217 C
 218       SUBROUTINE COUP4X(SW2,ZMASS,FMASS , GCHF)
 219 C
 220 C This subroutine sets up the coupling constant for the fermion-fermion-
 221 C Higgs vertex in the STANDARD MODEL.  The coupling is COMPLEX and the
 222 C array of the coupling specifies the chirality of the flowing-IN
 223 C fermion.  GCHF(1) denotes a left-handed coupling, and GCHF(2) a right-
 224 C handed coupling.
 225 C
 226 C INPUT:
 227 C       real    SW2            : square of sine of the weak angle
 228 C       real    ZMASS          : Z       mass
 229 C       real    FMASS          : fermion mass
 230 C
 231 C OUTPUT:
 232 C       complex GCHF(2)        : coupling of fermion and Higgs
 233 C
 234 #if defined(CERNLIB_IMPNONE)
 235       IMPLICIT NONE
 236 #endif
 237       COMPLEX*16 GCHF(2)
 238       REAL*8    SW2,ZMASS,FMASS,ALPHA,FOURPI,EZ,G
 239 C
 240       ALPHA=1.D0/128.D0
 241 C      ALPHA=1./REAL(137.0359895)
 242       FOURPI=4.D0*3.14159265358979323846D0
 243       EZ=SQRT(ALPHA*FOURPI)/SQRT(SW2*(1.D0-SW2))
 244       G=EZ*FMASS*0.5D0/ZMASS
 245 C
 246       GCHF(1) = DCMPLX( -G )
 247       GCHF(2) = DCMPLX( -G )
 248 C
 249       RETURN
 250       END
 251 C
 252 C ======================================================================
 253 C
 254       SUBROUTINE EAIXXX(EB,EA,SHLF,CHLF,PHI,NHE,NHA , EAI)
 255 C
 256 C This subroutine computes an off-shell electron wavefunction after
 257 C emitting a photon from the electron beam, with a special care for the
 258 C small angle region.  The momenta are measured in the laboratory frame,
 259 C where the e- beam is along the positive z axis.
 260 C
 261 C INPUT:
 262 C       real    EB             : energy (GeV)    of beam  e-
 263 C       real    EA             : energy (GeV)    of final photon
 264 C       real    SHLF           : sin(theta/2)    of final photon
 265 C       real    CHLF           : cos(theta/2)    of final photon
 266 C       real    PHI            : azimuthal angle of final photon
 267 C       integer NHE  = -1 or 1 : helicity        of beam  e-
 268 C       integer NHA  = -1 or 1 : helicity        of final photon
 269 C
 270 C OUTPUT:
 271 C       complex EAI(6)         : off-shell electron             |e',A,e>
 272 C
 273 #if defined(CERNLIB_IMPNONE)
 274       IMPLICIT NONE
 275 #endif
 276       COMPLEX*16 EAI(6),PHS
 277       REAL*8  EB,EA,SHLF,CHLF,PHI,ME,ALPHA,GAL,RNHE,X,C,S,D,COEFF,
 278      &        XNNP,XNNM,SNP,CSP
 279       INTEGER NHE,NHA,NN
 280 C
 281       ME   = 0.51099906D-3
 282       ALPHA=1./128.
 283       GAL  =SQRT(ALPHA*4.*3.14159265D0)
 284 C
 285       NN=NHA*NHE
 286       RNHE=NHE
 287       X=EA/EB
 288       C=(CHLF+SHLF)*(CHLF-SHLF)
 289       S=2.*CHLF*SHLF
 290       D=-1./(EA*EB*(4.*SHLF**2+(ME/EB)**2*C))
 291       COEFF=-NN*GAL*SQRT(EB)*D
 292       XNNP=X*(1+NN)
 293       XNNM=X*(1-NN)
 294       SNP=SIN(PHI)
 295       CSP=COS(PHI)
 296       PHS=DCMPLX( CSP , RNHE*SNP )
 297 C
 298       EAI((5-3*NHE)/2) = -RNHE*COEFF*ME*S*(1.+XNNP*.5)
 299       EAI((5-NHE)/2)   =  XNNP*COEFF*ME*CHLF**2*PHS
 300       EAI((5+NHE)/2)   =  RNHE*COEFF*EB*S*(-2.+XNNM)
 301       EAI((5+3*NHE)/2) =  XNNM*COEFF*EB*SHLF**2*PHS*2.
 302 C
 303       EAI(5) =  EB*DCMPLX( 1.-X , 1.-X*C )
 304       EAI(6) = -EB*X*S*DCMPLX( CSP , SNP )
 305 C
 306       RETURN
 307       END
 308 C
 309 C ----------------------------------------------------------------------
 310 C
 311       SUBROUTINE EAOXXX(EB,EA,SHLF,CHLF,PHI,NHE,NHA , EAO)
 312 C
 313 C This subroutine computes an off-shell positron wavefunction after
 314 C emitting a photon from the positron beam, with a special care for the
 315 C small angle region.  The momenta are measured in the laboratory frame,
 316 C where the e+ beam is along the negative z axis.
 317 C
 318 C INPUT:
 319 C       real    EB             : energy (GeV)    of beam  e+
 320 C       real    EA             : energy (GeV)    of final photon
 321 C       real    SHLF           : sin(theta/2)    of final photon
 322 C       real    CHLF           : cos(theta/2)    of final photon
 323 C       real    PHI            : azimuthal angle of final photon
 324 C       integer NHE  = -1 or 1 : helicity        of beam  e+
 325 C       integer NHA  = -1 or 1 : helicity        of final photon
 326 C
 327 C OUTPUT:
 328 C       complex EAO(6)         : off-shell positron             <e,A,e'|
 329 C
 330 #if defined(CERNLIB_IMPNONE)
 331       IMPLICIT NONE
 332 #endif
 333       COMPLEX*16 EAO(6),PHS
 334       REAL*8  EB,EA,SHLF,CHLF,PHI,ME,ALPHA,GAL,RNHE,X,C,S,D,COEFF,
 335      &        XNNP,XNNM,SNP,CSP
 336       INTEGER NHE,NHA,NN
 337 C
 338       ME   = 0.51099906D-3
 339       ALPHA=1./128.
 340       GAL  =SQRT(ALPHA*4.*3.14159265D0)
 341 C
 342       NN=NHA*NHE
 343       RNHE=NHE
 344       X=EA/EB
 345       C=(CHLF+SHLF)*(CHLF-SHLF)
 346       S=2.*CHLF*SHLF
 347       D=-1./(EA*EB*(4.*CHLF**2-(ME/EB)**2*C))
 348       COEFF=NN*GAL*SQRT(EB)*D
 349       XNNP=X*(1+NN)
 350       XNNM=X*(1-NN)
 351       SNP=SIN(PHI)
 352       CSP=COS(PHI)
 353       PHS=DCMPLX( CSP ,-RNHE*SNP )
 354 C
 355       EAO((5-3*NHE)/2) =              COEFF*ME*S*(1.+XNNP*.5)
 356       EAO((5-NHE)/2)   = RNHE*XNNP    *COEFF*ME*SHLF**2*PHS
 357       EAO((5+NHE)/2)   =              COEFF*EB*S*(-2.+XNNM)
 358       EAO((5+3*NHE)/2) = REAL(NHA-NHE)*COEFF*EB*X*CHLF**2*PHS*2.
 359 C
 360       EAO(5) = EB*DCMPLX( X-1. , X*C+1. )
 361       EAO(6) = EB*X*S*DCMPLX( CSP , SNP )
 362 C
 363       RETURN
 364       END
 365 C
 366 C ----------------------------------------------------------------------
 367 C
 368       SUBROUTINE FSIXXX(FI,SC,GC,FMASS,FWIDTH , FSI)
 369 C
 370 C this subroutine computes an off-shell fermion wavefunction from a
 371 C flowing-in external fermion and a vector boson.
 372 C
 373 C input:
 374 C       complex*16 fi(6)          : flow-in  fermion                   |fi>
 375 C       complex*16 sc(3)          : input    scalar                      s
 376 C       complex*16 gc(2)          : coupling constants                 gchf
 377 C       real*8    fmass          : mass  of output fermion f'
 378 C       real*8    fwidth         : width of output fermion f'
 379 C
 380 C output:
 381 C       complex fsi(6)         : off-shell fermion             |f',s,fi>
 382 C
 383 #if defined(CERNLIB_IMPNONE)
 384       IMPLICIT NONE
 385 #endif
 386       COMPLEX*16 FI(6),SC(3),FSI(6),GC(2),SL1,SL2,SR1,SR2,DS
 387       REAL*8     PF(0:3),FMASS,FWIDTH,PF2,P0P3,P0M3
 388 C
 389       FSI(5) = FI(5)-SC(2)
 390       FSI(6) = FI(6)-SC(3)
 391 C
 392       PF(0)=DBLE( FSI(5))
 393       PF(1)=DBLE( FSI(6))
 394       PF(2)=DIMAG(FSI(6))
 395       PF(3)=DIMAG(FSI(5))
 396       PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
 397 C
 398       DS=-SC(1)/DCMPLX(PF2-FMASS**2,MAX(DSIGN(FMASS*FWIDTH ,PF2),0D0))
 399       P0P3=PF(0)+PF(3)
 400       P0M3=PF(0)-PF(3)
 401       SL1=GC(1)*(P0P3*FI(1)+DCONJG(FSI(6))*FI(2))
 402       SL2=GC(1)*(P0M3*FI(2)      +FSI(6) *FI(1))
 403       SR1=GC(2)*(P0M3*FI(3)-DCONJG(FSI(6))*FI(4))
 404       SR2=GC(2)*(P0P3*FI(4)      -FSI(6) *FI(3))
 405 C
 406       FSI(1) = ( GC(1)*FMASS*FI(1) + SR1 )*DS
 407       FSI(2) = ( GC(1)*FMASS*FI(2) + SR2 )*DS
 408       FSI(3) = ( GC(2)*FMASS*FI(3) + SL1 )*DS
 409       FSI(4) = ( GC(2)*FMASS*FI(4) + SL2 )*DS
 410 C
 411       RETURN
 412       END
 413 C
 414 C ----------------------------------------------------------------------
 415 C
 416       SUBROUTINE FSOXXX(FO,SC,GC,FMASS,FWIDTH , FSO)
 417 C
 418 C this subroutine computes an off-shell fermion wavefunction from a
 419 C flowing-out external fermion and a vector boson.
 420 C
 421 C input:
 422 C       complex*16 fo(6)          : flow-out fermion                   <fo|
 423 C       complex*16 sc(6)          : input    scalar                      s
 424 C       complex*16 gc(2)          : coupling constants                 gchf
 425 C       real*8     fmass          : mass  of output fermion f'
 426 C       real*8     fwidth         : width of output fermion f'
 427 C
 428 C output:
 429 C       complex fso(6)         : off-shell fermion             <fo,s,f'|
 430 C
 431 #if defined(CERNLIB_IMPNONE)
 432       IMPLICIT NONE
 433 #endif
 434       COMPLEX*16 FO(6),SC(6),FSO(6),GC(2),SL1,SL2,SR1,SR2,DS
 435       REAL*8     PF(0:3),FMASS,FWIDTH,PF2,P0P3,P0M3
 436 C
 437       FSO(5) = FO(5)+SC(2)
 438       FSO(6) = FO(6)+SC(3)
 439 C
 440       PF(0)=DBLE( FSO(5))
 441       PF(1)=DBLE( FSO(6))
 442       PF(2)=DIMAG(FSO(6))
 443       PF(3)=DIMAG(FSO(5))
 444       PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
 445 C
 446       DS=-SC(1)/DCMPLX(PF2-FMASS**2,MAX(DSIGN(FMASS*FWIDTH ,PF2),0D0))
 447       P0P3=PF(0)+PF(3)
 448       P0M3=PF(0)-PF(3)
 449       SL1=GC(2)*(P0P3*FO(3)      +FSO(6) *FO(4))
 450       SL2=GC(2)*(P0M3*FO(4)+DCONJG(FSO(6))*FO(3))
 451       SR1=GC(1)*(P0M3*FO(1)      -FSO(6) *FO(2))
 452       SR2=GC(1)*(P0P3*FO(2)-DCONJG(FSO(6))*FO(1))
 453 C
 454       FSO(1) = ( GC(1)*FMASS*FO(1) + SL1 )*DS
 455       FSO(2) = ( GC(1)*FMASS*FO(2) + SL2 )*DS
 456       FSO(3) = ( GC(2)*FMASS*FO(3) + SR1 )*DS
 457       FSO(4) = ( GC(2)*FMASS*FO(4) + SR2 )*DS
 458 C
 459       RETURN
 460       END
 461 C
 462 C ----------------------------------------------------------------------
 463 C
 464       SUBROUTINE FVIXXX(FI,VC,G,FMASS,FWIDTH , FVI)
 465 C
 466 C this subroutine computes an off-shell fermion wavefunction from a
 467 C flowing-in external fermion and a vector boson.
 468 C
 469 C input:
 470 C       complex fi(6)          : flow-in  fermion                   |fi>
 471 C       complex vc(6)          : input    vector                      v
 472 C       real    g(2)           : coupling constants                  gvf
 473 C       real    fmass          : mass  of output fermion f'
 474 C       real    fwidth         : width of output fermion f'
 475 C
 476 C output:
 477 C       complex fvi(6)         : off-shell fermion             |f',v,fi>
 478 C
 479 #if defined(CERNLIB_IMPNONE)
 480       IMPLICIT NONE
 481 #endif
 482       COMPLEX*16 FI(6),VC(6),FVI(6),SL1,SL2,SR1,SR2,D
 483       REAL*8    G(2),PF(0:3),FMASS,FWIDTH,PF2
 484       REAL*8 RXZERO, RXONE
 485       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
 486       COMPLEX*16 CXIMAG
 487 C
 488       LOGICAL FIRST
 489       SAVE CXIMAG,FIRST
 490       DATA FIRST/.TRUE./
 491 C
 492 C          Fix compilation with g77
 493       IF(FIRST) THEN
 494         FIRST=.FALSE.
 495         CXIMAG=DCMPLX( RXZERO, RXONE )
 496       ENDIF
 497 C
 498       FVI(5) = FI(5)-VC(5)
 499       FVI(6) = FI(6)-VC(6)
 500 C
 501       PF(0)=DBLE( FVI(5))
 502       PF(1)=DBLE( FVI(6))
 503       PF(2)=DIMAG(FVI(6))
 504       PF(3)=DIMAG(FVI(5))
 505       PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
 506 C
 507       D=-RXONE/DCMPLX( PF2-FMASS**2,MAX(SIGN(FMASS*FWIDTH,PF2),RXZERO))
 508       SL1= (VC(1)+       VC(4))*FI(1)
 509      &    +(VC(2)-CXIMAG*VC(3))*FI(2)
 510       SL2= (VC(2)+CXIMAG*VC(3))*FI(1)
 511      &    +(VC(1)-       VC(4))*FI(2)
 512 C
 513       IF ( G(2) .NE. RXZERO ) THEN
 514          SR1= (VC(1)-       VC(4))*FI(3)
 515      &       -(VC(2)-CXIMAG*VC(3))*FI(4)
 516          SR2=-(VC(2)+CXIMAG*VC(3))*FI(3)
 517      &       +(VC(1)+       VC(4))*FI(4)
 518 C
 519          FVI(1) = ( G(1)*((PF(0)-PF(3))*SL1 -DCONJG(FVI(6))*SL2)
 520      &             +G(2)*FMASS*SR1)*D
 521          FVI(2) = ( G(1)*(      -FVI(6)*SL1 +(PF(0)+PF(3))*SL2)
 522      &             +G(2)*FMASS*SR2)*D
 523          FVI(3) = ( G(2)*((PF(0)+PF(3))*SR1 +DCONJG(FVI(6))*SR2)
 524      &             +G(1)*FMASS*SL1)*D
 525          FVI(4) = ( G(2)*(       FVI(6)*SR1 +(PF(0)-PF(3))*SR2)
 526      &             +G(1)*FMASS*SL2)*D
 527 C
 528       ELSE
 529          FVI(1) = G(1)*((PF(0)-PF(3))*SL1 -DCONJG(FVI(6))*SL2)*D
 530          FVI(2) = G(1)*(      -FVI(6)*SL1 +(PF(0)+PF(3))*SL2)*D
 531          FVI(3) = G(1)*FMASS*SL1*D
 532          FVI(4) = G(1)*FMASS*SL2*D
 533       END IF
 534 C
 535       RETURN
 536       END
 537 C
 538 C ----------------------------------------------------------------------
 539 C
 540       SUBROUTINE FVOXXX(FO,VC,G,FMASS,FWIDTH , FVO)
 541 C
 542 C this subroutine computes an off-shell fermion wavefunction from a
 543 C flowing-out external fermion and a vector boson.
 544 C
 545 C input:
 546 C       complex fo(6)          : flow-out fermion                   <fo|
 547 C       complex vc(6)          : input    vector                      v
 548 C       real    g(2)           : coupling constants                  gvf
 549 C       real    fmass          : mass  of output fermion f'
 550 C       real    fwidth         : width of output fermion f'
 551 C
 552 C output:
 553 C       complex fvo(6)         : off-shell fermion             <fo,v,f'|
 554 C
 555 #if defined(CERNLIB_IMPNONE)
 556       IMPLICIT NONE
 557 #endif
 558       COMPLEX*16 FO(6),VC(6),FVO(6),SL1,SL2,SR1,SR2,D
 559       REAL*8    G(2),PF(0:3),FMASS,FWIDTH,PF2
 560       REAL*8 RXZERO, RXONE
 561       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
 562       COMPLEX*16 CXIMAG
 563       LOGICAL FIRST
 564       SAVE CXIMAG,FIRST
 565       DATA FIRST/.TRUE./
 566 C
 567 C          Fix compilation with g77
 568       IF(FIRST) THEN
 569         FIRST=.FALSE.
 570         CXIMAG=DCMPLX( RXZERO, RXONE )
 571       ENDIF
 572 C
 573       FVO(5) = FO(5)+VC(5)
 574       FVO(6) = FO(6)+VC(6)
 575 C
 576       PF(0)=DBLE( FVO(5))
 577       PF(1)=DBLE( FVO(6))
 578       PF(2)=DIMAG(FVO(6))
 579       PF(3)=DIMAG(FVO(5))
 580       PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
 581 C
 582       D=-RXONE/DCMPLX( PF2-FMASS**2,MAX(SIGN(FMASS*FWIDTH,PF2),RXZERO))
 583       SL1= (VC(1)+       VC(4))*FO(3)
 584      &    +(VC(2)+CXIMAG*VC(3))*FO(4)
 585       SL2= (VC(2)-CXIMAG*VC(3))*FO(3)
 586      &    +(VC(1)-       VC(4))*FO(4)
 587 C
 588       IF ( G(2) .NE. RXZERO ) THEN
 589          SR1= (VC(1)-       VC(4))*FO(1)
 590      &       -(VC(2)+CXIMAG*VC(3))*FO(2)
 591          SR2=-(VC(2)-CXIMAG*VC(3))*FO(1)
 592      &       +(VC(1)+       VC(4))*FO(2)
 593 C
 594          FVO(1) = ( G(2)*( (PF(0)+PF(3))*SR1        +FVO(6)*SR2)
 595      &             +G(1)*FMASS*SL1)*D
 596          FVO(2) = ( G(2)*( DCONJG(FVO(6))*SR1 +(PF(0)-PF(3))*SR2)
 597      &             +G(1)*FMASS*SL2)*D
 598          FVO(3) = ( G(1)*( (PF(0)-PF(3))*SL1        -FVO(6)*SL2)
 599      &             +G(2)*FMASS*SR1)*D
 600          FVO(4) = ( G(1)*(-DCONJG(FVO(6))*SL1 +(PF(0)+PF(3))*SL2)
 601      &             +G(2)*FMASS*SR2)*D
 602 C
 603       ELSE
 604          FVO(1) = G(1)*FMASS*SL1*D
 605          FVO(2) = G(1)*FMASS*SL2*D
 606          FVO(3) = G(1)*( (PF(0)-PF(3))*SL1        -FVO(6)*SL2)*D
 607          FVO(4) = G(1)*(-DCONJG(FVO(6))*SL1 +(PF(0)+PF(3))*SL2)*D
 608       END IF
 609 C
 610       RETURN
 611       END
 612 C
 613 C ----------------------------------------------------------------------
 614 C
 615       SUBROUTINE GGGGXX(WM,W31,WP,W32,G, VERTEX)
 616 C
 617 C this subroutine computes an amplitude of the four-point coupling of
 618 C the w-, w+ and two w3/z/a.  the amplitude includes the contributions
 619 C of w exchange diagrams.  the internal w propagator is given in unitary
 620 C gauge.  if one sets wmass=0.0, then the gggg vertex is given (see sect
 621 C 2.9.1 of the manual).
 622 C
 623 C input:
 624 C       complex wm(0:3)        : flow-out w-                         wm
 625 C       complex w31(0:3)       : first    w3/z/a                     w31
 626 C       complex wp(0:3)        : flow-out w+                         wp
 627 C       complex w32(0:3)       : second   w3/z/a                     w32
 628 C       real    g              : coupling of w31 with w-/w+
 629 C                                                  (see the table below)
 630 C
 631 C the possible sets of the inputs are as follows:
 632 C   -------------------------------------------
 633 C   |  wm  |  w31 |  wp  |  w32 |  g31 |  g32 |
 634 C   -------------------------------------------
 635 C   |  w-  |  w3  |  w+  |  w3  |  gw  |  gw  |
 636 C   |  w-  |  w3  |  w+  |  z   |  gw  | gwwz |
 637 C   |  w-  |  w3  |  w+  |  a   |  gw  | gwwa |
 638 C   |  w-  |  z   |  w+  |  z   | gwwz | gwwz |
 639 C   |  w-  |  z   |  w+  |  a   | gwwz | gwwa |
 640 C   |  w-  |  a   |  w+  |  a   | gwwa | gwwa |
 641 C   -------------------------------------------
 642 C where all the bosons are defined by the flowing-out quantum number.
 643 C
 644 C output:
 645 C       complex vertex         : amplitude          gamma(wm,w31,wp,w32)
 646 C
 647 #if defined(CERNLIB_IMPNONE)
 648       IMPLICIT NONE
 649 #endif
 650       COMPLEX*16    WM(6),W31(6),WP(6),W32(6),VERTEX
 651       COMPLEX*16 DV1(0:3),DV2(0:3),DV3(0:3),DV4(0:3),
 652      &           DVERTX,V12,V13,V14,V23,V24,V34
 653       REAL*8       PWM(0:3),PW31(0:3),PWP(0:3),PW32(0:3),G
 654       REAL*8     DP1(0:3),DP2(0:3),DP3(0:3),DP4(0:3)
 655       REAL*8 RXZERO, RXONE
 656       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
 657 C
 658       PWM(0)=DBLE( WM(5))
 659       PWM(1)=DBLE( WM(6))
 660       PWM(2)=DIMAG(WM(6))
 661       PWM(3)=DIMAG(WM(5))
 662       PWP(0)=DBLE( WP(5))
 663       PWP(1)=DBLE( WP(6))
 664       PWP(2)=DIMAG(WP(6))
 665       PWP(3)=DIMAG(WP(5))
 666       PW31(0)=DBLE( W31(5))
 667       PW31(1)=DBLE( W31(6))
 668       PW31(2)=DIMAG(W31(6))
 669       PW31(3)=DIMAG(W31(5))
 670       PW32(0)=DBLE( W32(5))
 671       PW32(1)=DBLE( W32(6))
 672       PW32(2)=DIMAG(W32(6))
 673       PW32(3)=DIMAG(W32(5))
 674 C
 675       DV1(0)=DCMPLX(WM(1))
 676       DV1(1)=DCMPLX(WM(2))
 677       DV1(2)=DCMPLX(WM(3))
 678       DV1(3)=DCMPLX(WM(4))
 679       DP1(0)=DBLE(PWM(0))
 680       DP1(1)=DBLE(PWM(1))
 681       DP1(2)=DBLE(PWM(2))
 682       DP1(3)=DBLE(PWM(3))
 683       DV2(0)=DCMPLX(W31(1))
 684       DV2(1)=DCMPLX(W31(2))
 685       DV2(2)=DCMPLX(W31(3))
 686       DV2(3)=DCMPLX(W31(4))
 687       DP2(0)=DBLE(PW31(0))
 688       DP2(1)=DBLE(PW31(1))
 689       DP2(2)=DBLE(PW31(2))
 690       DP2(3)=DBLE(PW31(3))
 691       DV3(0)=DCMPLX(WP(1))
 692       DV3(1)=DCMPLX(WP(2))
 693       DV3(2)=DCMPLX(WP(3))
 694       DV3(3)=DCMPLX(WP(4))
 695       DP3(0)=DBLE(PWP(0))
 696       DP3(1)=DBLE(PWP(1))
 697       DP3(2)=DBLE(PWP(2))
 698       DP3(3)=DBLE(PWP(3))
 699       DV4(0)=DCMPLX(W32(1))
 700       DV4(1)=DCMPLX(W32(2))
 701       DV4(2)=DCMPLX(W32(3))
 702       DV4(3)=DCMPLX(W32(4))
 703       DP4(0)=DBLE(PW32(0))
 704       DP4(1)=DBLE(PW32(1))
 705       DP4(2)=DBLE(PW32(2))
 706       DP4(3)=DBLE(PW32(3))
 707 C
 708       V12= DV1(0)*DV2(0)-DV1(1)*DV2(1)-DV1(2)*DV2(2)-DV1(3)*DV2(3)
 709       V13= DV1(0)*DV3(0)-DV1(1)*DV3(1)-DV1(2)*DV3(2)-DV1(3)*DV3(3)
 710       V14= DV1(0)*DV4(0)-DV1(1)*DV4(1)-DV1(2)*DV4(2)-DV1(3)*DV4(3)
 711       V23= DV2(0)*DV3(0)-DV2(1)*DV3(1)-DV2(2)*DV3(2)-DV2(3)*DV3(3)
 712       V24= DV2(0)*DV4(0)-DV2(1)*DV4(1)-DV2(2)*DV4(2)-DV2(3)*DV4(3)
 713       V34= DV3(0)*DV4(0)-DV3(1)*DV4(1)-DV3(2)*DV4(2)-DV3(3)*DV4(3)
 714
 715          DVERTX = V14*V23 -V13*V24
 716 C
 717          VERTEX = DCMPLX( DVERTX ) * (G*G)
 718 C
 719       RETURN
 720       END
 721 C
 722 C ======================================================================
 723 C
 724       SUBROUTINE GGGXXX(WM,WP,W3,G , VERTEX)
 725 C
 726 C this subroutine computes an amplitude of the three-point coupling of
 727 C the gauge bosons.
 728 C
 729 C input:
 730 C       complex wm(6)          : vector               flow-out w-
 731 C       complex wp(6)          : vector               flow-out w+
 732 C       complex w3(6)          : vector               j3 or a    or z
 733 C       real    g              : coupling constant    gw or gwwa or gwwz
 734 C
 735 C output:
 736 C       complex vertex         : amplitude               gamma(wm,wp,w3)
 737 C
 738 #if defined(CERNLIB_IMPNONE)
 739       IMPLICIT NONE
 740 #endif
 741       COMPLEX*16 WM(6),WP(6),W3(6),VERTEX,
 742      &        XV1,XV2,XV3,V12,V23,V31,P12,P13,P21,P23,P31,P32
 743       REAL*8    PWM(0:3),PWP(0:3),PW3(0:3),G
 744       REAL*8 RXZERO, RTENTH
 745       PARAMETER( RXZERO=0.0D0, RTENTH=0.1D0 )
 746 C
 747       PWM(0)=DBLE( WM(5))
 748       PWM(1)=DBLE( WM(6))
 749       PWM(2)=DIMAG(WM(6))
 750       PWM(3)=DIMAG(WM(5))
 751       PWP(0)=DBLE( WP(5))
 752       PWP(1)=DBLE( WP(6))
 753       PWP(2)=DIMAG(WP(6))
 754       PWP(3)=DIMAG(WP(5))
 755       PW3(0)=DBLE( W3(5))
 756       PW3(1)=DBLE( W3(6))
 757       PW3(2)=DIMAG(W3(6))
 758       PW3(3)=DIMAG(W3(5))
 759 C
 760       V12=WM(1)*WP(1)-WM(2)*WP(2)-WM(3)*WP(3)-WM(4)*WP(4)
 761       V23=WP(1)*W3(1)-WP(2)*W3(2)-WP(3)*W3(3)-WP(4)*W3(4)
 762       V31=W3(1)*WM(1)-W3(2)*WM(2)-W3(3)*WM(3)-W3(4)*WM(4)
 763       XV1=RXZERO
 764       XV2=RXZERO
 765       XV3=RXZERO
 766       IF ( ABS(WM(1)) .NE. RXZERO ) THEN
 767          IF (ABS(WM(1)).GE.MAX(ABS(WM(2)),ABS(WM(3)),ABS(WM(4)))
 768      $        *RTENTH)
 769      &      XV1=PWM(0)/WM(1)
 770       ENDIF
 771       IF ( ABS(WP(1)) .NE. RXZERO) THEN
 772          IF (ABS(WP(1)).GE.MAX(ABS(WP(2)),ABS(WP(3)),ABS(WP(4)))
 773      $        *RTENTH)
 774      &      XV2=PWP(0)/WP(1)
 775       ENDIF
 776       IF ( ABS(W3(1)) .NE. RXZERO) THEN
 777          IF ( ABS(W3(1)).GE.MAX(ABS(W3(2)),ABS(W3(3)),ABS(W3(4)))
 778      $        *RTENTH)
 779      &      XV3=PW3(0)/W3(1)
 780       ENDIF
 781       P12= (PWM(0)-XV1*WM(1))*WP(1)-(PWM(1)-XV1*WM(2))*WP(2)
 782      &    -(PWM(2)-XV1*WM(3))*WP(3)-(PWM(3)-XV1*WM(4))*WP(4)
 783       P13= (PWM(0)-XV1*WM(1))*W3(1)-(PWM(1)-XV1*WM(2))*W3(2)
 784      &    -(PWM(2)-XV1*WM(3))*W3(3)-(PWM(3)-XV1*WM(4))*W3(4)
 785       P21= (PWP(0)-XV2*WP(1))*WM(1)-(PWP(1)-XV2*WP(2))*WM(2)
 786      &    -(PWP(2)-XV2*WP(3))*WM(3)-(PWP(3)-XV2*WP(4))*WM(4)
 787       P23= (PWP(0)-XV2*WP(1))*W3(1)-(PWP(1)-XV2*WP(2))*W3(2)
 788      &    -(PWP(2)-XV2*WP(3))*W3(3)-(PWP(3)-XV2*WP(4))*W3(4)
 789       P31= (PW3(0)-XV3*W3(1))*WM(1)-(PW3(1)-XV3*W3(2))*WM(2)
 790      &    -(PW3(2)-XV3*W3(3))*WM(3)-(PW3(3)-XV3*W3(4))*WM(4)
 791       P32= (PW3(0)-XV3*W3(1))*WP(1)-(PW3(1)-XV3*W3(2))*WP(2)
 792      &    -(PW3(2)-XV3*W3(3))*WP(3)-(PW3(3)-XV3*W3(4))*WP(4)
 793 C
 794       VERTEX = -(V12*(P13-P23)+V23*(P21-P31)+V31*(P32-P12))*G
 795 C
 796       RETURN
 797       END
 798       SUBROUTINE HIOXXX(FI,FO,GC,SMASS,SWIDTH , HIO)
 799 C
 800 C this subroutine computes an off-shell scalar current from an external
 801 C fermion pair.
 802 C
 803 C input:
 804 C       complex fi(6)          : flow-in  fermion                   |fi>
 805 C       complex fo(6)          : flow-out fermion                   <fo|
 806 C       complex gc(2)          : coupling constants                 gchf
 807 C       real    smass          : mass  of output scalar s
 808 C       real    swidth         : width of output scalar s
 809 C
 810 C output:
 811 C       complex hio(3)         : scalar current             j(<fi|s|fo>)
 812 C
 813 #if defined(CERNLIB_IMPNONE)
 814       IMPLICIT NONE
 815 #endif
 816       COMPLEX*16 FI(6),FO(6),HIO(3),GC(2),DN
 817       REAL*8  Q(0:3),SMASS,SWIDTH,Q2
 818 C
 819       HIO(2) = FO(5)-FI(5)
 820       HIO(3) = FO(6)-FI(6)
 821 C
 822       Q(0)=DBLE( HIO(2))
 823       Q(1)=DBLE( HIO(3))
 824       Q(2)=DIMAG(HIO(3))
 825       Q(3)=DIMAG(HIO(2))
 826       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
 827 C
 828       DN=-DCMPLX(Q2-SMASS**2,DMAX1(DSIGN(SMASS*SWIDTH,Q2),0.D0))
 829 C
 830       HIO(1) = ( GC(1)*(FO(1)*FI(1)+FO(2)*FI(2))
 831      &          +GC(2)*(FO(3)*FI(3)+FO(4)*FI(4)) )/DN
 832 C
 833       RETURN
 834       END
 835 C ----------------------------------------------------------------------
 836 C
 837       SUBROUTINE HSSSXX(S1,S2,S3,G,SMASS,SWIDTH , HSSS)
 838 C
 839 C This subroutine computes an off-shell scalar current from the four-
 840 C scalar coupling.
 841 C
 842 C INPUT:
 843 C       complex S1(3)          : first  scalar                        S1
 844 C       complex S2(3)          : second scalar                        S2
 845 C       complex S3(3)          : third  scalar                        S3
 846 C       real    G              : coupling constant                 GHHHH
 847 C       real    SMASS          : mass  of OUTPUT scalar S'
 848 C       real    SWIDTH         : width of OUTPUT scalar S'
 849 C
 850 C OUTPUT:
 851 C       complex HSSS(3)        : scalar current           J(S':S1,S2,S3)
 852 C
 853 #if defined(CERNLIB_IMPNONE)
 854       IMPLICIT NONE
 855 #endif
 856       COMPLEX*16 S1(3),S2(3),S3(3),HSSS(3),DG
 857       REAL*8     Q(0:3),G,SMASS,SWIDTH,Q2
 858 C
 859       HSSS(2) = S1(2)+S2(2)+S3(2)
 860       HSSS(3) = S1(3)+S2(3)+S3(3)
 861 C
 862       Q(0)=DBLE( HSSS(2))
 863       Q(1)=DBLE( HSSS(3))
 864       Q(2)=DIMAG(HSSS(3))
 865       Q(3)=DIMAG(HSSS(2))
 866       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
 867 C
 868       DG=-G/DCMPLX( Q2-SMASS**2,MAX(SIGN(SMASS*SWIDTH ,Q2),0.D0))
 869 C
 870       HSSS(1) = DG * S1(1)*S2(1)*S3(1)
 871 C
 872       RETURN
 873       END
 874 C ----------------------------------------------------------------------
 875 C
 876       SUBROUTINE HSSXXX(S1,S2,G,SMASS,SWIDTH , HSS)
 877 C
 878 C This subroutine computes an off-shell scalar current from the three-
 879 C scalar coupling.
 880 C
 881 C INPUT:
 882 C       complex S1(3)          : first  scalar                        S1
 883 C       complex S2(3)          : second scalar                        S2
 884 C       real    G              : coupling constant                  GHHH
 885 C       real    SMASS          : mass  of OUTPUT scalar S'
 886 C       real    SWIDTH         : width of OUTPUT scalar S'
 887 C
 888 C OUTPUT:
 889 C       complex HSS(3)         : scalar current              J(S':S1,S2)
 890 C
 891 #if defined(CERNLIB_IMPNONE)
 892       IMPLICIT NONE
 893 #endif
 894       COMPLEX*16 S1(3),S2(3),HSS(3),DG
 895       REAL*8  Q(0:3),G,SMASS,SWIDTH,Q2
 896 C
 897       HSS(2) = S1(2)+S2(2)
 898       HSS(3) = S1(3)+S2(3)
 899 C
 900       Q(0)=DBLE( HSS(2))
 901       Q(1)=DBLE( HSS(3))
 902       Q(2)=DIMAG(HSS(3))
 903       Q(3)=DIMAG(HSS(2))
 904       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
 905 C
 906       DG=-G/DCMPLX( Q2-SMASS**2, MAX(SIGN(SMASS*SWIDTH ,Q2),0.D0))
 907 C
 908       HSS(1) = DG*S1(1)*S2(1)
 909 C
 910       RETURN
 911       END
 912 C
 913 C ======================================================================
 914 C ----------------------------------------------------------------------
 915 C
 916       SUBROUTINE HVSXXX(VC,SC,G,SMASS,SWIDTH , HVS)
 917 C
 918 C this subroutine computes an off-shell scalar current from the vector-
 919 C scalar-scalar coupling.  the coupling is absent in the minimal sm in
 920 C unitary gauge.
 921 C
 922 C input:
 923 C       complex vc(6)          : input vector                          v
 924 C       complex sc(3)          : input scalar                          s
 925 C       complex g              : coupling constant (s charge)
 926 C       real    smass          : mass  of output scalar s'
 927 C       real    swidth         : width of output scalar s'
 928 C
 929 C examples of the coupling constant g for susy particles are as follows:
 930 C   -----------------------------------------------------------
 931 C   |    s1    | (q,i3) of s1  ||   v=a   |   v=z   |   v=w   |
 932 C   -----------------------------------------------------------
 933 C   | nu~_l    | (  0  , +1/2) ||   ---   |  gzn(1) |  gwf(1) |
 934 C   | e~_l     | ( -1  , -1/2) ||  gal(1) |  gzl(1) |  gwf(1) |
 935 C   | u~_l     | (+2/3 , +1/2) ||  gau(1) |  gzu(1) |  gwf(1) |
 936 C   | d~_l     | (-1/3 , -1/2) ||  gad(1) |  gzd(1) |  gwf(1) |
 937 C   -----------------------------------------------------------
 938 C   | e~_r-bar | ( +1  ,  0  ) || -gal(2) | -gzl(2) | -gwf(2) |
 939 C   | u~_r-bar | (-2/3 ,  0  ) || -gau(2) | -gzu(2) | -gwf(2) |
 940 C   | d~_r-bar | (+1/3 ,  0  ) || -gad(2) | -gzd(2) | -gwf(2) |
 941 C   -----------------------------------------------------------
 942 C where the sc charge is defined by the flowing-out quantum number.
 943 C
 944 C output:
 945 C       complex hvs(3)         : scalar current                j(s':v,s)
 946 C
 947 #if defined(CERNLIB_IMPNONE)
 948       IMPLICIT NONE
 949 #endif
 950       COMPLEX*16 VC(6),SC(3),HVS(3),DG,QVV,QPV,G
 951       REAL*8    QV(0:3),QP(0:3),QA(0:3),SMASS,SWIDTH,Q2
 952 C
 953       HVS(2) = VC(5)+SC(2)
 954       HVS(3) = VC(6)+SC(3)
 955 C
 956       QV(0)=DBLE(  VC(5))
 957       QV(1)=DBLE(  VC(6))
 958       QV(2)=DIMAG( VC(6))
 959       QV(3)=DIMAG( VC(5))
 960       QP(0)=DBLE(  SC(2))
 961       QP(1)=DBLE(  SC(3))
 962       QP(2)=DIMAG( SC(3))
 963       QP(3)=DIMAG( SC(2))
 964       QA(0)=DBLE( HVS(2))
 965       QA(1)=DBLE( HVS(3))
 966       QA(2)=DIMAG(HVS(3))
 967       QA(3)=DIMAG(HVS(2))
 968       Q2=QA(0)**2-(QA(1)**2+QA(2)**2+QA(3)**2)
 969 C
 970       DG=-G/DCMPLX( Q2-SMASS**2 , MAX(DSIGN( SMASS*SWIDTH ,Q2),0D0) )
 971       QVV=QV(0)*VC(1)-QV(1)*VC(2)-QV(2)*VC(3)-QV(3)*VC(4)
 972       QPV=QP(0)*VC(1)-QP(1)*VC(2)-QP(2)*VC(3)-QP(3)*VC(4)
 973 C
 974       HVS(1) = DG*(2D0*QPV+QVV)*SC(1)
 975 C
 976       RETURN
 977       END
 978 C
 979 C ----------------------------------------------------------------------
 980 C
 981       SUBROUTINE HVVXXX(V1,V2,G,SMASS,SWIDTH , HVV)
 982 C
 983 C this subroutine computes an off-shell scalar current from the vector-
 984 C vector-scalar coupling.
 985 C
 986 C input:
 987 C       complex v1(6)          : first  vector                        v1
 988 C       complex v2(6)          : second vector                        v2
 989 C       real    g              : coupling constant                  gvvh
 990 C       real    smass          : mass  of output scalar s
 991 C       real    swidth         : width of output scalar s
 992 C
 993 C output:
 994 C       complex hvv(3)         : off-shell scalar current     j(s:v1,v2)
 995 C
 996 #if defined(CERNLIB_IMPNONE)
 997       IMPLICIT NONE
 998 #endif
 999       COMPLEX*16 V1(6),V2(6),HVV(3),DG
1000       REAL*8    Q(0:3),G,SMASS,SWIDTH,Q2
1001       REAL*8 RXZERO
1002       PARAMETER( RXZERO=0.0D0 )
1003 C
1004       HVV(2) = V1(5)+V2(5)
1005       HVV(3) = V1(6)+V2(6)
1006 C
1007       Q(0)=DBLE( HVV(2))
1008       Q(1)=DBLE( HVV(3))
1009       Q(2)=DIMAG(HVV(3))
1010       Q(3)=DIMAG(HVV(2))
1011       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1012 C
1013       DG=-G/DCMPLX( Q2-SMASS**2 , MAX(SIGN( SMASS*SWIDTH ,Q2),RXZERO) )
1014 C
1015       HVV(1) = DG*(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
1016 C
1017       RETURN
1018       END
1019 C
1020 C ======================================================================
1021 C
1022       SUBROUTINE IOSXXX(FI,FO,SC,GC , VERTEX)
1023 C
1024 C This subroutine computes an amplitude of the fermion-fermion-scalar
1025 C coupling.
1026 C
1027 C INPUT:
1028 C       complex FI(6)          : flow-in  fermion                   |FI>
1029 C       complex FO(6)          : flow-out fermion                   <FO|
1030 C       complex SC(3)          : input    scalar                      S
1031 C       complex GC(2)          : coupling constants                 GCHF
1032 C
1033 C OUTPUT:
1034 C       complex VERTEX         : amplitude                     <FO|S|FI>
1035 C
1036 #if defined(CERNLIB_IMPNONE)
1037       IMPLICIT NONE
1038 #endif
1039       COMPLEX*16 FI(6),FO(6),SC(3),GC(2),VERTEX
1040 C
1041       VERTEX = SC(1)*( GC(1)*(FI(1)*FO(1)+FI(2)*FO(2))
1042      &                +GC(2)*(FI(3)*FO(3)+FI(4)*FO(4)) )
1043 C
1044       RETURN
1045       END
1046 C
1047 C ======================================================================
1048 C
1049       SUBROUTINE IOVXXX(FI,FO,VC,G , VERTEX)
1050 C
1051 C this subroutine computes an amplitude of the fermion-fermion-vector
1052 C coupling.
1053 C
1054 C input:
1055 C       complex fi(6)          : flow-in  fermion                   |fi>
1056 C       complex fo(6)          : flow-out fermion                   <fo|
1057 C       complex vc(6)          : input    vector                      v
1058 C       real    g(2)           : coupling constants                  gvf
1059 C
1060 C output:
1061 C       complex vertex         : amplitude                     <fo|v|fi>
1062 C
1063 #if defined(CERNLIB_IMPNONE)
1064       IMPLICIT NONE
1065 #endif
1066       COMPLEX*16 FI(6),FO(6),VC(6),VERTEX
1067       REAL*8    G(2)
1068       REAL*8 RXZERO, RXONE
1069       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1070       COMPLEX*16 CXIMAG
1071       LOGICAL FIRST
1072       SAVE CXIMAG,FIRST
1073       DATA FIRST/.TRUE./
1074 C
1075 C          Fix compilation with g77
1076       IF(FIRST) THEN
1077         FIRST=.FALSE.
1078         CXIMAG=DCMPLX( RXZERO, RXONE )
1079       ENDIF
1080 C
1081
1082       VERTEX =  G(1)*( (FO(3)*FI(1)+FO(4)*FI(2))*VC(1)
1083      &                +(FO(3)*FI(2)+FO(4)*FI(1))*VC(2)
1084      &                -(FO(3)*FI(2)-FO(4)*FI(1))*VC(3)*CXIMAG
1085      &                +(FO(3)*FI(1)-FO(4)*FI(2))*VC(4)        )
1086 C
1087       IF ( G(2) .NE. RXZERO ) THEN
1088          VERTEX = VERTEX
1089      &          + G(2)*( (FO(1)*FI(3)+FO(2)*FI(4))*VC(1)
1090      &                  -(FO(1)*FI(4)+FO(2)*FI(3))*VC(2)
1091      &                  +(FO(1)*FI(4)-FO(2)*FI(3))*VC(3)*CXIMAG
1092      &                  -(FO(1)*FI(3)-FO(2)*FI(4))*VC(4)        )
1093       END IF
1094 C
1095       RETURN
1096       END
1097 C
1098 C       Subroutine returns the desired fermion or
1099 C       anti-fermion spinor. ie., |f>
1100 C       A replacement for the HELAS routine IXXXXX
1101 C
1102 C       Adam Duff,  1992 August 31
1103 C       <duff@phenom.physics.wisc.edu>
1104 C
1105       SUBROUTINE IXXXXX(P,FMASS,NHEL,NSF,FI)
1106 C          P        IN: FOUR VECTOR MOMENTUM
1107 C          FMASS    IN: FERMION MASS
1108 C          NHEL     IN: SPINOR HELICITY, -1 OR 1
1109 C          NSF      IN: -1=ANTIFERMION, 1=FERMION
1110 C          FI       OUT: FERMION WAVEFUNCTION
1111 C
1112 C declare input/output variables
1113 C
1114 #if defined(CERNLIB_IMPNONE)
1115       IMPLICIT NONE
1116 #endif
1117       COMPLEX*16 FI(6)
1118       INTEGER*4 NHEL, NSF
1119       REAL*8 P(0:3), FMASS
1120       REAL*8 RXZERO, RXONE, RXTWO
1121       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
1122       REAL*8 PLAT, PABS, OMEGAP, OMEGAM, RS2PA, SPAZ
1123       COMPLEX*16 CXZERO
1124 C
1125 C declare local variables
1126 C
1127       LOGICAL FIRST
1128       SAVE CXZERO,FIRST
1129       DATA FIRST/.TRUE./
1130 C
1131 C          Fix compilation with g77
1132       IF(FIRST) THEN
1133         FIRST=.FALSE.
1134         CXZERO=DCMPLX( RXZERO, RXZERO )
1135       ENDIF
1136 C
1137 C define kinematic parameters
1138 C
1139       FI(5) = DCMPLX( P(0), P(3) ) * NSF
1140       FI(6) = DCMPLX( P(1), P(2) ) * NSF
1141       PLAT = SQRT( P(1)**2 + P(2)**2 )
1142       PABS = SQRT( P(1)**2 + P(2)**2 + P(3)**2 )
1143       OMEGAP = SQRT( P(0) + PABS )
1144 C
1145 C do massive fermion case
1146 C
1147       IF ( FMASS .NE. RXZERO ) THEN
1148          OMEGAM = FMASS / OMEGAP
1149          IF ( NSF .EQ. 1 ) THEN
1150             IF ( NHEL .EQ. 1 ) THEN
1151                IF ( P(3) .GE. RXZERO ) THEN
1152                   IF ( PLAT .EQ. RXZERO ) THEN
1153                      FI(1) = DCMPLX( OMEGAM, RXZERO )
1154                      FI(2) = CXZERO
1155                      FI(3) = DCMPLX( OMEGAP, RXZERO )
1156                      FI(4) = CXZERO
1157                   ELSE
1158                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1159                      SPAZ = SQRT( PABS + P(3) )
1160                      FI(1) = OMEGAM * RS2PA
1161      &                     * DCMPLX( SPAZ, RXZERO )
1162                      FI(2) = OMEGAM * RS2PA / SPAZ
1163      &                     * DCMPLX( P(1), P(2) )
1164                      FI(3) = OMEGAP * RS2PA
1165      &                     * DCMPLX( SPAZ, RXZERO )
1166                      FI(4) = OMEGAP * RS2PA / SPAZ
1167      &                     * DCMPLX( P(1), P(2) )
1168                   END IF
1169                ELSE
1170                   IF ( PLAT .EQ. RXZERO ) THEN
1171                      FI(1) = CXZERO
1172                      FI(2) = DCMPLX( OMEGAM, RXZERO )
1173                      FI(3) = CXZERO
1174                      FI(4) = DCMPLX( OMEGAP, RXZERO )
1175                   ELSE
1176                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1177                      SPAZ = SQRT( PABS - P(3) )
1178                      FI(1) = OMEGAM * RS2PA / SPAZ
1179      &                     * DCMPLX( PLAT, RXZERO )
1180                      FI(2) = OMEGAM * RS2PA * SPAZ / PLAT
1181      &                     * DCMPLX( P(1), P(2) )
1182                      FI(3) = OMEGAP * RS2PA / SPAZ
1183      &                     * DCMPLX( PLAT, RXZERO )
1184                      FI(4) = OMEGAP * RS2PA * SPAZ / PLAT
1185      &                     * DCMPLX( P(1), P(2) )
1186                   END IF
1187                END IF
1188             ELSE IF ( NHEL .EQ. -1 ) THEN
1189                IF ( P(3) .GE. RXZERO ) THEN
1190                   IF ( PLAT .EQ. RXZERO ) THEN
1191                      FI(1) = CXZERO
1192                      FI(2) = DCMPLX( OMEGAP, RXZERO )
1193                      FI(3) = CXZERO
1194                      FI(4) = DCMPLX( OMEGAM, RXZERO )
1195                   ELSE
1196                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1197                      SPAZ = SQRT( PABS + P(3) )
1198                      FI(1) = OMEGAP * RS2PA / SPAZ
1199      &                     * DCMPLX( -P(1), P(2) )
1200                      FI(2) = OMEGAP * RS2PA
1201      &                     * DCMPLX( SPAZ, RXZERO )
1202                      FI(3) = OMEGAM * RS2PA / SPAZ
1203      &                     * DCMPLX( -P(1), P(2) )
1204                      FI(4) = OMEGAM * RS2PA
1205      &                     * DCMPLX( SPAZ, RXZERO )
1206                   END IF
1207                ELSE
1208                   IF ( PLAT .EQ. RXZERO ) THEN
1209                      FI(1) = DCMPLX( -OMEGAP, RXZERO )
1210                      FI(2) = CXZERO
1211                      FI(3) = DCMPLX( -OMEGAM, RXZERO )
1212                      FI(4) = CXZERO
1213                   ELSE
1214                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1215                      SPAZ = SQRT( PABS - P(3) )
1216                      FI(1) = OMEGAP * RS2PA * SPAZ / PLAT
1217      &                     * DCMPLX( -P(1), P(2) )
1218                      FI(2) = OMEGAP * RS2PA / SPAZ
1219      &                     * DCMPLX( PLAT, RXZERO )
1220                      FI(3) = OMEGAM * RS2PA * SPAZ / PLAT
1221      &                     * DCMPLX( -P(1), P(2) )
1222                      FI(4) = OMEGAM * RS2PA / SPAZ
1223      &                     * DCMPLX( PLAT, RXZERO )
1224                   END IF
1225                END IF
1226             ELSE
1227                STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
1228             END IF
1229          ELSE IF ( NSF .EQ. -1 ) THEN
1230             IF ( NHEL .EQ. 1 ) THEN
1231                IF ( P(3) .GE. RXZERO ) THEN
1232                   IF ( PLAT .EQ. RXZERO ) THEN
1233                      FI(1) = CXZERO
1234                      FI(2) = DCMPLX( -OMEGAP, RXZERO )
1235                      FI(3) = CXZERO
1236                      FI(4) = DCMPLX( OMEGAM, RXZERO )
1237                   ELSE
1238                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1239                      SPAZ = SQRT( PABS + P(3) )
1240                      FI(1) = -OMEGAP * RS2PA / SPAZ
1241      &                     * DCMPLX( -P(1), P(2) )
1242                      FI(2) = -OMEGAP * RS2PA
1243      &                     * DCMPLX( SPAZ, RXZERO )
1244                      FI(3) = OMEGAM * RS2PA / SPAZ
1245      &                     * DCMPLX( -P(1), P(2) )
1246                      FI(4) = OMEGAM * RS2PA
1247      &                     * DCMPLX( SPAZ, RXZERO )
1248                   END IF
1249                ELSE
1250                   IF ( PLAT .EQ. RXZERO ) THEN
1251                      FI(1) = DCMPLX( OMEGAP, RXZERO )
1252                      FI(2) = CXZERO
1253                      FI(3) = DCMPLX( -OMEGAM, RXZERO )
1254                      FI(4) = CXZERO
1255                   ELSE
1256                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1257                      SPAZ = SQRT( PABS - P(3) )
1258                      FI(1) = -OMEGAP * RS2PA * SPAZ / PLAT
1259      &                     * DCMPLX( -P(1), P(2) )
1260                      FI(2) = -OMEGAP * RS2PA / SPAZ
1261      &                     * DCMPLX( PLAT, RXZERO )
1262                      FI(3) = OMEGAM * RS2PA * SPAZ / PLAT
1263      &                     * DCMPLX( -P(1), P(2) )
1264                      FI(4) = OMEGAM * RS2PA / SPAZ
1265      &                     * DCMPLX( PLAT, RXZERO )
1266                   END IF
1267                END IF
1268             ELSE IF ( NHEL .EQ. -1 ) THEN
1269                IF ( P(3) .GE. RXZERO ) THEN
1270                   IF ( PLAT .EQ. RXZERO ) THEN
1271                      FI(1) = DCMPLX( OMEGAM, RXZERO )
1272                      FI(2) = CXZERO
1273                      FI(3) = DCMPLX( -OMEGAP, RXZERO )
1274                      FI(4) = CXZERO
1275                   ELSE
1276                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1277                      SPAZ = SQRT( PABS + P(3) )
1278                      FI(1) = OMEGAM * RS2PA
1279      &                     * DCMPLX( SPAZ, RXZERO )
1280                      FI(2) = OMEGAM * RS2PA / SPAZ
1281      &                     * DCMPLX( P(1), P(2) )
1282                      FI(3) = -OMEGAP * RS2PA
1283      &                     * DCMPLX( SPAZ, RXZERO )
1284                      FI(4) = -OMEGAP * RS2PA / SPAZ
1285      &                     * DCMPLX( P(1), P(2) )
1286                   END IF
1287                ELSE
1288                   IF ( PLAT .EQ. RXZERO ) THEN
1289                      FI(1) = CXZERO
1290                      FI(2) = DCMPLX( OMEGAM, RXZERO )
1291                      FI(3) = CXZERO
1292                      FI(4) = DCMPLX( -OMEGAP, RXZERO )
1293                   ELSE
1294                      RS2PA = RXONE / SQRT( RXTWO * PABS )
1295                      SPAZ = SQRT( PABS - P(3) )
1296                      FI(1) = OMEGAM * RS2PA / SPAZ
1297      &                     * DCMPLX( PLAT, RXZERO )
1298                      FI(2) = OMEGAM * RS2PA * SPAZ / PLAT
1299      &                     * DCMPLX( P(1), P(2) )
1300                      FI(3) = -OMEGAP * RS2PA / SPAZ
1301      &                     * DCMPLX( PLAT, RXZERO )
1302                      FI(4) = -OMEGAP * RS2PA * SPAZ / PLAT
1303      &                     * DCMPLX( P(1), P(2) )
1304                   END IF
1305                END IF
1306             ELSE
1307                STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
1308             END IF
1309          ELSE
1310             STOP 'IXXXXX:  FERMION TYPE MUST BE +1,-1'
1311          END IF
1312 C
1313 C do massless fermion case
1314 C
1315       ELSE
1316          IF ( NSF .EQ. 1 ) THEN
1317             IF ( NHEL .EQ. 1 ) THEN
1318                IF ( P(3) .GE. RXZERO ) THEN
1319                   IF ( PLAT .EQ. RXZERO ) THEN
1320                      FI(1) = CXZERO
1321                      FI(2) = CXZERO
1322                      FI(3) = DCMPLX( OMEGAP, RXZERO )
1323                      FI(4) = CXZERO
1324                   ELSE
1325                      SPAZ = SQRT( PABS + P(3) )
1326                      FI(1) = CXZERO
1327                      FI(2) = CXZERO
1328                      FI(3) = DCMPLX( SPAZ, RXZERO )
1329                      FI(4) = RXONE / SPAZ
1330      &                     * DCMPLX( P(1), P(2) )
1331                   END IF
1332                ELSE
1333                   IF ( PLAT .EQ. RXZERO ) THEN
1334                      FI(1) = CXZERO
1335                      FI(2) = CXZERO
1336                      FI(3) = CXZERO
1337                      FI(4) = DCMPLX( OMEGAP, RXZERO )
1338                   ELSE
1339                      SPAZ = SQRT( PABS - P(3) )
1340                      FI(1) = CXZERO
1341                      FI(2) = CXZERO
1342                      FI(3) = RXONE / SPAZ
1343      &                     * DCMPLX( PLAT, RXZERO )
1344                      FI(4) = SPAZ / PLAT
1345      &                     * DCMPLX( P(1), P(2) )
1346                   END IF
1347                END IF
1348             ELSE IF ( NHEL .EQ. -1 ) THEN
1349                IF ( P(3) .GE. RXZERO ) THEN
1350                   IF ( PLAT .EQ. RXZERO ) THEN
1351                      FI(1) = CXZERO
1352                      FI(2) = DCMPLX( OMEGAP, RXZERO )
1353                      FI(3) = CXZERO
1354                      FI(4) = CXZERO
1355                   ELSE
1356                      SPAZ = SQRT( PABS + P(3) )
1357                      FI(1) = RXONE / SPAZ
1358      &                     * DCMPLX( -P(1), P(2) )
1359                      FI(2) = DCMPLX( SPAZ, RXZERO )
1360                      FI(3) = CXZERO
1361                      FI(4) = CXZERO
1362                   END IF
1363                ELSE
1364                   IF ( PLAT .EQ. RXZERO ) THEN
1365                      FI(1) = DCMPLX( -OMEGAP, RXZERO )
1366                      FI(2) = CXZERO
1367                      FI(3) = CXZERO
1368                      FI(4) = CXZERO
1369                   ELSE
1370                      SPAZ = SQRT( PABS - P(3) )
1371                      FI(1) = SPAZ / PLAT
1372      &                     * DCMPLX( -P(1), P(2) )
1373                      FI(2) = RXONE / SPAZ
1374      &                     * DCMPLX( PLAT, RXZERO )
1375                      FI(3) = CXZERO
1376                      FI(4) = CXZERO
1377                   END IF
1378                END IF
1379             ELSE
1380                STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
1381             END IF
1382          ELSE IF ( NSF .EQ. -1 ) THEN
1383             IF ( NHEL .EQ. 1 ) THEN
1384                IF ( P(3) .GE. RXZERO ) THEN
1385                   IF ( PLAT .EQ. RXZERO ) THEN
1386                      FI(1) = CXZERO
1387                      FI(2) = DCMPLX( -OMEGAP, RXZERO )
1388                      FI(3) = CXZERO
1389                      FI(4) = CXZERO
1390                   ELSE
1391                      SPAZ = SQRT( PABS + P(3) )
1392                      FI(1) = -RXONE / SPAZ
1393      &                     * DCMPLX( -P(1), P(2) )
1394                      FI(2) = DCMPLX( -SPAZ, RXZERO )
1395                      FI(3) = CXZERO
1396                      FI(4) = CXZERO
1397                   END IF
1398                ELSE
1399                   IF ( PLAT .EQ. RXZERO ) THEN
1400                      FI(1) = DCMPLX( OMEGAP, RXZERO )
1401                      FI(2) = CXZERO
1402                      FI(3) = CXZERO
1403                      FI(4) = CXZERO
1404                   ELSE
1405                      SPAZ = SQRT( PABS - P(3) )
1406                      FI(1) = -SPAZ / PLAT
1407      &                     * DCMPLX( -P(1), P(2) )
1408                      FI(2) = -RXONE / SPAZ
1409      &                     * DCMPLX( PLAT, RXZERO )
1410                      FI(3) = CXZERO
1411                      FI(4) = CXZERO
1412                   END IF
1413                END IF
1414             ELSE IF ( NHEL .EQ. -1 ) THEN
1415                IF ( P(3) .GE. RXZERO ) THEN
1416                   IF ( PLAT .EQ. RXZERO ) THEN
1417                      FI(1) = CXZERO
1418                      FI(2) = CXZERO
1419                      FI(3) = DCMPLX( -OMEGAP, RXZERO )
1420                      FI(4) = CXZERO
1421                   ELSE
1422                      SPAZ = SQRT( PABS + P(3) )
1423                      FI(1) = CXZERO
1424                      FI(2) = CXZERO
1425                      FI(3) = DCMPLX( -SPAZ, RXZERO )
1426                      FI(4) = -RXONE / SPAZ
1427      &                     * DCMPLX( P(1), P(2) )
1428                   END IF
1429                ELSE
1430                   IF ( PLAT .EQ. RXZERO ) THEN
1431                      FI(1) = CXZERO
1432                      FI(2) = CXZERO
1433                      FI(3) = CXZERO
1434                      FI(4) = DCMPLX( -OMEGAP, RXZERO )
1435                   ELSE
1436                      SPAZ = SQRT( PABS - P(3) )
1437                      FI(1) = CXZERO
1438                      FI(2) = CXZERO
1439                      FI(3) = -RXONE / SPAZ
1440      &                     * DCMPLX( PLAT, RXZERO )
1441                      FI(4) = -SPAZ / PLAT
1442      &                     * DCMPLX( P(1), P(2) )
1443                   END IF
1444                END IF
1445             ELSE
1446                STOP 'IXXXXX:  FERMION HELICITY MUST BE +1,-1'
1447             END IF
1448          ELSE
1449             STOP 'IXXXXX:  FERMION TYPE MUST BE +1,-1'
1450          END IF
1451       END IF
1452 C
1453 C done
1454 C
1455       RETURN
1456       END
1457 C
1458 C ----------------------------------------------------------------------
1459 C
1460       SUBROUTINE J3XXXX(FI,FO,GAF,GZF,ZMASS,ZWIDTH , J3)
1461 C
1462 C this subroutine computes the sum of photon and z currents with the
1463 C suitable weights ( j(w3) = cos(theta_w) j(z) + sin(theta_w) j(a) ).
1464 C the output j3 is useful as an input of vvvxxx, jvvxxx or w3w3xx.
1465 C the photon propagator is given in feynman gauge, and the z propagator
1466 C is given in unitary gauge.
1467 C
1468 C input:
1469 C       complex fi(6)          : flow-in  fermion                   |fi>
1470 C       complex fo(6)          : flow-out fermion                   <fo|
1471 C       real    gaf(2)         : fi couplings with a                 gaf
1472 C       real    gzf(2)         : fi couplings with z                 gzf
1473 C       real    zmass          : mass  of z
1474 C       real    zwidth         : width of z
1475 C
1476 C output:
1477 C       complex j3(6)          : w3 current             j^mu(<fo|w3|fi>)
1478 C
1479 #if defined(CERNLIB_IMPNONE)
1480       IMPLICIT NONE
1481 #endif
1482       COMPLEX*16 FI(6),FO(6),J3(6),
1483      &        C0L,C1L,C2L,C3L,CSL,C0R,C1R,C2R,C3R,CSR,DZ,DDIF
1484       REAL*8    GAF(2),GZF(2),Q(0:3),ZMASS,ZWIDTH,ZM2,ZMW,Q2,DA,WW,
1485      &        CW,SW,GN,GZ3L,GA3L
1486 C
1487       REAL*8 RXZERO, RXONE
1488       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1489       COMPLEX*16 CXIMAG
1490       LOGICAL FIRST
1491       SAVE CXIMAG,FIRST
1492       DATA FIRST/.TRUE./
1493 C
1494 C          Fix compilation with g77
1495       IF(FIRST) THEN
1496         FIRST=.FALSE.
1497         CXIMAG=DCMPLX( RXZERO, RXONE )
1498       ENDIF
1499 C
1500       J3(5) = FO(5)-FI(5)
1501       J3(6) = FO(6)-FI(6)
1502 C
1503       Q(0)=-DBLE( J3(5))
1504       Q(1)=-DBLE( J3(6))
1505       Q(2)=-DIMAG(J3(6))
1506       Q(3)=-DIMAG(J3(5))
1507       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1508       ZM2=ZMASS**2
1509       ZMW=ZMASS*ZWIDTH
1510 C
1511       DA=RXONE/Q2
1512       WW=MAX(DSIGN( ZMW ,Q2),RXZERO)
1513       DZ=RXONE/DCMPLX( Q2-ZM2 , WW )
1514       DDIF=DCMPLX( -ZM2 , WW )*DA*DZ
1515 C
1516 C ddif is the difference : ddif=da-dz
1517 C  for the running width, use below instead of the above ww,dz and ddif.
1518 C      ww=max( zwidth*q2/zmass ,r_zero)
1519 C      dz=r_one/dcmplx( q2-zm2 , ww )
1520 C      ddif=dcmplx( -zm2 , ww )*da*dz
1521 C
1522       CW=RXONE/SQRT(RXONE+(GZF(2)/GAF(2))**2)
1523       SW=SQRT((RXONE-CW)*(RXONE+CW))
1524       GN=GAF(2)*SW
1525       GZ3L=GZF(1)*CW
1526       GA3L=GAF(1)*SW
1527       C0L=  FO(3)*FI(1)+FO(4)*FI(2)
1528       C0R=  FO(1)*FI(3)+FO(2)*FI(4)
1529       C1L=-(FO(3)*FI(2)+FO(4)*FI(1))
1530       C1R=  FO(1)*FI(4)+FO(2)*FI(3)
1531       C2L= (FO(3)*FI(2)-FO(4)*FI(1))*CXIMAG
1532       C2R=(-FO(1)*FI(4)+FO(2)*FI(3))*CXIMAG
1533       C3L= -FO(3)*FI(1)+FO(4)*FI(2)
1534       C3R=  FO(1)*FI(3)-FO(2)*FI(4)
1535       CSL=(Q(0)*C0L-Q(1)*C1L-Q(2)*C2L-Q(3)*C3L)/ZM2
1536       CSR=(Q(0)*C0R-Q(1)*C1R-Q(2)*C2R-Q(3)*C3R)/ZM2
1537 C
1538       J3(1) =  GZ3L*DZ*(C0L-CSL*Q(0))+GA3L*C0L*DA
1539      &       + GN*(C0R*DDIF-CSR*Q(0)*DZ)
1540       J3(2) =  GZ3L*DZ*(C1L-CSL*Q(1))+GA3L*C1L*DA
1541      &       + GN*(C1R*DDIF-CSR*Q(1)*DZ)
1542       J3(3) =  GZ3L*DZ*(C2L-CSL*Q(2))+GA3L*C2L*DA
1543      &       + GN*(C2R*DDIF-CSR*Q(2)*DZ)
1544       J3(4) =  GZ3L*DZ*(C3L-CSL*Q(3))+GA3L*C3L*DA
1545      &       + GN*(C3R*DDIF-CSR*Q(3)*DZ)
1546 C
1547       RETURN
1548       END
1549 C
1550 C ----------------------------------------------------------------------
1551 C
1552       SUBROUTINE JEEXXX(EB,EF,SHLF,CHLF,PHI,NHB,NHF,NSF , JEE)
1553 C
1554 C This subroutine computes an off-shell photon wavefunction emitted from
1555 C the electron or positron beam, with a special care for the small angle
1556 C region.  The momenta are measured in the laboratory frame, where the
1557 C e- (e+) beam is along the positive (negative) z axis.
1558 C
1559 C INPUT:
1560 C       real    EB             : energy (GeV)    of beam  e-/e+
1561 C       real    EF             : energy (GeV)    of final e-/e+
1562 C       real    SHLF           : sin(theta/2)    of final e-/e+
1563 C       real    CHLF           : cos(theta/2)    of final e-/e+
1564 C       real    PHI            : azimuthal angle of final e-/e+
1565 C       integer NHB  = -1 or 1 : helicity        of beam  e-/e+
1566 C       integer NHF  = -1 or 1 : helicity        of final e-/e+
1567 C       integer NSF  = -1 or 1 : +1 for electron, -1 for positron
1568 C
1569 C OUTPUT:
1570 C       complex JEE(6)         : off-shell photon          J^mu(<e|A|e>)
1571 C
1572 #if defined(CERNLIB_IMPNONE)
1573       IMPLICIT NONE
1574 #endif
1575       COMPLEX*16 JEE(6),COEFF
1576       REAL*8  CS(2),EB,EF,SHLF,CHLF,PHI,ME,ALPHA,GAL,HI,SF,SFH,X,ME2,Q2,
1577      &        RFP,RFM,SNP,CSP,RXC,C,S
1578       INTEGER NHB,NHF,NSF
1579 C
1580       ME   =0.51099906D-3
1581       ALPHA=1./128.
1582       GAL  =SQRT(ALPHA*4.*3.14159265D0)
1583 C
1584       HI =NHB
1585       SF =NSF
1586       SFH=NHB*NSF
1587       CS((3+NSF)/2)=SHLF
1588       CS((3-NSF)/2)=CHLF
1589 C CS(1)=CHLF and CS(2)=SHLF for electron
1590 C CS(1)=SHLF and CS(2)=CHLF for positron
1591       X=EF/EB
1592       ME2=ME**2
1593       Q2=-4.*CS(2)**2*(EF*EB-ME2)
1594      &   +SF*(1.-X)**2/X*(SHLF+CHLF)*(SHLF-CHLF)*ME2
1595       RFP=(1+NSF)
1596       RFM=(1-NSF)
1597       SNP=SIN(PHI)
1598       CSP=COS(PHI)
1599 C
1600       IF (NHB.EQ.NHF) THEN
1601          RXC=2.*X/(1.-X)*CS(1)**2
1602          COEFF= GAL*2.*EB*SQRT(X)*CS(2)/Q2
1603      &         *(DCMPLX( RFP )-RFM*DCMPLX( CSP ,-SNP*HI ))*.5
1604          JEE(1) =  DCMPLX( 0.D0 )
1605          JEE(2) =  COEFF*DCMPLX( (1.+RXC)*CSP ,-SFH*SNP )
1606          JEE(3) =  COEFF*DCMPLX( (1.+RXC)*SNP , SFH*CSP )
1607          JEE(4) =  COEFF*(-SF*RXC/CS(1)*CS(2))
1608       ELSE
1609          COEFF= GAL*ME/Q2/SQRT(X)
1610      &         *(DCMPLX( RFP )+RFM*DCMPLX( CSP , SNP*HI ))*.5*HI
1611          JEE(1) = -COEFF*(1.+X)*CS(2)*DCMPLX( CSP , SFH*SNP )
1612          JEE(2) =  COEFF*(1.-X)*CS(1)
1613          JEE(3) =  JEE(2)*DCMPLX( 0.D0 , SFH )
1614          JEE(4) =  JEE(1)*SF*(1.-X)/(1.+X)
1615       ENDIF
1616 C
1617       C=(CHLF+SHLF)*(CHLF-SHLF)
1618       S=2.*CHLF*SHLF
1619 C
1620       JEE(5) = -EB*DCMPLX( 1.-X , SF-X*C )
1621       JEE(6) =  EB*X*S*DCMPLX( CSP , SNP )
1622 C
1623       RETURN
1624       END
1625 C
1626 C
1627 C ----------------------------------------------------------------------
1628 C
1629       SUBROUTINE JGGGXX(W1,W2,W3,G, JW3W)
1630 C
1631 C this subroutine computes an off-shell w+, w-, w3, z or photon current
1632 C from the four-point gauge boson coupling, including the contributions
1633 C of w exchange diagrams.  the vector propagator is given in feynman
1634 C gauge for a photon and in unitary gauge for w and z bosons.  if one
1635 C sets wmass=0.0, then the ggg-->g current is given (see sect 2.9.1 of
1636 C the manual).
1637 C
1638 C input:
1639 C       complex w1(6)          : first  vector                        w1
1640 C       complex w2(6)          : second vector                        w2
1641 C       complex w3(6)          : third  vector                        w3
1642 C       real    g             : first  coupling constant
1643 C                                                  (see the table below)
1644 C
1645 C output:
1646 C       complex jw3w(6)        : w current             j^mu(w':w1,w2,w3)
1647 C
1648 #if defined(CERNLIB_IMPNONE)
1649       IMPLICIT NONE
1650 #endif
1651       COMPLEX*16  W1(6),W2(6),W3(6),JW3W(6)
1652       COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
1653      &           JJ(0:3),DV,W32,W13
1654       REAL*8     P1(0:3),P2(0:3),P3(0:3),Q(0:3),G,DG2,Q2
1655 C
1656       REAL*8 RXZERO
1657       PARAMETER( RXZERO=0.0D0 )
1658 C
1659       JW3W(5) = W1(5)+W2(5)+W3(5)
1660       JW3W(6) = W1(6)+W2(6)+W3(6)
1661 C
1662       DW1(0)=DCMPLX(W1(1))
1663       DW1(1)=DCMPLX(W1(2))
1664       DW1(2)=DCMPLX(W1(3))
1665       DW1(3)=DCMPLX(W1(4))
1666       DW2(0)=DCMPLX(W2(1))
1667       DW2(1)=DCMPLX(W2(2))
1668       DW2(2)=DCMPLX(W2(3))
1669       DW2(3)=DCMPLX(W2(4))
1670       DW3(0)=DCMPLX(W3(1))
1671       DW3(1)=DCMPLX(W3(2))
1672       DW3(2)=DCMPLX(W3(3))
1673       DW3(3)=DCMPLX(W3(4))
1674       P1(0)=DBLE(      W1(5))
1675       P1(1)=DBLE(      W1(6))
1676       P1(2)=DBLE(DIMAG(W1(6)))
1677       P1(3)=DBLE(DIMAG(W1(5)))
1678       P2(0)=DBLE(      W2(5))
1679       P2(1)=DBLE(      W2(6))
1680       P2(2)=DBLE(DIMAG(W2(6)))
1681       P2(3)=DBLE(DIMAG(W2(5)))
1682       P3(0)=DBLE(      W3(5))
1683       P3(1)=DBLE(      W3(6))
1684       P3(2)=DBLE(DIMAG(W3(6)))
1685       P3(3)=DBLE(DIMAG(W3(5)))
1686       Q(0)=-(P1(0)+P2(0)+P3(0))
1687       Q(1)=-(P1(1)+P2(1)+P3(1))
1688       Q(2)=-(P1(2)+P2(2)+P3(2))
1689       Q(3)=-(P1(3)+P2(3)+P3(3))
1690
1691       Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
1692
1693       DG2=DBLE(G)*DBLE(G)
1694 C
1695       DV = 1.0D0/DCMPLX( Q2 )
1696
1697 C  for the running width, use below instead of the above dv.
1698 C      dv = 1.0d0/dcmplx( q2 -mv2 , dmax1(dwv*q2/dmv,0.d0) )
1699 C
1700       W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
1701 C
1702 C
1703       W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
1704 C
1705       JJ(0)=DG2*( DW1(0)*W32 - DW2(0)*W13 )
1706       JJ(1)=DG2*( DW1(1)*W32 - DW2(1)*W13 )
1707       JJ(2)=DG2*( DW1(2)*W32 - DW2(2)*W13 )
1708       JJ(3)=DG2*( DW1(3)*W32 - DW2(3)*W13 )
1709 C
1710       JW3W(1) = DCMPLX( JJ(0)*DV )
1711       JW3W(2) = DCMPLX( JJ(1)*DV )
1712       JW3W(3) = DCMPLX( JJ(2)*DV )
1713       JW3W(4) = DCMPLX( JJ(3)*DV )
1714 C
1715       RETURN
1716       END
1717 C
1718 C ----------------------------------------------------------------------
1719 C
1720       SUBROUTINE JGGXXX(V1,V2,G, JVV)
1721 C
1722 C this subroutine computes an off-shell vector current from the three-
1723 C point gauge boson coupling.  the vector propagator is given in feynman
1724 C gauge for a massless vector and in unitary gauge for a massive vector.
1725 C
1726 C input:
1727 C       complex v1(6)          : first  vector                        v1
1728 C       complex v2(6)          : second vector                        v2
1729 C       real    g              : coupling constant (see the table below)
1730 C
1731 C output:
1732 C       complex jvv(6)         : vector current            j^mu(v:v1,v2)
1733 C
1734 #if defined(CERNLIB_IMPNONE)
1735       IMPLICIT NONE
1736 #endif
1737       COMPLEX*16 V1(6),V2(6),JVV(6),J12(0:3),
1738      &        SV1,SV2,V12
1739       REAL*8    P1(0:3),P2(0:3),Q(0:3),G,GS,S
1740 C
1741       REAL*8 RXZERO
1742       PARAMETER( RXZERO=0.0D0 )
1743 C
1744       JVV(5) = V1(5)+V2(5)
1745       JVV(6) = V1(6)+V2(6)
1746 C
1747       P1(0)=DBLE( V1(5))
1748       P1(1)=DBLE( V1(6))
1749       P1(2)=DIMAG(V1(6))
1750       P1(3)=DIMAG(V1(5))
1751       P2(0)=DBLE( V2(5))
1752       P2(1)=DBLE( V2(6))
1753       P2(2)=DIMAG(V2(6))
1754       P2(3)=DIMAG(V2(5))
1755       Q(0)=-DBLE( JVV(5))
1756       Q(1)=-DBLE( JVV(6))
1757       Q(2)=-DIMAG(JVV(6))
1758       Q(3)=-DIMAG(JVV(5))
1759       S=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1760 C
1761       V12=V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4)
1762       SV1= (P2(0)-Q(0))*V1(1) -(P2(1)-Q(1))*V1(2)
1763      &    -(P2(2)-Q(2))*V1(3) -(P2(3)-Q(3))*V1(4)
1764       SV2=-(P1(0)-Q(0))*V2(1) +(P1(1)-Q(1))*V2(2)
1765      &    +(P1(2)-Q(2))*V2(3) +(P1(3)-Q(3))*V2(4)
1766       J12(0)=(P1(0)-P2(0))*V12 +SV1*V2(1) +SV2*V1(1)
1767       J12(1)=(P1(1)-P2(1))*V12 +SV1*V2(2) +SV2*V1(2)
1768       J12(2)=(P1(2)-P2(2))*V12 +SV1*V2(3) +SV2*V1(3)
1769       J12(3)=(P1(3)-P2(3))*V12 +SV1*V2(4) +SV2*V1(4)
1770 C
1771       GS=-G/S
1772 C
1773       JVV(1) = GS*J12(0)
1774       JVV(2) = GS*J12(1)
1775       JVV(3) = GS*J12(2)
1776       JVV(4) = GS*J12(3)
1777 C
1778       RETURN
1779       END
1780 C
1781 C ----------------------------------------------------------------------
1782 C
1783       SUBROUTINE JIOXXX(FI,FO,G,VMASS,VWIDTH , JIO)
1784 C
1785 C this subroutine computes an off-shell vector current from an external
1786 C fermion pair.  the vector boson propagator is given in feynman gauge
1787 C for a massless vector and in unitary gauge for a massive vector.
1788 C
1789 C input:
1790 C       complex fi(6)          : flow-in  fermion                   |fi>
1791 C       complex fo(6)          : flow-out fermion                   <fo|
1792 C       real    g(2)           : coupling constants                  gvf
1793 C       real    vmass          : mass  of output vector v
1794 C       real    vwidth         : width of output vector v
1795 C
1796 C output:
1797 C       complex jio(6)         : vector current          j^mu(<fo|v|fi>)
1798 C
1799 #if defined(CERNLIB_IMPNONE)
1800       IMPLICIT NONE
1801 #endif
1802       COMPLEX*16 FI(6),FO(6),JIO(6),C0,C1,C2,C3,CS,D
1803       REAL*8    G(2),Q(0:3),VMASS,VWIDTH,Q2,VM2,DD
1804 C
1805       REAL*8 RXZERO, RXONE
1806       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1807       COMPLEX*16 CXIMAG
1808       LOGICAL FIRST
1809       SAVE CXIMAG,FIRST
1810       DATA FIRST/.TRUE./
1811 C
1812 C          Fix compilation with g77
1813       IF(FIRST) THEN
1814         FIRST=.FALSE.
1815         CXIMAG=DCMPLX( RXZERO, RXONE )
1816       ENDIF
1817 C
1818       JIO(5) = FO(5)-FI(5)
1819       JIO(6) = FO(6)-FI(6)
1820 C
1821       Q(0)=DBLE( JIO(5))
1822       Q(1)=DBLE( JIO(6))
1823       Q(2)=DIMAG(JIO(6))
1824       Q(3)=DIMAG(JIO(5))
1825       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1826       VM2=VMASS**2
1827 C
1828       IF (VMASS.NE.RXZERO) THEN
1829 C
1830          D=RXONE/DCMPLX( Q2-VM2 , MAX(SIGN( VMASS*VWIDTH ,Q2),RXZERO) )
1831 C  for the running width, use below instead of the above d.
1832 C      d=r_one/dcmplx( q2-vm2 , max( vwidth*q2/vmass ,r_zero) )
1833 C
1834          IF (G(2).NE.RXZERO) THEN
1835 C
1836             C0=  G(1)*( FO(3)*FI(1)+FO(4)*FI(2))
1837      &          +G(2)*( FO(1)*FI(3)+FO(2)*FI(4))
1838             C1= -G(1)*( FO(3)*FI(2)+FO(4)*FI(1))
1839      &          +G(2)*( FO(1)*FI(4)+FO(2)*FI(3))
1840             C2=( G(1)*( FO(3)*FI(2)-FO(4)*FI(1))
1841      &          +G(2)*(-FO(1)*FI(4)+FO(2)*FI(3)))*CXIMAG
1842             C3=  G(1)*(-FO(3)*FI(1)+FO(4)*FI(2))
1843      &          +G(2)*( FO(1)*FI(3)-FO(2)*FI(4))
1844          ELSE
1845 C
1846             D=D*G(1)
1847             C0=  FO(3)*FI(1)+FO(4)*FI(2)
1848             C1= -FO(3)*FI(2)-FO(4)*FI(1)
1849             C2=( FO(3)*FI(2)-FO(4)*FI(1))*CXIMAG
1850             C3= -FO(3)*FI(1)+FO(4)*FI(2)
1851          END IF
1852 C
1853          CS=(Q(0)*C0-Q(1)*C1-Q(2)*C2-Q(3)*C3)/VM2
1854 C
1855          JIO(1) = (C0-CS*Q(0))*D
1856          JIO(2) = (C1-CS*Q(1))*D
1857          JIO(3) = (C2-CS*Q(2))*D
1858          JIO(4) = (C3-CS*Q(3))*D
1859 C
1860       ELSE
1861          DD=RXONE/Q2
1862 C
1863          IF (G(2).NE.RXZERO) THEN
1864             JIO(1) = ( G(1)*( FO(3)*FI(1)+FO(4)*FI(2))
1865      &                +G(2)*( FO(1)*FI(3)+FO(2)*FI(4)) )*DD
1866             JIO(2) = (-G(1)*( FO(3)*FI(2)+FO(4)*FI(1))
1867      &                +G(2)*( FO(1)*FI(4)+FO(2)*FI(3)) )*DD
1868             JIO(3) = ( G(1)*( FO(3)*FI(2)-FO(4)*FI(1))
1869      &                +G(2)*(-FO(1)*FI(4)+FO(2)*FI(3)))
1870      $           *DCMPLX(RXZERO,DD)
1871             JIO(4) = ( G(1)*(-FO(3)*FI(1)+FO(4)*FI(2))
1872      &                +G(2)*( FO(1)*FI(3)-FO(2)*FI(4)) )*DD
1873 C
1874          ELSE
1875             DD=DD*G(1)
1876 C
1877             JIO(1) =  ( FO(3)*FI(1)+FO(4)*FI(2))*DD
1878             JIO(2) = -( FO(3)*FI(2)+FO(4)*FI(1))*DD
1879             JIO(3) =  ( FO(3)*FI(2)-FO(4)*FI(1))*DCMPLX(RXZERO,DD)
1880             JIO(4) =  (-FO(3)*FI(1)+FO(4)*FI(2))*DD
1881          END IF
1882       END IF
1883 C
1884       RETURN
1885       END
1886 C ----------------------------------------------------------------------
1887 C
1888       SUBROUTINE JSSXXX(S1,S2,G,VMASS,VWIDTH , JSS)
1889 C
1890 C This subroutine computes an off-shell vector current from the vector-
1891 C scalar-scalar coupling.  The coupling is absent in the minimal SM in
1892 C unitary gauge.  The propagator is given in Feynman gauge for a
1893 C massless vector and in unitary gauge for a massive vector.
1894 C
1895 C INPUT:
1896 C       complex S1(3)          : first  scalar                        S1
1897 C       complex S2(3)          : second scalar                        S2
1898 C       real    G              : coupling constant (S1 charge)
1899 C       real    VMASS          : mass  of OUTPUT vector V
1900 C       real    VWIDTH         : width of OUTPUT vector V
1901 C
1902 C Examples of the coupling constant G for SUSY particles are as follows:
1903 C   -----------------------------------------------------------
1904 C   |    S1    | (Q,I3) of S1  ||   V=A   |   V=Z   |   V=W   |
1905 C   -----------------------------------------------------------
1906 C   | nu~_L    | (  0  , +1/2) ||   ---   |  GZN(1) |  GWF(1) |
1907 C   | e~_L     | ( -1  , -1/2) ||  GAL(1) |  GZL(1) |  GWF(1) |
1908 C   | u~_L     | (+2/3 , +1/2) ||  GAU(1) |  GZU(1) |  GWF(1) |
1909 C   | d~_L     | (-1/3 , -1/2) ||  GAD(1) |  GZD(1) |  GWF(1) |
1910 C   -----------------------------------------------------------
1911 C   | e~_R-bar | ( +1  ,  0  ) || -GAL(2) | -GZL(2) | -GWF(2) |
1912 C   | u~_R-bar | (-2/3 ,  0  ) || -GAU(2) | -GZU(2) | -GWF(2) |
1913 C   | d~_R-bar | (+1/3 ,  0  ) || -GAD(2) | -GZD(2) | -GWF(2) |
1914 C   -----------------------------------------------------------
1915 C where the S1 charge is defined by the flowing-OUT quantum number.
1916 C
1917 C OUTPUT:
1918 C       complex JSS(6)         : vector current            J^mu(V:S1,S2)
1919 C
1920 #if defined(CERNLIB_IMPNONE)
1921       IMPLICIT NONE
1922 #endif
1923       COMPLEX*16 S1(3),S2(3),JSS(6),DG,ADG
1924       REAL*8  PP(0:3),PA(0:3),Q(0:3),G,VMASS,VWIDTH,Q2,VM2,MP2,MA2,M2D
1925 C
1926       JSS(5) = S1(2)+S2(2)
1927       JSS(6) = S1(3)+S2(3)
1928 C
1929       Q(0)=DBLE( JSS(5))
1930       Q(1)=DBLE( JSS(6))
1931       Q(2)=DIMAG(JSS(6))
1932       Q(3)=DIMAG(JSS(5))
1933       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1934       VM2=VMASS**2
1935 C
1936       IF (VMASS.EQ.0.) GOTO 10
1937 C
1938       DG=G/DCMPLX( Q2-VM2, MAX(SIGN( VMASS*VWIDTH ,Q2),0.D0))
1939 C  For the running width, use below instead of the above DG.
1940 C      DG=G/dCMPLX( Q2-VM2 , MAX( VWIDTH*Q2/VMASS ,0.) )
1941 C
1942       ADG=DG*S1(1)*S2(1)
1943 C
1944       PP(0)=DBLE( S1(2))
1945       PP(1)=DBLE( S1(3))
1946       PP(2)=DIMAG(S1(3))
1947       PP(3)=DIMAG(S1(2))
1948       PA(0)=DBLE( S2(2))
1949       PA(1)=DBLE( S2(3))
1950       PA(2)=DIMAG(S2(3))
1951       PA(3)=DIMAG(S2(2))
1952       MP2=PP(0)**2-(PP(1)**2+PP(2)**2+PP(3)**2)
1953       MA2=PA(0)**2-(PA(1)**2+PA(2)**2+PA(3)**2)
1954       M2D=MP2-MA2
1955 C
1956       JSS(1) = ADG*( (PP(0)-PA(0)) - Q(0)*M2D/VM2)
1957       JSS(2) = ADG*( (PP(1)-PA(1)) - Q(1)*M2D/VM2)
1958       JSS(3) = ADG*( (PP(2)-PA(2)) - Q(2)*M2D/VM2)
1959       JSS(4) = ADG*( (PP(3)-PA(3)) - Q(3)*M2D/VM2)
1960 C
1961       RETURN
1962 C
1963   10  ADG=G*S1(1)*S2(1)/Q2
1964 C
1965       JSS(1) = ADG*DBLE( S1(2)-S2(2))
1966       JSS(2) = ADG*DBLE( S1(3)-S2(3))
1967       JSS(3) = ADG*DIMAG(S1(3)-S2(3))
1968       JSS(4) = ADG*DIMAG(S1(2)-S2(2))
1969 C
1970       RETURN
1971       END
1972 C
1973 C
1974 C ----------------------------------------------------------------------
1975 C
1976       SUBROUTINE JTIOXX(FI,FO,G , JIO)
1977 C
1978 C this subroutine computes an off-shell vector current from an external
1979 C fermion pair.  the vector boson propagator is not included in this
1980 C routine.
1981 C
1982 C input:
1983 C       complex fi(6)          : flow-in  fermion                   |fi>
1984 C       complex fo(6)          : flow-out fermion                   <fo|
1985 C       real    g(2)           : coupling constants                  gvf
1986 C
1987 C output:
1988 C       complex jio(6)         : vector current          j^mu(<fo|v|fi>)
1989 C
1990 #if defined(CERNLIB_IMPNONE)
1991       IMPLICIT NONE
1992 #endif
1993       COMPLEX*16 FI(6),FO(6),JIO(6)
1994       REAL*8    G(2)
1995 C
1996       REAL*8 RXZERO, RXONE
1997       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1998       COMPLEX*16 CXIMAG
1999       LOGICAL FIRST
2000       SAVE CXIMAG,FIRST
2001       DATA FIRST/.TRUE./
2002 C
2003 C          Fix compilation with g77
2004       IF(FIRST) THEN
2005         FIRST=.FALSE.
2006         CXIMAG=DCMPLX( RXZERO, RXONE )
2007       ENDIF
2008 C
2009       JIO(5) = FO(5)-FI(5)
2010       JIO(6) = FO(6)-FI(6)
2011 C
2012       IF ( G(2) .NE. RXZERO ) THEN
2013          JIO(1) = ( G(1)*( FO(3)*FI(1)+FO(4)*FI(2))
2014      &             +G(2)*( FO(1)*FI(3)+FO(2)*FI(4)) )
2015          JIO(2) = (-G(1)*( FO(3)*FI(2)+FO(4)*FI(1))
2016      &             +G(2)*( FO(1)*FI(4)+FO(2)*FI(3)) )
2017          JIO(3) = ( G(1)*( FO(3)*FI(2)-FO(4)*FI(1))
2018      &             +G(2)*(-FO(1)*FI(4)+FO(2)*FI(3)) )*CXIMAG
2019          JIO(4) = ( G(1)*(-FO(3)*FI(1)+FO(4)*FI(2))
2020      &             +G(2)*( FO(1)*FI(3)-FO(2)*FI(4)) )
2021 C
2022       ELSE
2023          JIO(1) =  ( FO(3)*FI(1)+FO(4)*FI(2))*G(1)
2024          JIO(2) = -( FO(3)*FI(2)+FO(4)*FI(1))*G(1)
2025          JIO(3) =  ( FO(3)*FI(2)-FO(4)*FI(1))*DCMPLX(RXZERO,G(1))
2026          JIO(4) =  (-FO(3)*FI(1)+FO(4)*FI(2))*G(1)
2027       END IF
2028 C
2029       RETURN
2030       END
2031 C ----------------------------------------------------------------------
2032 C
2033       SUBROUTINE JVSSXX(VC,S1,S2,G,VMASS,VWIDTH , JVSS)
2034 C
2035 C This subroutine computes an off-shell vector current from the vector-
2036 C vector-scalar-scalar coupling.  The vector propagator is given in
2037 C Feynman gauge for a massless vector and in unitary gauge for a massive
2038 C vector.
2039 C
2040 C INPUT:
2041 C       complex VC(6)          : input  vector                        V
2042 C       complex S1(3)          : first  scalar                        S1
2043 C       complex S2(3)          : second scalar                        S2
2044 C       real    G              : coupling constant                 GVVHH
2045 C       real    VMASS          : mass  of OUTPUT vector V'
2046 C       real    VWIDTH         : width of OUTPUT vector V'
2047 C
2048 C OUTPUT:
2049 C       complex JVSS(6)        : vector current         J^mu(V':V,S1,S2)
2050 C
2051 #if defined(CERNLIB_IMPNONE)
2052       IMPLICIT NONE
2053 #endif
2054       COMPLEX*16 VC(6),S1(3),S2(3),JVSS(6),DG
2055       REAL*8    Q(0:3),G,VMASS,VWIDTH,Q2,VK,VM2
2056 C
2057       JVSS(5) = VC(5)+S1(2)+S2(2)
2058       JVSS(6) = VC(6)+S1(3)+S2(3)
2059 C
2060       Q(0)=DBLE( JVSS(5))
2061       Q(1)=DBLE( JVSS(6))
2062       Q(2)=DIMAG(JVSS(6))
2063       Q(3)=DIMAG(JVSS(5))
2064       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
2065       VM2=VMASS**2
2066 C
2067       IF (VMASS.EQ.0.) GOTO 10
2068 C
2069       DG=G*S1(1)*S2(1)/DCMPLX( Q2-VM2,MAX(SIGN( VMASS*VWIDTH,Q2),0.D0))
2070 C  For the running width, use below instead of the above DG.
2071 C      DG=G*S1(1)*S2(1)/CMPLX( Q2-VM2 , MAX( VWIDTH*Q2/VMASS ,0.))
2072 C
2073       VK=(Q(0)*VC(1)-Q(1)*VC(2)-Q(2)*VC(3)-Q(3)*VC(4))/VM2
2074 C
2075       JVSS(1) = DG*(VC(1)-VK*Q(0))
2076       JVSS(2) = DG*(VC(2)-VK*Q(1))
2077       JVSS(3) = DG*(VC(3)-VK*Q(2))
2078       JVSS(4) = DG*(VC(4)-VK*Q(3))
2079 C
2080       RETURN
2081 C
2082   10  DG= G*S1(1)*S2(1)/Q2
2083 C
2084       JVSS(1) = DG*VC(1)
2085       JVSS(2) = DG*VC(2)
2086       JVSS(3) = DG*VC(3)
2087       JVSS(4) = DG*VC(4)
2088 C
2089       RETURN
2090       END
2091 C
2092 C
2093 C ----------------------------------------------------------------------
2094 C
2095       SUBROUTINE JVSXXX(VC,SC,G,VMASS,VWIDTH , JVS)
2096 C
2097 C this subroutine computes an off-shell vector current from the vector-
2098 C vector-scalar coupling.  the vector propagator is given in feynman
2099 C gauge for a massless vector and in unitary gauge for a massive vector.
2100 C
2101 C input:
2102 C       complex vc(6)          : input vector                          v
2103 C       complex sc(3)          : input scalar                          s
2104 C       real    g              : coupling constant                  gvvh
2105 C       real    vmass          : mass  of output vector v'
2106 C       real    vwidth         : width of output vector v'
2107 C
2108 C output:
2109 C       complex jvs(6)         : vector current             j^mu(v':v,s)
2110 C
2111 #if defined(CERNLIB_IMPNONE)
2112       IMPLICIT NONE
2113 #endif
2114       COMPLEX*16 VC(6),SC(3),JVS(6),DG,VK
2115       REAL*8    Q(0:3),VMASS,VWIDTH,Q2,VM2,G
2116 C
2117       JVS(5) = VC(5)+SC(2)
2118       JVS(6) = VC(6)+SC(3)
2119 C
2120       Q(0)=DBLE( JVS(5))
2121       Q(1)=DBLE( JVS(6))
2122       Q(2)=DIMAG(JVS(6))
2123       Q(3)=DIMAG(JVS(5))
2124       Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
2125       VM2=VMASS**2
2126 C
2127       IF (VMASS.EQ.0.) GOTO 10
2128 C
2129       DG=G*SC(1)/DCMPLX( Q2-VM2 , MAX(DSIGN( VMASS*VWIDTH ,Q2),0.D0) )
2130 C  for the running width, use below instead of the above dg.
2131 C      dg=g*sc(1)/dcmplx( q2-vm2 , max( vwidth*q2/vmass ,0.) )
2132 C
2133       VK=(-Q(0)*VC(1)+Q(1)*VC(2)+Q(2)*VC(3)+Q(3)*VC(4))/VM2
2134 C
2135       JVS(1) = DG*(Q(0)*VK+VC(1))
2136       JVS(2) = DG*(Q(1)*VK+VC(2))
2137       JVS(3) = DG*(Q(2)*VK+VC(3))
2138       JVS(4) = DG*(Q(3)*VK+VC(4))
2139 C
2140       RETURN
2141 C
2142   10  DG=G*SC(1)/Q2
2143 C
2144       JVS(1) = DG*VC(1)
2145       JVS(2) = DG*VC(2)
2146       JVS(3) = DG*VC(3)
2147       JVS(4) = DG*VC(4)
2148 C
2149       RETURN
2150       END
2151
2152
2153 C
2154 C ----------------------------------------------------------------------
2155 C
2156       SUBROUTINE JVVXXX(V1,V2,G,VMASS,VWIDTH , JVV)
2157 C
2158 C this subroutine computes an off-shell vector current from the three-
2159 C point gauge boson coupling.  the vector propagator is given in feynman
2160 C gauge for a massless vector and in unitary gauge for a massive vector.
2161 C
2162 C input:
2163 C       complex v1(6)          : first  vector                        v1
2164 C       complex v2(6)          : second vector                        v2
2165 C       real    g              : coupling constant (see the table below)
2166 C       real    vmass          : mass  of output vector v
2167 C       real    vwidth         : width of output vector v
2168 C
2169 C the possible sets of the inputs are as follows:
2170 C    ------------------------------------------------------------------
2171 C    |   v1   |   v2   |  jvv   |      g       |   vmass  |  vwidth   |
2172 C    ------------------------------------------------------------------
2173 C    |   w-   |   w+   |  a/z   |  gwwa/gwwz   | 0./zmass | 0./zwidth |
2174 C    | w3/a/z |   w-   |  w+    | gw/gwwa/gwwz |   wmass  |  wwidth   |
2175 C    |   w+   | w3/a/z |  w-    | gw/gwwa/gwwz |   wmass  |  wwidth   |
2176 C    ------------------------------------------------------------------
2177 C where all the bosons are defined by the flowing-out quantum number.
2178 C
2179 C output:
2180 C       complex jvv(6)         : vector current            j^mu(v:v1,v2)
2181 C
2182 #if defined(CERNLIB_IMPNONE)
2183       IMPLICIT NONE
2184 #endif
2185       COMPLEX*16 V1(6),V2(6),JVV(6),J12(0:3),JS,DG,
2186      &        SV1,SV2,S11,S12,S21,S22,V12
2187       REAL*8    P1(0:3),P2(0:3),Q(0:3),G,VMASS,VWIDTH,GS,S,VM2,M1,M2
2188 C
2189       REAL*8 RXZERO
2190       PARAMETER( RXZERO=0.0D0 )
2191 C
2192       JVV(5) = V1(5)+V2(5)
2193       JVV(6) = V1(6)+V2(6)
2194 C
2195       P1(0)=DBLE( V1(5))
2196       P1(1)=DBLE( V1(6))
2197       P1(2)=DIMAG(V1(6))
2198       P1(3)=DIMAG(V1(5))
2199       P2(0)=DBLE( V2(5))
2200       P2(1)=DBLE( V2(6))
2201       P2(2)=DIMAG(V2(6))
2202       P2(3)=DIMAG(V2(5))
2203       Q(0)=-DBLE( JVV(5))
2204       Q(1)=-DBLE( JVV(6))
2205       Q(2)=-DIMAG(JVV(6))
2206       Q(3)=-DIMAG(JVV(5))
2207       S=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
2208 C
2209       V12=V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4)
2210       SV1= (P2(0)-Q(0))*V1(1) -(P2(1)-Q(1))*V1(2)
2211      &    -(P2(2)-Q(2))*V1(3) -(P2(3)-Q(3))*V1(4)
2212       SV2=-(P1(0)-Q(0))*V2(1) +(P1(1)-Q(1))*V2(2)
2213      &    +(P1(2)-Q(2))*V2(3) +(P1(3)-Q(3))*V2(4)
2214       J12(0)=(P1(0)-P2(0))*V12 +SV1*V2(1) +SV2*V1(1)
2215       J12(1)=(P1(1)-P2(1))*V12 +SV1*V2(2) +SV2*V1(2)
2216       J12(2)=(P1(2)-P2(2))*V12 +SV1*V2(3) +SV2*V1(3)
2217       J12(3)=(P1(3)-P2(3))*V12 +SV1*V2(4) +SV2*V1(4)
2218 C
2219       IF ( VMASS .NE. RXZERO ) THEN
2220          VM2=VMASS**2
2221          M1=P1(0)**2-(P1(1)**2+P1(2)**2+P1(3)**2)
2222          M2=P2(0)**2-(P2(1)**2+P2(2)**2+P2(3)**2)
2223          S11=P1(0)*V1(1)-P1(1)*V1(2)-P1(2)*V1(3)-P1(3)*V1(4)
2224          S12=P1(0)*V2(1)-P1(1)*V2(2)-P1(2)*V2(3)-P1(3)*V2(4)
2225          S21=P2(0)*V1(1)-P2(1)*V1(2)-P2(2)*V1(3)-P2(3)*V1(4)
2226          S22=P2(0)*V2(1)-P2(1)*V2(2)-P2(2)*V2(3)-P2(3)*V2(4)
2227          JS=(V12*(-M1+M2) +S11*S12 -S21*S22)/VM2
2228 C
2229          DG=-G/DCMPLX( S-VM2 , MAX(SIGN( VMASS*VWIDTH ,S),RXZERO) )
2230 C
2231 C  for the running width, use below instead of the above dg.
2232 C         dg=-g/dcmplx( s-vm2 , max( vwidth*s/vmass ,r_zero) )
2233 C
2234          JVV(1) = DG*(J12(0)-Q(0)*JS)
2235          JVV(2) = DG*(J12(1)-Q(1)*JS)
2236          JVV(3) = DG*(J12(2)-Q(2)*JS)
2237          JVV(4) = DG*(J12(3)-Q(3)*JS)
2238 C
2239       ELSE
2240          GS=-G/S
2241 C
2242          JVV(1) = GS*J12(0)
2243          JVV(2) = GS*J12(1)
2244          JVV(3) = GS*J12(2)
2245          JVV(4) = GS*J12(3)
2246       END IF
2247 C
2248       RETURN
2249       END
2250 C
2251 C ----------------------------------------------------------------------
2252 C
2253       SUBROUTINE JW3WXX(W1,W2,W3,G1,G2,WMASS,WWIDTH,VMASS,VWIDTH , JW3W)
2254 C
2255 C this subroutine computes an off-shell w+, w-, w3, z or photon current
2256 C from the four-point gauge boson coupling, including the contributions
2257 C of w exchange diagrams.  the vector propagator is given in feynman
2258 C gauge for a photon and in unitary gauge for w and z bosons.  if one
2259 C sets wmass=0.0, then the ggg-->g current is given (see sect 2.9.1 of
2260 C the manual).
2261 C
2262 C input:
2263 C       complex w1(6)          : first  vector                        w1
2264 C       complex w2(6)          : second vector                        w2
2265 C       complex w3(6)          : third  vector                        w3
2266 C       real    g1             : first  coupling constant
2267 C       real    g2             : second coupling constant
2268 C                                                  (see the table below)
2269 C       real    wmass          : mass  of internal w
2270 C       real    wwidth         : width of internal w
2271 C       real    vmass          : mass  of output w'
2272 C       real    vwidth         : width of output w'
2273 C
2274 C the possible sets of the inputs are as follows:
2275 C   -------------------------------------------------------------------
2276 C   |  w1  |  w2  |  w3  | g1 | g2 |wmass|wwidth|vmass|vwidth || jw3w |
2277 C   -------------------------------------------------------------------
2278 C   |  w-  |  w3  |  w+  | gw |gwwz|wmass|wwidth|zmass|zwidth ||  z   |
2279 C   |  w-  |  w3  |  w+  | gw |gwwa|wmass|wwidth|  0. |  0.   ||  a   |
2280 C   |  w-  |  z   |  w+  |gwwz|gwwz|wmass|wwidth|zmass|zwidth ||  z   |
2281 C   |  w-  |  z   |  w+  |gwwz|gwwa|wmass|wwidth|  0. |  0.   ||  a   |
2282 C   |  w-  |  a   |  w+  |gwwa|gwwz|wmass|wwidth|zmass|zwidth ||  z   |
2283 C   |  w-  |  a   |  w+  |gwwa|gwwa|wmass|wwidth|  0. |  0.   ||  a   |
2284 C   -------------------------------------------------------------------
2285 C   |  w3  |  w-  |  w3  | gw | gw |wmass|wwidth|wmass|wwidth ||  w+  |
2286 C   |  w3  |  w+  |  w3  | gw | gw |wmass|wwidth|wmass|wwidth ||  w-  |
2287 C   |  w3  |  w-  |  z   | gw |gwwz|wmass|wwidth|wmass|wwidth ||  w+  |
2288 C   |  w3  |  w+  |  z   | gw |gwwz|wmass|wwidth|wmass|wwidth ||  w-  |
2289 C   |  w3  |  w-  |  a   | gw |gwwa|wmass|wwidth|wmass|wwidth ||  w+  |
2290 C   |  w3  |  w+  |  a   | gw |gwwa|wmass|wwidth|wmass|wwidth ||  w-  |
2291 C   |  z   |  w-  |  z   |gwwz|gwwz|wmass|wwidth|wmass|wwidth ||  w+  |
2292 C   |  z   |  w+  |  z   |gwwz|gwwz|wmass|wwidth|wmass|wwidth ||  w-  |
2293 C   |  z   |  w-  |  a   |gwwz|gwwa|wmass|wwidth|wmass|wwidth ||  w+  |
2294 C   |  z   |  w+  |  a   |gwwz|gwwa|wmass|wwidth|wmass|wwidth ||  w-  |
2295 C   |  a   |  w-  |  a   |gwwa|gwwa|wmass|wwidth|wmass|wwidth ||  w+  |
2296 C   |  a   |  w+  |  a   |gwwa|gwwa|wmass|wwidth|wmass|wwidth ||  w-  |
2297 C   -------------------------------------------------------------------
2298 C where all the bosons are defined by the flowing-out quantum number.
2299 C
2300 C output:
2301 C       complex jw3w(6)        : w current             j^mu(w':w1,w2,w3)
2302 C
2303 #if defined(CERNLIB_IMPNONE)
2304       IMPLICIT NONE
2305 #endif
2306       COMPLEX*16  W1(6),W2(6),W3(6),JW3W(6)
2307       COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
2308      &           JJ(0:3),J4(0:3),
2309      &           DV,W12,W32,W13,
2310      &           JQ
2311       REAL*8     G1,G2,WMASS,WWIDTH,VMASS,VWIDTH
2312       REAL*8     P1(0:3),P2(0:3),P3(0:3),Q(0:3),
2313      &           DG2,DMV,DWV,MV2,Q2
2314 C
2315       REAL*8 RXZERO
2316       PARAMETER( RXZERO=0.0D0 )
2317 C
2318       JW3W(5) = W1(5)+W2(5)+W3(5)
2319       JW3W(6) = W1(6)+W2(6)+W3(6)
2320 C
2321       DW1(0)=DCMPLX(W1(1))
2322       DW1(1)=DCMPLX(W1(2))
2323       DW1(2)=DCMPLX(W1(3))
2324       DW1(3)=DCMPLX(W1(4))
2325       DW2(0)=DCMPLX(W2(1))
2326       DW2(1)=DCMPLX(W2(2))
2327       DW2(2)=DCMPLX(W2(3))
2328       DW2(3)=DCMPLX(W2(4))
2329       DW3(0)=DCMPLX(W3(1))
2330       DW3(1)=DCMPLX(W3(2))
2331       DW3(2)=DCMPLX(W3(3))
2332       DW3(3)=DCMPLX(W3(4))
2333       P1(0)=DBLE(      W1(5))
2334       P1(1)=DBLE(      W1(6))
2335       P1(2)=DBLE(DIMAG(W1(6)))
2336       P1(3)=DBLE(DIMAG(W1(5)))
2337       P2(0)=DBLE(      W2(5))
2338       P2(1)=DBLE(      W2(6))
2339       P2(2)=DBLE(DIMAG(W2(6)))
2340       P2(3)=DBLE(DIMAG(W2(5)))
2341       P3(0)=DBLE(      W3(5))
2342       P3(1)=DBLE(      W3(6))
2343       P3(2)=DBLE(DIMAG(W3(6)))
2344       P3(3)=DBLE(DIMAG(W3(5)))
2345       Q(0)=-(P1(0)+P2(0)+P3(0))
2346       Q(1)=-(P1(1)+P2(1)+P3(1))
2347       Q(2)=-(P1(2)+P2(2)+P3(2))
2348       Q(3)=-(P1(3)+P2(3)+P3(3))
2349
2350
2351       Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
2352       DG2=DBLE(G1)*DBLE(G2)
2353       DMV=DBLE(VMASS)
2354       DWV=DBLE(VWIDTH)
2355       MV2=DMV**2
2356       IF (VMASS.EQ. RXZERO) THEN
2357       DV = 1.0D0/DCMPLX( Q2 )
2358       ELSE
2359       DV = 1.0D0/DCMPLX( Q2 -MV2 , DMAX1(DSIGN(DMV*DWV,Q2 ),0.D0) )
2360       ENDIF
2361 C  for the running width, use below instead of the above dv.
2362 C      dv = 1.0d0/dcmplx( q2 -mv2 , dmax1(dwv*q2/dmv,0.d0) )
2363 C
2364       W12=DW1(0)*DW2(0)-DW1(1)*DW2(1)-DW1(2)*DW2(2)-DW1(3)*DW2(3)
2365       W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
2366 C
2367       IF ( WMASS .NE. RXZERO ) THEN
2368          W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
2369 C
2370          J4(0)=DG2*( DW1(0)*W32 + DW3(0)*W12 - 2.D0*DW2(0)*W13 )
2371          J4(1)=DG2*( DW1(1)*W32 + DW3(1)*W12 - 2.D0*DW2(1)*W13 )
2372          J4(2)=DG2*( DW1(2)*W32 + DW3(2)*W12 - 2.D0*DW2(2)*W13 )
2373          J4(3)=DG2*( DW1(3)*W32 + DW3(3)*W12 - 2.D0*DW2(3)*W13 )
2374 C
2375          JJ(0)=J4(0)
2376          JJ(1)=J4(1)
2377          JJ(2)=J4(2)
2378          JJ(3)=J4(3)
2379
2380       ELSE
2381 C
2382          W12=DW1(0)*DW2(0)-DW1(1)*DW2(1)-DW1(2)*DW2(2)-DW1(3)*DW2(3)
2383          W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
2384          W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
2385 C
2386          J4(0)=DG2*( DW1(0)*W32 - DW2(0)*W13 )
2387          J4(1)=DG2*( DW1(1)*W32 - DW2(1)*W13 )
2388          J4(2)=DG2*( DW1(2)*W32 - DW2(2)*W13 )
2389          J4(3)=DG2*( DW1(3)*W32 - DW2(3)*W13 )
2390 C
2391          JJ(0)=J4(0)
2392          JJ(1)=J4(1)
2393          JJ(2)=J4(2)
2394          JJ(3)=J4(3)
2395
2396       END IF
2397 C
2398       IF ( VMASS .NE. RXZERO ) THEN
2399 C
2400          JQ=(JJ(0)*Q(0)-JJ(1)*Q(1)-JJ(2)*Q(2)-JJ(3)*Q(3))/MV2
2401 C
2402          JW3W(1) = DCMPLX( (JJ(0)-JQ*Q(0))*DV )
2403          JW3W(2) = DCMPLX( (JJ(1)-JQ*Q(1))*DV )
2404          JW3W(3) = DCMPLX( (JJ(2)-JQ*Q(2))*DV )
2405          JW3W(4) = DCMPLX( (JJ(3)-JQ*Q(3))*DV )
2406 C
2407       ELSE
2408 C
2409          JW3W(1) = DCMPLX( JJ(0)*DV )
2410          JW3W(2) = DCMPLX( JJ(1)*DV )
2411          JW3W(3) = DCMPLX( JJ(2)*DV )
2412          JW3W(4) = DCMPLX( JJ(3)*DV )
2413       END IF
2414 C
2415       RETURN
2416       END
2417 C
2418 C ----------------------------------------------------------------------
2419 C
2420       SUBROUTINE JWWWXX(W1,W2,W3,GWWA,GWWZ,ZMASS,ZWIDTH,WMASS,WWIDTH ,
2421      &                  JWWW)
2422 C
2423 C this subroutine computes an off-shell w+/w- current from the four-
2424 C point gauge boson coupling, including the contributions of photon and
2425 C z exchanges.  the vector propagators for the output w and the internal
2426 C z bosons are given in unitary gauge, and that of the internal photon
2427 C is given in feynman gauge.
2428 C
2429 C input:
2430 C       complex w1(6)          : first  vector                        w1
2431 C       complex w2(6)          : second vector                        w2
2432 C       complex w3(6)          : third  vector                        w3
2433 C       real    gwwa           : coupling constant of w and a       gwwa
2434 C       real    gwwz           : coupling constant of w and z       gwwz
2435 C       real    zmass          : mass  of internal z
2436 C       real    zwidth         : width of internal z
2437 C       real    wmass          : mass  of output w
2438 C       real    wwidth         : width of output w
2439 C
2440 C the possible sets of the inputs are as follows:
2441 C   -------------------------------------------------------------------
2442 C   |  w1  |  w2  |  w3  |gwwa|gwwz|zmass|zwidth|wmass|wwidth || jwww |
2443 C   -------------------------------------------------------------------
2444 C   |  w-  |  w+  |  w-  |gwwa|gwwz|zmass|zwidth|wmass|wwidth ||  w+  |
2445 C   |  w+  |  w-  |  w+  |gwwa|gwwz|zmass|zwidth|wmass|wwidth ||  w-  |
2446 C   -------------------------------------------------------------------
2447 C where all the bosons are defined by the flowing-out quantum number.
2448 C
2449 C output:
2450 C       complex jwww(6)        : w current             j^mu(w':w1,w2,w3)
2451 C
2452 #if defined(CERNLIB_IMPNONE)
2453       IMPLICIT NONE
2454 #endif
2455       COMPLEX*16  W1(6),W2(6),W3(6),JWWW(6)
2456       COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
2457      &           JJ(0:3),JS(0:3),JT(0:3),J4(0:3),
2458      &           JT12(0:3),JT32(0:3),J12(0:3),J32(0:3),
2459      &           DZS,DZT,DW,W12,W32,W13,P1W2,P2W1,P3W2,P2W3,
2460      &           JK12,JK32,JSW3,JTW1,P3JS,KSW3,P1JT,KTW1,JQ
2461       REAL*8     GWWA,GWWZ,ZMASS,ZWIDTH,WMASS,WWIDTH
2462       REAL*8     P1(0:3),P2(0:3),P3(0:3),Q(0:3),KS(0:3),KT(0:3),
2463      &           DGWWA2,DGWWZ2,DGW2,DMZ,DWZ,DMW,DWW,MZ2,MW2,Q2,KS2,KT2,
2464      &           DAS,DAT
2465 C
2466       JWWW(5) = W1(5)+W2(5)+W3(5)
2467       JWWW(6) = W1(6)+W2(6)+W3(6)
2468 C
2469       DW1(0)=DCMPLX(W1(1))
2470       DW1(1)=DCMPLX(W1(2))
2471       DW1(2)=DCMPLX(W1(3))
2472       DW1(3)=DCMPLX(W1(4))
2473       DW2(0)=DCMPLX(W2(1))
2474       DW2(1)=DCMPLX(W2(2))
2475       DW2(2)=DCMPLX(W2(3))
2476       DW2(3)=DCMPLX(W2(4))
2477       DW3(0)=DCMPLX(W3(1))
2478       DW3(1)=DCMPLX(W3(2))
2479       DW3(2)=DCMPLX(W3(3))
2480       DW3(3)=DCMPLX(W3(4))
2481       P1(0)=DBLE(      W1(5))
2482       P1(1)=DBLE(      W1(6))
2483       P1(2)=DBLE(DIMAG(W1(6)))
2484       P1(3)=DBLE(DIMAG(W1(5)))
2485       P2(0)=DBLE(      W2(5))
2486       P2(1)=DBLE(      W2(6))
2487       P2(2)=DBLE(DIMAG(W2(6)))
2488       P2(3)=DBLE(DIMAG(W2(5)))
2489       P3(0)=DBLE(      W3(5))
2490       P3(1)=DBLE(      W3(6))
2491       P3(2)=DBLE(DIMAG(W3(6)))
2492       P3(3)=DBLE(DIMAG(W3(5)))
2493       Q(0)=-(P1(0)+P2(0)+P3(0))
2494       Q(1)=-(P1(1)+P2(1)+P3(1))
2495       Q(2)=-(P1(2)+P2(2)+P3(2))
2496       Q(3)=-(P1(3)+P2(3)+P3(3))
2497       KS(0)=P1(0)+P2(0)
2498       KS(1)=P1(1)+P2(1)
2499       KS(2)=P1(2)+P2(2)
2500       KS(3)=P1(3)+P2(3)
2501       KT(0)=P2(0)+P3(0)
2502       KT(1)=P2(1)+P3(1)
2503       KT(2)=P2(2)+P3(2)
2504       KT(3)=P2(3)+P3(3)
2505       Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
2506       KS2=KS(0)**2-(KS(1)**2+KS(2)**2+KS(3)**2)
2507       KT2=KT(0)**2-(KT(1)**2+KT(2)**2+KT(3)**2)
2508       DGWWA2=DBLE(GWWA)**2
2509       DGWWZ2=DBLE(GWWZ)**2
2510       DGW2  =DGWWA2+DGWWZ2
2511       DMZ=DBLE(ZMASS)
2512       DWZ=DBLE(ZWIDTH)
2513       DMW=DBLE(WMASS)
2514       DWW=DBLE(WWIDTH)
2515       MZ2=DMZ**2
2516       MW2=DMW**2
2517 C
2518       DAS=-DGWWA2/KS2
2519       DAT=-DGWWA2/KT2
2520       DZS=-DGWWZ2/DCMPLX( KS2-MZ2 , DMAX1(DSIGN(DMZ*DWZ,KS2),0.D0) )
2521       DZT=-DGWWZ2/DCMPLX( KT2-MZ2 , DMAX1(DSIGN(DMZ*DWZ,KT2),0.D0) )
2522       DW =-1.0D0/DCMPLX( Q2 -MW2 , DMAX1(DSIGN(DMW*DWW,Q2 ),0.D0) )
2523 C  for the running width, use below instead of the above dw.
2524 C      dw =-1.0d0/dcmplx( q2 -mw2 , dmax1(dww*q2/dmw,0.d0) )
2525 C
2526       W12=DW1(0)*DW2(0)-DW1(1)*DW2(1)-DW1(2)*DW2(2)-DW1(3)*DW2(3)
2527       W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
2528 C
2529       P1W2= (P1(0)+KS(0))*DW2(0)-(P1(1)+KS(1))*DW2(1)
2530      &     -(P1(2)+KS(2))*DW2(2)-(P1(3)+KS(3))*DW2(3)
2531       P2W1= (P2(0)+KS(0))*DW1(0)-(P2(1)+KS(1))*DW1(1)
2532      &     -(P2(2)+KS(2))*DW1(2)-(P2(3)+KS(3))*DW1(3)
2533       P3W2= (P3(0)+KT(0))*DW2(0)-(P3(1)+KT(1))*DW2(1)
2534      &     -(P3(2)+KT(2))*DW2(2)-(P3(3)+KT(3))*DW2(3)
2535       P2W3= (P2(0)+KT(0))*DW3(0)-(P2(1)+KT(1))*DW3(1)
2536      &     -(P2(2)+KT(2))*DW3(2)-(P2(3)+KT(3))*DW3(3)
2537 C
2538       JT12(0)= (P1(0)-P2(0))*W12 + P2W1*DW2(0) - P1W2*DW1(0)
2539       JT12(1)= (P1(1)-P2(1))*W12 + P2W1*DW2(1) - P1W2*DW1(1)
2540       JT12(2)= (P1(2)-P2(2))*W12 + P2W1*DW2(2) - P1W2*DW1(2)
2541       JT12(3)= (P1(3)-P2(3))*W12 + P2W1*DW2(3) - P1W2*DW1(3)
2542       JT32(0)= (P3(0)-P2(0))*W32 + P2W3*DW2(0) - P3W2*DW3(0)
2543       JT32(1)= (P3(1)-P2(1))*W32 + P2W3*DW2(1) - P3W2*DW3(1)
2544       JT32(2)= (P3(2)-P2(2))*W32 + P2W3*DW2(2) - P3W2*DW3(2)
2545       JT32(3)= (P3(3)-P2(3))*W32 + P2W3*DW2(3) - P3W2*DW3(3)
2546 C
2547       JK12=(JT12(0)*KS(0)-JT12(1)*KS(1)-JT12(2)*KS(2)-JT12(3)*KS(3))/MZ2
2548       JK32=(JT32(0)*KT(0)-JT32(1)*KT(1)-JT32(2)*KT(2)-JT32(3)*KT(3))/MZ2
2549 C
2550       J12(0)=JT12(0)*(DAS+DZS)-KS(0)*JK12*DZS
2551       J12(1)=JT12(1)*(DAS+DZS)-KS(1)*JK12*DZS
2552       J12(2)=JT12(2)*(DAS+DZS)-KS(2)*JK12*DZS
2553       J12(3)=JT12(3)*(DAS+DZS)-KS(3)*JK12*DZS
2554       J32(0)=JT32(0)*(DAT+DZT)-KT(0)*JK32*DZT
2555       J32(1)=JT32(1)*(DAT+DZT)-KT(1)*JK32*DZT
2556       J32(2)=JT32(2)*(DAT+DZT)-KT(2)*JK32*DZT
2557       J32(3)=JT32(3)*(DAT+DZT)-KT(3)*JK32*DZT
2558 C
2559       JSW3=J12(0)*DW3(0)-J12(1)*DW3(1)-J12(2)*DW3(2)-J12(3)*DW3(3)
2560       JTW1=J32(0)*DW1(0)-J32(1)*DW1(1)-J32(2)*DW1(2)-J32(3)*DW1(3)
2561 C
2562       P3JS= (P3(0)-Q(0))*J12(0)-(P3(1)-Q(1))*J12(1)
2563      &     -(P3(2)-Q(2))*J12(2)-(P3(3)-Q(3))*J12(3)
2564       KSW3= (KS(0)-Q(0))*DW3(0)-(KS(1)-Q(1))*DW3(1)
2565      &     -(KS(2)-Q(2))*DW3(2)-(KS(3)-Q(3))*DW3(3)
2566       P1JT= (P1(0)-Q(0))*J32(0)-(P1(1)-Q(1))*J32(1)
2567      &     -(P1(2)-Q(2))*J32(2)-(P1(3)-Q(3))*J32(3)
2568       KTW1= (KT(0)-Q(0))*DW1(0)-(KT(1)-Q(1))*DW1(1)
2569      &     -(KT(2)-Q(2))*DW1(2)-(KT(3)-Q(3))*DW1(3)
2570 C
2571       JS(0)= (KS(0)-P3(0))*JSW3 + P3JS*DW3(0) - KSW3*J12(0)
2572       JS(1)= (KS(1)-P3(1))*JSW3 + P3JS*DW3(1) - KSW3*J12(1)
2573       JS(2)= (KS(2)-P3(2))*JSW3 + P3JS*DW3(2) - KSW3*J12(2)
2574       JS(3)= (KS(3)-P3(3))*JSW3 + P3JS*DW3(3) - KSW3*J12(3)
2575       JT(0)= (KT(0)-P1(0))*JTW1 + P1JT*DW1(0) - KTW1*J32(0)
2576       JT(1)= (KT(1)-P1(1))*JTW1 + P1JT*DW1(1) - KTW1*J32(1)
2577       JT(2)= (KT(2)-P1(2))*JTW1 + P1JT*DW1(2) - KTW1*J32(2)
2578       JT(3)= (KT(3)-P1(3))*JTW1 + P1JT*DW1(3) - KTW1*J32(3)
2579 C
2580       W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
2581 C
2582       J4(0)=DGW2*( DW1(0)*W32 + DW3(0)*W12 - 2.D0*DW2(0)*W13 )
2583       J4(1)=DGW2*( DW1(1)*W32 + DW3(1)*W12 - 2.D0*DW2(1)*W13 )
2584       J4(2)=DGW2*( DW1(2)*W32 + DW3(2)*W12 - 2.D0*DW2(2)*W13 )
2585       J4(3)=DGW2*( DW1(3)*W32 + DW3(3)*W12 - 2.D0*DW2(3)*W13 )
2586 C
2587 C      jj(0)=js(0)+jt(0)+j4(0)
2588 C      jj(1)=js(1)+jt(1)+j4(1)
2589 C      jj(2)=js(2)+jt(2)+j4(2)
2590 C      jj(3)=js(3)+jt(3)+j4(3)
2591
2592       JJ(0)=J4(0)
2593       JJ(1)=J4(1)
2594       JJ(2)=J4(2)
2595       JJ(3)=J4(3)
2596 C
2597       JQ=(JJ(0)*Q(0)-JJ(1)*Q(1)-JJ(2)*Q(2)-JJ(3)*Q(3))/MW2
2598 C
2599
2600       JWWW(1) = DCMPLX( (JJ(0)-JQ*Q(0))*DW )
2601       JWWW(2) = DCMPLX( (JJ(1)-JQ*Q(1))*DW )
2602       JWWW(3) = DCMPLX( (JJ(2)-JQ*Q(2))*DW )
2603       JWWW(4) = DCMPLX( (JJ(3)-JQ*Q(3))*DW )
2604 C
2605       RETURN
2606       END
2607
2608 C
2609 C ----------------------------------------------------------------------
2610 C
2611       SUBROUTINE MOM2CX(ESUM,MASS1,MASS2,COSTH1,PHI1 , P1,P2)
2612 C
2613 C This subroutine sets up two four-momenta in the two particle rest
2614 C frame.
2615 C
2616 C INPUT:
2617 C       real    ESUM           : energy sum of particle 1 and 2
2618 C       real    MASS1          : mass            of particle 1
2619 C       real    MASS2          : mass            of particle 2
2620 C       real    COSTH1         : cos(theta)      of particle 1
2621 C       real    PHI1           : azimuthal angle of particle 1
2622 C
2623 C OUTPUT:
2624 C       real    P1(0:3)        : four-momentum of particle 1
2625 C       real    P2(0:3)        : four-momentum of particle 2
2626 C
2627 #if defined(CERNLIB_IMPNONE)
2628       IMPLICIT NONE
2629 #endif
2630       REAL*8    P1(0:3),P2(0:3),
2631      &        ESUM,MASS1,MASS2,COSTH1,PHI1,MD2,ED,PP,SINTH1
2632 C
2633       MD2=(MASS1-MASS2)*(MASS1+MASS2)
2634       ED=MD2/ESUM
2635       IF (MASS1*MASS2.EQ.0.) THEN
2636       PP=(ESUM-ABS(ED))*0.5D0
2637 C
2638       ELSE
2639       PP=SQRT((MD2/ESUM)**2-2.0D0*(MASS1**2+MASS2**2)+ESUM**2)*0.5D0
2640       ENDIF
2641       SINTH1=SQRT((1.0D0-COSTH1)*(1.0D0+COSTH1))
2642 C
2643       P1(0) = MAX((ESUM+ED)*0.5D0,0.D0)
2644       P1(1) = PP*SINTH1*COS(PHI1)
2645       P1(2) = PP*SINTH1*SIN(PHI1)
2646       P1(3) = PP*COSTH1
2647 C
2648       P2(0) = MAX((ESUM-ED)*0.5D0,0.D0)
2649       P2(1) = -P1(1)
2650       P2(2) = -P1(2)
2651       P2(3) = -P1(3)
2652 C
2653       RETURN
2654       END
2655 C **********************************************************************
2656 C
2657       SUBROUTINE MOMNTX(ENERGY,MASS,COSTH,PHI , P)
2658 C
2659 C This subroutine sets up a four-momentum from the four inputs.
2660 C
2661 C INPUT:
2662 C       real    ENERGY         : energy
2663 C       real    MASS           : mass
2664 C       real    COSTH          : cos(theta)
2665 C       real    PHI            : azimuthal angle
2666 C
2667 C OUTPUT:
2668 C       real    P(0:3)         : four-momentum
2669 C
2670 #if defined(CERNLIB_IMPNONE)
2671       IMPLICIT NONE
2672 #endif
2673       REAL*8    P(0:3),ENERGY,MASS,COSTH,PHI,PP,SINTH
2674 C
2675       P(0) = ENERGY
2676       IF (ENERGY.EQ.MASS) THEN
2677          P(1) = 0.
2678          P(2) = 0.
2679          P(3) = 0.
2680       ELSE
2681          PP=SQRT((ENERGY-MASS)*(ENERGY+MASS))
2682          SINTH=SQRT((1.-COSTH)*(1.+COSTH))
2683          P(3) = PP*COSTH
2684          IF (PHI.EQ.0.) THEN
2685             P(1) = PP*SINTH
2686             P(2) = 0.
2687          ELSE
2688             P(1) = PP*SINTH*COS(PHI)
2689             P(2) = PP*SINTH*SIN(PHI)
2690          ENDIF
2691       ENDIF
2692       RETURN
2693       END
2694 C
2695 C
2696 C
2697 C       Subroutine returns the desired fermion or
2698 C       anti-fermion anti-spinor. ie., <f|
2699 C       A replacement for the HELAS routine OXXXXX
2700 C
2701 C       Adam Duff,  1992 August 31
2702 C       <duff@phenom.physics.wisc.edu>
2703 C
2704       SUBROUTINE OXXXXX(P,FMASS,NHEL,NSF,FO)
2705 C
2706 C          P            IN: FOUR VECTOR MOMENTUM
2707 C          FMASS        IN: FERMION MASS
2708 C          NHEL         IN: ANTI-SPINOR HELICITY, -1 OR 1
2709 C          NSF          IN: -1=ANTIFERMION, 1=FERMION
2710 C          FO           OUT: FERMION WAVEFUNCTION
2711 C
2712 C declare input/output variables
2713 C
2714 #if defined(CERNLIB_IMPNONE)
2715       IMPLICIT NONE
2716 #endif
2717       COMPLEX*16 FO(6)
2718       INTEGER*4 NHEL, NSF
2719       REAL*8 P(0:3), FMASS
2720 C
2721 C declare local variables
2722 C
2723       REAL*8 RXZERO, RXONE, RXTWO
2724       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
2725       REAL*8 PLAT, PABS, OMEGAP, OMEGAM, RS2PA, SPAZ
2726       COMPLEX*16 CXZERO
2727       LOGICAL FIRST
2728       SAVE CXZERO,FIRST
2729       DATA FIRST/.TRUE./
2730 C
2731 C          Fix compilation with g77
2732       IF(FIRST) THEN
2733         FIRST=.FALSE.
2734         CXZERO=DCMPLX( RXZERO, RXZERO )
2735       ENDIF
2736 C
2737 C define kinematic parameters
2738 C
2739       FO(5) = DCMPLX( P(0), P(3) ) * NSF
2740       FO(6) = DCMPLX( P(1), P(2) ) * NSF
2741       PLAT = SQRT( P(1)**2 + P(2)**2 )
2742       PABS = SQRT( P(1)**2 + P(2)**2 + P(3)**2 )
2743       OMEGAP = SQRT( P(0) + PABS )
2744 C
2745 C do massive fermion case
2746 C
2747       IF ( FMASS .NE. RXZERO ) THEN
2748          OMEGAM = FMASS / OMEGAP
2749          IF ( NSF .EQ. 1 ) THEN
2750             IF ( NHEL .EQ. 1 ) THEN
2751                IF ( P(3) .GE. RXZERO ) THEN
2752                   IF ( PLAT .EQ. RXZERO ) THEN
2753                      FO(1) = DCMPLX( OMEGAP, RXZERO )
2754                      FO(2) = CXZERO
2755                      FO(3) = DCMPLX( OMEGAM, RXZERO )
2756                      FO(4) = CXZERO
2757                   ELSE
2758                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2759                      SPAZ = SQRT( PABS + P(3) )
2760                      FO(1) = OMEGAP * RS2PA
2761      &                     * DCMPLX( SPAZ, RXZERO )
2762                      FO(2) = OMEGAP * RS2PA / SPAZ
2763      &                     * DCMPLX( P(1), -P(2) )
2764                      FO(3) = OMEGAM * RS2PA
2765      &                     * DCMPLX( SPAZ, RXZERO )
2766                      FO(4) = OMEGAM * RS2PA / SPAZ
2767      &                     * DCMPLX( P(1), -P(2) )
2768                   END IF
2769                ELSE
2770                   IF ( PLAT .EQ. RXZERO ) THEN
2771                      FO(1) = CXZERO
2772                      FO(2) = DCMPLX( OMEGAP, RXZERO )
2773                      FO(3) = CXZERO
2774                      FO(4) = DCMPLX( OMEGAM, RXZERO )
2775                   ELSE
2776                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2777                      SPAZ = SQRT( PABS - P(3) )
2778                      FO(1) = OMEGAP * RS2PA / SPAZ
2779      &                     * DCMPLX( PLAT, RXZERO )
2780                      FO(2) = OMEGAP * RS2PA * SPAZ / PLAT
2781      &                     * DCMPLX( P(1), -P(2) )
2782                      FO(3) = OMEGAM * RS2PA / SPAZ
2783      &                     * DCMPLX( PLAT, RXZERO )
2784                      FO(4) = OMEGAM * RS2PA * SPAZ / PLAT
2785      &                     * DCMPLX( P(1), -P(2) )
2786                   END IF
2787                END IF
2788             ELSE IF ( NHEL .EQ. -1 ) THEN
2789                IF ( P(3) .GE. RXZERO ) THEN
2790                   IF ( PLAT .EQ. RXZERO ) THEN
2791                      FO(1) = CXZERO
2792                      FO(2) = DCMPLX( OMEGAM, RXZERO )
2793                      FO(3) = CXZERO
2794                      FO(4) = DCMPLX( OMEGAP, RXZERO )
2795                   ELSE
2796                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2797                      SPAZ = SQRT( PABS + P(3) )
2798                      FO(1) = OMEGAM * RS2PA / SPAZ
2799      &                     * DCMPLX( -P(1), -P(2) )
2800                      FO(2) = OMEGAM * RS2PA
2801      &                     * DCMPLX( SPAZ, RXZERO )
2802                      FO(3) = OMEGAP * RS2PA / SPAZ
2803      &                     * DCMPLX( -P(1), -P(2) )
2804                      FO(4) = OMEGAP * RS2PA
2805      &                     * DCMPLX( SPAZ, RXZERO )
2806                   END IF
2807                ELSE
2808                   IF ( PLAT .EQ. RXZERO ) THEN
2809                      FO(1) = DCMPLX( -OMEGAM, RXZERO )
2810                      FO(2) = CXZERO
2811                      FO(3) = DCMPLX( -OMEGAP, RXZERO )
2812                      FO(4) = CXZERO
2813                   ELSE
2814                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2815                      SPAZ = SQRT( PABS - P(3) )
2816                      FO(1) = OMEGAM * RS2PA * SPAZ / PLAT
2817      &                     * DCMPLX( -P(1), -P(2) )
2818                      FO(2) = OMEGAM * RS2PA / SPAZ
2819      &                     * DCMPLX( PLAT, RXZERO )
2820                      FO(3) = OMEGAP * RS2PA * SPAZ / PLAT
2821      &                     * DCMPLX( -P(1), -P(2) )
2822                      FO(4) = OMEGAP * RS2PA / SPAZ
2823      &                     * DCMPLX( PLAT, RXZERO )
2824                   END IF
2825                END IF
2826             ELSE
2827                STOP 'OXXXXX:  FERMION HELICITY MUST BE +1,-1'
2828             END IF
2829          ELSE IF ( NSF .EQ. -1 ) THEN
2830             IF ( NHEL .EQ. 1 ) THEN
2831                IF ( P(3) .GE. RXZERO ) THEN
2832                   IF ( PLAT .EQ. RXZERO ) THEN
2833                      FO(1) = CXZERO
2834                      FO(2) = DCMPLX( OMEGAM, RXZERO )
2835                      FO(3) = CXZERO
2836                      FO(4) = DCMPLX( -OMEGAP, RXZERO )
2837                   ELSE
2838                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2839                      SPAZ = SQRT( PABS + P(3) )
2840                      FO(1) = OMEGAM * RS2PA / SPAZ
2841      &                     * DCMPLX( -P(1), -P(2) )
2842                      FO(2) = OMEGAM * RS2PA
2843      &                     * DCMPLX( SPAZ, RXZERO )
2844                      FO(3) = -OMEGAP * RS2PA / SPAZ
2845      &                     * DCMPLX( -P(1), -P(2) )
2846                      FO(4) = -OMEGAP * RS2PA
2847      &                     * DCMPLX( SPAZ, RXZERO )
2848                   END IF
2849                ELSE
2850                   IF ( PLAT .EQ. RXZERO ) THEN
2851                      FO(1) = DCMPLX( -OMEGAM, RXZERO )
2852                      FO(2) = CXZERO
2853                      FO(3) = DCMPLX( OMEGAP, RXZERO )
2854                      FO(4) = CXZERO
2855                   ELSE
2856                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2857                      SPAZ = SQRT( PABS - P(3) )
2858                      FO(1) = OMEGAM * RS2PA * SPAZ / PLAT
2859      &                     * DCMPLX( -P(1), -P(2) )
2860                      FO(2) = OMEGAM * RS2PA / SPAZ
2861      &                     * DCMPLX( PLAT, RXZERO )
2862                      FO(3) = -OMEGAP * RS2PA * SPAZ / PLAT
2863      &                     * DCMPLX( -P(1), -P(2) )
2864                      FO(4) = -OMEGAP * RS2PA / SPAZ
2865      &                     * DCMPLX( PLAT, RXZERO )
2866                   END IF
2867                END IF
2868             ELSE IF ( NHEL .EQ. -1 ) THEN
2869                IF ( P(3) .GE. RXZERO ) THEN
2870                   IF ( PLAT .EQ. RXZERO ) THEN
2871                      FO(1) = DCMPLX( -OMEGAP, RXZERO )
2872                      FO(2) = CXZERO
2873                      FO(3) = DCMPLX( OMEGAM, RXZERO )
2874                      FO(4) = CXZERO
2875                   ELSE
2876                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2877                      SPAZ = SQRT( PABS + P(3) )
2878                      FO(1) = -OMEGAP * RS2PA
2879      &                     * DCMPLX( SPAZ, RXZERO )
2880                      FO(2) = -OMEGAP * RS2PA / SPAZ
2881      &                     * DCMPLX( P(1), -P(2) )
2882                      FO(3) = OMEGAM * RS2PA
2883      &                     * DCMPLX( SPAZ, RXZERO )
2884                      FO(4) = OMEGAM * RS2PA / SPAZ
2885      &                     * DCMPLX( P(1), -P(2) )
2886                   END IF
2887                ELSE
2888                   IF ( PLAT .EQ. RXZERO ) THEN
2889                      FO(1) = CXZERO
2890                      FO(2) = DCMPLX( -OMEGAP, RXZERO )
2891                      FO(3) = CXZERO
2892                      FO(4) = DCMPLX( OMEGAM, RXZERO )
2893                   ELSE
2894                      RS2PA = RXONE / SQRT( RXTWO * PABS )
2895                      SPAZ = SQRT( PABS - P(3) )
2896                      FO(1) = -OMEGAP * RS2PA / SPAZ
2897      &                     * DCMPLX( PLAT, RXZERO )
2898                      FO(2) = -OMEGAP * RS2PA * SPAZ / PLAT
2899      &                     * DCMPLX( P(1), -P(2) )
2900                      FO(3) = OMEGAM * RS2PA / SPAZ
2901      &                     * DCMPLX( PLAT, RXZERO )
2902                      FO(4) = OMEGAM * RS2PA * SPAZ / PLAT
2903      &                     * DCMPLX( P(1), -P(2) )
2904                   END IF
2905                END IF
2906             ELSE
2907                STOP 'OXXXXX:  FERMION HELICITY MUST BE +1,-1'
2908             END IF
2909          ELSE
2910             STOP 'OXXXXX:  FERMION TYPE MUST BE +1,-1'
2911          END IF
2912 C
2913 C do massless case
2914 C
2915       ELSE
2916          IF ( NSF .EQ. 1 ) THEN
2917             IF ( NHEL .EQ. 1 ) THEN
2918                IF ( P(3) .GE. RXZERO ) THEN
2919                   IF ( PLAT .EQ. RXZERO ) THEN
2920                      FO(1) = DCMPLX( OMEGAP, RXZERO )
2921                      FO(2) = CXZERO
2922                      FO(3) = CXZERO
2923                      FO(4) = CXZERO
2924                   ELSE
2925                      SPAZ = SQRT( PABS + P(3) )
2926                      FO(1) = DCMPLX( SPAZ, RXZERO )
2927                      FO(2) = RXONE / SPAZ
2928      &                     * DCMPLX( P(1), -P(2) )
2929                      FO(3) = CXZERO
2930                      FO(4) = CXZERO
2931                   END IF
2932                ELSE
2933                   IF ( PLAT .EQ. RXZERO ) THEN
2934                      FO(1) = CXZERO
2935                      FO(2) = DCMPLX( OMEGAP, RXZERO )
2936                      FO(3) = CXZERO
2937                      FO(4) = CXZERO
2938                   ELSE
2939                      SPAZ = SQRT( PABS - P(3) )
2940                      FO(1) = RXONE / SPAZ
2941      &                     * DCMPLX( PLAT, RXZERO )
2942                      FO(2) = SPAZ / PLAT
2943      &                     * DCMPLX( P(1), -P(2) )
2944                      FO(3) = CXZERO
2945                      FO(4) = CXZERO
2946                   END IF
2947                END IF
2948             ELSE IF ( NHEL .EQ. -1 ) THEN
2949                IF ( P(3) .GE. RXZERO ) THEN
2950                   IF ( PLAT .EQ. RXZERO ) THEN
2951                      FO(1) = CXZERO
2952                      FO(2) = CXZERO
2953                      FO(3) = CXZERO
2954                      FO(4) = DCMPLX( OMEGAP, RXZERO )
2955                   ELSE
2956                      SPAZ = SQRT( PABS + P(3) )
2957                      FO(1) = CXZERO
2958                      FO(2) = CXZERO
2959                      FO(3) = RXONE / SPAZ
2960      &                     * DCMPLX( -P(1), -P(2) )
2961                      FO(4) = DCMPLX( SPAZ, RXZERO )
2962                   END IF
2963                ELSE
2964                   IF ( PLAT .EQ. RXZERO ) THEN
2965                      FO(1) = CXZERO
2966                      FO(2) = CXZERO
2967                      FO(3) = DCMPLX( -OMEGAP, RXZERO )
2968                      FO(4) = CXZERO
2969                   ELSE
2970                      SPAZ = SQRT( PABS - P(3) )
2971                      FO(1) = CXZERO
2972                      FO(2) = CXZERO
2973                      FO(3) = SPAZ / PLAT
2974      &                     * DCMPLX( -P(1), -P(2) )
2975                      FO(4) = RXONE / SPAZ
2976      &                     * DCMPLX( PLAT, RXZERO )
2977                   END IF
2978                END IF
2979             ELSE
2980                STOP 'OXXXXX:  FERMION HELICITY MUST BE +1,-1'
2981             END IF
2982          ELSE IF ( NSF .EQ. -1 ) THEN
2983             IF ( NHEL .EQ. 1 ) THEN
2984                IF ( P(3) .GE. RXZERO ) THEN
2985                   IF ( PLAT .EQ. RXZERO ) THEN
2986                      FO(1) = CXZERO
2987                      FO(2) = CXZERO
2988                      FO(3) = CXZERO
2989                      FO(4) = DCMPLX( -OMEGAP, RXZERO )
2990                   ELSE
2991                      SPAZ = SQRT( PABS + P(3) )
2992                      FO(1) = CXZERO
2993                      FO(2) = CXZERO
2994                      FO(3) = -RXONE / SPAZ
2995      &                     * DCMPLX( -P(1), -P(2) )
2996                      FO(4) = DCMPLX( -SPAZ, RXZERO )
2997                   END IF
2998                ELSE
2999                   IF ( PLAT .EQ. RXZERO ) THEN
3000                      FO(1) = CXZERO
3001                      FO(2) = CXZERO
3002                      FO(3) = DCMPLX( OMEGAP, RXZERO )
3003                      FO(4) = CXZERO
3004                   ELSE
3005                      SPAZ = SQRT( PABS - P(3) )
3006                      FO(1) = CXZERO
3007                      FO(2) = CXZERO
3008                      FO(3) = -SPAZ / PLAT
3009      &                     * DCMPLX( -P(1), -P(2) )
3010                      FO(4) = -RXONE / SPAZ
3011      &                     * DCMPLX( PLAT, RXZERO )
3012                   END IF
3013                END IF
3014             ELSE IF ( NHEL .EQ. -1 ) THEN
3015                IF ( P(3) .GE. RXZERO ) THEN
3016                   IF ( PLAT .EQ. RXZERO ) THEN
3017                      FO(1) = DCMPLX( -OMEGAP, RXZERO )
3018                      FO(2) = CXZERO
3019                      FO(3) = CXZERO
3020                      FO(4) = CXZERO
3021                   ELSE
3022                      SPAZ = SQRT( PABS + P(3) )
3023                      FO(1) = DCMPLX( -SPAZ, RXZERO )
3024                      FO(2) = -RXONE / SPAZ
3025      &                     * DCMPLX( P(1), -P(2) )
3026                      FO(3) = CXZERO
3027                      FO(4) = CXZERO
3028                   END IF
3029                ELSE
3030                   IF ( PLAT .EQ. RXZERO ) THEN
3031                      FO(1) = CXZERO
3032                      FO(2) = DCMPLX( -OMEGAP, RXZERO )
3033                      FO(3) = CXZERO
3034                      FO(4) = CXZERO
3035                   ELSE
3036                      SPAZ = SQRT( PABS - P(3) )
3037                      FO(1) = -RXONE / SPAZ
3038      &                     * DCMPLX( PLAT, RXZERO )
3039                      FO(2) = -SPAZ / PLAT
3040      &                     * DCMPLX( P(1), -P(2) )
3041                      FO(3) = CXZERO
3042                      FO(4) = CXZERO
3043                   END IF
3044                END IF
3045             ELSE
3046                STOP 'OXXXXX:  FERMION HELICITY MUST BE +1,-1'
3047             END IF
3048          ELSE
3049             STOP 'OXXXXX:  FERMION TYPE MUST BE +1,-1'
3050          END IF
3051       END IF
3052 C
3053 C done
3054 C
3055       RETURN
3056       END
3057 C
3058 C ----------------------------------------------------------------------
3059 C
3060       SUBROUTINE ROTXXX(P,Q , PROT)
3061 C
3062 C this subroutine performs the spacial rotation of a four-momentum.
3063 C the momentum p is assumed to be given in the frame where the spacial
3064 C component of q points the positive z-axis.  prot is the momentum p
3065 C rotated to the frame where q is given.
3066 C
3067 C input:
3068 C       real    p(0:3)         : four-momentum p in q(1)=q(2)=0 frame
3069 C       real    q(0:3)         : four-momentum q in the rotated frame
3070 C
3071 C output:
3072 C       real    prot(0:3)      : four-momentum p in the rotated frame
3073 C
3074 #if defined(CERNLIB_IMPNONE)
3075       IMPLICIT NONE
3076 #endif
3077       REAL*8    P(0:3),Q(0:3),PROT(0:3),QT2,QT,PSGN,QQ,P1
3078 C
3079       REAL*8 RXZERO, RXONE
3080       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
3081 C
3082       PROT(0) = P(0)
3083 C
3084       QT2=Q(1)**2+Q(2)**2
3085 C
3086       IF ( QT2 .EQ. RXZERO ) THEN
3087           IF ( Q(3) .EQ. RXZERO ) THEN
3088              PROT(1) = P(1)
3089              PROT(2) = P(2)
3090              PROT(3) = P(3)
3091           ELSE
3092              PSGN=DSIGN(RXONE,Q(3))
3093              PROT(1) = P(1)*PSGN
3094              PROT(2) = P(2)*PSGN
3095              PROT(3) = P(3)*PSGN
3096           ENDIF
3097       ELSE
3098           QQ=SQRT(QT2+Q(3)**2)
3099           QT=SQRT(QT2)
3100           P1=P(1)
3101           PROT(1) = Q(1)*Q(3)/QQ/QT*P1 -Q(2)/QT*P(2) +Q(1)/QQ*P(3)
3102           PROT(2) = Q(2)*Q(3)/QQ/QT*P1 +Q(1)/QT*P(2) +Q(2)/QQ*P(3)
3103           PROT(3) =          -QT/QQ*P1               +Q(3)/QQ*P(3)
3104       ENDIF
3105 C
3106       RETURN
3107       END
3108 C ======================================================================
3109 C
3110       SUBROUTINE SSSSXX(S1,S2,S3,S4,G , VERTEX)
3111 C
3112 C This subroutine computes an amplitude of the four-scalar coupling.
3113 C
3114 C INPUT:
3115 C       complex S1(3)          : first  scalar                        S1
3116 C       complex S2(3)          : second scalar                        S2
3117 C       complex S3(3)          : third  scalar                        S3
3118 C       complex S4(3)          : fourth scalar                        S4
3119 C       real    G              : coupling constant                 GHHHH
3120 C
3121 C OUTPUT:
3122 C       complex VERTEX         : amplitude            Gamma(S1,S2,S3,S4)
3123 C
3124 #if defined(CERNLIB_IMPNONE)
3125       IMPLICIT NONE
3126 #endif
3127       COMPLEX*16 S1(3),S2(3),S3(3),S4(3),VERTEX
3128       REAL*8     G
3129 C
3130       VERTEX = G*S1(1)*S2(1)*S3(1)*S4(1)
3131 C
3132       RETURN
3133       END
3134 C
3135 C ======================================================================
3136 C
3137       SUBROUTINE SSSXXX(S1,S2,S3,G , VERTEX)
3138 C
3139 C This subroutine computes an amplitude of the three-scalar coupling.
3140 C
3141 C INPUT:
3142 C       complex S1(3)          : first  scalar                        S1
3143 C       complex S2(3)          : second scalar                        S2
3144 C       complex S3(3)          : third  scalar                        S3
3145 C       real    G              : coupling constant                  GHHH
3146 C
3147 C OUTPUT:
3148 C       complex VERTEX         : amplitude               Gamma(S1,S2,S3)
3149 C
3150 #if defined(CERNLIB_IMPNONE)
3151       IMPLICIT NONE
3152 #endif
3153       COMPLEX*16 S1(3),S2(3),S3(3),VERTEX
3154       REAL*8    G
3155 C
3156       VERTEX = G*S1(1)*S2(1)*S3(1)
3157 C
3158       RETURN
3159       END
3160 C
3161 C
3162 C ----------------------------------------------------------------------
3163 C
3164       SUBROUTINE SXXXXX(P,NSS , SC)
3165 C
3166 C This subroutine computes a complex SCALAR wavefunction.
3167 C
3168 C INPUT:
3169 C       real    P(0:3)         : four-momentum of scalar boson
3170 C       integer NSS  = -1 or 1 : +1 for final, -1 for initial
3171 C
3172 C OUTPUT:
3173 C       complex SC(3)          : scalar wavefunction                   S
3174 C
3175 #if defined(CERNLIB_IMPNONE)
3176       IMPLICIT NONE
3177 #endif
3178       COMPLEX*16 SC(3)
3179       REAL*8    P(0:3)
3180       INTEGER NSS
3181 C
3182       SC(1) = DCMPLX( 1.0 )
3183       SC(2) = DCMPLX(P(0),P(3))*NSS
3184       SC(3) = DCMPLX(P(1),P(2))*NSS
3185 C
3186       RETURN
3187       END
3188 C
3189 C ======================================================================
3190 C
3191       SUBROUTINE VSSXXX(VC,S1,S2,G , VERTEX)
3192 C
3193 C this subroutine computes an amplitude from the vector-scalar-scalar
3194 C coupling.  the coupling is absent in the minimal sm in unitary gauge.
3195 C
3196 C       complex vc(6)          : input  vector                        v
3197 C       complex s1(3)          : first  scalar                        s1
3198 C       complex s2(3)          : second scalar                        s2
3199 C       complex g              : coupling constant (s1 charge)
3200 C
3201 C examples of the coupling constant g for susy particles are as follows:
3202 C   -----------------------------------------------------------
3203 C   |    s1    | (q,i3) of s1  ||   v=a   |   v=z   |   v=w   |
3204 C   -----------------------------------------------------------
3205 C   | nu~_l    | (  0  , +1/2) ||   ---   |  gzn(1) |  gwf(1) |
3206 C   | e~_l     | ( -1  , -1/2) ||  gal(1) |  gzl(1) |  gwf(1) |
3207 C   | u~_l     | (+2/3 , +1/2) ||  gau(1) |  gzu(1) |  gwf(1) |
3208 C   | d~_l     | (-1/3 , -1/2) ||  gad(1) |  gzd(1) |  gwf(1) |
3209 C   -----------------------------------------------------------
3210 C   | e~_r-bar | ( +1  ,  0  ) || -gal(2) | -gzl(2) | -gwf(2) |
3211 C   | u~_r-bar | (-2/3 ,  0  ) || -gau(2) | -gzu(2) | -gwf(2) |
3212 C   | d~_r-bar | (+1/3 ,  0  ) || -gad(2) | -gzd(2) | -gwf(2) |
3213 C   -----------------------------------------------------------
3214 C where the s1 charge is defined by the flowing-out quantum number.
3215 C
3216 C output:
3217 C       complex vertex         : amplitude                gamma(v,s1,s2)
3218 C
3219 #if defined(CERNLIB_IMPNONE)
3220       IMPLICIT NONE
3221 #endif
3222       COMPLEX*16 VC(6),S1(3),S2(3),VERTEX,G
3223       REAL*8    P(0:3)
3224 C
3225       P(0)=DBLE( S1(2)-S2(2))
3226       P(1)=DBLE( S1(3)-S2(3))
3227       P(2)=DIMAG(S1(3)-S2(3))
3228       P(3)=DIMAG(S1(2)-S2(2))
3229 C
3230       VERTEX = G*S1(1)*S2(1)
3231      &        *(VC(1)*P(0)-VC(2)*P(1)-VC(3)*P(2)-VC(4)*P(3))
3232 C
3233       RETURN
3234       END
3235 C
3236       SUBROUTINE VVSSXX(V1,V2,S1,S2,G , VERTEX)
3237 C
3238 C This subroutine computes an amplitude of the vector-vector-scalar-
3239 C scalar coupling.
3240 C
3241 C INPUT:
3242 C       complex V1(6)          : first  vector                        V1
3243 C       complex V2(6)          : second vector                        V2
3244 C       complex S1(3)          : first  scalar                        S1
3245 C       complex S2(3)          : second scalar                        S2
3246 C       real    G              : coupling constant                 GVVHH
3247 C
3248 C OUTPUT:
3249 C       complex VERTEX         : amplitude            Gamma(V1,V2,S1,S2)
3250 C
3251 #if defined(CERNLIB_IMPNONE)
3252       IMPLICIT NONE
3253 #endif
3254       COMPLEX*16 V1(6),V2(6),S1(3),S2(3),VERTEX
3255       REAL*8    G
3256 C
3257       VERTEX = G*S1(1)*S2(1)
3258      &        *(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
3259 C
3260       RETURN
3261       END
3262 C
3263 C
3264 C ======================================================================
3265 C
3266       SUBROUTINE VVSXXX(V1,V2,SC,G , VERTEX)
3267 C
3268 C this subroutine computes an amplitude of the vector-vector-scalar
3269 C coupling.
3270 C
3271 C input:
3272 C       complex v1(6)          : first  vector                        v1
3273 C       complex v2(6)          : second vector                        v2
3274 C       complex sc(3)          : input  scalar                        s
3275 C       real    g              : coupling constant                  gvvh
3276 C
3277 C output:
3278 C       complex vertex         : amplitude                gamma(v1,v2,s)
3279 C
3280 #if defined(CERNLIB_IMPNONE)
3281       IMPLICIT NONE
3282 #endif
3283       COMPLEX*16 V1(6),V2(6),SC(3),VERTEX
3284       REAL*8    G
3285 C
3286       VERTEX = G*SC(1)*(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
3287 C
3288       RETURN
3289       END
3290 C
3291 C ======================================================================
3292 C
3293       SUBROUTINE VVVXXX(WM,WP,W3,G , VERTEX)
3294 C
3295 C this subroutine computes an amplitude of the three-point coupling of
3296 C the gauge bosons.
3297 C
3298 C input:
3299 C       complex wm(6)          : vector               flow-out w-
3300 C       complex wp(6)          : vector               flow-out w+
3301 C       complex w3(6)          : vector               j3 or a    or z
3302 C       real    g              : coupling constant    gw or gwwa or gwwz
3303 C
3304 C output:
3305 C       complex vertex         : amplitude               gamma(wm,wp,w3)
3306 C
3307 #if defined(CERNLIB_IMPNONE)
3308       IMPLICIT NONE
3309 #endif
3310       COMPLEX*16 WM(6),WP(6),W3(6),VERTEX,
3311      &        XV1,XV2,XV3,V12,V23,V31,P12,P13,P21,P23,P31,P32
3312       REAL*8    PWM(0:3),PWP(0:3),PW3(0:3),G
3313 C
3314       REAL*8 RXZERO, RTENTH
3315       PARAMETER( RXZERO=0.0D0, RTENTH=0.1D0 )
3316 C
3317       PWM(0)=DBLE( WM(5))
3318       PWM(1)=DBLE( WM(6))
3319       PWM(2)=DIMAG(WM(6))
3320       PWM(3)=DIMAG(WM(5))
3321       PWP(0)=DBLE( WP(5))
3322       PWP(1)=DBLE( WP(6))
3323       PWP(2)=DIMAG(WP(6))
3324       PWP(3)=DIMAG(WP(5))
3325       PW3(0)=DBLE( W3(5))
3326       PW3(1)=DBLE( W3(6))
3327       PW3(2)=DIMAG(W3(6))
3328       PW3(3)=DIMAG(W3(5))
3329 C
3330       V12=WM(1)*WP(1)-WM(2)*WP(2)-WM(3)*WP(3)-WM(4)*WP(4)
3331       V23=WP(1)*W3(1)-WP(2)*W3(2)-WP(3)*W3(3)-WP(4)*W3(4)
3332       V31=W3(1)*WM(1)-W3(2)*WM(2)-W3(3)*WM(3)-W3(4)*WM(4)
3333       XV1=RXZERO
3334       XV2=RXZERO
3335       XV3=RXZERO
3336       IF ( ABS(WM(1)) .NE. RXZERO ) THEN
3337          IF (ABS(WM(1)).GE.MAX(ABS(WM(2)),ABS(WM(3)),ABS(WM(4)))
3338      $        *RTENTH)
3339      &      XV1=PWM(0)/WM(1)
3340       ENDIF
3341       IF ( ABS(WP(1)) .NE. RXZERO) THEN
3342          IF (ABS(WP(1)).GE.MAX(ABS(WP(2)),ABS(WP(3)),ABS(WP(4)))
3343      $        *RTENTH)
3344      &      XV2=PWP(0)/WP(1)
3345       ENDIF
3346       IF ( ABS(W3(1)) .NE. RXZERO) THEN
3347          IF ( ABS(W3(1)).GE.MAX(ABS(W3(2)),ABS(W3(3)),ABS(W3(4)))
3348      $        *RTENTH)
3349      &      XV3=PW3(0)/W3(1)
3350       ENDIF
3351       P12= (PWM(0)-XV1*WM(1))*WP(1)-(PWM(1)-XV1*WM(2))*WP(2)
3352      &    -(PWM(2)-XV1*WM(3))*WP(3)-(PWM(3)-XV1*WM(4))*WP(4)
3353       P13= (PWM(0)-XV1*WM(1))*W3(1)-(PWM(1)-XV1*WM(2))*W3(2)
3354      &    -(PWM(2)-XV1*WM(3))*W3(3)-(PWM(3)-XV1*WM(4))*W3(4)
3355       P21= (PWP(0)-XV2*WP(1))*WM(1)-(PWP(1)-XV2*WP(2))*WM(2)
3356      &    -(PWP(2)-XV2*WP(3))*WM(3)-(PWP(3)-XV2*WP(4))*WM(4)
3357       P23= (PWP(0)-XV2*WP(1))*W3(1)-(PWP(1)-XV2*WP(2))*W3(2)
3358      &    -(PWP(2)-XV2*WP(3))*W3(3)-(PWP(3)-XV2*WP(4))*W3(4)
3359       P31= (PW3(0)-XV3*W3(1))*WM(1)-(PW3(1)-XV3*W3(2))*WM(2)
3360      &    -(PW3(2)-XV3*W3(3))*WM(3)-(PW3(3)-XV3*W3(4))*WM(4)
3361       P32= (PW3(0)-XV3*W3(1))*WP(1)-(PW3(1)-XV3*W3(2))*WP(2)
3362      &    -(PW3(2)-XV3*W3(3))*WP(3)-(PW3(3)-XV3*W3(4))*WP(4)
3363 C
3364       VERTEX = -(V12*(P13-P23)+V23*(P21-P31)+V31*(P32-P12))*G
3365 C
3366       RETURN
3367       END
3368 C
3369 C
3370 C       Subroutine returns the value of evaluated
3371 C       helicity basis boson polarisation wavefunction.
3372 C       Replaces the HELAS routine VXXXXX
3373 C
3374 C       Adam Duff,  1992 September 3
3375 C       <duff@phenom.physics.wisc.edu>
3376 C
3377       SUBROUTINE VXXXXX(P,VMASS,NHEL,NSV,VC)
3378 C
3379 C          P            IN: BOSON FOUR MOMENTUM
3380 C          VMASS        IN: BOSON MASS
3381 C          NHEL         IN: BOSON HELICITY
3382 C          NSV          IN: INCOMING (-1) OR OUTGOING (+1)
3383 C          VC           OUT: BOSON WAVEFUNCTION
3384 C
3385 C declare input/output variables
3386 C
3387 #if defined(CERNLIB_IMPNONE)
3388       IMPLICIT NONE
3389 #endif
3390       COMPLEX*16 VC(6)
3391       INTEGER*4 NHEL, NSV
3392       REAL*8 P(0:3), VMASS
3393 C
3394 C declare local variables
3395 C
3396       REAL*8 RXZERO, RXONE, RXTWO
3397       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
3398       REAL*8 PLAT, PABS, RS2, RPLAT, RPABS, RDEN
3399       COMPLEX*16 CXZERO
3400       LOGICAL FIRST
3401       SAVE CXZERO,FIRST
3402       DATA FIRST/.TRUE./
3403 C
3404 C          Fix compilation with g77
3405       IF(FIRST) THEN
3406         FIRST=.FALSE.
3407         CXZERO=DCMPLX( RXZERO, RXZERO )
3408       ENDIF
3409 C
3410 C define internal/external momenta
3411 C
3412       IF ( NSV**2 .NE. 1 ) THEN
3413          STOP 'VXXXXX:  NSV IS NOT ONE OF -1, +1'
3414       END IF
3415 C
3416       RS2 = SQRT( RXONE / RXTWO )
3417       VC(5) = DCMPLX( P(0), P(3) ) * NSV
3418       VC(6) = DCMPLX( P(1), P(2) ) * NSV
3419       PLAT = SQRT( P(1)**2 + P(2)**2 )
3420       PABS = SQRT( P(1)**2 + P(2)**2 + P(3)**2 )
3421 C
3422 C calculate polarisation four vectors
3423 C
3424       IF ( NHEL**2 .EQ. 1 ) THEN
3425          IF ( (PABS .EQ. RXZERO) .OR. (PLAT .EQ. RXZERO) ) THEN
3426             VC(1) = CXZERO
3427             VC(2) = DCMPLX( -NHEL * RS2 * DSIGN( RXONE, P(3) ), RXZERO )
3428             VC(3) = DCMPLX( RXZERO, NSV * RS2 )
3429             VC(4) = CXZERO
3430          ELSE
3431             RPLAT = RXONE / PLAT
3432             RPABS = RXONE / PABS
3433             VC(1) = CXZERO
3434             VC(2) = DCMPLX( -NHEL * RS2 * RPABS * RPLAT * P(1) * P(3),
3435      &                     -NSV * RS2 * RPLAT * P(2) )
3436             VC(3) = DCMPLX( -NHEL * RS2 * RPABS * RPLAT * P(2) * P(3),
3437      &                     NSV * RS2 * RPLAT * P(1) )
3438             VC(4) = DCMPLX( NHEL * RS2 * RPABS * PLAT,
3439      &                     RXZERO )
3440          END IF
3441       ELSE IF ( NHEL .EQ. 0 ) THEN
3442          IF ( VMASS .GT. RXZERO ) THEN
3443             IF ( PABS .EQ. RXZERO ) THEN
3444                VC(1) = CXZERO
3445                VC(2) = CXZERO
3446                VC(3) = CXZERO
3447                VC(4) = DCMPLX( RXONE, RXZERO )
3448             ELSE
3449                RDEN = P(0) / ( VMASS * PABS )
3450                VC(1) = DCMPLX( PABS / VMASS, RXZERO )
3451                VC(2) = DCMPLX( RDEN * P(1), RXZERO )
3452                VC(3) = DCMPLX( RDEN * P(2), RXZERO )
3453                VC(4) = DCMPLX( RDEN * P(3), RXZERO )
3454             END IF
3455          ELSE
3456             STOP  'VXXXXX: NHEL = 0 IS ONLY FOR MASSIVE BOSONS'
3457          END IF
3458       ELSE IF ( NHEL .EQ. 4 ) THEN
3459          IF ( VMASS .GT. RXZERO ) THEN
3460             RDEN = RXONE / VMASS
3461             VC(1) = DCMPLX( RDEN * P(0), RXZERO )
3462             VC(2) = DCMPLX( RDEN * P(1), RXZERO )
3463             VC(3) = DCMPLX( RDEN * P(2), RXZERO )
3464             VC(4) = DCMPLX( RDEN * P(3), RXZERO )
3465          ELSEIF (VMASS .EQ. RXZERO) THEN
3466             RDEN = RXONE / P(0)
3467             VC(1) = DCMPLX( RDEN * P(0), RXZERO )
3468             VC(2) = DCMPLX( RDEN * P(1), RXZERO )
3469             VC(3) = DCMPLX( RDEN * P(2), RXZERO )
3470             VC(4) = DCMPLX( RDEN * P(3), RXZERO )
3471          ELSE
3472             STOP 'VXXXXX: NHEL = 4 IS ONLY FOR M>=0'
3473          END IF
3474       ELSE
3475          STOP 'VXXXXX:  NHEL IS NOT ONE OF -1, 0, 1 OR 4'
3476       END IF
3477 C
3478 C done
3479 C
3480       RETURN
3481       END
3482 C
3483 C ----------------------------------------------------------------------
3484 C
3485       SUBROUTINE W3W3XX(WM,W31,WP,W32,G31,G32,WMASS,WWIDTH , VERTEX)
3486 C
3487 C this subroutine computes an amplitude of the four-point coupling of
3488 C the w-, w+ and two w3/z/a.  the amplitude includes the contributions
3489 C of w exchange diagrams.  the internal w propagator is given in unitary
3490 C gauge.  if one sets wmass=0.0, then the gggg vertex is given (see sect
3491 C 2.9.1 of the manual).
3492 C
3493 C input:
3494 C       complex wm(0:3)        : flow-out w-                         wm
3495 C       complex w31(0:3)       : first    w3/z/a                     w31
3496 C       complex wp(0:3)        : flow-out w+                         wp
3497 C       complex w32(0:3)       : second   w3/z/a                     w32
3498 C       real    g31            : coupling of w31 with w-/w+
3499 C       real    g32            : coupling of w32 with w-/w+
3500 C                                                  (see the table below)
3501 C       real    wmass          : mass  of w
3502 C       real    wwidth         : width of w
3503 C
3504 C the possible sets of the inputs are as follows:
3505 C   -------------------------------------------
3506 C   |  wm  |  w31 |  wp  |  w32 |  g31 |  g32 |
3507 C   -------------------------------------------
3508 C   |  w-  |  w3  |  w+  |  w3  |  gw  |  gw  |
3509 C   |  w-  |  w3  |  w+  |  z   |  gw  | gwwz |
3510 C   |  w-  |  w3  |  w+  |  a   |  gw  | gwwa |
3511 C   |  w-  |  z   |  w+  |  z   | gwwz | gwwz |
3512 C   |  w-  |  z   |  w+  |  a   | gwwz | gwwa |
3513 C   |  w-  |  a   |  w+  |  a   | gwwa | gwwa |
3514 C   -------------------------------------------
3515 C where all the bosons are defined by the flowing-out quantum number.
3516 C
3517 C output:
3518 C       complex vertex         : amplitude          gamma(wm,w31,wp,w32)
3519 C
3520 #if defined(CERNLIB_IMPNONE)
3521       IMPLICIT NONE
3522 #endif
3523       COMPLEX*16    WM(6),W31(6),WP(6),W32(6),VERTEX
3524       COMPLEX*16 DV1(0:3),DV2(0:3),DV3(0:3),DV4(0:3),DVERTX,
3525      &           V12,V13,V14,V23,V24,V34
3526       REAL*8     G31,G32,WMASS,WWIDTH
3527 C
3528       REAL*8 RXZERO, RXONE
3529       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
3530
3531       DV1(0)=DCMPLX(WM(1))
3532       DV1(1)=DCMPLX(WM(2))
3533       DV1(2)=DCMPLX(WM(3))
3534       DV1(3)=DCMPLX(WM(4))
3535       DV2(0)=DCMPLX(W31(1))
3536       DV2(1)=DCMPLX(W31(2))
3537       DV2(2)=DCMPLX(W31(3))
3538       DV2(3)=DCMPLX(W31(4))
3539       DV3(0)=DCMPLX(WP(1))
3540       DV3(1)=DCMPLX(WP(2))
3541       DV3(2)=DCMPLX(WP(3))
3542       DV3(3)=DCMPLX(WP(4))
3543       DV4(0)=DCMPLX(W32(1))
3544       DV4(1)=DCMPLX(W32(2))
3545       DV4(2)=DCMPLX(W32(3))
3546       DV4(3)=DCMPLX(W32(4))
3547 C
3548       IF ( DBLE(WMASS) .NE. RXZERO ) THEN
3549 C         dm2inv = r_one / dmw2
3550 C
3551          V12= DV1(0)*DV2(0)-DV1(1)*DV2(1)-DV1(2)*DV2(2)-DV1(3)*DV2(3)
3552          V13= DV1(0)*DV3(0)-DV1(1)*DV3(1)-DV1(2)*DV3(2)-DV1(3)*DV3(3)
3553          V14= DV1(0)*DV4(0)-DV1(1)*DV4(1)-DV1(2)*DV4(2)-DV1(3)*DV4(3)
3554          V23= DV2(0)*DV3(0)-DV2(1)*DV3(1)-DV2(2)*DV3(2)-DV2(3)*DV3(3)
3555          V24= DV2(0)*DV4(0)-DV2(1)*DV4(1)-DV2(2)*DV4(2)-DV2(3)*DV4(3)
3556          V34= DV3(0)*DV4(0)-DV3(1)*DV4(1)-DV3(2)*DV4(2)-DV3(3)*DV4(3)
3557 C
3558          DVERTX = V12*V34 +V14*V23 -2.D0*V13*V24
3559 C
3560          VERTEX = DCMPLX( DVERTX ) * (G31*G32)
3561 C
3562       ELSE
3563          V12= DV1(0)*DV2(0)-DV1(1)*DV2(1)-DV1(2)*DV2(2)-DV1(3)*DV2(3)
3564          V13= DV1(0)*DV3(0)-DV1(1)*DV3(1)-DV1(2)*DV3(2)-DV1(3)*DV3(3)
3565          V14= DV1(0)*DV4(0)-DV1(1)*DV4(1)-DV1(2)*DV4(2)-DV1(3)*DV4(3)
3566          V23= DV2(0)*DV3(0)-DV2(1)*DV3(1)-DV2(2)*DV3(2)-DV2(3)*DV3(3)
3567          V24= DV2(0)*DV4(0)-DV2(1)*DV4(1)-DV2(2)*DV4(2)-DV2(3)*DV4(3)
3568          V34= DV3(0)*DV4(0)-DV3(1)*DV4(1)-DV3(2)*DV4(2)-DV3(3)*DV4(3)
3569 C
3570
3571          DVERTX = V14*V23 -V13*V24
3572 C
3573          VERTEX = DCMPLX( DVERTX ) * (G31*G32)
3574       END IF
3575 C
3576       RETURN
3577       END
3578 C
3579 C ======================================================================
3580 C
3581       SUBROUTINE WWWWXX(WM1,WP1,WM2,WP2,GWWA,GWWZ,ZMASS,ZWIDTH , VERTEX)
3582 C
3583 C this subroutine computes an amplitude of the four-point w-/w+
3584 C coupling, including the contributions of photon and z exchanges.  the
3585 C photon propagator is given in feynman gauge and the z propagator is
3586 C given in unitary gauge.
3587 C
3588 C input:
3589 C       complex wm1(0:3)       : first  flow-out w-                  wm1
3590 C       complex wp1(0:3)       : first  flow-out w+                  wp1
3591 C       complex wm2(0:3)       : second flow-out w-                  wm2
3592 C       complex wp2(0:3)       : second flow-out w+                  wp2
3593 C       real    gwwa           : coupling constant of w and a       gwwa
3594 C       real    gwwz           : coupling constant of w and z       gwwz
3595 C       real    zmass          : mass  of z
3596 C       real    zwidth         : width of z
3597 C
3598 C output:
3599 C       complex vertex         : amplitude        gamma(wm1,wp1,wm2,wp2)
3600 C
3601 #if defined(CERNLIB_IMPNONE)
3602       IMPLICIT NONE
3603 #endif
3604       COMPLEX*16    WM1(6),WP1(6),WM2(6),WP2(6),VERTEX
3605       COMPLEX*16 DV1(0:3),DV2(0:3),DV3(0:3),DV4(0:3),
3606      &           J12(0:3),J34(0:3),J14(0:3),J32(0:3),DVERTX,
3607      &           SV1,SV2,SV3,SV4,TV1,TV2,TV3,TV4,DZS,DZT,
3608      &           V12,V13,V14,V23,V24,V34,JS12,JS34,JS14,JS32,JS,JT
3609       REAL*8       PWM1(0:3),PWP1(0:3),PWM2(0:3),PWP2(0:3),
3610      &           GWWA,GWWZ,ZMASS,ZWIDTH
3611       REAL*8     Q(0:3),K(0:3),DP1(0:3),DP2(0:3),DP3(0:3),DP4(0:3),
3612      &           DGWWA2,DGWWZ2,DGW2,DMZ,DWIDTH,S,T,DAS,DAT
3613 C
3614       REAL*8 RXZERO, RXONE, RXTWO
3615       PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
3616 C
3617       PWM1(0)=DBLE( WM1(5))
3618       PWM1(1)=DBLE( WM1(6))
3619       PWM1(2)=DIMAG(WM1(6))
3620       PWM1(3)=DIMAG(WM1(5))
3621       PWP1(0)=DBLE( WP1(5))
3622       PWP1(1)=DBLE( WP1(6))
3623       PWP1(2)=DIMAG(WP1(6))
3624       PWP1(3)=DIMAG(WP1(5))
3625       PWM2(0)=DBLE( WM2(5))
3626       PWM2(1)=DBLE( WM2(6))
3627       PWM2(2)=DIMAG(WM2(6))
3628       PWM2(3)=DIMAG(WM2(5))
3629       PWP2(0)=DBLE( WP2(5))
3630       PWP2(1)=DBLE( WP2(6))
3631       PWP2(2)=DIMAG(WP2(6))
3632       PWP2(3)=DIMAG(WP2(5))
3633 C
3634       DV1(0)=DCMPLX(WM1(1))
3635       DV1(1)=DCMPLX(WM1(2))
3636       DV1(2)=DCMPLX(WM1(3))
3637       DV1(3)=DCMPLX(WM1(4))
3638       DP1(0)=DBLE(PWM1(0))
3639       DP1(1)=DBLE(PWM1(1))
3640       DP1(2)=DBLE(PWM1(2))
3641       DP1(3)=DBLE(PWM1(3))
3642       DV2(0)=DCMPLX(WP1(1))
3643       DV2(1)=DCMPLX(WP1(2))
3644       DV2(2)=DCMPLX(WP1(3))
3645       DV2(3)=DCMPLX(WP1(4))
3646       DP2(0)=DBLE(PWP1(0))
3647       DP2(1)=DBLE(PWP1(1))
3648       DP2(2)=DBLE(PWP1(2))
3649       DP2(3)=DBLE(PWP1(3))
3650       DV3(0)=DCMPLX(WM2(1))
3651       DV3(1)=DCMPLX(WM2(2))
3652       DV3(2)=DCMPLX(WM2(3))
3653       DV3(3)=DCMPLX(WM2(4))
3654       DP3(0)=DBLE(PWM2(0))
3655       DP3(1)=DBLE(PWM2(1))
3656       DP3(2)=DBLE(PWM2(2))
3657       DP3(3)=DBLE(PWM2(3))
3658       DV4(0)=DCMPLX(WP2(1))
3659       DV4(1)=DCMPLX(WP2(2))
3660       DV4(2)=DCMPLX(WP2(3))
3661       DV4(3)=DCMPLX(WP2(4))
3662       DP4(0)=DBLE(PWP2(0))
3663       DP4(1)=DBLE(PWP2(1))
3664       DP4(2)=DBLE(PWP2(2))
3665       DP4(3)=DBLE(PWP2(3))
3666       DGWWA2=DBLE(GWWA)**2
3667       DGWWZ2=DBLE(GWWZ)**2
3668       DGW2  =DGWWA2+DGWWZ2
3669       DMZ   =DBLE(ZMASS)
3670       DWIDTH=DBLE(ZWIDTH)
3671 C
3672       V12= DV1(0)*DV2(0)-DV1(1)*DV2(1)-DV1(2)*DV2(2)-DV1(3)*DV2(3)
3673       V13= DV1(0)*DV3(0)-DV1(1)*DV3(1)-DV1(2)*DV3(2)-DV1(3)*DV3(3)
3674       V14= DV1(0)*DV4(0)-DV1(1)*DV4(1)-DV1(2)*DV4(2)-DV1(3)*DV4(3)
3675       V23= DV2(0)*DV3(0)-DV2(1)*DV3(1)-DV2(2)*DV3(2)-DV2(3)*DV3(3)
3676       V24= DV2(0)*DV4(0)-DV2(1)*DV4(1)-DV2(2)*DV4(2)-DV2(3)*DV4(3)
3677       V34= DV3(0)*DV4(0)-DV3(1)*DV4(1)-DV3(2)*DV4(2)-DV3(3)*DV4(3)
3678 C
3679       Q(0)=DP1(0)+DP2(0)
3680       Q(1)=DP1(1)+DP2(1)
3681       Q(2)=DP1(2)+DP2(2)
3682       Q(3)=DP1(3)+DP2(3)
3683       K(0)=DP1(0)+DP4(0)
3684       K(1)=DP1(1)+DP4(1)
3685       K(2)=DP1(2)+DP4(2)
3686       K(3)=DP1(3)+DP4(3)
3687 C
3688       S=Q(0)**2-Q(1)**2-Q(2)**2-Q(3)**2
3689       T=K(0)**2-K(1)**2-K(2)**2-K(3)**2
3690 C
3691       DAS=-RXONE/S
3692       DAT=-RXONE/T
3693       DZS=-RXONE/DCMPLX( S-DMZ**2 , DMAX1(DSIGN(DMZ*DWIDTH,S),RXZERO) )
3694       DZT=-RXONE/DCMPLX( T-DMZ**2 , DMAX1(DSIGN(DMZ*DWIDTH,T),RXZERO) )
3695 C
3696       SV1= (DP2(0)+Q(0))*DV1(0) -(DP2(1)+Q(1))*DV1(1)
3697      &    -(DP2(2)+Q(2))*DV1(2) -(DP2(3)+Q(3))*DV1(3)
3698       SV2=-(DP1(0)+Q(0))*DV2(0) +(DP1(1)+Q(1))*DV2(1)
3699      &    +(DP1(2)+Q(2))*DV2(2) +(DP1(3)+Q(3))*DV2(3)
3700       SV3= (DP4(0)-Q(0))*DV3(0) -(DP4(1)-Q(1))*DV3(1)
3701      &    -(DP4(2)-Q(2))*DV3(2) -(DP4(3)-Q(3))*DV3(3)
3702       SV4=-(DP3(0)-Q(0))*DV4(0) +(DP3(1)-Q(1))*DV4(1)
3703      &    +(DP3(2)-Q(2))*DV4(2) +(DP3(3)-Q(3))*DV4(3)
3704 C
3705       TV1= (DP4(0)+K(0))*DV1(0) -(DP4(1)+K(1))*DV1(1)
3706      &    -(DP4(2)+K(2))*DV1(2) -(DP4(3)+K(3))*DV1(3)
3707       TV2=-(DP3(0)-K(0))*DV2(0) +(DP3(1)-K(1))*DV2(1)
3708      &    +(DP3(2)-K(2))*DV2(2) +(DP3(3)-K(3))*DV2(3)
3709       TV3= (DP2(0)-K(0))*DV3(0) -(DP2(1)-K(1))*DV3(1)
3710      &    -(DP2(2)-K(2))*DV3(2) -(DP2(3)-K(3))*DV3(3)
3711       TV4=-(DP1(0)+K(0))*DV4(0) +(DP1(1)+K(1))*DV4(1)
3712      &    +(DP1(2)+K(2))*DV4(2) +(DP1(3)+K(3))*DV4(3)
3713 C
3714       J12(0)=(DP1(0)-DP2(0))*V12 +SV1*DV2(0) +SV2*DV1(0)
3715       J12(1)=(DP1(1)-DP2(1))*V12 +SV1*DV2(1) +SV2*DV1(1)
3716       J12(2)=(DP1(2)-DP2(2))*V12 +SV1*DV2(2) +SV2*DV1(2)
3717       J12(3)=(DP1(3)-DP2(3))*V12 +SV1*DV2(3) +SV2*DV1(3)
3718       J34(0)=(DP3(0)-DP4(0))*V34 +SV3*DV4(0) +SV4*DV3(0)
3719       J34(1)=(DP3(1)-DP4(1))*V34 +SV3*DV4(1) +SV4*DV3(1)
3720       J34(2)=(DP3(2)-DP4(2))*V34 +SV3*DV4(2) +SV4*DV3(2)
3721       J34(3)=(DP3(3)-DP4(3))*V34 +SV3*DV4(3) +SV4*DV3(3)
3722 C
3723       J14(0)=(DP1(0)-DP4(0))*V14 +TV1*DV4(0) +TV4*DV1(0)
3724       J14(1)=(DP1(1)-DP4(1))*V14 +TV1*DV4(1) +TV4*DV1(1)
3725       J14(2)=(DP1(2)-DP4(2))*V14 +TV1*DV4(2) +TV4*DV1(2)
3726       J14(3)=(DP1(3)-DP4(3))*V14 +TV1*DV4(3) +TV4*DV1(3)
3727       J32(0)=(DP3(0)-DP2(0))*V23 +TV3*DV2(0) +TV2*DV3(0)
3728       J32(1)=(DP3(1)-DP2(1))*V23 +TV3*DV2(1) +TV2*DV3(1)
3729       J32(2)=(DP3(2)-DP2(2))*V23 +TV3*DV2(2) +TV2*DV3(2)
3730       J32(3)=(DP3(3)-DP2(3))*V23 +TV3*DV2(3) +TV2*DV3(3)
3731 C
3732       JS12=Q(0)*J12(0)-Q(1)*J12(1)-Q(2)*J12(2)-Q(3)*J12(3)
3733       JS34=Q(0)*J34(0)-Q(1)*J34(1)-Q(2)*J34(2)-Q(3)*J34(3)
3734       JS14=K(0)*J14(0)-K(1)*J14(1)-K(2)*J14(2)-K(3)*J14(3)
3735       JS32=K(0)*J32(0)-K(1)*J32(1)-K(2)*J32(2)-K(3)*J32(3)
3736 C
3737       JS=J12(0)*J34(0)-J12(1)*J34(1)-J12(2)*J34(2)-J12(3)*J34(3)
3738       JT=J14(0)*J32(0)-J14(1)*J32(1)-J14(2)*J32(2)-J14(3)*J32(3)
3739 C
3740       DVERTX = (V12*V34 +V14*V23 -RXTWO*V13*V24)*DGW2
3741
3742 C     &        +(dzs*dgwwz2+das*dgwwa2)*js -dzs*dgwwz2*js12*js34/dmz**2
3743 C     &        +(dzt*dgwwz2+dat*dgwwa2)*jt -dzt*dgwwz2*js14*js32/dmz**2
3744 C
3745       VERTEX = -DCMPLX( DVERTX )
3746 C
3747       RETURN
3748       END