C++InterfacetoTauola
initwksw.f
1 
2  SUBROUTINE initwkswdelt(mode,IDEX,IDFX,SVAR,SWSQEFF, DELTSQ, DeltV, GMU, ALPHAINV, AMZi, GAMMZi, KEYGSW,
3  &regsw1,cimgsw1,regsw2,cimgsw2,regsw3,cimgsw3,regsw4,cimgsw4,regsw6,cimgsw6 )
4 
5 
6 ! initialization routine coupling masses etc., explicitly varying SWSQ
7  IMPLICIT REAL*8 (a-h,o-z)
8  COMMON / t_beampm / ene ,amin,amfin,ide,idf
9  REAL*8 ene ,amin,amfin
10  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
11  & ,xupgi0 ,xupzi0 ,xupgf0 ,xupzf0
12  & ,ndiag0,ndiaga,keya,keyz
13  & ,itce,jtce,itcf,jtcf,kolor
14  REAL*8 ss,poln,t3e,qe,t3f,qf
15  & ,xupgi0(2),xupzi0(2),xupgf0(2),xupzf0(2)
16  COMMON / t_gauspm1/vvcor, zetvpi, gamvpi
17  & ,xupgi ,xupzi ,xupgf ,xupzf
18  COMPLEX*16 vvcor, zetvpi, gamvpi
19  COMPLEX*16 xupgi(2),xupzi(2),xupgf(2),xupzf(2)
20 
21  COMMON / t_gswprmn /swsq,amw,amz,amh,amtop,gammz
22  REAL*8 swsq,amw,amz,amh,amtop,gammz
23  COMMON / t_ewn / gmun, alphainvn
24  REAL*8 gmun, alphainvn
25  COMPLEX *16 gsw(10)
26  REAL*8 pi
27  DATA pi /3.141592653589793238462643d0/
28  gsw(1) = dcmplx(regsw1,cimgsw1)
29  gsw(2) = dcmplx(regsw2,cimgsw2)
30  gsw(3) = dcmplx(regsw3,cimgsw3)
31  gsw(4) = dcmplx(regsw4,cimgsw4)
32  ! GSW(5) out
33  gsw(6) = dcmplx(regsw6,cimgsw6)
34 
35 C PRINT *, ' initwksw GSW = ', SWSQEFF, ReGSW1, CImGSW1, ReGSW2, CImGSW2, ReGSW6, CImGSW6
36 
37 C SWSQ = sin2 (theta Weinberg)
38 C AMW,AMZ = W & Z boson masses respectively
39 C AMH = the Higgs mass
40 C AMTOP = the top mass
41 C GAMMZ = Z0 width
42 C
43  ene=sqrt(svar)/2
44  amin=0.511d-3
45  swsq=swsqeff
46  amz=amzi !91.1887
47  gammz=gammzi !2.4952
48  gmun=gmu
49  alphainvn=alphainv
50 
51 
52 C Gfermi=1.16639d-5
53  gfermi=gmu
54 
55  zetvpi = gfermi *amz**2 *alphainv /(dsqrt(2.d0)*8.d0*pi)
56  $ *(swsq*(1d0-swsq)) *16d0
57  $ * gsw(1)
58 C updated following KK2f_defaults
59 C IF( KEYGSW.NE.0) THEN
60 C GAMMZ=2.50072032
61 C ENDIF
62 
63 
64 
65 
66  gamvpi = 1d0 /(2d0-gsw(6))
67 
68 C PRINT *, ' initwksw ZetVPi, GamVPi = ', GSW(1), ZetVPi, GamVPi
69 
70 
71  IF (idfx.EQ. 11) then
72  idf=2 ! denotes tau +2 tau-
73  amfin=0.511d-3 !this mass is irrelevant if small, used in ME only
74  ELSEIF (idfx.EQ.-11) then
75  idf=-2 ! denotes tau -2 tau-
76  amfin=0.511d-3 !this mass is irrelevant if small, used in ME only
77  ELSEIF (idfx.EQ. 15) then
78  idf=2 ! denotes tau +2 tau-
79  amfin=1.77703 !this mass is irrelevant if small, used in ME only
80  ELSEIF (idfx.EQ.-15) then
81  idf=-2 ! denotes tau -2 tau-
82  amfin=1.77703 !this mass is irrelevant if small, used in ME only
83  ELSE
84  WRITE(*,*) 'INITWKSW: WRONG IDFX'
85  stop
86  ENDIF
87 
88  IF (idex.EQ. 11) then !electron
89  ide= 2
90  amin=0.511d-3
91  ELSEIF (idex.EQ.-11) then !positron
92  ide=-2
93  amin=0.511d-3
94  ELSEIF (idex.EQ. 13) then !mu+
95  ide= 2
96  amin=0.105659
97  ELSEIF (idex.EQ.-13) then !mu-
98  ide=-2
99  amin=0.105659
100  ELSEIF (idex.EQ. 1) then !d
101  ide= 4
102  amin=0.05
103  ELSEIF (idex.EQ.- 1) then !d~
104  ide=-4
105  amin=0.05
106  ELSEIF (idex.EQ. 2) then !u
107  ide= 3
108  amin=0.02
109  ELSEIF (idex.EQ.- 2) then !u~
110  ide=-3
111  amin=0.02
112  ELSEIF (idex.EQ. 3) then !s
113  ide= 4
114  amin=0.3
115  ELSEIF (idex.EQ.- 3) then !s~
116  ide=-4
117  amin=0.3
118  ELSEIF (idex.EQ. 4) then !c
119  ide= 3
120  amin=1.3
121  ELSEIF (idex.EQ.- 4) then !c~
122  ide=-3
123  amin=1.3
124  ELSEIF (idex.EQ. 5) then !b
125  ide= 4
126  amin=4.5
127  ELSEIF (idex.EQ.- 5) then !b~
128  ide=-4
129  amin=4.5
130  ELSEIF (idex.EQ. 12) then !nu_e
131  ide= 1
132  amin=0.1d-3
133  ELSEIF (idex.EQ.- 12) then !nu_e~
134  ide=-1
135  amin=0.1d-3
136  ELSEIF (idex.EQ. 14) then !nu_mu
137  ide= 1
138  amin=0.1d-3
139  ELSEIF (idex.EQ.- 14) then !nu_mu~
140  ide=-1
141  amin=0.1d-3
142  ELSEIF (idex.EQ. 16) then !nu_tau
143  ide= 1
144  amin=0.1d-3
145  ELSEIF (idex.EQ.- 16) then !nu_tau~
146  ide=-1
147  amin=0.1d-3
148 
149  ELSE
150  WRITE(*,*) 'INITWKSW: WRONG IDEX'
151  stop
152  ENDIF
153 
154 C ----------------------------------------------------------------------
155 C
156 C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
157 C
158 C called by : KORALZ
159 C ----------------------------------------------------------------------
160  itce=ide/iabs(ide)
161  jtce=(1-itce)/2
162  itcf=idf/iabs(idf)
163  jtcf=(1-itcf)/2
164  CALL t_givizo( ide, 1,aizor,qe,kdumm)
165  CALL t_givizo( ide,-1,aizol,qe,kdumm)
166  xupgi(1)=qe
167  xupgi(2)=qe
168  t3e = (aizol+aizor)/2.
169  xupzi(1)=(aizor-qe*(swsq+deltsq)*gsw(3)-qe*deltv)/sqrt(swsq*(1-swsq))
170  xupzi(2)=(aizol-qe*(swsq+deltsq)*gsw(3)-qe*deltv)/sqrt(swsq*(1-swsq))
171  ve =(xupzi(1)+xupzi(2))/2.
172  CALL t_givizo( idf, 1,aizor,qf,kolor)
173  CALL t_givizo( idf,-1,aizol,qf,kolor)
174  xupgf(1)=qf
175  xupgf(2)=qf
176  t3f = (aizol+aizor)/2.
177  xupzf(1)=(aizor-qf*(swsq+deltsq)*gsw(2)-qf*deltv)/sqrt(swsq*(1-swsq))
178  xupzf(2)=(aizol-qf*(swsq+deltsq)*gsw(2)-qf*deltv)/sqrt(swsq*(1-swsq))
179  vf =(xupzf(1)+xupzf(2))/2.
180 
181 * Coupling costants times EW form-factors
182  deno = dsqrt(swsq*(1d0-swsq))
183  ! Ve = (2*T3e -4*Qe*m_Sw2*CorEle)/Deno
184  ! Vf = (2*T3f -4*Qf*m_Sw2*CorFin)/Deno
185  ! Ae = 2*T3e /Deno
186  ! Af = 2*T3f /Deno
187 * Angle dependent double-vector extra-correction
188  vvcef = ( (t3e) *(t3f)
189  $ -(qe*swsq+deltsq) *(t3f) *gsw(3) -qe*(t3f)*deltv
190  $ -(qf*swsq+deltsq) *(t3e) *gsw(2) -qf*(t3e)*deltv
191  $ + (qe*swsq) *(qf*swsq) *gsw(4)
192  $ + 2*qe*qf*deltsq*swsq + 2*qe*qf*deltv*swsq )/deno**2
193 
194  vvcor = 1d0
195  IF(keygsw.EQ.1) THEN
196  vvcor = vvcef/(ve*vf)
197  ENDIF
198 C
199 C PRINT *,' initwksw VVCor = ', VVCor
200  ndiag0=2
201  ndiaga=11
202  keya = 1
203  keyz = 1
204 C
205 C
206  RETURN
207  END
208  FUNCTION t_bornew(MODE,KEYGSW,SVAR,COSTHE,TA,TB)
209 C ----------------------------------------------------------------------
210 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
211 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
212 C EXAMPLE OF THE METHOD APPLIED THERE
213 C INPUT PARAMETERS ARE: SVAR -- transfer
214 C COSTHE -- cosine of angle between tau+ and 1st beam
215 C TA,TB -- helicity states of tau+ tau-
216 C mode -- parameter for mass terms; 1 means mass terms are on.
217 C keyGSW -- keyGSW=0 gamma propagator is off
218 C keyGSW=10 running Z width
219 C
220 C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
221 C ----------------------------------------------------------------------
222  IMPLICIT REAL*8(a-h,o-z)
223  COMMON / t_beampm / ene ,amin,amfin,ide,idf
224  REAL*8 ene ,amin,amfin
225  COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
226  & ,xupgi0 ,xupzi0 ,xupgf0 ,xupzf0
227  & ,ndiag0,ndiaga,keya,keyz
228  & ,itce,jtce,itcf,jtcf,kolor
229  REAL*8 ss,poln,t3e,qe,t3f,qf
230  & ,xupgi0(2),xupzi0(2),xupgf0(2),xupzf0(2)
231  COMMON / t_gauspm1/vvcor, zetvpi, gamvpi
232  & ,xupgi ,xupzi ,xupgf ,xupzf
233  COMPLEX*16 vvcor, zetvpi, gamvpi
234  COMPLEX*16 xupgi(2),xupzi(2),xupgf(2),xupzf(2)
235  COMMON / t_ewn / gmun, alphainvn
236  REAL*8 gmun, alphainvn
237 
238 
239  REAL*8 seps1,seps2
240 C=====================================================================
241  COMMON / t_gswprmn /swsq,amw,amz,amh,amtop,gammz
242  REAL*8 swsq,amw,amz,amh,amtop,gammz
243 C SWSQ = sin2 (theta Weinberg)
244 C AMW,AMZ = W & Z boson masses respectively
245 C AMH = the Higgs mass
246 C AMTOP = the top mass
247 C GAMMZ = Z0 width
248  COMPLEX*16 aborn(2,2),aphot(2,2),azett(2,2)
249  COMPLEX*16 xupzfp(2),xupzip(2),xupzif(2,2)
250  COMPLEX*16 abornm(2,2),aphotm(2,2),azettm(2,2)
251  COMPLEX*16 propa,propz
252  COMPLEX*16 xr,xi
253  COMPLEX*16 xupf,xupi
254  COMPLEX*16 xthing
255  DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
256  DATA mode0 /-5/
257  DATA ide0 /-55/
258  DATA svar0,cost0 /-5.d0,-6.d0/
259  DATA pi /3.141592653589793238462643d0/
260  DATA seps1,seps2 /0d0,0d0/
261 
262 C
263 C MEMORIZATION =========================================================
264  IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
265  $ .OR.ide0.NE.ide)THEN
266 C
267 
268  ! PRINT *,' T_BORN EW loop ( ',sqrt(svar),XUPGI(1),')= ', VVcor, ZetVPi!, GamVPi
269  ! PRINT *,' T_BORN new( ',mode,')= ',SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
270 C ** SWITCH OF MEMORISATION
271 C IDE0=IDE
272 C MODE0=MODE
273 C SVAR0=SVAR
274 C COST0=COSTHE
275 C ** PROPAGATORS
276  sinthe=sqrt(1.d0-costhe**2)
277  beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
278 ! BETA=1.D0! Dec 10, 2019 mass term may need to be killed for EW tests
279 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
280  xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5d0*beta*(xupzf(1)-xupzf(2))
281  xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5d0*beta*(xupzf(1)-xupzf(2))
282  xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5d0*(xupzi(1)-xupzi(2))
283  xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5d0*(xupzi(1)-xupzi(2))
284  xupzif(1,1)=(0.5d0*(xupzi(1)+xupzi(2))+0.5d0*(xupzi(1)-xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2))+0.5d0*beta*(xupzf(1)-xupzf(2)))
285  $ +(0.5d0*(xupzi(1)+xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2)))*(vvcor-1)
286  xupzif(1,2)=(0.5d0*(xupzi(1)+xupzi(2))+0.5d0*(xupzi(1)-xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2))-0.5d0*beta*(xupzf(1)-xupzf(2)))
287  $ +(0.5d0*(xupzi(1)+xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2)))*(vvcor-1)
288  xupzif(2,1)=(0.5d0*(xupzi(1)+xupzi(2))-0.5d0*(xupzi(1)-xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2))+0.5d0*beta*(xupzf(1)-xupzf(2)))
289  $ +(0.5d0*(xupzi(1)+xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2)))*(vvcor-1)
290  xupzif(2,2)=(0.5d0*(xupzi(1)+xupzi(2))-0.5d0*(xupzi(1)-xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2))-0.5d0*beta*(xupzf(1)-xupzf(2)))
291  $ +(0.5d0*(xupzi(1)+xupzi(2)))*(0.5d0*(xupzf(1)+xupzf(2)))*(vvcor-1)
292 
293 C FINAL STATE VECTOR COUPLING
294  xupf =0.5d0*(xupzf(1)+xupzf(2))
295  xupi =0.5d0*(xupzi(1)+xupzi(2))
296  xthing =0d0
297 
298 
299  propa =1d0/svar*gamvpi
300 C use running width
301  propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)*zetvpi
302  alphainv=alphainvn
303 
304  IF( keygsw. eq. 2.OR.keygsw.EQ.0) THEN
305  propz =1d0/dcmplx(svar-amz**2,amz*gammz)*zetvpi
306  ELSEIF( keygsw. eq. 10) THEN
307 ! PROPZ =1D0/DCMPLX(SVAR-AMZ**2,AMZ*GAMMZ)*ZetV ! this form need
308 ! redefined M_Z and G_Z.
309 ! Below variant with this rescaling implemented
310  propz =1d0/dcmplx(svar-amz**2/(1+gammz**2/amz**2), ! alternative to
311  $ amz*gammz /(1+gammz**2/amz**2)) ! running width
312  $ *zetvpi
313  propz =propz*dcmplx(1,-gammz/amz/(1+gammz**2/amz**2))
314  ENDIF
315 
316 C needed for variants v1 and v2 of effective Born. Factors change normalization of
317 C t_born by (alpha(M_Z)/alpha(0))^2
318 
319  propz =propz/alphainv
320  propa =propa/alphainv
321 C to achieve in QED point-like low energy limit
322 C q_i^2q_f^2 (1+cos_theta^2) for Born
323  propz =propz*137.03604
324  propa =propa*137.03604
325 
326 
327  IF (keygsw.EQ.0) propa=0.d0
328  DO 50 i=1,2
329  DO 50 j=1,2
330  regula= (3-2*i)*(3-2*j) + costhe
331  regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
332  aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
333  azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
334  azett(i,j)=propz*(xupzif(i,j)+xthing)*regula ! with electroweak effects in.
335  aborn(i,j)=aphot(i,j)+azett(i,j)
336  aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
337  azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
338  abornm(i,j)=aphotm(i,j)+azettm(i,j)
339  50 CONTINUE
340  ENDIF
341 C
342 C******************
343 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
344 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
345 C* HELICITY CONSERVATION EXPLICITLY OBEYED
346  polar1= (seps1)
347  polar2= (-seps2)
348  born=0d0
349  DO 150 i=1,2
350  helic= 3-2*i
351  DO 150 j=1,2
352  helit=3-2*j
353  factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
354  factom=factor*(1+helit*ta)*(1-helit*tb)
355  factor=factor*(1+helit*ta)*(1+helit*tb)
356 
357  born=born+cdabs(aborn(i,j))**2*factor
358 C MASS TERM IN BORN
359  IF (mode.GE.1) THEN
360  born=born+cdabs(abornm(i,j))**2*factom
361  ENDIF
362 
363  150 CONTINUE
364 C************
365  funt=born
366  IF(funt.LT.0.d0) funt=born
367 
368 C
369  IF (svar.GT.4d0*amfin**2) THEN
370 C PHASE SPACE THRESHOLD FACTOR and SVAR**2
371  thresh=sqrt(1-4d0*amfin**2/svar)
372  t_bornew= funt*thresh *svar**2
373  ELSE
374  thresh=0.d0
375  t_bornew=0.d0
376  ENDIF
377  END
378 
379  FUNCTION t_gamm(MODE,SVAR,COSTHE,TA,TB)
380  IMPLICIT REAL*8(a-h,o-z)
381  real*8 r(1:4, 1:4)
382  COMMON / t_beampm / ene ,amin,amfin,ide,idf
383  REAL*8 ene ,amin,amfin
384  common /dipolegam/ adip,bdip,iqeddip,ifgammold
385  ifgammold=1
386  a=0.0
387  b=0.0
388  iqed=iqeddip ! on/off for SM dipole in gamma gamma to tau tau
389  theta=acos(costhe)
390  e=sqrt(svar)/2.
391  IF (ifgammold.eq.0) then
392  t_gamm=0.d0
393  return
394  endif
395 
396  call dipolgammarij(iqed, e, theta, a, b, r)
397  funt=r(4,4)+ ta*r(4,3)+ tb*r(3,4)+ ta*tb*r(3,3)
398  IF (svar.GT.4d0*amfin**2) THEN
399 C PHASE SPACE THRESHOLD FACTOR and SVAR**2
400  thresh=sqrt(1-4d0*amfin**2/svar)
401  t_gamm= funt*thresh *svar**2
402  ELSE
403  thresh=0.d0
404  t_gamm=0.d0
405  ENDIF
406  write(*,*) 't_gamm stary',r(4,4), ta,r(1,1),tb,r(2,2),r(3,3)
407  END
408  FUNCTION t_gammnew(MODE,SVAR,COSTHE,TA,TB)
409  IMPLICIT REAL*8(a-h,o-z)
410  real*8 r(1:4, 1:4)
411  COMMON / t_beampm / ene ,amin,amfin,ide,idf
412  REAL*8 ene ,amin,amfin
413  common /dipolegam/ adip,bdip,iqeddip,ifgammold
414 
415  a=adip
416  b=bdip
417  iqed=1 ! SM dipole always on in SM+NP gamma gamma to tau tau
418  theta=acos(costhe)
419  e=sqrt(svar)/2.
420  call dipolgammarij(iqed, e, theta, a, b, r)
421  funt=r(4,4)+ ta*r(4,3)+ tb*r(3,4)+ ta*tb*r(3,3)
422  IF (svar.GT.4d0*amfin**2) THEN
423 C PHASE SPACE THRESHOLD FACTOR and SVAR**2
424  thresh=sqrt(1-4d0*amfin**2/svar)
425  t_gammnew= funt*thresh *svar**2
426  ELSE
427  thresh=0.d0
428  t_gammnew=0.d0
429  ENDIF
430  t_gammnew=0.d0
431  write(*,*) 't_gamm nowy',t_gammnew
432  END
433 
434  Subroutine dipolgammarij (iqed, E, theta, A, B, R)
435 
436 c For gamma gamma -> tau- tau+ reaction
437 c
438 c Calculates spin-correlation coefficients R(i,j) as functions of beam
439 c energy E = sqrt(s)/2 (in GeV) and scattaring angle theta.
440 c V is velocity of tau lepton, gam is Lorentz factor, alpha is
441 c fine-structure constant.
442 c Anomalous magnetic moment A1 = ASM + A and electric dipole moment B are
443 c real constants.
444 c Parameters A, B describe 'New Physics'.
445 c ASM is anomalous magnetic moment of tau in Standard Model (SM):
446 c S.Eidelman and M.Passera. Mod.Phys.Lett. A22, 159-179, 2007.
447 c
448 c Order of coefficients:
449 c i = 1,2,3 correspond to S_x, S_y, S_z for tau-,
450 c j = 1,2,3 correspond to S'_x, S'_y, S'_z for tau+
451 c i,j = 4 (equivalent to tt) corresponds to 1 (no spin dependence).
452 
453  Implicit none
454  integer iqed
455 C integer i, j ! not used here
456  real*8 e, theta, a, b
457  real*8 r(1:4, 1:4)
458  real*8 asm, a1, gam, e4, v,d2
459  real*8 fnorm
460  real*8 pi/3.141592653589793238d0/,alpha/7.2973525693d-3/
461  real*8 mtau ! mass of tau in GeV
462  COMMON / t_beampm / ene ,amin,amfin,ide,idf
463  integer ide,idf
464  REAL*8 ene ,amin,amfin
465  integer k,l
466  mtau=amfin
467  v = dsqrt(1.d0 -(mtau/e)**2) ! tau velocity
468  e4 = (4.0d0 *pi *alpha)**2 ! (electric charge)^4
469  gam = e/mtau ! Lorentz factor for tau
470  fnorm= 2*(2*pi**2*alpha**2) ! for normalization as in tauspinner.
471 c Cotribution to anomalous magnetic moment in Standard Model:
472  asm = 1.17721d-3
473  asm = asm *iqed ! switch for dipole SM part
474  a1 = asm + a ! Standard Model + New Physics
475 
476  d2 = (v**2 *dcos(theta)**2 -1.d0)**2 ! factor in denominator
477 
478  r(1,1) = e4 *(-11.d0 *v**4 +28.d0 *a1 *v**2 +
479  $ 4.d0 *(v**2 -2.d0) *dcos(2.d0*theta) *v**2 -
480  $ (v**2 -4.d0 *a1 -2.d0) *dcos(4.d0*theta) *v**2 +
481  $ 22.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)
482 
483  r(1,2) = e4 *v *b *(dcos(4.d0*theta) *v**2 +15.d0 *v**2 +
484  $ 4.d0 *cos(2.d0*theta) -20.d0) /(4.d0 *d2)
485 
486  r(1,3) = e4 *v**2 *gam *((a1 -1.d0) *v**2 +
487  $ (v**2 +a1 *(v**2 -2.d0) -1.d0) *dcos(2.d0*theta) +1.d0) *
488  $ dsin(2.d0*theta) /(2.d0 *d2)
489 
490  r(1,4) = 0.d0
491 
492  r(2,1) = -r(1,2)
493 
494  r(2,2) = e4 *(-dcos(4.d0*theta) *v**4 -11.d0 *v**4 +
495  $ 16.d0 *a1 *v**2 +4.d0 *(v**2 +4.d0 *a1) *dcos(2.d0*theta) *v**2 +
496  $ 16.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)
497 
498  r(2,3) = e4 *v *gam *b *(dcos(2.d0*theta)*v**2 -3.d0*v**2 +2.d0) *
499  $ dsin(2.d0*theta) /(2.d0 *d2)
500 
501  r(2,4) = 0.d0
502 
503  r(3,1) = r(1,3)
504 
505  r(3,2) = -r(2,3)
506 
507  r(3,3) = e4 *(-4.d0 *dcos(2.d0*theta) *v**4 +11.d0 *v**4 +
508  $ 36.d0 *a1 *v**2 +(v**2 -4.d0 *a1 -2.d0) *dcos(4.d0*theta) *v**2 +
509  $ 2.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)
510 
511  r(3,4) = 0.d0
512 
513  r(4,1) = 0.d0
514 
515  r(4,2) = 0.d0
516 
517  r(4,3) = 0.d0
518 
519  r(4,4) = e4 *(-dcos(4.d0*theta) *v**4 -11.d0 *v**4 -16.d0 *a1 *v**2 +
520  $ 4.d0 *(v**2 -4.d0 *a1 -2.d0) *dcos(2.d0*theta) *v**2 +
521  $ 8.d0 * v**2 +32.d0 *a1 +8.d0) /(8.d0 *d2)
522 
523 ! to normalize as in tauspinner
524  ! FNORM=1.0
525  do k=1,4
526  do l=1,4
527  r(k,l)=r(k,l)/fnorm
528  enddo
529  enddo
530 
531  return
532  end
533 
534  Subroutine dipolqq(iqed,Energy,theta,channel,Amz0,Gamz0,sin2W0,alphaQED,ReA0,ImA0, ReB,ImB, ReX,ImX,ReY,ImY,GSWr,GSWi,R)
535  implicit none
536  real*8 amz0,gamz0,sin2w0,alphaqed,gswr(7),gswi(7)
537  integer iqed,channel,k,j,k0,j0
538  real*8 energy,theta,rea,ima,rea0,ima0,reb,imb,rex,imx,rey,imy,xnor,thet, buf
539  real*8 r(1:4, 1:4),r0(1:4, 1:4),sign(4)
540 ! ReA0=0
541 ! ImA0=0
542 ! ReB=0
543 ! ImB=0
544 ! ReX=0
545 ! ImX=0
546 ! ReY=0
547 ! ImY=0
548  sign(1)=-1.0
549  sign(2)= 1.0
550  sign(3)=-1.0
551  sign(4)= 1.0
552  thet=acos(cos(theta))
553 
554  call dipolqqrijradcor(iqed,energy, thet,amz0,gamz0,sin2w0,alphaqed,
555  # gswr,gswi, rea0, ima0, reb, imb, rex, imx, rey, imy, r0, channel)
556  ! we normalize here, but we passresult for x-section as well.
557  xnor=r0(4,4)
558  do k=1,4
559  do j=1,4
560  r0(k,j)=r0(k,j)/xnor*sign(k)*sign(j)
561  enddo
562  enddo
563  r0(4,4)=xnor
564 
565  do k=1,4
566  do j=1,4
567  if(j.eq.1) then
568  buf=r0(k,2)
569  r0(k,2)=r0(k,j)
570  r0(k,1)=-buf
571  endif
572 
573  enddo
574  enddo
575 
576 
577  do j=1,4
578  do k=1,4
579  if(k.eq.1) then
580  buf=r0(2,j)
581  r0(2,j)=r0(k,j)
582  r0(1,j)=-buf
583  endif
584 
585  enddo
586  enddo
587 
588 
589 
590  do k0=1,4
591  do j0=1,4
592  k=k0+1
593  j=j0+1
594  if(j.eq.5) j=1
595  if(k.eq.5) k=1
596  r(k,j)=r0(k0,j0)
597  enddo
598  enddo
599  ! if(R0(3,3)/R0(4,4).gt.1d0) then
600  ! write(*,*) "z fortranu thet=", thet," energy= ",energy, " Rzz=", R0(3,3)/R0(4,4)," ", R0(3,3)," ",R0(4,4)
601 
602  ! write(*,*) "iqed,Energy, thet, channel ",iqed,Energy, thet, channel
603  ! write(*,*) "Amz0,Gamz0,sin2W0,alphaQED ", Amz0,Gamz0,sin2W0,alphaQED
604  ! write(*,*) "ReA0, ImA0, ReB, ImB ", ReA0, ImA0, ReB, ImB
605  ! write(*,*) "ReX, ImX, ReY, ImY ", ReX, ImX, ReY, ImY
606  ! write(*,*) "GSWr ", GSWr
607  ! write(*,*) "GSWi ", GSWi
608  ! endif
609  ! IF(R(4,4)/R(1,1).gt.1.0) then
610  ! write(*,*) "==============",R(4,4)/R(1,1)
611  ! write(*,*) 'energuy,theta,channel= ',Energy,' ',cos(theta),' ',channel
612  ! do k=1,4
613  ! write(*,7) k, R(k,1)/R(1,1),R(k,2)/R(1,1),R(k,3)/R(1,1),R(k,4)/R(1,1)
614  ! enddo
615  ! endif
616 
617 
618  7 format(1x,i4, 4(2x,f12.6))
619 
620  end
621 
622  Subroutine dipolqqrijradcor (iqed,Energy,theta,Amz0,Gamz0,sin2W0,alphaQED,
623  # gswr,gswi,rea0,ima0,reb,imb,rex,imx,rey,imy,r,channel)
624 C version of dipole routine with electroweak form factors included,
625 C If ported to different applications, beware of electroweak schemes, constant
626 C and form-factors initialization
627 C qed anomalous magnetic moment coded, but commented out
628 C (comments) #####################################
629 C 1) code is adopted to run with KKMC
630 C 2) and EW for channel=1 (initial state e, final state tau or mu.
631 C 3) internal activation of electroweak KKMC libs with IfGSW=KeyElw, call BornV_GetKeyElw(KeyElw)
632 C 4) note distinct conventions, papers: e-Print: 2307.03526 versus Eur.Phys.J.C 79 (2019) 6, 480
633 C that is about photon/Z normalization too
634 C ########################################################
635 
636 c PROCESS fermion_i + fermionbar_i -> fermion_f- + fermionbar_f
637 c
638 c IMPROVED BORN APPROXIMATION (IBA) with radiative corrections
639 c Integer parameter channel = 1 for LEPTONS, 2 -for UP QUARK, 3 - for DOWN QUARK
640 c
641 c Photon and Z-boson exchanges are included.
642 c Calculates spin-correlation coefficients R(i,j) as functions of beam
643 c energy: Energy = sqrt(s)/2 (in GeV), and scattaring angle theta.
644 c
645 c Anomalous magnetic formfactor A = ReA +i*ImA,
646 c and electric dipole formfactor B = ReB + i*ImB are complex.
647 c Also weak anomalous magnetic formfactor X = ReX + i*ImX and
648 c weak electric dipole formfactor Y = ReY + i*ImY are complex.
649 c
650 c V is velocity of tau, gam is Lorentz factor, alpha is fine-structure constant,
651 c G_F is Fermi constant of weak interaction.
652 c RSM(i,j) - spin-correlation matrix in IBA without anomalous moments.
653 c RDM(i,j) - spin-correlation matrix, linear in anomalous dipole moments A, B, X, Y.
654 c R(i,j) - total spin-correlation matrix.
655 c Order of coefficients R(ij): i = 1,2,3 correspond to S_x, S_y, S_z for tau-,
656 c j = 1,2,3 correspond to S'_x, S'_y, S'_z for tau+,
657 c i = j = 4 corresponds to 1 (no spin dependence).
658 c
659 c Functions below and their notation are from the paper:
660 c E. Richter-Was, Z. Was Eur. Phys.J. C (2019) 79, 480
661 c
662 c indices i = 1, 2, 3 correspond to lepton, up quark, down quark, respectively
663 c Gamma_vp - vacuum polarization factor (function of s)
664 c Rho_11, Rho_12, Rho_13 =rho_{lepton,i}(s,t) for i= L,Up,Down
665 c K_1,K_2,K_3 = K_i(s,t) for i= L, Up, Down
666 c K_11,K_12,K_13 = K_{lepton,i}(s,t) for i = L,Up,Down
667 
668  Implicit none
669 ! INCLUDE '../../basf2/KK2f/GPS.fi'
670  real*8 m_sw2,m_mz,m_gammz,m_gmu,m_alfinv,m_swsq,m_pi
671  ! INCLUDE '../../basf2/bornv/BornV.fi'
672 cc not needed external function K_1, K_2, K_3, K_11, K_12, K_13
673 cc not needed external function Rho_11,Rho_12,Rho_13,Gamma_vp
674  complex*16 a,b,x,y, pg, pz, den, vi, vf
675  complex*16 k1,k2,k3,k11,k12,k13,gammavp,rho11,rho12,rho13
676  DOUBLE COMPLEX gsw(100)
677  real*8 amz0,gamz0,sin2w0,alphaqed,gswr(7),gswi(7)
678  DOUBLE PRECISION svar,costhetd
679  DOUBLE COMPLEX rhoew, vpgamma, corele, corfin, corelefin, vvcef
680  integer i, j, channel,iqed,kff,ifgsw,keyelw,ifprint
681  real*8 energy,theta,rea,ima,rea0,ima0,reb,imb,rex,imx,rey,imy,gam,v
682  real*8 arqed, aiqed, ar1, ai1
683  real*8 e, f, m, qi, qf, ai, af, sw2, sw, cw2,cw, s, t, fvivf
684  real*8 r(1:4, 1:4), rsm(1:4, 1:4), rdm(1:4, 1:4)
685  real*8 pi/3.141592653589793238d0/,alpha/7.2973525693d-3/
686  real*8 mz/91.1876d0/,gz/2.4952d0/ ! Mass and width of Z in GeV
687 cc real*8 v0/-0.03783d0/,a0/-0.50123d0/ ! vector and axial couplings of leptons
688  real*8 mtau/1.77686d0/ ! mass of tau in GeV
689 
690  real*8 g_f/1.1663788d-5/ ! from PDG Fermi constant, GeV^{-2}
691  real*8 g_mu/1.166389d-5/,mw/80.353d0/ ! G_mu and M_W from Eur.Phys.J. C (2019) 79, 480
692 
693 C ==================================================================
694 C definition of form-factors as <in-line> functions. Complex form-factors
695 C initialized from KKMC libraries will be used. Note that electroweak
696 C scheme need to be reconsidered if function is used for other purposes.
697  k3(s,t) = corele
698  k1(s,t) = corele
699  k13(s,t) = corelefin
700  k2(s,t) = corele
701  k12(s,t) = corelefin
702  k11(s,t) = corelefin
703 
704  gammavp(s) = vpgamma
705  rho11(s,t) = rhoew
706  rho12(s,t) = rhoew
707  rho13(s,t) = rhoew
708 C ===================================================================
709 
710 
711  kff=15 ! 415-400
712  svar=4*energy**2
713  costhetd=cos(theta)
714 ! RhoEW = GSW(1)
715 ! VPgamma = GSW(6)
716 ! CorEle = GSW(2)
717 ! CorFin = GSW(3)
718 ! CorEleFin = GSW(4)
719  ifgsw=1 ! switch to activate electroweak formfactor
720 ! call BornV_GetKeyElw(KeyElw)
721  ifgsw=keyelw
722  ifgsw=1
723  m_pi=pi
724  m_gmu=g_mu
725  ifprint=0 ! switch to activate prints
726  If(ifgsw.eq.1) then
727 ! CALL BornV_InterpoGSW(KFf,Svar,CosThetD)
728 ! CALL BornV_GetGSW(GSW)
729  rhoew = dcmplx(gswr(1),gswi(1))
730  vpgamma = dcmplx(gswr(7),gswi(7))
731  vpgamma = 1.d0 /(2.d0-dcmplx(gswr(7),gswi(7)))
732  corele = dcmplx(gswr(2),gswi(2))
733  corfin = dcmplx(gswr(3),gswi(3))
734  corelefin = dcmplx(gswr(4),gswi(4))
735  sw2 = sin2w0 ! m_Sw2 ! 0.22351946d0 ! sin(theta_w)^2 from Eur.Phys.J. C (2019) 79, 480
736  m_swsq=sw2
737  m_sw2=sw2
738 ! alpha=1.d0/m_AlfInv
739  alpha=alphaqed
740  m_alfinv=1.d0/alpha
741  mz=amz0 ! m_MZ
742  m_mz=mz
743  gz=m_gammz*svar/mz**2 ! running width
744  m_gammz=gamz0
745  gz=gamz0*svar/mz**2 ! running width
746  if (ifprint.eq.1) then
747  write(*,*) ' m_KeyElw= ', keyelw
748  write(*,*) 'sw2,alpha,Mz,Gz,svar=', sw2,alpha,mz,gz,svar
749  write(*,*) 'G_mu, m_gmu= ', g_mu, m_gmu
750  write(*,*) 'RhoEW = ',rhoew!, GSWr(1)
751  write(*,*) 'VPgamma = ',vpgamma!, GSWr(6)
752  write(*,*) 'gsw(6) (7)= ', gswr(6), gswr(7)
753  write(*,*) 'CorEle = ', corele!, GSWr(2)
754  write(*,*) 'CorFin = ',corfin !, GSWr(3)
755  write(*,*) 'CorEleFin = ',corelefin!, GSWr(5)
756  write(*,*) 'CorEleFin -CorEle*CorFin = ',corelefin-corele*corfin
757 ! stop
758  endif
759  else
760  rhoew = 1.d0
761  vpgamma = 1.d0
762  corele = 1.d0
763  corfin = 1.d0
764  corelefin = 1.d0
765  sw2 = 0.22351946d0 ! sin(theta_w)^2 from Eur.Phys.J. C (2019) 79, 480
766  m_swsq=sw2
767  endif
768 
769  v= dsqrt(1.d0 -(mtau/energy)**2) ! tau velocity
770 
771 c Contributions to ReA(s) and ImA(s) from QED in one loop
772  if(iqed.eq.10) then ! Warning: temporarily this contribution is blocked
773  arqed = -alpha*mtau**2 /(pi*v*4.d0*energy**2)*dlog((1.d0+v)/(1.d0-v))
774  aiqed = alpha *mtau**2 /(v *4.d0 *energy**2)
775  else
776  arqed = 0.d0 !switch for dipole QED part
777  aiqed = 0.d0
778  endif
779 
780 
781 
782  rea = arqed + rea0 ! QED + new physics
783  ima = aiqed + ima0 ! QED + new physics
784 
785 
786  m= mtau
787  gam= energy/mtau ! Lorentz factor of tau
788  e= dsqrt(4.0d0 *pi *alpha) ! electric charge
789 c Below coupling f=g_w/(2*cos(theta_w))=e/(2*cos(theta_w)*sin(theta_w)),
790 c which is expressed via Fermi constant:
791  f = mz *dsqrt(g_mu *dsqrt(2.d0))
792 
793 
794  sw = dsqrt(sw2) ! sin(theta_w)
795  cw2 = 1.d0 - sw2
796  cw = dsqrt(cw2)
797 
798 c Mandelstam variables s and t (mass of tau is included in t)
799  s = 4.d0 *energy**2
800  t = m**2 -0.5d0 *s *(1.d0 -v*dcos(theta))
801 
802 c Denominator of the Z propagator with running width
803  den = s - mz**2 + (0.d0, 1.d0) *gz *mz
804 c Constant for effective photon like couplings correction coming from Z boson
805  fvivf = 4.d0 *(f**2/e**2) *sw**4 ! WARNING: the *sw**4 factor come from vector couplings numerators Fvivf
806 c Note: Fvivf can also be written as Fvivf= sw**2/cw**2, or through the constant G_mu:
807 c Fvivf= sqrt(2) *G_mu *Mz**2 *sw**4 /(pi*alpha) This agrees with definition below.
808 
809 c FINAL FERMION IS TAU LEPTON (f=lepton) with couplings:
810  qf= -1.d0
811  vf= -0.03783d0 ! from PDG, used in previous code (v0 coupling)
812  af= -0.50123d0 ! from PDG, used in previous code (a0 coupling)
813  if(ifgsw.eq.1) then
814  vf= -0.5d0 - 2.d0 *qf *sw2 * k1(s,t) ! is function for lepto
815  af= -0.5d0
816  endif
817 
818 c
819 c COUPLINGS FOR INITIAL FERMIONS e u or d
820  IF (channel.EQ.1) THEN ! LEPTONS
821  qi= -1.d0
822  vi= -0.03783d0 ! from PDG, used in the previous code
823  ai= -0.50123d0 ! from PDG, used in the previous code
824 c
825 c Effective Z propagator (reduces to 1/Den without rad. corr.):
826  pz = rho11(s,t) /den
827  if(ifgsw.eq.1) then
828  vi= -0.5d0 - 2.d0 *qi *sw2 *k1(s,t)
829  ai= -0.5d0
830  if (ifprint.eq.1) write(*,*) 'no EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
831  fvivf=
832  $ m_gmu *m_mz**2 *m_alfinv /(dsqrt(2.d0)*8.d0*m_pi)
833  $ *(m_sw2*(1.d0-m_sw2)) *16.d0
834  fvivf=fvivf *m_sw2/(1.d0-m_sw2) ! why this factor? Because vi ai couplings in KKMC
835  ! are divided by deno= 4 sqrt(m_swsq*(1d0-m_swsq)) in addition multiplied by 2 and we are not using it here
836  ! Ve = (2*T3e -4*Qe*m_swsq)/deno and Ae = 2*T3e/deno.
837  if (ifprint.eq.1) then
838  write(*,*) 'with EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
839  stop
840  endif
841  endif
842 
843 C Effective photon propagator with v_{if}-v_i*v_f correction to Z exchange
844 C (NOTE: Pz reduces to 1/s if rad. corr. are off):
845  pg = 1.d0/s *(gammavp(s) + fvivf *s/den*
846  $ rho11(s,t)*(k11(s,t)-k1(s,t)*k1(s,t)))
847 
848 c
849  ELSE IF (channel.EQ.2) THEN ! UP QUARK
850  qi= +2.d0/3.d0
851  vi= +0.266d0 ! from PDG, used in the previous code
852  ai= +0.519d0 ! from PDG, used in the previous code
853  if(ifgsw.eq.1) then
854  vi= +0.5d0 - 2.d0 *qi *sw2 *k2(s,t)
855  ai= +0.5d0
856  if (ifprint.eq.1) write(*,*) 'no EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
857  fvivf=
858  $ m_gmu *m_mz**2 *m_alfinv /(dsqrt(2.d0)*8.d0*m_pi)
859  $ *(m_sw2*(1.d0-m_sw2)) *16.d0
860  fvivf=fvivf *m_sw2/(1.d0-m_sw2) ! why this factor? Because vi ai couplings in KKMC
861  ! are divided by deno= 4 sqrt(m_swsq*(1d0-m_swsq)) in addition multiplied by 2.
862  ! Ve = (2*T3e -4*Qe*m_swsq)/deno and Ae = 2*T3e/deno. Note that Fvivf is used only for calculation of
863  ! (vv_if-v_i*v_f)
864  if (ifprint.eq.1) then
865  write(*,*) 'with EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
866  stop
867  endif
868  endif
869 
870 c
871 c Effective Z propagator (reduces to 1/Den without rad. corr.):
872  pz = rho12(s,t) /den
873 c Effective photon propagator with v_{if}-v_i*v_f correction to Z exchange
874 C (NOTE: Pz reduces to 1/s if rad. corr. are off):
875  pg = 1.d0/s *(gammavp(s) + fvivf *s/den*
876  $ rho12(s,t)*(k12(s,t)-k1(s,t)*k2(s,t)))
877 c
878  ELSE IF (channel.EQ.3) THEN ! DOWN QUARK
879  qi= -1.d0/3.d0
880  vi= -0.38d0 ! from PDG, used in the previous code
881  ai= -0.527d0 ! from PDG, used in the previous code
882  if(ifgsw.eq.1) then
883  vi= -0.5d0 - 2.d0 *qi *sw2 *k3(s,t)
884  ai= -0.5d0
885  if (ifprint.eq.1) write(*,*) 'no EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
886  fvivf=
887  $ m_gmu *m_mz**2 *m_alfinv /(dsqrt(2.d0)*8.d0*m_pi)
888  $ *(m_sw2*(1.d0-m_sw2)) *16.d0
889  fvivf=fvivf *m_sw2/(1.d0-m_sw2) ! why this factor? Because vi ai couplings in KKMC
890  ! are divided by deno= 4 sqrt(m_swsq*(1d0-m_swsq)) in addition multiplied by 2.
891  ! Ve = (2*T3e -4*Qe*m_swsq)/deno and Ae = 2*T3e/deno. Note that Fvivf is used only for calculation of
892  ! (vv_if-v_i*v_f)
893  if (ifprint.eq.1) then
894  write(*,*) 'with EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
895  stop
896  endif
897  endif
898 
899 c
900 c Effective Z propagator (reduces to 1/Den without rad. corr.):
901  pz = rho13(s,t) /den
902 c Effective photon propagator (reduces to 1/s without rad. corr.):
903  pg = 1.d0/s *(gammavp(s) + fvivf *s/den*
904  $ rho13(s,t) *(k13(s,t)-k1(s,t)*k3(s,t)))
905 c
906  ELSE ! OUTPUT WILL BE ZERO
907  qi= 0.d0
908  vi= (0.d0, 0.d0)
909  ai= 0.d0
910  pz= (0.d0, 0.d0)
911  pg= (0.d0, 0.d0)
912  END IF
913 
914 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
915 c Complex magnetic and electric form-factors
916  a = dcmplx(rea, ima)
917  b = dcmplx(reb, imb)
918  x = dcmplx(rex, imx)
919  y = dcmplx(rey, imy)
920 
921 c Contributions to A(s), X(s) from QED in one loop: NOT INCLUDED
922 cc ArQED = -alpha *mtau**2 /(PI *V *4.d0 *E**2) *dlog((1.d0+V)/(1.d0-V))
923 cc AiQED = alpha *mtau**2 /(V *4.d0 *E**2)
924 cc Ar1 = ArQED + Ar ! QED + new physics
925 cc Ai1 = AiQED + Ai ! QED + new physics
926 
927 
928 
929 c ----------------------------------------------------------------
930 c CONTRIBUTIONS to R(i,j) IN STANDARD MODEL NO DIPOLE MOMENTS
931 c ----------------------------------------------------------------
932 c
933 
934  rsm(1,1)= 4.d0*gam**2*m**4*(e**2*(1.d0 + gam**2)*qf*qi*
935  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(pg) +
936  - f**2*dconjg(pz)*
937  - (-(af**2*f**2*(-1.d0 + gam**2)*pz*
938  - (ai**2 + vi*dconjg(vi))) +
939  - (1.d0 + gam**2)*dconjg(vf)*
940  - (ai**2*f**2*pz*vf +
941  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))))*
942  - dsin(theta)**2
943 
944  rsm(1,2)= (0.d0, -4.d0)*af*f**2*gam**4*m**4*v*
945  - (e**2*pz*qf*qi*vi*dconjg(pg) +
946  - dconjg(pz)*(-(ai**2*f**2*pz*vf) -
947  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi) +
948  - f**2*pz*dconjg(vf)*(ai**2 + vi*dconjg(vi))))*
949  - dsin(theta)**2
950 
951  rsm(2,1)= rsm(1,2)
952 
953  rsm(2,2)= 4.d0*gam**2*(-1.d0 + gam**2)*m**4*
954  - (-(e**2*qf*qi*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
955  - dconjg(pg)) +
956  - dconjg(pz)*(af**2*f**4*pz*
957  - (ai**2 + vi*dconjg(vi)) -
958  - f**2*dconjg(vf)*
959  - (ai**2*f**2*pz*vf +
960  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))))*
961  - dsin(theta)**2
962 
963  rsm(1,3)= 4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
964  - (af*ai*f**2*pz*v -
965  - 2.d0*gam**2*(-1.d0 + v**2)*
966  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dcos(theta)) +
967  - f**2*dconjg(pz)*
968  - (af*(ai*v*(e**2*pg*qf*qi +
969  - f**2*pz*vf*(vi + dconjg(vi))) +
970  - 2.d0*af*f**2*pz*(1.d0 + gam**2*(-1.d0 + v**2))*
971  - (ai**2 + vi*dconjg(vi))*dcos(theta)) +
972  - dconjg(vf)*
973  - (ai*f**2*pz*(af*v*vi -
974  - 2.d0*ai*gam**2*(-1.d0 + v**2)*vf*dcos(theta)) +
975  - dconjg(vi)*
976  - (af*ai*f**2*pz*v -
977  - 2.d0*gam**2*(-1.d0 + v**2)*
978  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dcos(theta))))
979  - )*dsin(theta)
980 
981  rsm(3,1)= rsm(1,3)
982 
983  rsm(2,3)= (0.d0, -2.d0)*af*f**2*gam**3*m**4*v*
984  - (e**2*pz*qf*qi*vi*dconjg(pg) +
985  - dconjg(pz)*(-(ai**2*f**2*pz*vf) -
986  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi) +
987  - f**2*pz*dconjg(vf)*(ai**2 + vi*dconjg(vi))))*
988  - dsin(2.d0*theta)
989 
990  rsm(3,2)= rsm(2,3)
991 
992  rsm(3,3)= 2.d0*gam**2*m**4*(-(e**4*pg*qf**2*qi**2*dconjg(pg)*
993  - (1.d0 - 3.d0*gam**2 +
994  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
995  - dcos(2.d0*theta))) -
996  - e**2*f**2*qf*qi*(pz*dconjg(pg)*
997  - (-4.d0*af*ai*gam**2*v*dcos(theta) +
998  - vf*vi*(1.d0 - 3.d0*gam**2 +
999  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1000  - dcos(2.d0*theta))) +
1001  - pg*dconjg(pz)*
1002  - (-4.d0*af*ai*gam**2*v*dcos(theta) +
1003  - dconjg(vf)*dconjg(vi)*
1004  - (1.d0 - 3.d0*gam**2 +
1005  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1006  - dcos(2.d0*theta)))) +
1007  - f**4*pz*dconjg(pz)*
1008  - (af*(4*ai*gam**2*v*vf*(vi + dconjg(vi))*
1009  - dcos(theta) +
1010  - af*(ai**2 + vi*dconjg(vi))*
1011  - (1.d0 + gam**2*(-1.d0 + 4.d0*v**2) +
1012  - (-1.d0 + 5.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1013  - dcos(2.d0*theta))) +
1014  - dconjg(vf)*
1015  - (dconjg(vi)*
1016  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1017  - vf*vi*(-1.d0 + 3.d0*gam**2 +
1018  - (1.d0 - 3.d0*gam**2 - 4.d0*gam**4*(-1.d0 + v**2))*
1019  - dcos(2.d0*theta))) +
1020  - ai*(4.d0*af*gam**2*v*vi*dcos(theta) -
1021  - ai*vf*(1.d0 - 3.d0*gam**2 +
1022  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1023  - dcos(2.d0*theta))))))
1024 
1025  rsm(1,4)= -4.d0*f**2*gam**3*m**4*
1026  - (e**2*pz*qf*qi*dconjg(pg)*
1027  - (2.d0*ai*vf + af*v*vi*dcos(theta)) +
1028  - dconjg(pz)*(af*v*
1029  - (ai**2*f**2*pz*vf +
1030  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1031  - dcos(theta) +
1032  - dconjg(vf)*
1033  - (2.d0*ai*(e**2*pg*qf*qi +
1034  - f**2*pz*vf*(vi + dconjg(vi))) +
1035  - af*f**2*pz*v*(ai**2 + vi*dconjg(vi))*
1036  - dcos(theta))))*dsin(theta)
1037 
1038  rsm(4,1)= rsm(1,4)
1039 
1040  rsm(2,4)= (0.d0, 4.d0)*af*ai*f**2*gam**3*m**4*v*
1041  - (e**2*pz*qf*qi*dconjg(pg) +
1042  - dconjg(pz)*(-(e**2*pg*qf*qi) -
1043  - f**2*pz*vf*(vi + dconjg(vi)) +
1044  - f**2*pz*dconjg(vf)*(vi + dconjg(vi))))*
1045  - dsin(theta)
1046 
1047  rsm(4,2)= rsm(2,4)
1048 
1049  rsm(3,4)= -4.d0*f**2*gam**2*m**4*
1050  - (e**2*gam**2*pz*qf*qi*dconjg(pg)*
1051  - (2.d0*ai*vf*dcos(theta) + af*v*vi*(1.d0 + dcos(theta)**2)) +
1052  - (dconjg(pz)*(dconjg(vi)*
1053  - (4.d0*af**2*ai*f**2*(-1.d0 + gam**2)*pz*dcos(theta) +
1054  - 4.d0*ai*f**2*gam**2*pz*vf*dconjg(vf)*
1055  - dcos(theta) +
1056  - af*gam**2*v*
1057  - (e**2*pg*qf*qi +
1058  - f**2*pz*vi*(vf + dconjg(vf)))*
1059  - (3.d0 + dcos(2.d0*theta))) +
1060  - ai*(gam**2*dconjg(vf)*
1061  - (4.d0*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1062  - dcos(theta) +
1063  - af*ai*f**2*pz*v*(3.d0 + dcos(2.d0*theta))) +
1064  - af*f**2*pz*
1065  - (4.d0*af*(-1.d0 + gam**2)*vi*dcos(theta) +
1066  - ai*gam**2*v*vf*(3.d0 + dcos(2.d0*theta))))))/2.d0)
1067 
1068  rsm(4,3)= rsm(3,4)
1069 
1070 
1071  rsm(4,4)= 2.d0*gam**2*m**4*(e**4*pg*qf**2*qi**2*dconjg(pg)*
1072  - (1.d0 + 3.d0*gam**2 + (-1.d0 + gam**2)*dcos(2.d0*theta)) +
1073  - e**2*f**2*qf*qi*(pz*dconjg(pg)*
1074  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1075  - vf*vi*(1.d0 + 3.d0*gam**2 + (-1.d0 + gam**2)*dcos(2.d0*theta))
1076  - ) + pg*dconjg(pz)*
1077  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1078  - dconjg(vf)*dconjg(vi)*
1079  - (1.d0 + 3.d0*gam**2 + (-1.d0 + gam**2)*dcos(2.d0*theta)))) -
1080  - f**4*pz*dconjg(pz)*
1081  - (-(af*(4.d0*ai*gam**2*v*vf*(vi + dconjg(vi))*
1082  - dcos(theta) +
1083  - af*(-1.d0 + gam**2)*(ai**2 + vi*dconjg(vi))*
1084  - (3.d0 + dcos(2.d0*theta)))) -
1085  - dconjg(vf)*
1086  - (ai*(4.d0*af*gam**2*v*vi*dcos(theta) +
1087  - ai*vf*(1.d0 + 3.d0*gam**2 +
1088  - (-1.d0 + gam**2)*dcos(2.d0*theta))) +
1089  - dconjg(vi)*
1090  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1091  - vf*vi*(1.d0 + 3.d0*gam**2 +
1092  - (-1.d0 + gam**2)*dcos(2.d0*theta))))))
1093 
1094 
1095 c -------------------------------------------------------------
1096 c CONTRIBUTION TO R(i,j), LINEAR IN DIPOLE MOMENTS A, B, X, Y
1097 c -------------------------------------------------------------
1098 
1099  rdm(1,1)= 8.d0*gam**4*m**4*(e**2*qf*qi*
1100  - (e**2*pg*qf*qi*(a + dconjg(a)) +
1101  - f**2*pz*vi*(x + vf*dconjg(a)))*dconjg(pg) +
1102  - f**2*dconjg(pz)*
1103  - (dconjg(vf)*(ai**2*f**2*pz*x +
1104  - (a*e**2*pg*qf*qi + f**2*pz*vi*x)*dconjg(vi)) +
1105  - (ai**2*f**2*pz*vf +
1106  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1107  - dconjg(x)))*dsin(theta)**2
1108 
1109  rdm(1,2)= 4.d0*gam**4*m**4*v*(e**2*qf*qi*
1110  - (e**2*pg*qf*qi*(b + dconjg(b)) +
1111  - f**2*pz*vi*(y - (0.d0, 1.d0)*af*dconjg(a) +
1112  - vf*dconjg(b)))*dconjg(pg) +
1113  - f**2*dconjg(pz)*
1114  - (dconjg(vf)*(ai**2*f**2*pz*y +
1115  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dconjg(vi)) +
1116  - (0.d0, 1.d0)*af*(dconjg(vi)*
1117  - (a*e**2*pg*qf*qi +
1118  - f**2*pz*vi*(x - dconjg(x))) +
1119  - ai**2*f**2*pz*(x - dconjg(x))) +
1120  - (ai**2*f**2*pz*vf +
1121  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1122  - dconjg(y)))*dsin(theta)**2
1123 
1124  rdm(2,1)= -4.d0*gam**4*m**4*v*(e**2*qf*qi*
1125  - (e**2*pg*qf*qi*(b + dconjg(b)) +
1126  - f**2*pz*vi*(y + (0,1)*af*dconjg(a) +
1127  - vf*dconjg(b)))*dconjg(pg) +
1128  - f**2*dconjg(pz)*
1129  - (dconjg(vf)*(ai**2*f**2*pz*y +
1130  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dconjg(vi)) -
1131  - (0.d0 ,1.d0)*af*(dconjg(vi)*
1132  - (a*e**2*pg*qf*qi +
1133  - f**2*pz*vi*(x - dconjg(x))) +
1134  - ai**2*f**2*pz*(x - dconjg(x))) +
1135  - (ai**2*f**2*pz*vf +
1136  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1137  - dconjg(y)))*dsin(theta)**2
1138 
1139  rdm(2,2)= 0.d0
1140 
1141  rdm(1,3)= 4.d0*gam**5*m**4*(-(e**2*qf*qi*dconjg(pg)*
1142  - ((0.d0, -1.d0)*ai*f**2*pz*v*
1143  - (y - (0.d0, 1.d0)*af*dconjg(a) - vf*dconjg(b)) +
1144  - (e**2*pg*qf*qi*(-2 + v**2)*(a + dconjg(a)) +
1145  - f**2*pz*vi*
1146  - ((-2.d0 + v**2)*(x + vf*dconjg(a)) +
1147  - (0.d0, 1.d0)*af*v**2*dconjg(b)))*dcos(theta))) +
1148  - f**2*dconjg(pz)*
1149  - (dconjg(vf)*(ai*
1150  - ((0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) -
1151  - ai*f**2*pz*(-2.d0 + v**2)*x*dcos(theta)) +
1152  - dconjg(vi)*
1153  - ((0.d0, 1.d0)*ai*f**2*pz*v*y -
1154  - (-2.d0 + v**2)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1155  - dcos(theta))) +
1156  - dconjg(vi)*
1157  - ((0.d0,-1.d0)*ai*f**2*pz*v*vf*dconjg(y) -
1158  - (-2.d0 + v**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1159  - dconjg(x)*dcos(theta) +
1160  - af*v*(ai*f**2*pz*(x + dconjg(x)) +
1161  - (0.d0, 1.d0)*v*
1162  - (b*e**2*pg*qf*qi +
1163  - f**2*pz*vi*(y - dconjg(y)))*dcos(theta)))
1164  - + ai*((0.d0, -1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1165  - dconjg(y) +
1166  - ai*f**2*pz*(2.d0 - v**2)*vf*dconjg(x)*
1167  - dcos(theta) +
1168  - af*v*(a*e**2*pg*qf*qi +
1169  - f**2*pz*vi*(x + dconjg(x)) +
1170  - (0.d0, 1.d0)*ai*f**2*pz*v*(y - dconjg(y))*
1171  - dcos(theta)))))*dsin(theta)
1172 
1173  rdm(3,1)= 4.d0*gam**5*m**4*(e**2*qf*qi*dconjg(pg)*
1174  - (ai*f**2*pz*v*(af*dconjg(a) -
1175  - (0.d0, 1.d0)*(y - vf*dconjg(b))) -
1176  - (e**2*pg*qf*qi*(-2.d0 + v**2)*(a + dconjg(a)) +
1177  - f**2*pz*vi*
1178  - ((-2.d0 + v**2)*(x + vf*dconjg(a)) -
1179  - (0.d0, 1.d0)*af*v**2*dconjg(b)))*dcos(theta)) +
1180  - f**2*dconjg(pz)*
1181  - (dconjg(vf)*(-(ai*
1182  - ((0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) +
1183  - ai*f**2*pz*(-2.d0 + v**2)*x*dcos(theta))) +
1184  - dconjg(vi)*
1185  - ((0.d0, -1.d0)*ai*f**2*pz*v*y -
1186  - (-2.d0 + v**2)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1187  - dcos(theta))) +
1188  - dconjg(vi)*
1189  - ((0.d0, 1.d0)*ai*f**2*pz*v*vf*dconjg(y) -
1190  - (-2.d0 + v**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1191  - dconjg(x)*dcos(theta) +
1192  - af*v*(ai*f**2*pz*(x + dconjg(x)) -
1193  - (0.d0, 1.d0)*v*
1194  - (b*e**2*pg*qf*qi +
1195  - f**2*pz*vi*(y - dconjg(y)))*dcos(theta)))
1196  - + ai*((0.d0,1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1197  - dconjg(y) +
1198  - ai*f**2*pz*(2.d0 - v**2)*vf*dconjg(x)*
1199  - dcos(theta) +
1200  - af*v*(a*e**2*pg*qf*qi +
1201  - f**2*pz*vi*(x + dconjg(x)) -
1202  - (0.d0, 1.d0)*ai*f**2*pz*v*(y - dconjg(y))*
1203  - dcos(theta)))))*dsin(theta)
1204 
1205  rdm(2,3)= (0.d0, 4.d0)*gam**5*m**4*v*
1206  - (e**2*qf*qi*dconjg(pg)*
1207  - (ai*f**2*pz*v*(x - vf*dconjg(a) +
1208  - (0.d0, 1.d0)*af*dconjg(b)) +
1209  - (0.d0, 1.d0)*(e**2*pg*qf*qi*(b + dconjg(b)) +
1210  - f**2*pz*vi*
1211  - (y + (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b)))*
1212  - dcos(theta)) +
1213  - f**2*dconjg(pz)*
1214  - (dconjg(vf)*(ai*
1215  - (a*e**2*pg*qf*qi*v + f**2*pz*v*vi*x +
1216  - (0.d0, 1.d0)*ai*f**2*pz*y*dcos(theta)) +
1217  - dconjg(vi)*
1218  - (ai*f**2*pz*v*x +
1219  - (0.d0, 1.d0)*(b*e**2*pg*qf*qi + f**2*pz*vi*y)*
1220  - dcos(theta))) +
1221  - (0.d0, 1.d0)*(dconjg(vi)*
1222  - ((0.d0, 1.d0)*ai*f**2*pz*v*vf*dconjg(x) +
1223  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(y)*
1224  - dcos(theta) +
1225  - af*(ai*f**2*pz*v*(y + dconjg(y)) -
1226  - (0.d0, 1.d0)*
1227  - (a*e**2*pg*qf*qi +
1228  - f**2*pz*vi*(x - dconjg(x)))*
1229  - dcos(theta))) +
1230  - ai*((0.d0, 1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1231  - dconjg(x) +
1232  - ai*f**2*pz*vf*dconjg(y)*dcos(theta) +
1233  - af*(v*(b*e**2*pg*qf*qi +
1234  - f**2*pz*vi*(y + dconjg(y))) -
1235  - (0.d0, 1.d0)*ai*f**2*pz*(x - dconjg(x))*
1236  - dcos(theta))))))*dsin(theta)
1237 
1238  rdm(3,2)= 4.d0*gam**5*m**4*v*(e**2*qf*qi*dconjg(pg)*
1239  - (ai*f**2*pz*v*((0.d0, 1.d0)*(x - vf*dconjg(a)) +
1240  - af*dconjg(b)) +
1241  - (e**2*pg*qf*qi*(b + dconjg(b)) +
1242  - f**2*pz*vi*
1243  - (y - (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b)))*
1244  - dcos(theta)) +
1245  - f**2*dconjg(pz)*
1246  - (dconjg(vf)*(ai*
1247  - ((0.d0, 1.d0)*v*(a*e**2*pg*qf*qi + f**2*pz*vi*x) +
1248  - ai*f**2*pz*y*dcos(theta)) +
1249  - dconjg(vi)*
1250  - ((0.d0, 1.d0)*ai*f**2*pz*v*x +
1251  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dcos(theta))) +
1252  - dconjg(vi)*
1253  - ((0.d0, -1.d0)*ai*f**2*pz*v*vf*dconjg(x) +
1254  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(y)*
1255  - dcos(theta) +
1256  - af*(ai*f**2*pz*v*(y + dconjg(y)) +
1257  - (0.d0, 1.d0)*(a*e**2*pg*qf*qi +
1258  - f**2*pz*vi*(x - dconjg(x)))*dcos(theta)))
1259  - + ai*((0.d0, -1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1260  - dconjg(x) +
1261  - ai*f**2*pz*vf*dconjg(y)*dcos(theta) +
1262  - af*(v*(b*e**2*pg*qf*qi +
1263  - f**2*pz*vi*(y + dconjg(y))) +
1264  - (0.d0, 1.d0)*ai*f**2*pz*(x - dconjg(x))*dcos(theta)
1265  - ))))*dsin(theta)
1266 
1267  rdm(3,3)= -8.d0*gam**6*m**4*(-1.d0 + v**2)*dcos(theta)*
1268  - (e**2*qf*qi*dconjg(pg)*
1269  - ((a*e**2*pg*qf*qi + f**2*pz*vi*x)*dcos(theta) +
1270  - dconjg(a)*(af*ai*f**2*pz*v +
1271  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dcos(theta))) +
1272  - f**2*dconjg(pz)*
1273  - (ai*(af*v*(a*e**2*pg*qf*qi +
1274  - f**2*pz*vi*(x + dconjg(x))) +
1275  - ai*f**2*pz*x*dconjg(vf)*dcos(theta) +
1276  - ai*f**2*pz*vf*dconjg(x)*dcos(theta)) +
1277  - dconjg(vi)*
1278  - (af*ai*f**2*pz*v*(x + dconjg(x)) +
1279  - ((a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1280  - dconjg(vf) +
1281  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(x))*
1282  - dcos(theta))))
1283 
1284  rdm(1,4)= -4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1285  - (ai*f**2*pz*((1 + gam**2)*(x + vf*dconjg(a)) -
1286  - (0.d0, 1.d0)*af*(-1.d0 + gam**2)*dconjg(b)) +
1287  - (0.d0, 1.d0)*gam**2*v*
1288  - (e**2*pg*qf*qi*(b - dconjg(b)) +
1289  - f**2*pz*vi*
1290  - (y - (0.d0, 1.d0)*af*dconjg(a) - vf*dconjg(b)))*
1291  - dcos(theta)) +
1292  - f**2*dconjg(pz)*
1293  - (dconjg(vf)*(ai*
1294  - ((1.d0 + gam**2)*
1295  - (a*e**2*pg*qf*qi + f**2*pz*vi*x) +
1296  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*y*dcos(theta)) +
1297  - dconjg(vi)*
1298  - (ai*f**2*(1.d0 + gam**2)*pz*x +
1299  - (0.d0, 1.d0)*gam**2*v*
1300  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dcos(theta)))
1301  - + ai*((1.d0 + gam**2)*
1302  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(x) -
1303  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*vf*dconjg(y)*
1304  - dcos(theta) +
1305  - af*((0.d0, 1.d0)*(-1.d0 + gam**2)*
1306  - (b*e**2*pg*qf*qi +
1307  - f**2*pz*vi*(y - dconjg(y))) +
1308  - ai*f**2*gam**2*pz*v*(x + dconjg(x))*
1309  - dcos(theta))) +
1310  - dconjg(vi)*
1311  - (ai*f**2*(1.d0 + gam**2)*pz*vf*dconjg(x) -
1312  - (0.d0, 1.d0)*gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1313  - dconjg(y)*dcos(theta) +
1314  - af*((0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*
1315  - (y - dconjg(y)) +
1316  - gam**2*v*
1317  - (a*e**2*pg*qf*qi +
1318  - f**2*pz*vi*(x + dconjg(x)))*dcos(theta)))
1319  - ))*dsin(theta)
1320 
1321  rdm(4,1)= -4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1322  - (ai*f**2*pz*((1.d0 + gam**2)*(x + vf*dconjg(a)) +
1323  - (0.d0, 1.d0)*af*(-1.d0 + gam**2)*dconjg(b)) -
1324  - (0.d0, 1.d0)*gam**2*v*
1325  - (e**2*pg*qf*qi*(b - dconjg(b)) +
1326  - f**2*pz*vi*
1327  - (y + (0.d0, 1.d0)*af*dconjg(a) - vf*dconjg(b)))*
1328  - dcos(theta)) +
1329  - f**2*dconjg(pz)*
1330  - (dconjg(vf)*(ai*
1331  - ((1.d0 + gam**2)*
1332  - (a*e**2*pg*qf*qi + f**2*pz*vi*x) -
1333  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*y*dcos(theta)) +
1334  - dconjg(vi)*
1335  - (ai*f**2*(1.d0 + gam**2)*pz*x -
1336  - (0.d0, 1.d0)*gam**2*v*
1337  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dcos(theta)))
1338  - + ai*((1.d0 + gam**2)*
1339  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(x) +
1340  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*vf*dconjg(y)*
1341  - dcos(theta) +
1342  - af*((0.d0, -1.d0)*(-1.d0 + gam**2)*
1343  - (b*e**2*pg*qf*qi +
1344  - f**2*pz*vi*(y - dconjg(y))) +
1345  - ai*f**2*gam**2*pz*v*(x + dconjg(x))*
1346  - dcos(theta))) +
1347  - dconjg(vi)*
1348  - (ai*f**2*(1.d0 + gam**2)*pz*vf*dconjg(x) +
1349  - (0.d0, 1.d0)*gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1350  - dconjg(y)*dcos(theta) +
1351  - af*((0.d0, -1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*
1352  - (y - dconjg(y)) +
1353  - gam**2*v*
1354  - (a*e**2*pg*qf*qi +
1355  - f**2*pz*vi*(x + dconjg(x)))*dcos(theta)))
1356  - ))*dsin(theta)
1357 
1358  rdm(2,4)= 4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1359  - (ai*f**2*gam**2*pz*v*
1360  - (y + (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b)) -
1361  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1362  - (e**2*pg*qf*qi*(a - dconjg(a)) +
1363  - f**2*pz*vi*
1364  - (x - vf*dconjg(a) + (0.d0, 1.d0)*af*dconjg(b)))*
1365  - dcos(theta)) +
1366  - f**2*dconjg(pz)*
1367  - (dconjg(vf)*(ai*
1368  - (gam**2*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) -
1369  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*x*dcos(theta)) +
1370  - dconjg(vi)*
1371  - (ai*f**2*gam**2*pz*v*y -
1372  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1373  - (a*e**2*pg*qf*qi + f**2*pz*vi*x)*dcos(theta)))
1374  - + ai*(gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1375  - dconjg(y) +
1376  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*vf*dconjg(x)*
1377  - dcos(theta) +
1378  - af*((0.d0, -1.d0)*gam**2*v*
1379  - (a*e**2*pg*qf*qi +
1380  - f**2*pz*vi*(x - dconjg(x))) +
1381  - ai*f**2*(-1.d0 + gam**2)*pz*(y + dconjg(y))*
1382  - dcos(theta))) +
1383  - dconjg(vi)*
1384  - ((0.d0, 1.d0)*((0.d0, -1.d0)*ai*f**2*gam**2*pz*v*vf*
1385  - dconjg(y) +
1386  - (-1.d0 + gam**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1387  - dconjg(x)*dcos(theta)) +
1388  - af*((0.d0, -1.d0)*ai*f**2*gam**2*pz*v*
1389  - (x - dconjg(x)) +
1390  - (-1.d0 + gam**2)*
1391  - (b*e**2*pg*qf*qi +
1392  - f**2*pz*vi*(y + dconjg(y)))*dcos(theta)))
1393  - ))*dsin(theta)
1394 
1395  rdm(4,2)= 4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1396  - (-(ai*f**2*gam**2*pz*v*
1397  - (y - (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b))) -
1398  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1399  - (e**2*pg*qf*qi*(a - dconjg(a)) +
1400  - f**2*pz*vi*
1401  - (x - vf*dconjg(a) - (0.d0, 1.d0)*af*dconjg(b)))*
1402  - dcos(theta)) -
1403  - f**2*dconjg(pz)*
1404  - (dconjg(vf)*(ai*
1405  - (gam**2*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) +
1406  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*x*dcos(theta)) +
1407  - dconjg(vi)*
1408  - (ai*f**2*gam**2*pz*v*y +
1409  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1410  - (a*e**2*pg*qf*qi + f**2*pz*vi*x)*dcos(theta)))
1411  - + ai*(gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1412  - dconjg(y) -
1413  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*vf*dconjg(x)*
1414  - dcos(theta) +
1415  - af*((0.d0, 1.d0)*gam**2*v*
1416  - (a*e**2*pg*qf*qi +
1417  - f**2*pz*vi*(x - dconjg(x))) +
1418  - ai*f**2*(-1.d0 + gam**2)*pz*(y + dconjg(y))*
1419  - dcos(theta))) +
1420  - dconjg(vi)*
1421  - ((0.d0, -1.d0)*((0.d0, 1.d0)*ai*f**2*gam**2*pz*v*vf*
1422  - dconjg(y) +
1423  - (-1.d0 + gam**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1424  - dconjg(x)*dcos(theta)) +
1425  - af*((0.d0, 1.d0)*ai*f**2*gam**2*pz*v*
1426  - (x - dconjg(x)) +
1427  - (-1.d0 + gam**2)*
1428  - (b*e**2*pg*qf*qi +
1429  - f**2*pz*vi*(y + dconjg(y)))*dcos(theta)))
1430  - ))*dsin(theta)
1431 
1432  rdm(3,4)= -2.d0*gam**4*m**4*(f**2*dconjg(pz)*
1433  - (-(af*v*(ai**2*f**2*pz*(x + dconjg(x)) +
1434  - dconjg(vi)*
1435  - (a*e**2*pg*qf*qi +
1436  - f**2*pz*vi*(x + dconjg(x))))*
1437  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1438  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1439  - (0.d0, 1.d0)*((0.d0, -4.d0)*ai*
1440  - (e**2*pg*qf*qi +
1441  - f**2*pz*vf*(vi + dconjg(vi)))*dconjg(x)*
1442  - dcos(theta) +
1443  - v*(ai**2*f**2*pz*vf +
1444  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))
1445  - *dconjg(y)*
1446  - (2.d0 + gam**2*(-1.d0 + v**2) +
1447  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1448  - dconjg(vf)*
1449  - ((0.d0, -1.d0)*ai*((0.d0, 4.d0)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1450  - dcos(theta) +
1451  - ai*f**2*pz*v*y*
1452  - (2.d0 + gam**2*(-1.d0 + v**2) +
1453  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1454  - dconjg(vi)*
1455  - (4.d0*ai*f**2*pz*x*dcos(theta) -
1456  - (0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y)*
1457  - (2.d0 + gam**2*(-1.d0 + v**2) +
1458  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))) +
1459  - e**2*qf*qi*dconjg(pg)*
1460  - (4.d0*ai*f**2*pz*(x + vf*dconjg(a))*dcos(theta) -
1461  - (0.d0, 1.d0)*v*(e**2*pg*qf*qi*(b - dconjg(b))*
1462  - (2.d0 + gam**2*(-1.d0 + v**2) +
1463  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1464  - f**2*pz*vi*
1465  - ((0.d0, -1.d0)*af*dconjg(a)*
1466  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1467  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1468  - (y - vf*dconjg(b))*
1469  - (2.d0 + gam**2*(-1.d0 + v**2) +
1470  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))))
1471 
1472  rdm(4,3)= -2.d0*gam**4*m**4*(f**2*dconjg(pz)*
1473  - (-(af*v*(ai**2*f**2*pz*(x + dconjg(x)) +
1474  - dconjg(vi)*
1475  - (a*e**2*pg*qf*qi +
1476  - f**2*pz*vi*(x + dconjg(x))))*
1477  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1478  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) -
1479  - (0.d0, 1.d0)*((0.d0, 4.d0)*ai*
1480  - (e**2*pg*qf*qi +
1481  - f**2*pz*vf*(vi + dconjg(vi)))*dconjg(x)*
1482  - dcos(theta) +
1483  - v*(ai**2*f**2*pz*vf +
1484  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))
1485  - *dconjg(y)*
1486  - (2.d0 + gam**2*(-1.d0 + v**2) +
1487  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1488  - dconjg(vf)*
1489  - ((0.d0, 1.d0)*ai*((0.d0, -4.d0)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1490  - dcos(theta) +
1491  - ai*f**2*pz*v*y*
1492  - (2.d0 + gam**2*(-1.d0 + v**2) +
1493  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1494  - dconjg(vi)*
1495  - (4.d0*ai*f**2*pz*x*dcos(theta) +
1496  - (0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y)*
1497  - (2.d0 + gam**2*(-1.d0 + v**2) +
1498  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))) +
1499  - e**2*qf*qi*dconjg(pg)*
1500  - (4.d0*ai*f**2*pz*(x + vf*dconjg(a))*dcos(theta) +
1501  - (0.d0, 1.d0)*v*(e**2*pg*qf*qi*(b - dconjg(b))*
1502  - (2.d0 + gam**2*(-1.d0 + v**2) +
1503  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1504  - f**2*pz*vi*
1505  - ((0.d0, 1.d0)*af*dconjg(a)*
1506  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1507  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1508  - (y - vf*dconjg(b))*
1509  - (2.d0 + gam**2*(-1.d0 + v**2) +
1510  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))))
1511 
1512  rdm(4,4)= 8.d0*gam**4*m**4*(e**4*pg*qf**2*qi**2*(a + dconjg(a))*
1513  - dconjg(pg) +
1514  - e**2*f**2*qf*qi*(pz*vi*x*dconjg(pg) +
1515  - a*pg*dconjg(pz)*dconjg(vf)*dconjg(vi) +
1516  - pg*dconjg(pz)*dconjg(vi)*dconjg(x) +
1517  - a*af*ai*pg*v*dconjg(pz)*dcos(theta) +
1518  - pz*dconjg(a)*dconjg(pg)*
1519  - (vf*vi + af*ai*v*dcos(theta))) +
1520  - f**4*pz*dconjg(pz)*
1521  - (ai*(ai*x*dconjg(vf) + ai*vf*dconjg(x) +
1522  - af*v*vi*x*dcos(theta) +
1523  - af*v*vi*dconjg(x)*dcos(theta)) +
1524  - dconjg(vi)*
1525  - (vi*(x*dconjg(vf) + vf*dconjg(x)) +
1526  - af*ai*v*(x + dconjg(x))*dcos(theta))))
1527 
1528 c-------------------------------------------------------------
1529 c TOTAL SPIN-CORRELATION MATRIX: STANDARD MODEL + DIPOLE MOMENTS
1530 c-------------------------------------------------------------
1531  do i=1,4
1532  do j=1,4
1533  r(i,j)= rsm(i,j) + rdm(i,j)
1534  enddo
1535  enddo
1536 
1537  return
1538  end
Generated on Tue Mar 12 2024 22:57:50 for C++InterfacetoTauola by   doxygen 1.8.7