
      real function WmassWeiKWCC03(Wmass,Wwidth)
***************************************************************************
* Function to reweight the CC03 KORELW MC as function of the W-mass and Wwidth.
* This function has to be called from your ALPHA job. The nominal Wmass and
* Wwidth, which
* was used to generate the KORELW MC sample is determined in this routine.
* The routine takes also care of Mass dependend SM-parameters.
*
* A standalone version is also available: 
*                         WmassWeiCC03(Wmass,Wnom,Pxq,Pyq,Pzq,Eq)
*
* input:  real Wmass
*         real Wwidth   (if <0 then Wwidth is calculated automatically)
* output: real function WmassWeiKWCC03    (Weight)
*
* Author: Marcel Vreeswijk, Jan. 1997
****************************************************************************
      IMPLICIT NONE
#include "qdecl.h"
#include "qcde.h"
#include "qhac.h"

      real Eq(4),Pxq(4),Pyq(4),Pzq(4),rpart(4)
      real mult,del
      INTEGER IKPAR,nlink
      integer iwp,iwm,id,i,idf,idl,
     +        icom2,icop2,icom,icop
      real Wmassnom,Wmass,Wwidth,Wwidthnom,xsec,xnom
      double precision born4m,bpp1(4),bpp2(4),bpp3(4),bpp4(4)
#include "qmacro.h"




      do i=1,4
        Eq(i)=0
        Pxq(i)=0
        Pyq(i)=0
        Pzq(i)=0
        rpart(i)=0
      enddo

      icom2=0
      icop2=0
      icom=0
      icop=0
      mult=1

      IWP = KPDIR('W+',KMONTE)
      if (iwp.le.0) goto 999

      do 20 i=1,KNDAU(IWP)
        id=kdau(IWP,i)
        if (KTPCOD(id).ne.1) then
          if (KTPCOD(id).lt.100) then ! leptonic
            if (QCH(id).ne.0.) then
              idf=id
              icop=KTPCOD(id)
            else
              idl=id
              icop2=KTPCOD(id)
            endif
          else                  ! hadronic
            del=abs(abs(QCH(id))-0.333333)
            if (del.lt.0.01) then ! charge=1/3= type
              idf=id
              icop=KTPCOD(id)
            else
              idl=id
              icop2=KTPCOD(id)
            endif
          endif
        else                    ! init. stat rad
          mult=-1
        endif
 20   continue
      icop=icop*mult
      icop2=icop2*mult
      Eq(1)=QE(idf)
      Pxq(1)=Qx(idf)
      Pyq(1)=Qy(idf)
      Pzq(1)=Qz(idf)
      Eq(2)=QE(idl)
      Pxq(2)=Qx(idl)
      Pyq(2)=Qy(idl)
      Pzq(2)=Qz(idl)
      rpart(1)=icop
      rpart(2)=icop2


 99   continue

      IWM = KPDIR('W-',KMONTE)
      if (iwm.le.0) goto 999

      mult=1
      do 10 i=1,KNDAU(IWM)
        id=kdau(IWM,i)
c        print *,'part id ',kndau(iwm),KTPCOD(id)
        if (KTPCOD(id).ne.1) then
          if (KTPCOD(id).lt.100) then ! leptonic
            if (QCH(id).ne.0.) then
              idf=id
              icom=KTPCOD(id)
c              print *,'icom lep',icom,mult
            else
              idl=id
              icom2=KTPCOD(id)
            endif
          else                  ! hadronic
            del=abs(abs(QCH(id))-0.333333)
            if (del.lt.0.01) then ! charge=1/3=down type
              idf=id
              icom=KTPCOD(id)
c              print *,'icom had',icom,mult
            else
              idl=id
              icom2=KTPCOD(id)
            endif
          endif
        else ! initial state rad
          mult=-1
        endif
 10   continue


      icom=icom*mult
      icom2=icom2*mult
c     print *,'icom ',icom
      Eq(3)=QE(idf)
      Pxq(3)=Qx(idf)
      Pyq(3)=Qy(idf)
      Pzq(3)=Qz(idf)
      Eq(4)=QE(idl)
      Pxq(4)=Qx(idl)
      Pyq(4)=Qy(idl)
      Pzq(4)=Qz(idl)
      rpart(3)=icom
      rpart(4)=icom2


*** end getting four fermions

      
C transfer to Maciek definition.
!
!   COMMON / MOMDEC / Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
!      W- ;W+ ;f (W-);\bar f (W-);f (W+);\bar f (W+);
! qqqq          d        ua         u        da
! en            e-       na         n        e+        
C
C i used
C qqqq          d        ua        da        u  
C en            e-       na        e+        n  
C and W+ first 
      do i=1,4
        bpp4(1)=Pxq(1)
        bpp4(2)=Pyq(1)
        bpp4(3)=Pzq(1)
        bpp4(4)=Eq(1)

        bpp3(1)=Pxq(2)
        bpp3(2)=Pyq(2)
        bpp3(3)=Pzq(2)
        bpp3(4)=Eq(2)
        

        bpp1(1)=Pxq(3)
        bpp1(2)=Pyq(3)
        bpp1(3)=Pzq(3)
        bpp1(4)=Eq(3)

        bpp2(1)=Pxq(4)
        bpp2(2)=Pyq(4)
        bpp2(3)=Pzq(4)
        bpp2(4)=Eq(4)
      enddo
      
      ikpar=nlink('KPAR',0)
      if (ikpar.ne.0) then
        Wmassnom=rw(ikpar+28)
        Wwidthnom=rw(ikpar+29)
      endif
      
c      print *,' Pxq Pyq Pzq Eq '
c      print *,Pxq(1),Pyq(1),Pzq(1),Eq(1),Eq(2),Eq(3),Eq(4)
      

      xsec=0.
      xnom=0.
      WmassWeiKWCC03=1.
      if (Eq(1).gt.0 .and. Eq(2).gt.0 .and.
     +    Eq(3).gt.0 .and. Eq(4).gt.0) then
        call marinitkw(Wmassnom,Wwidthnom)
        xnom=born4m(bpp1,bpp2,bpp3,bpp4)
        if (Wwidth.lt.0.) then
          Wwidth=(Wmass/Wmassnom)**3 * Wwidthnom
        endif
        call marinitkw(Wmass,Wwidth)
        xsec=born4m(bpp1,bpp2,bpp3,bpp4)
      endif
      
C      print *,'xsec ',xsec,Wmass,Wmassnom

      if (xnom.gt.0 .and. xsec.gt.0) then
        WmassWeiKWCC03=xsec/xnom
      else
        print *,'WmassWeiKWCC03: ---> Xsec zero ',xsec,xnom
      endif


      return
 999  continue
      print *,'WmassWeiKWCC03: ---> NO W FOUND'
      end


***********************************************************************
* the function returns


      subroutine marinitkw(Wmass,Wwidth)
*     ******************

      implicit real*8 (a-h,o-z)
#include "qdecl.h"
#include "qcde.h"
#include "qhac.h"

      common / matpar / dpi,ceuler     

!.. general purpose common of parameters
      common / weking / ene,amaz,gammz,amel,amfin,xk0,sinw2,ide,idf 
      common / wekin2 / amaw,gammw,gmu,alphaw   

!.. necessary for kineXX to construct 4momenta, should be set to 0
      common / momset / qeff1(4),qeff2(4),sphum(4),sphot(100,4),nphot  

! This common can be everywhere, contains various switches
      COMMON / KeyKey/  KeyRad,KeyPhy,KeyTek,KeyMis,KeyDwm,KeyDwp

!   --Anomalous WWV Couplings, for WWamgc only
      COMMON / ancoco / g1(2),kap(2),lam(2),g4(2),g5(2),kapt(2),lamt(2)
      COMPLEX*16 g1,kap,lam,g4,g5,kapt,lamt

      real Wmass,Wwidth

#include "qmacro.h"




!.. for clarity set to zero ALL params in common blocks, 
!.. then reset the important ones.
!.. parameters not reset are DUMMY for born calculation

      DPI=3.1415926535897932D0

      do i=1,4
        qeff1(i)=0d0
        qeff2(i)=0d0
        sphum(i)=0d0
        do j=1,100
          sphot(j,i)=0d0
        enddo
      enddo
!.. nphot MUST be always kept 0 if one uses kineXX routines
      nphot=0

      ene=0d0
      amaz=0d0
      gammz=0d0
      amel=0d0
      amfin=0d0
      xk0=0d0
      sinw2=0d0
      ide=0
      idf=0

      amaw=0d0
      gammw=0d0
      gmu=0d0
      alphaw=0d0

      KeyRad=0
      KeyPhy=0
      KeyTek=0
      KeyMis=0
      KeyDwm=0
      KeyDwp=0
******************************************************************
*
* NOMINAL W mass
      amaw  = 1d0*(Wmass)
      gammw = 1d0*(Wwidth)
******************************************************************

!.. alive parameters in common blocks. all others are dummy for Born.
      ikpar=nlink('KPAR',0)

      KeyMix = MOD(int(rw(ikpar+2+4)),10)
      Key4f = MOD(int(rw(ikpar+2+4)),100)/10
      KeyISR = MOD(int(rw(ikpar+2+1)),10)

c      print *,'keyrad ',keyisr

      amaz  = rw(ikpar+26)
C      gammz =  2.4974d0
      gammz = rw(ikpar+27)

C      gammw =  2.086092d0

C      alphaw=  1d0/128.07D0

      alphaw=  1d0/rw(ikpar+25)
c      print *,'alphaw '
      alphaw=1./
     +     (1./alphaw+2./3./DPI*dlog(amaw/dble(rw(ikpar+28))))     
c      print *,'alphaw myself ',alphaw
      ALFWIN=1D0/alphaw

      GMU=rw(ikpar+24)

c.. sinw2 is calculated from gmu,alfwin,amaz
      A2 = DPI / ( ALFWIN*SQRT(2D0)*GMU )
      SINW2 = ( 1-SQRT( 1-(4*A2/AMAZ**2) ) )/2D0  

c.. alpha_w
      ALPHAW = 1D0/ ALFWIN
!-----------------------------------------------------------------------
!-- sinw2 electroweak mixing angle
      IF(KeyMix .EQ. 1) THEN
         A2 = DPI / ( ALFWIN*SQRT(2D0)*GMU )
         SINW2 = ( 1-SQRT( 1-(4*A2/AMAZ**2) ) )/2D0 
      ELSEIF(KeyMix .EQ. 2) THEN
         SINW2 = 1D0 -AMAW**2/AMAZ**2
      ELSE 
!-- new definition for LEP200 workshop
         sinw2 = Dpi * alphaw /( sqrt(2d0) * amaw**2 * gmu )
      ENDIF
C
C      sinw2 =  0.231030912451578502d0 ! 1D0 -AMAW**2/AMAZ**2

c      print *,'Keymix Sinw2 ',Keymix, sinw2
c      print *,'mass wi Z ',amaz,gammz
c      print *,'mass wi W ',amaw,gammw
c      print *,'alpha W   ',alphaW

C 80.25     sinw2=0.233
C  Pxq Pyq Pzq Eq 
C  -2.748567      -27.51628       10.17260       29.46561       53.24711    
C   59.30864       29.97610    
C keyrad            0
C Keymix Sinw2            1  0.231030912451579     
C mass wi Z    91.1887969970703        2.49740028381348     
C mass wi W    80.2500000000000        2.08609199523926     
C alpha W     7.808229427741232E-003
C  xsec   8.6385608E-06
C
C 70.75
c   8.259672       12.68537      -28.42242       32.20210       53.11419    
c   45.17196       40.67379    
c keyrad            0
c Keymix Sinw2            1  0.231030912451579     
c mass wi Z    91.1887969970703        2.49740028381348     
c mass wi W    79.7500000000000        2.04734230041504     
c alpha W     7.813720223127877E-003
c keyrad            0
c Keymix Sinw2            1  0.231030912451579     
c mass wi Z    91.1887969970703        2.49740028381348     
c mass wi W    80.2500000000000        2.08609223365784     
c alpha W     7.813720223127877E-003
c xsec   6.2464317E-04
c--------------- fout adressering keyrad, keymix
c keyrad            1
c Keymix Sinw2            0  0.230915770158055     
c mass wi Z    91.1887969970703        2.49740028381348     
c mass wi W    80.2500000000000        2.08609199523926     
c alpha W     7.808229427741232E-003
c xsec   8.6385608E-06


 
!... various switches, only relevant are 
!    keyred, keyacc, keyspn, keywu, keyzet
!    keyred must be set to 5 in order to desactivate overwriting of
!                                      externally supplied 4momenta
!      KeyRed = MOD(KeyPhy,100000)/10000
!      KeyAcc = MOD(KeyMis,1000)/100
!      KeySpn = MOD(KeyPhy,10000)/1000
!      KeyWu  = MOD(KeyPhy,1000000)/100000
!      KeyZet = MOD(KeyPhy,1000)/100
      keyphy = 51011  ! KeyWu, KeyRed, KeySpn,KeyZet, KeyMas, KeyBra
      keymis = 100   ! 100*KeyACC +10*Key4f +KeyMix 

C.. math. constants


! DONE in xsec
!----------------------------------------------------------------------!
! Setting up the anomalous couplings as given in the paper:            !
!     K. Hagiwara, R.D. Peccei, D. Zeppenfeld and K. Hikasa,           !
!                 Nucl. Phys. B282 (1987) 253.                         !
! The variables used in this routine correspond to the following       !
! contants defined in the above paper:                                 !
!           constant name     corresponding variable                   ! 
!                g_1^V                g1(2)                            !
!                kappa_V              kap(2)                           !
!                lambda_V             lam(2)                           !
!                g_4^V                g4(2)                            !
!                g_5^V                g5(2)                            !
!                kappa-tilde_V        kapt(2)                          !
!                lambda-tilde_V       lamt(2)                          ! 
!----------------------------------------------------------------------!      

! Set up the constants (within SM: g1=kap=1, lam=g4=g5=kapt=lamt=0):
! 1) for WWgamma vertex
      g1(1)   = (1.0d0, 0.0d0)
      kap(1)  = (1.0d0, 0.0d0)
      lam(1)  = (0.0d0, 0.0d0)
      g4(1)   = (0.0d0, 0.0d0)
      g5(1)   = (0.0d0, 0.0d0)
      kapt(1) = (0.0d0, 0.0d0)
      lamt(1) = (0.0d0, 0.0d0)
! 2 for WWZ vertex
      g1(2)   = (1.0d0, 0.0d0)
      kap(2)  = (1.0d0, 0.0d0)
      lam(2)  = (0.0d0, 0.0d0)
      g4(2)   = (0.0d0, 0.0d0)
      g5(2)   = (0.0d0, 0.0d0)
      kapt(2) = (0.0d0, 0.0d0)
      lamt(2) = (0.0d0, 0.0d0)
C now transform the 5 relevant couplings:
C      DZ, XG, XZ, YG, YZ

      sinw2 =  0.231030912451578502d0 ! 1D0 -AMAW**2/AMAZ**2
      ctgw=1.d0/tan(asin(sqrt(sinw2)))
      
      DZ=0.d0
      XG=0.d0
      XZ=0.d0
      YG=0.d0
      YZ=0.d0
C--
C g1
C--
      gzww=DZ+ctgw
      c=gzww/ctgw
      g1(2)=c
c      print *,g1(2)
C--
C KG
C--
      c=XG+1.d0
      kap(1)=c
C--
C LG
C--
      c=YG
      lam(1)=c
C--
C KZ
C--
      c=XZ/gzww+1.d0
      kap(2)=c
C--
C LZ
C--
      c=YZ/ctgw
      lam(2)=c

      end

      double precision function born4m(bpp1,bpp2,bpp3,bpp4)
*     ************************************************
      implicit real*8 (a-h,o-z)
!.. wwborn reads its 4momenta from this common
      common / bormom / bq1(4),bq2(4),bp1(4),bp2(4),bp3(4),bp4(4)
      dimension bb(4),bpp1(4),bpp2(4),bpp3(4),bpp4(4)
      
      data init / 0 /

cmv      if(init.eq.0) then
cmv        init=1
cmv        call marinit
cmv        write(6,*)'initialization'
cmv      endif

      do i=1,4
        bq1(i)=bpp1(i)+bpp2(i)
        bq2(i)=bpp3(i)+bpp4(i)
        bp1(i)=bpp1(i)
        bp2(i)=bpp2(i)
        bp3(i)=bpp3(i)
        bp4(i)=bpp4(i)
        bb(i)=bq1(i)+bq2(i)
      enddo

      sprim=prodm(bb,bb)
      s1=prodm(bq1,bq1)
      s2=prodm(bq2,bq2)

! distinction svar-sprim relevant only in the presence of photons
      svar=sprim

!      write(6,*) 'ma: svar,sprim,s1,s2=',svar,sprim,s1,s2
      born4m = wwborn(svar,sprim,s1,s2)

      end

      function 
     @  born6a(svar,sprim,s1,s2,costhe,phi,cosde1,phi1,cosde2,phi2)
*     *************************************************************
! distinction svar-sprim relevant only in the presence of photons
      implicit real*8 (a-h,o-z)
      dimension amdec0(4),dq1(4),dq2(4),dp1(4),dp2(4),dp3(4),dp4(4)
      
!.. wwborn reads its 4momenta from this common
      common / bormom / bq1(4),bq2(4),bp1(4),bp2(4),bp3(4),bp4(4)

      data init / 0 /

cmv      if(init.eq.0) then
cmv        init=1
cmv        call marinit
cmv        write(6,*)'initialization'
cmv      endif

!.. final state fermions must be massless
        do 11,i=1,4
 11     amdec0(i)=0d0

!.. construct 4momenta out of angles.
!.. two constructions are given here, based on different 
!.. definition of angles.
!.. note, that kineww/kinett send output 4momenta in common / bormom /
!.. in the absence of photons they are identical to arguments dqI, dpI.
      itype = 1
cmv      if(itype.eq.1) then
cmv        call kineww(svar,sprim,costhe,phi,cosde1,phi1,cosde2,phi2,
cmv     $          sqrt(s1),sqrt(s2),amdec0,   dq1,dq2,dp1,dp2,dp3,dp4)
cmv      elseif(itype.eq.2) then
cmv        call kinett(svar,sprim,costhe,phi,cosde1,phi1,cosde2,phi2,
cmv     $          sqrt(s1),sqrt(s2),amdec0,   dq1,dq2,dp1,dp2,dp3,dp4)
cmv      endif

!.. to test rotational invariance around the beams
!..       do ij=1,10
!..       ph1=0.5d0
!..       call RXTOD3(PH1,bq1,bq1)
!..       call RXTOD3(PH1,bq2,bq2)
!..       call RXTOD3(PH1,bp1,bp1)
!..       call RXTOD3(PH1,bp2,bp2)
!..       call RXTOD3(PH1,bp3,bp3)
!..       call RXTOD3(PH1,bp4,bp4)
!..
!..       born6a = wwborn(svar,sprim,s1,s2)
!.. 
!..c       call dumpb(6)
!..c       write(6,*)'masses',dmas2(bp1),dmas2(bp2),dmas2(bp3),dmas2(bp4)
!..       write(6,*)'born=',born6a
!..       enddo

       born6a = wwborn(svar,sprim,s1,s2)

      end

c......==================== BORN ==========================
      function wwborn(svar,sprim,s1,s2)
*     **************************************
c This function provides a value of a differential born cross section
c for W+W- pair production and decay in e+e- scattering. 
c INPUT: svar   - CMS energy squared (in GeV**2)
c        costhe - cosine of the W- polar angle in the CMS
c                 of the incoming e+e- with z-axis pointing 
c                 in the e- direction
c Written by: Wieslaw Placzek            date: 20.07.1994
c Last update: 27.07.1994                by: W.P.
c
      implicit real*8 (a-h,o-z)
      common / matpar / pi,ceuler     
      common / weking / ene,amaz,gammz,amel,amfin,xk0,sinw2,ide,idf 
      common / wekin2 / amaw,gammw,gmu,alphaw   
      common / bormom / bq1(4),bq2(4),bp1(4),bp2(4),bp3(4),bp4(4)
      common / angles / costhe,phi,cosde1,phi1,cosde2,phi2
      common / decays / iflav(4), amdec(4), br(2), brel
C      save   / momset /,/ momdec /,/ decays /,/ weking /,/ wekin2 /
C on retire momdec ( S.J. 3/2/97)
      save   / momset /,/ decays /,/ weking /,/ wekin2 /
! This common can be everywhere, contains various switches
      COMMON / KeyKey/  KeyRad,KeyPhy,KeyTek,KeyMis,KeyDwm,KeyDwp
! polarization amplitudes for WW production for left and right 
! handed electrons
      complex*16 awwel(3,3),awwer(3,3)
! polarization amplitudes for W-decays
      complex*16 adwm(3),adwp(3)
! auxilliary
c      complex*16 amww(3,3)
      complex*16 amwwr, amwwl,BWEXA1,BWEXA2
      dimension amdec0(4),dq1(4),dq2(4),dp1(4),dp2(4),dp3(4),dp4(4)
      save

      KeySpn = MOD(KeyPhy,10000)/1000
      KeyRed = MOD(KeyPhy,100000)/10000
      KeyAcc = MOD(KeyMis,1000)/100
      KeyWu  = MOD(KeyPhy,1000000)/100000
!
c.. first make decay products massless
cc        write(6,*)'WWborn dumps before'
cc        call dumpb(6)
!=======================================================================
! for born fed with external angles from common / bormom / (25.11.96 MS)
      keyred=5
! for born fed with external angles from common / bormom / (25.11.96 MS)
!=======================================================================
      if(KeyRed.eq.0) then
        do 11,i=1,4
 11     amdec0(i)=0d0
        stop
cmv        call kineww(svar,sprim,costhe,phi,cosde1,phi1,cosde2,phi2,
cmv     $          sqrt(s1),sqrt(s2),amdec0,   dq1,dq2,dp1,dp2,dp3,dp4)
       elseif(KeyRed.eq.1) then
c.. brute force, no 4 mom conserv.
        bp1mod=bp1(4)/dsqrt(bp1(1)**2+bp1(2)**2+bp1(3)**2)
        bp2mod=bp2(4)/dsqrt(bp2(1)**2+bp2(2)**2+bp2(3)**2)
        bp3mod=bp3(4)/dsqrt(bp3(1)**2+bp3(2)**2+bp3(3)**2)
        bp4mod=bp4(4)/dsqrt(bp4(1)**2+bp4(2)**2+bp4(3)**2)
        do 15,i=1,3
        bp1(i)=bp1(i)*bp1mod
        bp2(i)=bp2(i)*bp2mod
        bp3(i)=bp3(i)*bp3mod
 15     bp4(i)=bp4(i)*bp4mod
!=======================================================================
! for born fed with external angles from common / bormom / (25.11.96 MS)
      elseif(KeyRed.eq.5) then
        continue
! for born fed with external angles from common / bormom / (25.11.96 MS)
!=======================================================================
      else
        write(6,*)'wwborn==>wrong KEYRED:',keyred
        stop
      endif
c        write(6,*)'WWborn dumps after'
c        call dumpb(6)
c      write(6,*)'masses',dmas2(bp1),dmas2(bp2),dmas2(bp3),dmas2(bp4)
      wwborn=0

!! here was a bug A. Vallasi 30.06.96
!ms      wlambd=abs(sprim**2+s1**2+s2**2-2*sprim*s1-2*sprim*s2-2*s1*s2)
       wlambd=max(0d0,abs((sprim-s1-s2)**2 -4*s1*s2))
!      tvar=-(sprim-s1-s2-dsqrt(wlambd)*costhe)/2
cmv tvar is recalculated at a later stage.
cmv      tvar=-1d0/2d0*(dsqrt(wlambd)*(1d0-costhe)
cmv     $              +4d0*s1*s2/(sprim-s1-s2+dsqrt(wlambd)) ) 
!====================================================================
!   tvar cannot be calculated from angles if external 4moms are given
!   (25.11.96 MS)
      dq1(3)=dsqrt(sprim)/2d0
      dq1(2)=0d0
      dq1(1)=0d0
      dq1(4)=dsqrt(sprim)/2d0
      tvar= s1 -2*prodm(dq1,bq1)
!   (25.11.96 MS)
!   tvar cannot be calculated from angles if external 4moms are given
!====================================================================
!... Amplitudes for WW production
      IF(keyacc.EQ.0) THEN 
        call wwprod(sprim,tvar,bq1,bq2,awwel,awwer) 
      ELSE
        call wwamgc(sprim,tvar,bq1,bq2,awwel,awwer) 
      ENDIF
!... Amplitudes for W-decays
      call wdecay(bq1,bp1,bp2,adwm)  
      call wdecay(bq2,bp3,bp4,adwp) 
!... Exact W-propagators:
      IF(KeyWu.EQ.0) THEN
        BWEXA1=1D0/dcmplx((S1-AMAW**2),(S1/AMAW*GAMMW))
        BWEXA2=1D0/dcmplx((S2-AMAW**2),(S2/AMAW*GAMMW))
      ELSEIF(KeyWu.EQ.1) THEN
        BWEXA1=1D0/dcmplx((S1-AMAW**2),(AMAW*GAMMW))
        BWEXA2=1D0/dcmplx((S2-AMAW**2),(AMAW*GAMMW))
      ELSEIF(KeyWu.EQ.2) THEN
        BWEXA1=1D0/dcmplx((S1-AMAW**2),0D0)
        BWEXA2=1D0/dcmplx((S2-AMAW**2),0D0)
      ELSE
        WRITE(6,*)'BWIGN==> Wrong KeyWu=',keywu
        STOP
      ENDIF


!... Polarization amplitudes for WW production and decay
      if(keyspn.eq.0)then
        xmatr=0
        do 10 l2=1,3
        do 10 l1=1,3

          amwwl=awwel(l1,l2)*adwm(l1)*adwp(l2)*BWEXA1*BWEXA2 !ms,zw
          amwwr=awwer(l1,l2)*adwm(l1)*adwp(l2)*BWEXA1*BWEXA2 !ms,zw
          xmatr=xmatr + amwwl*dconjg(amwwl) !ms
          xmatr=xmatr + amwwr*dconjg(amwwr) !ms
 10     continue
      elseif(keyspn.eq.1)then
        xmatr=0
        amwwl=(0,0)
        amwwr=(0,0)
        do 20 l2=1,3
        do 20 l1=1,3
          amwwl=amwwl+awwel(l1,l2)*adwm(l1)*adwp(l2)*BWEXA1*BWEXA2 !ms,zw
          amwwr=amwwr+awwer(l1,l2)*adwm(l1)*adwp(l2)*BWEXA1*BWEXA2 !ms,zw
 20     continue
        xmatr=xmatr + amwwl*dconjg(amwwl) !ms
        xmatr=xmatr + amwwr*dconjg(amwwr) !ms
      endif
      fkin=1D0
      wwborn=fkin*xmatr 
      end

      subroutine wwprod(s,t,q1,q2,awwel,awwer)
*     ****************************************
c This routine calculates polarization amplitudes for the process
c e+e- --> W+W-, for on-shell W's in Born approximation. Calculation
c is done in the CMS of e+e- with z-axis pointing along e- direction. 
c It is based on the paper: 
c K. Kolodziej, M. Zralek, Phys. Rev. D43 (1991) 43;
c INPUT: s   - center mass energy squared (in GeV**2)
c        t   - transfer (in GeV**2)
c        q1(4) - four-momentum of W-  
c        q2(4) - four-momentum of W+
c OUTPUT: awwel(3,3) - complex array containing polarization amplitudes 
c                     for left-handed electron 
c                     {M_0(-,+,l1,l2) in eq. (31)}
c         awwel(3,3) - complex array containing polarization amplitudes 
c                     for right-handed electron 
c                     {M_0(+,-,l1,l2) in eq. (31)}
c
c Written by: Wieslaw Placzek            date: 01.07.1994
c Last update: 02.08.1994                by: W.P.
c
      implicit real*8 (a-h,o-z)
      real*8 q1(4),q2(4)
      complex*16 awwel(3,3),awwer(3,3)
      common / matpar / pi,ceuler     
      common / weking / ene,amaz,gammz,amel,amfin,xk0,sinw2,ide,idf 
      common / wekin2 / amaw,gammw,gmu,alphaw   
! This common can be everywhere, contains various switches
      COMMON / KeyKey/  KeyRad,KeyPhy,KeyTek,KeyMis,KeyDwm,KeyDwp
! polarization vectors of W-bosons, eq. (9)
      real*8 eps1(4,3),eps2(4,3)
      real*8 e1(4),e2(4),p1(4),pmq(4)
      complex*16 f1e1,f1e2,f1q2,fkz3
      data init /0/
      save init,zetl,etal,etar
      save
!
!...initialization
!--------------------------------------------------------------------- 
      if (init.eq.0) then
        init=1 
        KeyZet = MOD(KeyPhy,1000)/100
!... electroweak coefficient factors, eq. (30)
        zetl=0.5d0/sinw2
        etal=1-zetl
        etar=1
      endif  
!... calculation
!=====================================================================
!... four-momentum of the incoming electron in CMS 
!    (+z axis along the electron beam)
      ecm=dsqrt(s)   
      p1(1)=0d0
      p1(2)=0d0
      p1(3)=ecm/2d0
      p1(4)=p1(3)
!... calculation of polarization vectors of W-bosons, eq. (9)
      call polvec(q1,eps1)
      call polvec(q2,eps2)
!... calculation of the polarization amplitudes, eq. (31)
!    {note: we use different normalization!}
!    (note: electric charge squared equals: 4*pi*alphaw}
      wsp=-4*pi*alphaw*ecm
      do 10 k=1,4
 10   pmq(k)=p1(k)-q1(k)
      do 20 l2=1,3
      do 20 l1=1,3
        do 25 k=1,4
          e1(k)=eps1(k,l1)
 25       e2(k)=eps2(k,l2)
        f1e1=dcmplx(e1(1),-e1(2))
        f1e2=dcmplx(e2(1),-e2(2))
        f1q2=dcmplx(q2(1),-q2(2))
        e1e2=prodm(e1,e2)
        e1q2=prodm(e1,q2)
        e2q1=prodm(e2,q1) 
        if(keyzet.eq.0)then
        awwel(l1,l2)=( 2*(1/s - etal/dcmplx(s-amaz**2,s/amaz*gammz))*
     &                   (e1q2*f1e2 - e1e2*f1q2 - e2q1*f1e1) + 
     &                 zetl/t*fkz3(e2,pmq,e1) ) *wsp   
        awwer(l1,l2)= -2*(1/s - etar/dcmplx(s-amaz**2,s/amaz*gammz))*
     &                   (e1q2*conjg(f1e2) - e1e2*conjg(f1q2) -
     &                    e2q1*conjg(f1e1)) *wsp
        elseif(keyzet.eq.1)then
        awwel(l1,l2)=( 2*(1/s - etal/dcmplx(s-amaz**2,amaz*gammz))*
     &                   (e1q2*f1e2 - e1e2*f1q2 - e2q1*f1e1) + 
     &                 zetl/t*fkz3(e2,pmq,e1) ) *wsp   
        awwer(l1,l2)= -2*(1/s - etar/dcmplx(s-amaz**2,amaz*gammz))*
     &                   (e1q2*conjg(f1e2) - e1e2*conjg(f1q2) -
     &                    e2q1*conjg(f1e1)) *wsp
        elseif(keyzet.eq.2)then
        awwel(l1,l2)=( 2*(1/s - etal/(s-amaz**2))*
     &                   (e1q2*f1e2 - e1e2*f1q2 - e2q1*f1e1) + 
     &                 zetl/t*fkz3(e2,pmq,e1) ) *wsp   
        awwer(l1,l2)= -2*(1/s - etar/(s-amaz**2))*
     &                   (e1q2*conjg(f1e2) - e1e2*conjg(f1q2) -
     &                    e2q1*conjg(f1e1)) *wsp
        else
          write(6,*)'wrong KEYZET:',keyzet
        endif
 20   continue
      end  

      function prodm(p,q)
*     *******************
c Scalar product of the four-vectors p and q in Minkowski space;
c note: p_0=p(4), q_0=q(4)
c 
      implicit real*8 (a-h,o-z)
      real*8 p(4),q(4)
      prodm=p(4)*q(4)-p(3)*q(3)-p(2)*q(2)-p(1)*q(1)
      end

      subroutine polvec(q,eps)
*     **************************
c Calculation of polarization vectors of a vector boson 
c in the rectangular basis, see eq. (9); see also K. Hagiwara 
c and D. Zeppenfeld, Nucl. Phys. B274 (1986) 1, eq. (3.47).
c     INPUT:  q(4)     - four-momentum of the vector boson
c     OUTPUT: eps(4,3) - three polarization four-vector
c
c Written by: Wieslaw Placzek            date: 01.07.1994
c Last update: 27.07.1994                by: W.P.
c
      implicit real*8 (a-h,o-z)
      real*8 q(4),eps(4,3)   
      save
!
      qt2=q(1)**2+q(2)**2
      qt =sqrt(qt2)
      aq2=qt2+q(3)**2
      aq =sqrt(aq2)
      am =sqrt(q(4)**2-aq2) 
      do 10 l=1,3
      do 10 k=1,4
   10 eps(k,l)=0
      if(aq.lt.1d-10) then
        eps(1,1)=1
        eps(2,2)=1
        if (am.gt.1d-10) eps(3,3)=1 
      elseif (qt2.lt.1d-10) then
        eps(1,1)=q(3)/aq
        eps(2,2)=q(3)/aq
        if (am.gt.1d-10) then
          eps(3,3)=q(4)/am/aq*q(3)
          eps(4,3)=aq/am
        endif
      else
        ws1=1/aq/qt
        eps(1,1)= ws1*q(1)*q(3)
        eps(2,1)= ws1*q(2)*q(3)
        eps(3,1)=-ws1*qt2
        eps(4,1)= 0
        ws2=1/qt 
        eps(1,2)=-ws2*q(2)
        eps(2,2)= ws2*q(1)
        eps(3,2)= 0
        eps(4,2)= 0
        if (am.gt.1d-10) then
          ws3=q(4)/am/aq 
          eps(1,3)=ws3*q(1)
          eps(2,3)=ws3*q(2)
          eps(3,3)=ws3*q(3)
          eps(4,3)=ws3*aq2/q(4)
        endif
      endif  
      end

      function fkz3(a,b,c)
*     *******************************
c Function F_3 of four-vectors contracted with Dirac matrices; 
c see eq. (19)
c
c Written by: Wieslaw Placzek            date: 01.07.1994
c
      real*8 a(4),b(4),c(4)
      complex*16 fkz3
      fkz3=(a(4) +a(3)) * ((b(4) -b(3))*dcmplx(c(1),-c(2)) -
     &                     dcmplx(b(1),-b(2))*(c(4) -c(3)))+ 
     &     dcmplx(a(1),-a(2)) * ((b(4) +b(3))*(c(4) -c(3)) -
     &                     dcmplx(b(1), b(2))*dcmplx(c(1),-c(2)))
      end        

      subroutine wdecay(q,p1,p2,adw)
*     ******************************
c This routine calculates polarization amplitudes for W decays 
c into massless fermions. It is based on the paper: 
c K. Hagiwara et al., Nucl. Phys. B282 (1987) 253; see Appendix C.
c No CKM-mixing matrix elements incuded here.
c INPUT: q(4)        - four-momentum of W  
c        p1(4),p2(4) - four-momenta of decay products
c OUTPUT: adw(3) - complex array containing W decay amplitudes 
c                   {M(lambda,sigma_1,sigma_2) in eq. (C.16)}
c
c Written by: Wieslaw Placzek            date: 20.07.1994
c Last update: 02.08.1994                by: W.P.
c
      implicit real*8 (a-h,o-z)
      common / matpar / pi,ceuler
      common / weking / ene,amaz,gammz,amel,amfin,xk0,sinw2,ide,idf 
      common / wekin2 / amaw,gammw,gmu,alphaw   
      save / matpar /, / weking /, / wekin2 /
      real*8 q(4),p1(4),p2(4)
      complex*16 adw(3)
      real*8 eps(4,3),e(4)
      complex*16 sfunhz
!
!... calculation of polarization vectors of W in rectangular basis
      call polvec(q,eps)
!... calculation of the W decay amplitudes
      do 10 l=1,3
        do 15 k=1,4
 15       e(k)=eps(k,l)
        adw(l)=  sqrt(4*pi*alphaw) /sqrt(2*sinw2)
     $        *2*sqrt(p1(4)*p2(4)) *sfunhz(p1,e,p2)
 10   continue
      end

      function sfunhz(p1,a,p2)
*     ********************************** *
c Spinorial string S(pi,a,pf) for massless spinors chi(pi), chi(pf);
c a(4) - given four-vector, 
c see K. Hagiwara et al., Nucl. Phys. B282 (1987) 253; Appendix C. 
c
c Written by: Wieslaw Placzek            date: 20.07.1994
c Last update: 01.08.1994                by: W.P.
c
      implicit real*8 (a-h,o-z)
      real*8 a(4),p1(4),p2(4)
      complex*16 z1p,z2m,zam,zap
      complex*16 sfunhz
!
      x1p=p1(3)+p1(4)
      x2p=p2(3)+p2(4)
      z1p=dcmplx(p1(1), p1(2))
      z2m=dcmplx(p2(1),-p2(2))
      xam=a(4)-a(3)
      xap=a(4)+a(3)
      zam=dcmplx(a(1),-a(2))
      zap=dcmplx(a(1), a(2))
      if (x1p.gt.1d-20 .and. x2p.gt.1d-20) then
        fac=0.5/sqrt(p1(4)*p2(4)*x1p*x2p)
        sfunhz=fac *( x1p*(x2p*xam-z2m*zap) + z1p*(z2m*xap-x2p*zam) )
      elseif (x1p.gt.1d-20) then
        sfunhz= (x1p*zap - z1p*xap)/sqrt(2*p1(4)*x1p)
      elseif (x2p.gt.1d-20) then
        sfunhz= (z2m*xap - x2p*zam)/sqrt(2*p2(4)*x2p)
      else
        sfunhz=xap
      endif
      end

!======================================================================
!============= Born Version with Anomalous Couplings ==================
!======================================================================

      subroutine WWamgc(s,t,q1,q2,awwel,awwer)
*     ****************************************
!----------------------------------------------------------------------!
! This routine calculates polarization amplitudes for the process:     !
!              e-(p1) e+(p2) ---> W-(q1) W+(q2)                        !
! Calculation  is done in the CMS of e+e- with z-axis pointing along   !
! the e- direction. These amplitudes include general type three boson  !
! coupling constant as given in the paper:                             !
!     K. Hagiwara, R.D. Peccei, D. Zeppenfeld and K. Hikasa,           !
!                 Nucl. Phys. B282 (1987) 253.                         !
! The same formalism as in the subroutine wwprod (where only the SM    !
! coupling constant are included) is used here.                        !
! INPUT: s   - center mass energy squared (in GeV**2)                  !
!        t   - transfer momentum squared (in GeV**2)                   !
!        q1(4) - four-momentum of W-                                   !
!        q2(4) - four-momentum of W+                                   !
! OUTPUT: awwel(3,3) - complex array containing polarization           !
!                      amplitudes for left-handed electron,            !
!                      M_0(-,+,l1,l2);                                 !
!         awwel(3,3) - complex array containing polarization           ! 
!                      amplitudes for right-handed electron,           ! 
!                      M_0(+,-,l1,l2).                                 !
! Note: Before first use of this routine general type coupling         !
!       constant have to be set up in the routine setacc.              !
!----------------------------------------------------------------------! 
! Written by: Wieslaw Placzek                 Knoxville, November 1995 !
! Last update: 16.11.1995            by: W.P.                          !
!----------------------------------------------------------------------!
      IMPLICIT REAL*8 (a-h,o-z)
      COMPLEX*16 zi
      PARAMETER ( zi = (0d0,1d0) )
      REAL*8 q1(4),q2(4)
      COMPLEX*16 awwel(3,3),awwer(3,3)
      COMMON / weking / ene,amaz,gammz,amel,amfin,xk0,sinw2,ide,idf 
      COMMON / wekin2 / amaw,gammw,gmu,alphaw 
      common / matpar / pi,ceuler  
! This common can be everywhere, contains various switches
      COMMON / KeyKey /  KeyRad,KeyPhy,KeyTek,KeyMis,KeyDwm,KeyDwp
      SAVE / weking /, / wekin2 /, / KeyKey /, / matpar /
! Polarization vectors of W-bosons          
      REAL*8 eps1(4,3),eps2(4,3)
      REAL*8 e1(4),e2(4),p1(4),pmq(4)
! General type 3-boson coupling constants
      COMPLEX*16 GamV(4,2)
      COMPLEX*16 dZ,Gagm,Gagp,GaZm,GaZp,fkz3
      DATA init /0/
      SAVE init,KeyZet,zetl,etal,etar
!--------------------------------------------------------------------- 
! Initialization ...
      IF (init.eq.0) THEN
        init=1 
        KeyZet = MOD(KeyPhy,1000)/100
! Electroweak coefficient factors
        zetl=0.5/sinw2
        etal=1-zetl
        etar=1
! Set up anomaluos couplings constants
!!        CALL setacc ! MOVED TO FILEXP
      ENDIF  
!--------------------------------------------------------------------- 
! Calculation ...
! Options for Z-boson width
! a) running Z-width
      IF (KeyZet.eq.0) THEN
        dZ = DCMPLX(s-amaz**2,s/amaz*gammz)
! b) constant Z-width
      ELSEIF (KeyZet.eq.1) THEN
        dZ = DCMPLX(s-amaz**2,amaz*gammz)
! c) zero Z-width
      ELSEIF (KeyZet.eq.2) THEN
        dZ = s-amaz**2
      ELSE
         WRITE(6,*)'>>> Wrong KeyZet:',KeyZet
      ENDIF
! Four-momentum of the incoming electron in CMS (+z axis along e-)
      p1(1) = 0
      p1(2) = 0
      p1(3) = SQRT(s)/2
      p1(4) = p1(3)
! Calculation of W polarization vectors
      CALL polvec(q1,eps1)
      CALL polvec(q2,eps2)
      DO k = 1,4
        pmq(k) = p1(k) - q1(k)
      ENDDO
      wsp =-4*pi*alphaw*SQRT(s)
! Calculation of the polarization amplitudes
      DO l2 = 1,3
        DO l1 = 1,3
          DO k = 1,4
            e1(k)=eps1(k,l1)
            e2(k)=eps2(k,l2)
          ENDDO
! Calculation of the 3-boson couplings
          CALL WWVgcc(s,amaw,q1,e1,q2,e2,GamV)
! WWgamma vertex
          Gagm = GamV(1,1) - zi*GamV(2,1)
          Gagp = GamV(1,1) + zi*GamV(2,1)          
! WWZ vertex
          GaZm = GamV(1,2) - zi*GamV(2,2)
          GaZp = GamV(1,2) + zi*GamV(2,2) 
! Polarization amplitudes         
          awwel(l1,l2) = wsp*( Gagm/s - GaZm*etal/dZ 
     &                       + zetl/t *fkz3(e2,pmq,e1) )
          awwer(l1,l2) = wsp*( Gagp/s - GaZp*etar/dZ )
        ENDDO
      ENDDO
      END

      subroutine setacc_UNUSED
*     *****************
!----------------------------------------------------------------------!
! Setting up the anomalous couplings constants as given in the paper:  !
!     K. Hagiwara, R.D. Peccei, D. Zeppenfeld and K. Hikasa,           !
!                 Nucl. Phys. B282 (1987) 253.                         !
! Note: This subroutine has to be called prior first use of the        !
! subroutine WWamgc.                                                   !
! The variables used in this routine correspond to the following       !
! contants defined in the above paper:                                 !
!           constant name     corresponding variable                   ! 
!                g_1^V                g1(2)                            !
!                kappa_V              kap(2)                           !
!                lambda_V             lam(2)                           !
!                g_4^V                g4(2)                            !
!                g_5^V                g5(2)                            !
!                kappa-tilde_V        kapt(2)                          !
!                lambda-tilde_V       lamt(2)                          ! 
!----------------------------------------------------------------------! 
! Written by: Wieslaw Placzek                 Knoxville, November 1995 !
! Last update: 14.11.1995            by: W.P.                          !
!----------------------------------------------------------------------!
      COMMON / ancoco / g1(2),kap(2),lam(2),g4(2),g5(2),kapt(2),lamt(2)
      COMPLEX*16 g1,kap,lam,g4,g5,kapt,lamt
      SAVE  / ancoco /
! Set up the constants (within SM: g1=kap=1, lam=g4=g5=kapt=lamt=0):
! 1) for WWgamma vertex
      g1(1)   = (1.0, 0.0)
      kap(1)  = (1.0, 0.0)
      lam(1)  = (0.0, 0.0)
      g4(1)   = (0.0, 0.0)
      g5(1)   = (0.0, 0.0)
      kapt(1) = (0.0, 0.0)
      lamt(1) = (0.0, 0.0)
! 2) for WWZ vertex
      g1(2)   = (1.0, 0.0)
      kap(2)  = (1.0, 0.0)
      lam(2)  = (0.0, 0.0)
      g4(2)   = (0.0, 0.0)
      g5(2)   = (0.0, 0.0)
      kapt(2) = (0.0, 0.0)
      lamt(2) = (0.0, 0.0)
      END

      subroutine WWVgcc(s,amW,q1,eps1,q2,eps2,GamV)
*     *********************************************
!----------------------------------------------------------------------!
! This routine calculates four-vector Gamma_V^mu including general     !
! type 3-boson WWV couplings (V=gamma,Z) as defined in the paper:      !
!     K. Hagiwara, R.D. Peccei, D. Zeppenfeld and K. Hikasa,           !
!                 Nucl. Phys. B282 (1987) 253.                         !
! INPUT: s       - center mass energy squared (in GeV**2)              !
!        amW     - W-boson mass                                        !
!        q1(4)   - four-momentum of W-                                 !
!        eps1(4) - polarizarion vector of W-                           !
!        q2(4)   - four-momentum of W+                                 !
!        eps1(4) - polarizarion vector of W+                           !
! OUTPUT: GamV(4,2) - 2 complex number four-vectors Gamma_V^mu:        !
!                     GamV(4,1) for WWgamma vertex,                    !
!                     GamV(4,2) for WWZ vertex.                        !
!----------------------------------------------------------------------! 
! Written by: Wieslaw Placzek                 Knoxville, November 1995 !
! Last update: 15.11.1995            by: W.P.                          !
!----------------------------------------------------------------------!
      IMPLICIT REAL*8 (a-h,o-z)
      COMPLEX*16 zi
      PARAMETER ( zi = (0d0,1d0) )
      REAL*8 q1(4),eps1(4),q2(4),eps2(4)
      COMPLEX*16 GamV(4,2)
      COMMON / ancoco / g1(2),kap(2),lam(2),g4(2),g5(2),kapt(2),lamt(2)
      COMPLEX*16 g1,kap,lam,g4,g5,kapt,lamt
      SAVE  / ancoco /
! 7 form factors f_i^V
      COMPLEX*16 f1(2),f2(2),f3(2),f4(2),f5(2),f6(2),f7(2)
      REAL*8 P(4),Q(4),Peps(4),Qeps(4)
!
      amW2 = amW**2
! Set up form factors
      DO i = 1,2
         f1(i) = g1(i) + s/(2*amW2) *lam(i)
         f2(i) = lam(i)
         f3(i) = g1(i) + kap(i) + lam(i)
         f4(i) = g4(i)
         f5(i) = g5(i)
         f6(i) = kapt(i) - lamt(i)
         f7(i) = -0.5*lamt(i)
      ENDDO
! Sum and difference of W's 4-momenta
      DO k = 1,4
         P(k) = q1(k) + q2(k)
         Q(k) = q1(k) - q2(k)
      ENDDO
! Scalar products of various 4-vectors
      e1e2 = prodm(eps1,eps2)
      Pe1  = prodm(P,eps1)
      Pe2  = prodm(P,eps2)
! Coefficients for f5 and f6 
      CALL epsabc(P,eps1,eps2,Peps)
      CALL epsabc(Q,eps1,eps2,Qeps)
! Coefficient for f7
      PQeps = prodm(P,Qeps)
! Calculate Gamma_V^mu
      DO i = 1,2
        DO k = 1,4
          GamV(k,i) = ( f1(i)*e1e2 - f2(i)/amW2*Pe1*Pe2 
     &                - f7(i)/amW2*PQeps )*Q(k)
     &              + (-f3(i) + zi*f4(i) )*Pe2*eps1(k)   
     &              + ( f3(i) + zi*f4(i) )*Pe1*eps2(k)   
     &              + zi*f5(i)*Qeps(k) - f6(i)*Peps(k)
        ENDDO
      ENDDO
      END

      subroutine epsabc(a,b,c,q)
*     **************************
!----------------------------------------------------------------------!
! This routine calculates four-vector q according to the formula:      !
!                                                                      !
!  q^mu = epsilon^{mu,alpha,beta,gamma} a_alpha b_beta c_gamma,        !
!                                                                      !
! where epsilon is a totally antisymmetric tensor in Bjorken & Drell   !
! convention, and a, b, c are four-vectors.                            !
!----------------------------------------------------------------------! 
! Written by: Wieslaw Placzek                 Knoxville, November 1995 !
! Last update: 16.11.1995            by: W.P.                          !
!----------------------------------------------------------------------!
      IMPLICIT REAL*8 (a-h,o-z)
      REAL*8 a(4),b(4),c(4),q(4)
!
      b1c2 = b(1)*c(2) - b(2)*c(1)
      b1c4 = b(1)*c(4) - b(4)*c(1)
      b2c3 = b(2)*c(3) - b(3)*c(2)
      b2c4 = b(2)*c(4) - b(4)*c(2)
      b3c1 = b(3)*c(1) - b(1)*c(3)
      b3c4 = b(3)*c(4) - b(4)*c(3)
! 4-vector q^mu
      q(1) = a(2)*b3c4 - a(3)*b2c4 + a(4)*b2c3
      q(2) =-a(1)*b3c4 + a(3)*b1c4 + a(4)*b3c1
      q(3) = a(1)*b2c4 - a(2)*b1c4 + a(4)*b1c2  
      q(4) = a(1)*b2c3 + a(2)*b3c1 + a(3)*b1c2
      END




