anomGamwt.f

      Subroutine dipolgammarij (iqed, E, theta, A, B, R)

c              For  gamma gamma -> tau- tau+  reaction
c
c  Calculates spin-correlation coefficients R(i,j) as functions of beam
c  energy E = sqrt(s)/2 (in GeV) and scattaring angle theta.
c  V is velocity of tau lepton, gam is Lorentz factor, alpha is
c  fine-structure constant.
c  Anomalous magnetic moment A1 = ASM + A and electric dipole moment B are
c  real constants.
c  Parameters A, B describe 'New Physics'.
c  ASM is anomalous magnetic moment of tau in Standard Model (SM):
c  S.Eidelman and M.Passera. Mod.Phys.Lett. A22, 159-179, 2007.
c
c  Order of coefficients:
c     i = 1,2,3 correspond to S_x, S_y, S_z for tau-,
c     j = 1,2,3 correspond to S'_x, S'_y, S'_z for tau+
c     i,j = 4 (equivalent to tt) corresponds to 1 (no spin dependence).

      Implicit none
      integer iqed
      integer i, j                         ! not used here
      real*8 e, theta, a, b, asm, a1, gam, e4, v
      real*8 r(1:4, 1:4)
      real*8 pi/3.141592653589793238d0/,alpha/7.2973525693d-3/
      real*8 mtau/1.77686d0/               ! mass of tau in GeV

      v = dsqrt(1.d0 -(mtau/e)**2)         ! tau velocity 
      e4 = (4.0d0 *pi *alpha)**2           ! (electric charge)^4
      gam = e/mtau                         ! Lorentz factor for tau

c  Cotribution to anomalous magnetic moment in Standard Model:
      asm = 1.17721d-3
      asm = asm *iqed                     ! switch for dipole SM part
      a1 = asm + a                        ! Standard Model + New Physics

      d2 = (v**2 *dcos(theta)**2 -1.d0)**2        ! factor in denominator

      r(1,1) = e4 *(-11.d0 *v**4 +28.d0 *a1 *v**2 +
     $ 4.d0 *(v**2 -2.d0) *dcos(2.d0*theta) *v**2 -
     $ (v**2 -4.d0 *a1 -2.d0) *dcos(4.d0*theta) *v**2 +
     $ 22.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)

      r(1,2) = e4 *v *b *(dcos(4.d0*theta) *v**2 +15.d0 *v**2 +
     $ 4.d0 *cos(2.d0*theta) -20.d0) /(4.d0 *d2)

      r(1,3) = e4 *v**2 *gam *((a1 -1.d0) *v**2 +
     $ (v**2 +a1 *(v**2 -2.d0) -1.d0) *dcos(2.d0*theta) +1.d0) *
     $ dsin(2.d0*theta) /(2.d0 *d2)

      r(1,4) = 0.d0

      r(2,1) = -r(1,2)

      r(2,2) = e4 *(-dcos(4.d0*theta) *v**4 -11.d0 *v**4 +
     $ 16.d0 *a1 *v**2 +4.d0 *(v**2 +4.d0 *a1) *dcos(2.d0*theta) *v**2 +
     $ 16.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)

      r(2,3) = e4 *v *gam *b *(dcos(2.d0*theta)*v**2 -3.d0*v**2 +2.d0) *
     $ dsin(2.d0*theta) /(2.d0 *d2)

      r(2,4) = 0.d0

      r(3,1) = r(1,3)

      r(3,2) = -r(2,3)

      r(3,3) = e4 *(-4.d0 *dcos(2.d0*theta) *v**4 +11.d0 *v**4 +
     $ 36.d0 *a1 *v**2 +(v**2 -4.d0 *a1 -2.d0) *dcos(4.d0*theta) *v**2 +
     $ 2.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)

      r(3,4) = 0.d0

      r(4,1) = 0.d0

      r(4,2) = 0.d0

      r(4,3) = 0.d0

      r(4,4) = e4 *(-dcos(4.d0*theta) *v**4 -11.d0 *v**4 -16.d0 *a1 *v**2 +
     $ 4.d0 *(v**2 -4.d0 *a1 -2.d0) *dcos(2.d0*theta) *v**2 +
     $ 8.d0 * v**2 +32.d0 *a1 +8.d0) /(8.d0 *d2)

      return
      end