1 #include "isajet/pilot.h"
2 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 ***
11 C *********************************************************************
13 C Converted to double precision by W. Long and T. Seltzer for MadGraph.
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
19 C ======================================================================
21 SUBROUTINE BOOSTX(P,Q , PBOOST)
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
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
33 C real pboost(0:3) : four-momentum p in the boosted frame
35 #if defined(CERNLIB_IMPNONE)
38 REAL*8 P(0:3),Q(0:3),PBOOST(0:3),PQ,QQ,M,LF
40 PARAMETER( RXZERO=0.0D0 )
42 QQ=Q(1)**2+Q(2)**2+Q(3)**2
44 IF ( QQ .NE. RXZERO ) THEN
45 PQ=P(1)*Q(1)+P(2)*Q(2)+P(3)*Q(3)
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
62 C **********************************************************************
64 SUBROUTINE COUP1X(SW2 , GW,GWWA,GWWZ)
66 C this subroutine sets up the coupling constants of the gauge bosons in
70 C real sw2 : square of sine of the weak angle
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
77 #if defined(CERNLIB_IMPNONE)
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 )
86 C alpha = r_one / r_ialph
87 FOURPI = RXFOUR * RXPI
88 EE=SQRT( ALPHA * FOURPI )
90 CW=SQRT( RXONE - SW2 )
99 C ----------------------------------------------------------------------
101 SUBROUTINE COUP2X(SW2 , GAL,GAU,GAD,GWF,GZN,GZL,GZU,GZD,G1)
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.
109 C real sw2 : square of sine of the weak angle
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
122 #if defined(CERNLIB_IMPNONE)
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
128 REAL*8 RXZERO, RXHALF, RXONE, RXTWO, RTHREE, RXFOUR, RXOTE
130 PARAMETER( RXZERO=0.0D0, RXHALF=0.5D0, RXONE=1.0D0, RXTWO=2.0D0,
132 PARAMETER( RXFOUR=4.0D0, RXOTE=128.0D0 )
133 PARAMETER( RXPI=3.14159265358979323846D0, RIALPH=137.0359895D0 )
135 ALPHA = RXONE / RXOTE
136 C alpha = r_one / r_ialph
137 FOURPI = RXFOUR * RXPI
138 EE=SQRT( ALPHA * FOURPI )
140 CW=SQRT( RXONE - SW2 )
146 GAU(1) = -EE*RXTWO/RTHREE
147 GAU(2) = -EE*RXTWO/RTHREE
150 GWF(1) = -EE/SQRT(RXTWO*SW2)
154 GZL(1) = -EZ*(-RXHALF+SW2)
156 GZU(1) = -EZ*( RXHALF-SW2*RXTWO/RTHREE)
157 GZU(2) = EY* RXTWO/RTHREE
158 GZD(1) = -EZ*(-RXHALF+SW2 /RTHREE)
166 C ----------------------------------------------------------------------
168 SUBROUTINE COUP3X(SW2,ZMASS,HMASS ,
169 & GWWH,GZZH,GHHH,GWWHH,GZZHH,GHHHH)
171 C this subroutine sets up the coupling constants of the gauge bosons and
172 C higgs boson in the standard model.
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
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
187 #if defined(CERNLIB_IMPNONE)
190 REAL*8 SW2,ZMASS,HMASS,GWWH,GZZH,GHHH,GWWHH,GZZHH,GHHHH,
191 & ALPHA,FOURPI,EE2,SC2,V
193 REAL*8 RXHALF, RXONE, RXTWO, RTHREE, RXFOUR, RXOTE
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 )
199 ALPHA = RXONE / RXOTE
200 C alpha = r_one / r_ialph
201 FOURPI = RXFOUR * RXPI
203 SC2=SW2*( RXONE - SW2 )
204 V = RXTWO * ZMASS*SQRT(SC2)/SQRT(EE2)
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
216 C ----------------------------------------------------------------------
218 SUBROUTINE COUP4X(SW2,ZMASS,FMASS , GCHF)
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-
227 C real SW2 : square of sine of the weak angle
228 C real ZMASS : Z mass
229 C real FMASS : fermion mass
232 C complex GCHF(2) : coupling of fermion and Higgs
234 #if defined(CERNLIB_IMPNONE)
238 REAL*8 SW2,ZMASS,FMASS,ALPHA,FOURPI,EZ,G
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
246 GCHF(1) = DCMPLX( -G )
247 GCHF(2) = DCMPLX( -G )
252 C ======================================================================
254 SUBROUTINE EAIXXX(EB,EA,SHLF,CHLF,PHI,NHE,NHA , EAI)
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.
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
271 C complex EAI(6) : off-shell electron |e',A,e>
273 #if defined(CERNLIB_IMPNONE)
276 COMPLEX*16 EAI(6),PHS
277 REAL*8 EB,EA,SHLF,CHLF,PHI,ME,ALPHA,GAL,RNHE,X,C,S,D,COEFF,
283 GAL =SQRT(ALPHA*4.*3.14159265D0)
288 C=(CHLF+SHLF)*(CHLF-SHLF)
290 D=-1./(EA*EB*(4.*SHLF**2+(ME/EB)**2*C))
291 COEFF=-NN*GAL*SQRT(EB)*D
296 PHS=DCMPLX( CSP , RNHE*SNP )
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.
303 EAI(5) = EB*DCMPLX( 1.-X , 1.-X*C )
304 EAI(6) = -EB*X*S*DCMPLX( CSP , SNP )
309 C ----------------------------------------------------------------------
311 SUBROUTINE EAOXXX(EB,EA,SHLF,CHLF,PHI,NHE,NHA , EAO)
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.
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
328 C complex EAO(6) : off-shell positron <e,A,e'|
330 #if defined(CERNLIB_IMPNONE)
333 COMPLEX*16 EAO(6),PHS
334 REAL*8 EB,EA,SHLF,CHLF,PHI,ME,ALPHA,GAL,RNHE,X,C,S,D,COEFF,
340 GAL =SQRT(ALPHA*4.*3.14159265D0)
345 C=(CHLF+SHLF)*(CHLF-SHLF)
347 D=-1./(EA*EB*(4.*CHLF**2-(ME/EB)**2*C))
348 COEFF=NN*GAL*SQRT(EB)*D
353 PHS=DCMPLX( CSP ,-RNHE*SNP )
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.
360 EAO(5) = EB*DCMPLX( X-1. , X*C+1. )
361 EAO(6) = EB*X*S*DCMPLX( CSP , SNP )
366 C ----------------------------------------------------------------------
368 SUBROUTINE FSIXXX(FI,SC,GC,FMASS,FWIDTH , FSI)
370 C this subroutine computes an off-shell fermion wavefunction from a
371 C flowing-in external fermion and a vector boson.
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'
381 C complex fsi(6) : off-shell fermion |f',s,fi>
383 #if defined(CERNLIB_IMPNONE)
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
396 PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
398 DS=-SC(1)/DCMPLX(PF2-FMASS**2,MAX(DSIGN(FMASS*FWIDTH ,PF2),0D0))
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))
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
414 C ----------------------------------------------------------------------
416 SUBROUTINE FSOXXX(FO,SC,GC,FMASS,FWIDTH , FSO)
418 C this subroutine computes an off-shell fermion wavefunction from a
419 C flowing-out external fermion and a vector boson.
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'
429 C complex fso(6) : off-shell fermion <fo,s,f'|
431 #if defined(CERNLIB_IMPNONE)
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
444 PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
446 DS=-SC(1)/DCMPLX(PF2-FMASS**2,MAX(DSIGN(FMASS*FWIDTH ,PF2),0D0))
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))
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
462 C ----------------------------------------------------------------------
464 SUBROUTINE FVIXXX(FI,VC,G,FMASS,FWIDTH , FVI)
466 C this subroutine computes an off-shell fermion wavefunction from a
467 C flowing-in external fermion and a vector boson.
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'
477 C complex fvi(6) : off-shell fermion |f',v,fi>
479 #if defined(CERNLIB_IMPNONE)
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
485 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
492 C Fix compilation with g77
495 CXIMAG=DCMPLX( RXZERO, RXONE )
505 PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
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)
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)
519 FVI(1) = ( G(1)*((PF(0)-PF(3))*SL1 -DCONJG(FVI(6))*SL2)
521 FVI(2) = ( G(1)*( -FVI(6)*SL1 +(PF(0)+PF(3))*SL2)
523 FVI(3) = ( G(2)*((PF(0)+PF(3))*SR1 +DCONJG(FVI(6))*SR2)
525 FVI(4) = ( G(2)*( FVI(6)*SR1 +(PF(0)-PF(3))*SR2)
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
538 C ----------------------------------------------------------------------
540 SUBROUTINE FVOXXX(FO,VC,G,FMASS,FWIDTH , FVO)
542 C this subroutine computes an off-shell fermion wavefunction from a
543 C flowing-out external fermion and a vector boson.
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'
553 C complex fvo(6) : off-shell fermion <fo,v,f'|
555 #if defined(CERNLIB_IMPNONE)
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
561 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
567 C Fix compilation with g77
570 CXIMAG=DCMPLX( RXZERO, RXONE )
580 PF2=PF(0)**2-(PF(1)**2+PF(2)**2+PF(3)**2)
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)
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)
594 FVO(1) = ( G(2)*( (PF(0)+PF(3))*SR1 +FVO(6)*SR2)
596 FVO(2) = ( G(2)*( DCONJG(FVO(6))*SR1 +(PF(0)-PF(3))*SR2)
598 FVO(3) = ( G(1)*( (PF(0)-PF(3))*SL1 -FVO(6)*SL2)
600 FVO(4) = ( G(1)*(-DCONJG(FVO(6))*SL1 +(PF(0)+PF(3))*SL2)
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
613 C ----------------------------------------------------------------------
615 SUBROUTINE GGGGXX(WM,W31,WP,W32,G, VERTEX)
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).
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)
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.
645 C complex vertex : amplitude gamma(wm,w31,wp,w32)
647 #if defined(CERNLIB_IMPNONE)
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)
656 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
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))
683 DV2(0)=DCMPLX(W31(1))
684 DV2(1)=DCMPLX(W31(2))
685 DV2(2)=DCMPLX(W31(3))
686 DV2(3)=DCMPLX(W31(4))
699 DV4(0)=DCMPLX(W32(1))
700 DV4(1)=DCMPLX(W32(2))
701 DV4(2)=DCMPLX(W32(3))
702 DV4(3)=DCMPLX(W32(4))
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)
715 DVERTX = V14*V23 -V13*V24
717 VERTEX = DCMPLX( DVERTX ) * (G*G)
722 C ======================================================================
724 SUBROUTINE GGGXXX(WM,WP,W3,G , VERTEX)
726 C this subroutine computes an amplitude of the three-point coupling of
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
736 C complex vertex : amplitude gamma(wm,wp,w3)
738 #if defined(CERNLIB_IMPNONE)
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 )
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)
766 IF ( ABS(WM(1)) .NE. RXZERO ) THEN
767 IF (ABS(WM(1)).GE.MAX(ABS(WM(2)),ABS(WM(3)),ABS(WM(4)))
771 IF ( ABS(WP(1)) .NE. RXZERO) THEN
772 IF (ABS(WP(1)).GE.MAX(ABS(WP(2)),ABS(WP(3)),ABS(WP(4)))
776 IF ( ABS(W3(1)) .NE. RXZERO) THEN
777 IF ( ABS(W3(1)).GE.MAX(ABS(W3(2)),ABS(W3(3)),ABS(W3(4)))
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)
794 VERTEX = -(V12*(P13-P23)+V23*(P21-P31)+V31*(P32-P12))*G
798 SUBROUTINE HIOXXX(FI,FO,GC,SMASS,SWIDTH , HIO)
800 C this subroutine computes an off-shell scalar current from an external
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
811 C complex hio(3) : scalar current j(<fi|s|fo>)
813 #if defined(CERNLIB_IMPNONE)
816 COMPLEX*16 FI(6),FO(6),HIO(3),GC(2),DN
817 REAL*8 Q(0:3),SMASS,SWIDTH,Q2
826 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
828 DN=-DCMPLX(Q2-SMASS**2,DMAX1(DSIGN(SMASS*SWIDTH,Q2),0.D0))
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
835 C ----------------------------------------------------------------------
837 SUBROUTINE HSSSXX(S1,S2,S3,G,SMASS,SWIDTH , HSSS)
839 C This subroutine computes an off-shell scalar current from the four-
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'
851 C complex HSSS(3) : scalar current J(S':S1,S2,S3)
853 #if defined(CERNLIB_IMPNONE)
856 COMPLEX*16 S1(3),S2(3),S3(3),HSSS(3),DG
857 REAL*8 Q(0:3),G,SMASS,SWIDTH,Q2
859 HSSS(2) = S1(2)+S2(2)+S3(2)
860 HSSS(3) = S1(3)+S2(3)+S3(3)
866 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
868 DG=-G/DCMPLX( Q2-SMASS**2,MAX(SIGN(SMASS*SWIDTH ,Q2),0.D0))
870 HSSS(1) = DG * S1(1)*S2(1)*S3(1)
874 C ----------------------------------------------------------------------
876 SUBROUTINE HSSXXX(S1,S2,G,SMASS,SWIDTH , HSS)
878 C This subroutine computes an off-shell scalar current from the three-
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'
889 C complex HSS(3) : scalar current J(S':S1,S2)
891 #if defined(CERNLIB_IMPNONE)
894 COMPLEX*16 S1(3),S2(3),HSS(3),DG
895 REAL*8 Q(0:3),G,SMASS,SWIDTH,Q2
904 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
906 DG=-G/DCMPLX( Q2-SMASS**2, MAX(SIGN(SMASS*SWIDTH ,Q2),0.D0))
908 HSS(1) = DG*S1(1)*S2(1)
913 C ======================================================================
914 C ----------------------------------------------------------------------
916 SUBROUTINE HVSXXX(VC,SC,G,SMASS,SWIDTH , HVS)
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
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'
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.
945 C complex hvs(3) : scalar current j(s':v,s)
947 #if defined(CERNLIB_IMPNONE)
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
968 Q2=QA(0)**2-(QA(1)**2+QA(2)**2+QA(3)**2)
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)
974 HVS(1) = DG*(2D0*QPV+QVV)*SC(1)
979 C ----------------------------------------------------------------------
981 SUBROUTINE HVVXXX(V1,V2,G,SMASS,SWIDTH , HVV)
983 C this subroutine computes an off-shell scalar current from the vector-
984 C vector-scalar coupling.
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
994 C complex hvv(3) : off-shell scalar current j(s:v1,v2)
996 #if defined(CERNLIB_IMPNONE)
999 COMPLEX*16 V1(6),V2(6),HVV(3),DG
1000 REAL*8 Q(0:3),G,SMASS,SWIDTH,Q2
1002 PARAMETER( RXZERO=0.0D0 )
1004 HVV(2) = V1(5)+V2(5)
1005 HVV(3) = V1(6)+V2(6)
1011 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1013 DG=-G/DCMPLX( Q2-SMASS**2 , MAX(SIGN( SMASS*SWIDTH ,Q2),RXZERO) )
1015 HVV(1) = DG*(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
1020 C ======================================================================
1022 SUBROUTINE IOSXXX(FI,FO,SC,GC , VERTEX)
1024 C This subroutine computes an amplitude of the fermion-fermion-scalar
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
1034 C complex VERTEX : amplitude <FO|S|FI>
1036 #if defined(CERNLIB_IMPNONE)
1039 COMPLEX*16 FI(6),FO(6),SC(3),GC(2),VERTEX
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)) )
1047 C ======================================================================
1049 SUBROUTINE IOVXXX(FI,FO,VC,G , VERTEX)
1051 C this subroutine computes an amplitude of the fermion-fermion-vector
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
1061 C complex vertex : amplitude <fo|v|fi>
1063 #if defined(CERNLIB_IMPNONE)
1066 COMPLEX*16 FI(6),FO(6),VC(6),VERTEX
1068 REAL*8 RXZERO, RXONE
1069 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1075 C Fix compilation with g77
1078 CXIMAG=DCMPLX( RXZERO, RXONE )
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) )
1087 IF ( G(2) .NE. RXZERO ) THEN
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) )
1098 C Subroutine returns the desired fermion or
1099 C anti-fermion spinor. ie., |f>
1100 C A replacement for the HELAS routine IXXXXX
1102 C Adam Duff, 1992 August 31
1103 C <duff@phenom.physics.wisc.edu>
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
1112 C declare input/output variables
1114 #if defined(CERNLIB_IMPNONE)
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
1125 C declare local variables
1131 C Fix compilation with g77
1134 CXZERO=DCMPLX( RXZERO, RXZERO )
1137 C define kinematic parameters
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 )
1145 C do massive fermion case
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 )
1155 FI(3) = DCMPLX( OMEGAP, RXZERO )
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) )
1170 IF ( PLAT .EQ. RXZERO ) THEN
1172 FI(2) = DCMPLX( OMEGAM, RXZERO )
1174 FI(4) = DCMPLX( OMEGAP, RXZERO )
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) )
1188 ELSE IF ( NHEL .EQ. -1 ) THEN
1189 IF ( P(3) .GE. RXZERO ) THEN
1190 IF ( PLAT .EQ. RXZERO ) THEN
1192 FI(2) = DCMPLX( OMEGAP, RXZERO )
1194 FI(4) = DCMPLX( OMEGAM, RXZERO )
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 )
1208 IF ( PLAT .EQ. RXZERO ) THEN
1209 FI(1) = DCMPLX( -OMEGAP, RXZERO )
1211 FI(3) = DCMPLX( -OMEGAM, RXZERO )
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 )
1227 STOP 'IXXXXX: FERMION HELICITY MUST BE +1,-1'
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
1234 FI(2) = DCMPLX( -OMEGAP, RXZERO )
1236 FI(4) = DCMPLX( OMEGAM, RXZERO )
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 )
1250 IF ( PLAT .EQ. RXZERO ) THEN
1251 FI(1) = DCMPLX( OMEGAP, RXZERO )
1253 FI(3) = DCMPLX( -OMEGAM, RXZERO )
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 )
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 )
1273 FI(3) = DCMPLX( -OMEGAP, RXZERO )
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) )
1288 IF ( PLAT .EQ. RXZERO ) THEN
1290 FI(2) = DCMPLX( OMEGAM, RXZERO )
1292 FI(4) = DCMPLX( -OMEGAP, RXZERO )
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) )
1307 STOP 'IXXXXX: FERMION HELICITY MUST BE +1,-1'
1310 STOP 'IXXXXX: FERMION TYPE MUST BE +1,-1'
1313 C do massless fermion case
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
1322 FI(3) = DCMPLX( OMEGAP, RXZERO )
1325 SPAZ = SQRT( PABS + P(3) )
1328 FI(3) = DCMPLX( SPAZ, RXZERO )
1329 FI(4) = RXONE / SPAZ
1330 & * DCMPLX( P(1), P(2) )
1333 IF ( PLAT .EQ. RXZERO ) THEN
1337 FI(4) = DCMPLX( OMEGAP, RXZERO )
1339 SPAZ = SQRT( PABS - P(3) )
1342 FI(3) = RXONE / SPAZ
1343 & * DCMPLX( PLAT, RXZERO )
1345 & * DCMPLX( P(1), P(2) )
1348 ELSE IF ( NHEL .EQ. -1 ) THEN
1349 IF ( P(3) .GE. RXZERO ) THEN
1350 IF ( PLAT .EQ. RXZERO ) THEN
1352 FI(2) = DCMPLX( OMEGAP, RXZERO )
1356 SPAZ = SQRT( PABS + P(3) )
1357 FI(1) = RXONE / SPAZ
1358 & * DCMPLX( -P(1), P(2) )
1359 FI(2) = DCMPLX( SPAZ, RXZERO )
1364 IF ( PLAT .EQ. RXZERO ) THEN
1365 FI(1) = DCMPLX( -OMEGAP, RXZERO )
1370 SPAZ = SQRT( PABS - P(3) )
1372 & * DCMPLX( -P(1), P(2) )
1373 FI(2) = RXONE / SPAZ
1374 & * DCMPLX( PLAT, RXZERO )
1380 STOP 'IXXXXX: FERMION HELICITY MUST BE +1,-1'
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
1387 FI(2) = DCMPLX( -OMEGAP, RXZERO )
1391 SPAZ = SQRT( PABS + P(3) )
1392 FI(1) = -RXONE / SPAZ
1393 & * DCMPLX( -P(1), P(2) )
1394 FI(2) = DCMPLX( -SPAZ, RXZERO )
1399 IF ( PLAT .EQ. RXZERO ) THEN
1400 FI(1) = DCMPLX( OMEGAP, RXZERO )
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 )
1414 ELSE IF ( NHEL .EQ. -1 ) THEN
1415 IF ( P(3) .GE. RXZERO ) THEN
1416 IF ( PLAT .EQ. RXZERO ) THEN
1419 FI(3) = DCMPLX( -OMEGAP, RXZERO )
1422 SPAZ = SQRT( PABS + P(3) )
1425 FI(3) = DCMPLX( -SPAZ, RXZERO )
1426 FI(4) = -RXONE / SPAZ
1427 & * DCMPLX( P(1), P(2) )
1430 IF ( PLAT .EQ. RXZERO ) THEN
1434 FI(4) = DCMPLX( -OMEGAP, RXZERO )
1436 SPAZ = SQRT( PABS - P(3) )
1439 FI(3) = -RXONE / SPAZ
1440 & * DCMPLX( PLAT, RXZERO )
1441 FI(4) = -SPAZ / PLAT
1442 & * DCMPLX( P(1), P(2) )
1446 STOP 'IXXXXX: FERMION HELICITY MUST BE +1,-1'
1449 STOP 'IXXXXX: FERMION TYPE MUST BE +1,-1'
1458 C ----------------------------------------------------------------------
1460 SUBROUTINE J3XXXX(FI,FO,GAF,GZF,ZMASS,ZWIDTH , J3)
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.
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
1477 C complex j3(6) : w3 current j^mu(<fo|w3|fi>)
1479 #if defined(CERNLIB_IMPNONE)
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
1487 REAL*8 RXZERO, RXONE
1488 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1494 C Fix compilation with g77
1497 CXIMAG=DCMPLX( RXZERO, RXONE )
1507 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1512 WW=MAX(DSIGN( ZMW ,Q2),RXZERO)
1513 DZ=RXONE/DCMPLX( Q2-ZM2 , WW )
1514 DDIF=DCMPLX( -ZM2 , WW )*DA*DZ
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
1522 CW=RXONE/SQRT(RXONE+(GZF(2)/GAF(2))**2)
1523 SW=SQRT((RXONE-CW)*(RXONE+CW))
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
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)
1550 C ----------------------------------------------------------------------
1552 SUBROUTINE JEEXXX(EB,EF,SHLF,CHLF,PHI,NHB,NHF,NSF , JEE)
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.
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
1570 C complex JEE(6) : off-shell photon J^mu(<e|A|e>)
1572 #if defined(CERNLIB_IMPNONE)
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
1582 GAL =SQRT(ALPHA*4.*3.14159265D0)
1589 C CS(1)=CHLF and CS(2)=SHLF for electron
1590 C CS(1)=SHLF and CS(2)=CHLF for positron
1593 Q2=-4.*CS(2)**2*(EF*EB-ME2)
1594 & +SF*(1.-X)**2/X*(SHLF+CHLF)*(SHLF-CHLF)*ME2
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))
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)
1617 C=(CHLF+SHLF)*(CHLF-SHLF)
1620 JEE(5) = -EB*DCMPLX( 1.-X , SF-X*C )
1621 JEE(6) = EB*X*S*DCMPLX( CSP , SNP )
1627 C ----------------------------------------------------------------------
1629 SUBROUTINE JGGGXX(W1,W2,W3,G, JW3W)
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
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)
1646 C complex jw3w(6) : w current j^mu(w':w1,w2,w3)
1648 #if defined(CERNLIB_IMPNONE)
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
1657 PARAMETER( RXZERO=0.0D0 )
1659 JW3W(5) = W1(5)+W2(5)+W3(5)
1660 JW3W(6) = W1(6)+W2(6)+W3(6)
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))
1676 P1(2)=DBLE(DIMAG(W1(6)))
1677 P1(3)=DBLE(DIMAG(W1(5)))
1680 P2(2)=DBLE(DIMAG(W2(6)))
1681 P2(3)=DBLE(DIMAG(W2(5)))
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))
1691 Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
1695 DV = 1.0D0/DCMPLX( Q2 )
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) )
1700 W32=DW3(0)*DW2(0)-DW3(1)*DW2(1)-DW3(2)*DW2(2)-DW3(3)*DW2(3)
1703 W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
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 )
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 )
1718 C ----------------------------------------------------------------------
1720 SUBROUTINE JGGXXX(V1,V2,G, JVV)
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.
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)
1732 C complex jvv(6) : vector current j^mu(v:v1,v2)
1734 #if defined(CERNLIB_IMPNONE)
1737 COMPLEX*16 V1(6),V2(6),JVV(6),J12(0:3),
1739 REAL*8 P1(0:3),P2(0:3),Q(0:3),G,GS,S
1742 PARAMETER( RXZERO=0.0D0 )
1744 JVV(5) = V1(5)+V2(5)
1745 JVV(6) = V1(6)+V2(6)
1759 S=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
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)
1781 C ----------------------------------------------------------------------
1783 SUBROUTINE JIOXXX(FI,FO,G,VMASS,VWIDTH , JIO)
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.
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
1797 C complex jio(6) : vector current j^mu(<fo|v|fi>)
1799 #if defined(CERNLIB_IMPNONE)
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
1805 REAL*8 RXZERO, RXONE
1806 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
1812 C Fix compilation with g77
1815 CXIMAG=DCMPLX( RXZERO, RXONE )
1818 JIO(5) = FO(5)-FI(5)
1819 JIO(6) = FO(6)-FI(6)
1825 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1828 IF (VMASS.NE.RXZERO) THEN
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) )
1834 IF (G(2).NE.RXZERO) THEN
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))
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)
1853 CS=(Q(0)*C0-Q(1)*C1-Q(2)*C2-Q(3)*C3)/VM2
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
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
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
1886 C ----------------------------------------------------------------------
1888 SUBROUTINE JSSXXX(S1,S2,G,VMASS,VWIDTH , JSS)
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.
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
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.
1918 C complex JSS(6) : vector current J^mu(V:S1,S2)
1920 #if defined(CERNLIB_IMPNONE)
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
1926 JSS(5) = S1(2)+S2(2)
1927 JSS(6) = S1(3)+S2(3)
1933 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
1936 IF (VMASS.EQ.0.) GOTO 10
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.) )
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)
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)
1963 10 ADG=G*S1(1)*S2(1)/Q2
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))
1974 C ----------------------------------------------------------------------
1976 SUBROUTINE JTIOXX(FI,FO,G , JIO)
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
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
1988 C complex jio(6) : vector current j^mu(<fo|v|fi>)
1990 #if defined(CERNLIB_IMPNONE)
1993 COMPLEX*16 FI(6),FO(6),JIO(6)
1996 REAL*8 RXZERO, RXONE
1997 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
2003 C Fix compilation with g77
2006 CXIMAG=DCMPLX( RXZERO, RXONE )
2009 JIO(5) = FO(5)-FI(5)
2010 JIO(6) = FO(6)-FI(6)
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)) )
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)
2031 C ----------------------------------------------------------------------
2033 SUBROUTINE JVSSXX(VC,S1,S2,G,VMASS,VWIDTH , JVSS)
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
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'
2049 C complex JVSS(6) : vector current J^mu(V':V,S1,S2)
2051 #if defined(CERNLIB_IMPNONE)
2054 COMPLEX*16 VC(6),S1(3),S2(3),JVSS(6),DG
2055 REAL*8 Q(0:3),G,VMASS,VWIDTH,Q2,VK,VM2
2057 JVSS(5) = VC(5)+S1(2)+S2(2)
2058 JVSS(6) = VC(6)+S1(3)+S2(3)
2064 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
2067 IF (VMASS.EQ.0.) GOTO 10
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.))
2073 VK=(Q(0)*VC(1)-Q(1)*VC(2)-Q(2)*VC(3)-Q(3)*VC(4))/VM2
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))
2082 10 DG= G*S1(1)*S2(1)/Q2
2093 C ----------------------------------------------------------------------
2095 SUBROUTINE JVSXXX(VC,SC,G,VMASS,VWIDTH , JVS)
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.
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'
2109 C complex jvs(6) : vector current j^mu(v':v,s)
2111 #if defined(CERNLIB_IMPNONE)
2114 COMPLEX*16 VC(6),SC(3),JVS(6),DG,VK
2115 REAL*8 Q(0:3),VMASS,VWIDTH,Q2,VM2,G
2117 JVS(5) = VC(5)+SC(2)
2118 JVS(6) = VC(6)+SC(3)
2124 Q2=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
2127 IF (VMASS.EQ.0.) GOTO 10
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.) )
2133 VK=(-Q(0)*VC(1)+Q(1)*VC(2)+Q(2)*VC(3)+Q(3)*VC(4))/VM2
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))
2154 C ----------------------------------------------------------------------
2156 SUBROUTINE JVVXXX(V1,V2,G,VMASS,VWIDTH , JVV)
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.
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
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.
2180 C complex jvv(6) : vector current j^mu(v:v1,v2)
2182 #if defined(CERNLIB_IMPNONE)
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
2190 PARAMETER( RXZERO=0.0D0 )
2192 JVV(5) = V1(5)+V2(5)
2193 JVV(6) = V1(6)+V2(6)
2207 S=Q(0)**2-(Q(1)**2+Q(2)**2+Q(3)**2)
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)
2219 IF ( VMASS .NE. RXZERO ) THEN
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
2229 DG=-G/DCMPLX( S-VM2 , MAX(SIGN( VMASS*VWIDTH ,S),RXZERO) )
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) )
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)
2251 C ----------------------------------------------------------------------
2253 SUBROUTINE JW3WXX(W1,W2,W3,G1,G2,WMASS,WWIDTH,VMASS,VWIDTH , JW3W)
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
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'
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.
2301 C complex jw3w(6) : w current j^mu(w':w1,w2,w3)
2303 #if defined(CERNLIB_IMPNONE)
2306 COMPLEX*16 W1(6),W2(6),W3(6),JW3W(6)
2307 COMPLEX*16 DW1(0:3),DW2(0:3),DW3(0:3),
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
2316 PARAMETER( RXZERO=0.0D0 )
2318 JW3W(5) = W1(5)+W2(5)+W3(5)
2319 JW3W(6) = W1(6)+W2(6)+W3(6)
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))
2335 P1(2)=DBLE(DIMAG(W1(6)))
2336 P1(3)=DBLE(DIMAG(W1(5)))
2339 P2(2)=DBLE(DIMAG(W2(6)))
2340 P2(3)=DBLE(DIMAG(W2(5)))
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))
2351 Q2 =Q(0)**2 -(Q(1)**2 +Q(2)**2 +Q(3)**2)
2352 DG2=DBLE(G1)*DBLE(G2)
2356 IF (VMASS.EQ. RXZERO) THEN
2357 DV = 1.0D0/DCMPLX( Q2 )
2359 DV = 1.0D0/DCMPLX( Q2 -MV2 , DMAX1(DSIGN(DMV*DWV,Q2 ),0.D0) )
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) )
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)
2367 IF ( WMASS .NE. RXZERO ) THEN
2368 W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
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 )
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)
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 )
2398 IF ( VMASS .NE. RXZERO ) THEN
2400 JQ=(JJ(0)*Q(0)-JJ(1)*Q(1)-JJ(2)*Q(2)-JJ(3)*Q(3))/MV2
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 )
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 )
2418 C ----------------------------------------------------------------------
2420 SUBROUTINE JWWWXX(W1,W2,W3,GWWA,GWWZ,ZMASS,ZWIDTH,WMASS,WWIDTH ,
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.
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
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.
2450 C complex jwww(6) : w current j^mu(w':w1,w2,w3)
2452 #if defined(CERNLIB_IMPNONE)
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,
2466 JWWW(5) = W1(5)+W2(5)+W3(5)
2467 JWWW(6) = W1(6)+W2(6)+W3(6)
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))
2483 P1(2)=DBLE(DIMAG(W1(6)))
2484 P1(3)=DBLE(DIMAG(W1(5)))
2487 P2(2)=DBLE(DIMAG(W2(6)))
2488 P2(3)=DBLE(DIMAG(W2(5)))
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))
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
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) )
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)
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)
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)
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
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
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)
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)
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)
2580 W13=DW1(0)*DW3(0)-DW1(1)*DW3(1)-DW1(2)*DW3(2)-DW1(3)*DW3(3)
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 )
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)
2597 JQ=(JJ(0)*Q(0)-JJ(1)*Q(1)-JJ(2)*Q(2)-JJ(3)*Q(3))/MW2
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 )
2609 C ----------------------------------------------------------------------
2611 SUBROUTINE MOM2CX(ESUM,MASS1,MASS2,COSTH1,PHI1 , P1,P2)
2613 C This subroutine sets up two four-momenta in the two particle rest
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
2624 C real P1(0:3) : four-momentum of particle 1
2625 C real P2(0:3) : four-momentum of particle 2
2627 #if defined(CERNLIB_IMPNONE)
2630 REAL*8 P1(0:3),P2(0:3),
2631 & ESUM,MASS1,MASS2,COSTH1,PHI1,MD2,ED,PP,SINTH1
2633 MD2=(MASS1-MASS2)*(MASS1+MASS2)
2635 IF (MASS1*MASS2.EQ.0.) THEN
2636 PP=(ESUM-ABS(ED))*0.5D0
2639 PP=SQRT((MD2/ESUM)**2-2.0D0*(MASS1**2+MASS2**2)+ESUM**2)*0.5D0
2641 SINTH1=SQRT((1.0D0-COSTH1)*(1.0D0+COSTH1))
2643 P1(0) = MAX((ESUM+ED)*0.5D0,0.D0)
2644 P1(1) = PP*SINTH1*COS(PHI1)
2645 P1(2) = PP*SINTH1*SIN(PHI1)
2648 P2(0) = MAX((ESUM-ED)*0.5D0,0.D0)
2655 C **********************************************************************
2657 SUBROUTINE MOMNTX(ENERGY,MASS,COSTH,PHI , P)
2659 C This subroutine sets up a four-momentum from the four inputs.
2662 C real ENERGY : energy
2664 C real COSTH : cos(theta)
2665 C real PHI : azimuthal angle
2668 C real P(0:3) : four-momentum
2670 #if defined(CERNLIB_IMPNONE)
2673 REAL*8 P(0:3),ENERGY,MASS,COSTH,PHI,PP,SINTH
2676 IF (ENERGY.EQ.MASS) THEN
2681 PP=SQRT((ENERGY-MASS)*(ENERGY+MASS))
2682 SINTH=SQRT((1.-COSTH)*(1.+COSTH))
2688 P(1) = PP*SINTH*COS(PHI)
2689 P(2) = PP*SINTH*SIN(PHI)
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
2701 C Adam Duff, 1992 August 31
2702 C <duff@phenom.physics.wisc.edu>
2704 SUBROUTINE OXXXXX(P,FMASS,NHEL,NSF,FO)
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
2712 C declare input/output variables
2714 #if defined(CERNLIB_IMPNONE)
2719 REAL*8 P(0:3), FMASS
2721 C declare local variables
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
2731 C Fix compilation with g77
2734 CXZERO=DCMPLX( RXZERO, RXZERO )
2737 C define kinematic parameters
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 )
2745 C do massive fermion case
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 )
2755 FO(3) = DCMPLX( OMEGAM, RXZERO )
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) )
2770 IF ( PLAT .EQ. RXZERO ) THEN
2772 FO(2) = DCMPLX( OMEGAP, RXZERO )
2774 FO(4) = DCMPLX( OMEGAM, RXZERO )
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) )
2788 ELSE IF ( NHEL .EQ. -1 ) THEN
2789 IF ( P(3) .GE. RXZERO ) THEN
2790 IF ( PLAT .EQ. RXZERO ) THEN
2792 FO(2) = DCMPLX( OMEGAM, RXZERO )
2794 FO(4) = DCMPLX( OMEGAP, RXZERO )
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 )
2808 IF ( PLAT .EQ. RXZERO ) THEN
2809 FO(1) = DCMPLX( -OMEGAM, RXZERO )
2811 FO(3) = DCMPLX( -OMEGAP, RXZERO )
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 )
2827 STOP 'OXXXXX: FERMION HELICITY MUST BE +1,-1'
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
2834 FO(2) = DCMPLX( OMEGAM, RXZERO )
2836 FO(4) = DCMPLX( -OMEGAP, RXZERO )
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 )
2850 IF ( PLAT .EQ. RXZERO ) THEN
2851 FO(1) = DCMPLX( -OMEGAM, RXZERO )
2853 FO(3) = DCMPLX( OMEGAP, RXZERO )
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 )
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 )
2873 FO(3) = DCMPLX( OMEGAM, RXZERO )
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) )
2888 IF ( PLAT .EQ. RXZERO ) THEN
2890 FO(2) = DCMPLX( -OMEGAP, RXZERO )
2892 FO(4) = DCMPLX( OMEGAM, RXZERO )
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) )
2907 STOP 'OXXXXX: FERMION HELICITY MUST BE +1,-1'
2910 STOP 'OXXXXX: FERMION TYPE MUST BE +1,-1'
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 )
2925 SPAZ = SQRT( PABS + P(3) )
2926 FO(1) = DCMPLX( SPAZ, RXZERO )
2927 FO(2) = RXONE / SPAZ
2928 & * DCMPLX( P(1), -P(2) )
2933 IF ( PLAT .EQ. RXZERO ) THEN
2935 FO(2) = DCMPLX( OMEGAP, RXZERO )
2939 SPAZ = SQRT( PABS - P(3) )
2940 FO(1) = RXONE / SPAZ
2941 & * DCMPLX( PLAT, RXZERO )
2943 & * DCMPLX( P(1), -P(2) )
2948 ELSE IF ( NHEL .EQ. -1 ) THEN
2949 IF ( P(3) .GE. RXZERO ) THEN
2950 IF ( PLAT .EQ. RXZERO ) THEN
2954 FO(4) = DCMPLX( OMEGAP, RXZERO )
2956 SPAZ = SQRT( PABS + P(3) )
2959 FO(3) = RXONE / SPAZ
2960 & * DCMPLX( -P(1), -P(2) )
2961 FO(4) = DCMPLX( SPAZ, RXZERO )
2964 IF ( PLAT .EQ. RXZERO ) THEN
2967 FO(3) = DCMPLX( -OMEGAP, RXZERO )
2970 SPAZ = SQRT( PABS - P(3) )
2974 & * DCMPLX( -P(1), -P(2) )
2975 FO(4) = RXONE / SPAZ
2976 & * DCMPLX( PLAT, RXZERO )
2980 STOP 'OXXXXX: FERMION HELICITY MUST BE +1,-1'
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
2989 FO(4) = DCMPLX( -OMEGAP, RXZERO )
2991 SPAZ = SQRT( PABS + P(3) )
2994 FO(3) = -RXONE / SPAZ
2995 & * DCMPLX( -P(1), -P(2) )
2996 FO(4) = DCMPLX( -SPAZ, RXZERO )
2999 IF ( PLAT .EQ. RXZERO ) THEN
3002 FO(3) = DCMPLX( OMEGAP, RXZERO )
3005 SPAZ = SQRT( PABS - P(3) )
3008 FO(3) = -SPAZ / PLAT
3009 & * DCMPLX( -P(1), -P(2) )
3010 FO(4) = -RXONE / SPAZ
3011 & * DCMPLX( PLAT, RXZERO )
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 )
3022 SPAZ = SQRT( PABS + P(3) )
3023 FO(1) = DCMPLX( -SPAZ, RXZERO )
3024 FO(2) = -RXONE / SPAZ
3025 & * DCMPLX( P(1), -P(2) )
3030 IF ( PLAT .EQ. RXZERO ) THEN
3032 FO(2) = DCMPLX( -OMEGAP, RXZERO )
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) )
3046 STOP 'OXXXXX: FERMION HELICITY MUST BE +1,-1'
3049 STOP 'OXXXXX: FERMION TYPE MUST BE +1,-1'
3058 C ----------------------------------------------------------------------
3060 SUBROUTINE ROTXXX(P,Q , PROT)
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.
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
3072 C real prot(0:3) : four-momentum p in the rotated frame
3074 #if defined(CERNLIB_IMPNONE)
3077 REAL*8 P(0:3),Q(0:3),PROT(0:3),QT2,QT,PSGN,QQ,P1
3079 REAL*8 RXZERO, RXONE
3080 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
3086 IF ( QT2 .EQ. RXZERO ) THEN
3087 IF ( Q(3) .EQ. RXZERO ) THEN
3092 PSGN=DSIGN(RXONE,Q(3))
3098 QQ=SQRT(QT2+Q(3)**2)
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)
3108 C ======================================================================
3110 SUBROUTINE SSSSXX(S1,S2,S3,S4,G , VERTEX)
3112 C This subroutine computes an amplitude of the four-scalar coupling.
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
3122 C complex VERTEX : amplitude Gamma(S1,S2,S3,S4)
3124 #if defined(CERNLIB_IMPNONE)
3127 COMPLEX*16 S1(3),S2(3),S3(3),S4(3),VERTEX
3130 VERTEX = G*S1(1)*S2(1)*S3(1)*S4(1)
3135 C ======================================================================
3137 SUBROUTINE SSSXXX(S1,S2,S3,G , VERTEX)
3139 C This subroutine computes an amplitude of the three-scalar coupling.
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
3148 C complex VERTEX : amplitude Gamma(S1,S2,S3)
3150 #if defined(CERNLIB_IMPNONE)
3153 COMPLEX*16 S1(3),S2(3),S3(3),VERTEX
3156 VERTEX = G*S1(1)*S2(1)*S3(1)
3162 C ----------------------------------------------------------------------
3164 SUBROUTINE SXXXXX(P,NSS , SC)
3166 C This subroutine computes a complex SCALAR wavefunction.
3169 C real P(0:3) : four-momentum of scalar boson
3170 C integer NSS = -1 or 1 : +1 for final, -1 for initial
3173 C complex SC(3) : scalar wavefunction S
3175 #if defined(CERNLIB_IMPNONE)
3182 SC(1) = DCMPLX( 1.0 )
3183 SC(2) = DCMPLX(P(0),P(3))*NSS
3184 SC(3) = DCMPLX(P(1),P(2))*NSS
3189 C ======================================================================
3191 SUBROUTINE VSSXXX(VC,S1,S2,G , VERTEX)
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.
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)
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.
3217 C complex vertex : amplitude gamma(v,s1,s2)
3219 #if defined(CERNLIB_IMPNONE)
3222 COMPLEX*16 VC(6),S1(3),S2(3),VERTEX,G
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))
3230 VERTEX = G*S1(1)*S2(1)
3231 & *(VC(1)*P(0)-VC(2)*P(1)-VC(3)*P(2)-VC(4)*P(3))
3236 SUBROUTINE VVSSXX(V1,V2,S1,S2,G , VERTEX)
3238 C This subroutine computes an amplitude of the vector-vector-scalar-
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
3249 C complex VERTEX : amplitude Gamma(V1,V2,S1,S2)
3251 #if defined(CERNLIB_IMPNONE)
3254 COMPLEX*16 V1(6),V2(6),S1(3),S2(3),VERTEX
3257 VERTEX = G*S1(1)*S2(1)
3258 & *(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
3264 C ======================================================================
3266 SUBROUTINE VVSXXX(V1,V2,SC,G , VERTEX)
3268 C this subroutine computes an amplitude of the vector-vector-scalar
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
3278 C complex vertex : amplitude gamma(v1,v2,s)
3280 #if defined(CERNLIB_IMPNONE)
3283 COMPLEX*16 V1(6),V2(6),SC(3),VERTEX
3286 VERTEX = G*SC(1)*(V1(1)*V2(1)-V1(2)*V2(2)-V1(3)*V2(3)-V1(4)*V2(4))
3291 C ======================================================================
3293 SUBROUTINE VVVXXX(WM,WP,W3,G , VERTEX)
3295 C this subroutine computes an amplitude of the three-point coupling of
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
3305 C complex vertex : amplitude gamma(wm,wp,w3)
3307 #if defined(CERNLIB_IMPNONE)
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
3314 REAL*8 RXZERO, RTENTH
3315 PARAMETER( RXZERO=0.0D0, RTENTH=0.1D0 )
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)
3336 IF ( ABS(WM(1)) .NE. RXZERO ) THEN
3337 IF (ABS(WM(1)).GE.MAX(ABS(WM(2)),ABS(WM(3)),ABS(WM(4)))
3341 IF ( ABS(WP(1)) .NE. RXZERO) THEN
3342 IF (ABS(WP(1)).GE.MAX(ABS(WP(2)),ABS(WP(3)),ABS(WP(4)))
3346 IF ( ABS(W3(1)) .NE. RXZERO) THEN
3347 IF ( ABS(W3(1)).GE.MAX(ABS(W3(2)),ABS(W3(3)),ABS(W3(4)))
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)
3364 VERTEX = -(V12*(P13-P23)+V23*(P21-P31)+V31*(P32-P12))*G
3370 C Subroutine returns the value of evaluated
3371 C helicity basis boson polarisation wavefunction.
3372 C Replaces the HELAS routine VXXXXX
3374 C Adam Duff, 1992 September 3
3375 C <duff@phenom.physics.wisc.edu>
3377 SUBROUTINE VXXXXX(P,VMASS,NHEL,NSV,VC)
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
3385 C declare input/output variables
3387 #if defined(CERNLIB_IMPNONE)
3392 REAL*8 P(0:3), VMASS
3394 C declare local variables
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
3404 C Fix compilation with g77
3407 CXZERO=DCMPLX( RXZERO, RXZERO )
3410 C define internal/external momenta
3412 IF ( NSV**2 .NE. 1 ) THEN
3413 STOP 'VXXXXX: NSV IS NOT ONE OF -1, +1'
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 )
3422 C calculate polarisation four vectors
3424 IF ( NHEL**2 .EQ. 1 ) THEN
3425 IF ( (PABS .EQ. RXZERO) .OR. (PLAT .EQ. RXZERO) ) THEN
3427 VC(2) = DCMPLX( -NHEL * RS2 * DSIGN( RXONE, P(3) ), RXZERO )
3428 VC(3) = DCMPLX( RXZERO, NSV * RS2 )
3431 RPLAT = RXONE / PLAT
3432 RPABS = RXONE / PABS
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,
3441 ELSE IF ( NHEL .EQ. 0 ) THEN
3442 IF ( VMASS .GT. RXZERO ) THEN
3443 IF ( PABS .EQ. RXZERO ) THEN
3447 VC(4) = DCMPLX( RXONE, RXZERO )
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 )
3456 STOP 'VXXXXX: NHEL = 0 IS ONLY FOR MASSIVE BOSONS'
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
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 )
3472 STOP 'VXXXXX: NHEL = 4 IS ONLY FOR M>=0'
3475 STOP 'VXXXXX: NHEL IS NOT ONE OF -1, 0, 1 OR 4'
3483 C ----------------------------------------------------------------------
3485 SUBROUTINE W3W3XX(WM,W31,WP,W32,G31,G32,WMASS,WWIDTH , VERTEX)
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).
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
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.
3518 C complex vertex : amplitude gamma(wm,w31,wp,w32)
3520 #if defined(CERNLIB_IMPNONE)
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
3528 REAL*8 RXZERO, RXONE
3529 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0 )
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))
3548 IF ( DBLE(WMASS) .NE. RXZERO ) THEN
3549 C dm2inv = r_one / dmw2
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)
3558 DVERTX = V12*V34 +V14*V23 -2.D0*V13*V24
3560 VERTEX = DCMPLX( DVERTX ) * (G31*G32)
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)
3571 DVERTX = V14*V23 -V13*V24
3573 VERTEX = DCMPLX( DVERTX ) * (G31*G32)
3579 C ======================================================================
3581 SUBROUTINE WWWWXX(WM1,WP1,WM2,WP2,GWWA,GWWZ,ZMASS,ZWIDTH , VERTEX)
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.
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
3599 C complex vertex : amplitude gamma(wm1,wp1,wm2,wp2)
3601 #if defined(CERNLIB_IMPNONE)
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
3614 REAL*8 RXZERO, RXONE, RXTWO
3615 PARAMETER( RXZERO=0.0D0, RXONE=1.0D0, RXTWO=2.0D0 )
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))
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
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)
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
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) )
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)
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)
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)
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)
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)
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)
3740 DVERTX = (V12*V34 +V14*V23 -RXTWO*V13*V24)*DGW2
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
3745 VERTEX = -DCMPLX( DVERTX )