C++ Interface to Tauola
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 dipolgamma(iqed, E, theta, A, B, R)
435 !
436 ! This routine performs adjustment of frame orientation between conventions used by Korchin and the one of TauSpinner
437 ! which is inherited from conventions of KORALB-KKMC: that consist of rotations
438 ! Also index shift is performed from xyzt (1234) to txyz [in c++ that will be seen as (0123)]
439 
440  implicit none
441  real*8 amz0,gamz0,sin2w0,alphaqed,gswr(7),gswi(7)
442  integer iqed,channel,k,j,k0,j0
443  real*8 e,theta,a,b,xnor,thet, buf
444  real*8 r(1:4, 1:4),r0(1:4, 1:4),sign(4)
445 ! ReA0=0
446 ! ImA0=0
447 ! ReB=0
448 ! ImB=0
449 ! ReX=0
450 ! ImX=0
451 ! ReY=0
452 ! ImY=0
453  sign(1)=-1.0
454  sign(2)= 1.0
455  sign(3)=-1.0
456  sign(4)= 1.0
457  thet=acos(cos(theta))
458 
459  call dipolgammarij(iqed, e, theta, a, b, r0)
460 
461 
462 
463 ! rotation from q ~q --> tau- tau+ to ~q q --> tau+ tau- frames.
464 ! WARNING: we normalize R0, but x-section is set into R0(4,4).
465  xnor=r0(4,4)
466  do k=1,4
467  do j=1,4
468  r0(k,j)=r0(k,j)/xnor*sign(k)*sign(j)
469  enddo
470  enddo
471  r0(4,4)=xnor
472 
473  ! rotation in second index by angle pi/2 around z-axis
474  do k=1,4
475  do j=1,4
476  if(j.eq.1) then
477  buf=r0(k,2)
478  r0(k,2)=r0(k,j)
479  r0(k,1)=-buf
480  endif
481 
482  enddo
483  enddo
484 
485 ! rotation in first index by angle pi/2 around z-axis
486  do j=1,4
487  do k=1,4
488  if(k.eq.1) then
489  buf=r0(2,j)
490  r0(2,j)=r0(k,j)
491  r0(1,j)=-buf
492  endif
493 
494  enddo
495  enddo
496 
497 
498 ! shift order of entries from x,y,z,t to t,x,y,z (necessary for C++ conventions).
499  do k0=1,4
500  do j0=1,4
501  k=k0+1
502  j=j0+1
503  if(j.eq.5) j=1
504  if(k.eq.5) k=1
505  r(k,j)=r0(k0,j0)
506  enddo
507  enddo
508  ! if(R0(3,3)/R0(4,4).gt.1d0) then
509  ! write(*,*) "z fortranu thet=", thet," energy= ",energy, " Rzz=", R0(3,3)/R0(4,4)," ", R0(3,3)," ",R0(4,4)
510 
511  ! write(*,*) "iqed,Energy, thet, channel ",iqed,Energy, thet, channel
512  ! write(*,*) "Amz0,Gamz0,sin2W0,alphaQED ", Amz0,Gamz0,sin2W0,alphaQED
513  ! write(*,*) "ReA0, ImA0, ReB, ImB ", ReA0, ImA0, ReB, ImB
514  ! write(*,*) "ReX, ImX, ReY, ImY ", ReX, ImX, ReY, ImY
515  ! write(*,*) "GSWr ", GSWr
516  ! write(*,*) "GSWi ", GSWi
517  ! endif
518  ! IF(R(4,4)/R(1,1).gt.1.0) then
519  ! write(*,*) "==============",R(4,4)/R(1,1)
520  ! write(*,*) 'energuy,theta,channel= ',Energy,' ',cos(theta),' ',channel
521  ! do k=1,4
522  ! 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)
523  ! enddo
524  ! endif
525 
526 
527  7 format(1x,i4, 4(2x,f12.6))
528 
529  end
530  Subroutine dipolgammarijold (iqed, E, theta, A, B, R)
531 
532 c For gamma gamma -> tau- tau+ reaction
533 c
534 c Calculates spin-correlation coefficients R(i,j) as functions of beam
535 c energy E = sqrt(s)/2 (in GeV) and scattaring angle theta.
536 c V is velocity of tau lepton, gam is Lorentz factor, alpha is
537 c fine-structure constant.
538 c Anomalous magnetic moment A1 = ASM + A and electric dipole moment B are
539 c real constants.
540 c Parameters A, B describe 'New Physics'.
541 c ASM is anomalous magnetic moment of tau in Standard Model (SM):
542 c S.Eidelman and M.Passera. Mod.Phys.Lett. A22, 159-179, 2007.
543 c
544 c Order of coefficients:
545 c i = 1,2,3 correspond to S_x, S_y, S_z for tau-,
546 c j = 1,2,3 correspond to S'_x, S'_y, S'_z for tau+
547 c i,j = 4 (equivalent to tt) corresponds to 1 (no spin dependence).
548 
549  Implicit none
550  integer iqed
551 C integer i, j ! not used here
552  real*8 e, theta, a, b
553  real*8 r(1:4, 1:4)
554  real*8 asm, a1, gam, e4, v,d2
555  real*8 fnorm
556  real*8 pi/3.141592653589793238d0/,alpha/7.2973525693d-3/
557  real*8 mtau ! mass of tau in GeV
558  COMMON / t_beampm / ene ,amin,amfin,ide,idf
559  integer IDE,IDF
560  real*8 ene ,amin,amfin
561  integer k,l
562  mtau=amfin
563  v = dsqrt(1.d0 -(mtau/e)**2) ! tau velocity
564  e4 = (4.0d0 *pi *alpha)**2 ! (electric charge)^4
565  gam = e/mtau ! Lorentz factor for tau
566  fnorm= 2*(2*pi**2*alpha**2) ! for normalization as in tauspinner.
567 c Cotribution to anomalous magnetic moment in Standard Model:
568  asm = 1.17721d-3
569  asm = asm *iqed ! switch for dipole SM part
570  a1 = asm + a ! Standard Model + New Physics
571 
572  d2 = (v**2 *dcos(theta)**2 -1.d0)**2 ! factor in denominator
573 
574  r(1,1) = e4 *(-11.d0 *v**4 +28.d0 *a1 *v**2 +
575  $ 4.d0 *(v**2 -2.d0) *dcos(2.d0*theta) *v**2 -
576  $ (v**2 -4.d0 *a1 -2.d0) *dcos(4.d0*theta) *v**2 +
577  $ 22.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)
578 
579  r(1,2) = e4 *v *b *(dcos(4.d0*theta) *v**2 +15.d0 *v**2 +
580  $ 4.d0 *cos(2.d0*theta) -20.d0) /(4.d0 *d2)
581 
582  r(1,3) = e4 *v**2 *gam *((a1 -1.d0) *v**2 +
583  $ (v**2 +a1 *(v**2 -2.d0) -1.d0) *dcos(2.d0*theta) +1.d0) *
584  $ dsin(2.d0*theta) /(2.d0 *d2)
585 
586  r(1,4) = 0.d0
587 
588  r(2,1) = -r(1,2)
589 
590  r(2,2) = e4 *(-dcos(4.d0*theta) *v**4 -11.d0 *v**4 +
591  $ 16.d0 *a1 *v**2 +4.d0 *(v**2 +4.d0 *a1) *dcos(2.d0*theta) *v**2 +
592  $ 16.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)
593 
594  r(2,3) = e4 *v *gam *b *(dcos(2.d0*theta)*v**2 -3.d0*v**2 +2.d0) *
595  $ dsin(2.d0*theta) /(2.d0 *d2)
596 
597  r(2,4) = 0.d0
598 
599  r(3,1) = r(1,3)
600 
601  r(3,2) = -r(2,3)
602 
603  r(3,3) = e4 *(-4.d0 *dcos(2.d0*theta) *v**4 +11.d0 *v**4 +
604  $ 36.d0 *a1 *v**2 +(v**2 -4.d0 *a1 -2.d0) *dcos(4.d0*theta) *v**2 +
605  $ 2.d0 *v**2 -32.d0 *a1 -8.d0) /(8.d0 *d2)
606 
607  r(3,4) = 0.d0
608 
609  r(4,1) = 0.d0
610 
611  r(4,2) = 0.d0
612 
613  r(4,3) = 0.d0
614 
615  r(4,4) = e4 *(-dcos(4.d0*theta) *v**4 -11.d0 *v**4 -16.d0 *a1 *v**2 +
616  $ 4.d0 *(v**2 -4.d0 *a1 -2.d0) *dcos(2.d0*theta) *v**2 +
617  $ 8.d0 * v**2 +32.d0 *a1 +8.d0) /(8.d0 *d2)
618 
619 ! to normalize as in tauspinner
620  ! FNORM=1.0
621  do k=1,4
622  do l=1,4
623  r(k,l)=r(k,l)/fnorm
624  enddo
625  enddo
626 
627  return
628  end
629 
630  Subroutine dipolgammarij (iqed, E, theta, A, B, R)
631 
632 c For gamma gamma -> tau- tau+ reaction
633 c
634 c 14 March 2025: updated with exact calculation, in which
635 c all orders in dipole moments are included up to the 4th order.
636 c earlier it was linearized part only. For discovery evaluation versions
637 c should be equivalent.
638 c
639 c Calculates spin-correlation coefficients R(i,j) as functions of beam
640 c energy E = sqrt(s)/2 (in GeV) and scattaring angle theta.
641 c V is velocity of tau lepton, gam is Lorentz factor, alpha is
642 c fine-structure constant.
643 c Anomalous magnetic moment A1 = ASM + A and electric dipole moment B are
644 c real constants.
645 c Parameters A, B describe 'New Physics'.
646 c ASM is anomalous magnetic moment of tau in Standard Model (SM):
647 c S.Eidelman and M.Passera. Mod.Phys.Lett. A22, 159-179, 2007.
648 c
649 c Order of coefficients:
650 c i = 1,2,3 correspond to S_x, S_y, S_z for tau-,
651 c j = 1,2,3 correspond to S'_x, S'_y, S'_z for tau+
652 c i,j = 4 (equivalent to tt) corresponds to 1 (no spin dependence).
653 
654  Implicit none
655  integer iqed
656  integer i, j ! not used here
657  real*8 e, theta, a, b
658  real*8 r(1:4, 1:4)
659  real*8 asm, a1, gam, e4, v, d2
660  real*8 fnorm
661  real*8 pi/3.141592653589793238d0/,alpha/7.2973525693d-3/
662  real*8 mtau ! mass of tau in GeV
663  COMMON / t_beampm / ene ,amin,amfin,ide,idf
664  integer IDE,IDF
665  real*8 ene ,amin,amfin
666  integer k,l
667 
668  mtau=amfin
669  v = dsqrt(1.d0 -(mtau/e)**2) ! tau velocity
670  e4 = (4.0d0 *pi *alpha)**2 ! (electric charge)^4
671  gam = e/mtau ! Lorentz factor for tau
672  fnorm= 2*(2*pi**2*alpha**2) ! for normalization as in tauspinner.
673 
674 c Cotribution to anomalous magnetic moment in Standard Model:
675  asm = 1.17721d-3
676  asm = asm *iqed ! switch for dipole SM part
677  a1 = asm + a ! Standard Model + New Physics
678 
679  d2 = (v**2 *dcos(theta)**2 -1.d0)**2 ! factor in denominator NOT used below
680 
681 
682  r(1,1) =
683  $ (e4*(-8.d0 + 8.d0*gam**2 - 16.d0*b**2*gam**2 - 4.d0*gam**4 -
684  $ 16.d0*a1*gam**4 - 20.d0*a1**2*gam**4 - 8.d0*a1**3*gam**4 -
685  $ 2.d0*a1**4*gam**4 + 12.d0*b**2*gam**4 - 8.d0*a1*b**2*gam**4 -
686  $ 4.d0*a1**2*b**2*gam**4 - 2.d0*b**4*gam**4 +
687  $ 2.d0*a1**4*gam**6 - 8.d0*b**2*gam**6 +
688  $ 4.d0*a1**2*b**2*gam**6 + 2.d0*b**4*gam**6 - a1**4*gam**8 -
689  $ 2.d0*a1**2*b**2*gam**8 - b**4*gam**8 +
690  $ gam**2*(-8.d0 + 2.d0*
691  $ (4.d0 + 4.d0*a1**3 + a1**4 - 6.d0*b**2 + b**4 +
692  $ 4.d0*a1*(2.d0 + b**2) + 2.d0*a1**2*(5.d0 + b**2))*gam**2 -
693  $ 4.d0*(a1**4 - 4.d0*b**2 + 2.d0*a1**2*b**2 + b**4)*gam**4 +
694  $ 3.d0*(a1**2 + b**2)**2*gam**6)*v**2*dcos(theta)**2 +
695  $ (a1**2 + b**2)**2*gam**8*v**6*dcos(theta)**6 -
696  $ gam**4*v**4*dcos(theta)**4*
697  $ (4.d0 + 8.d0*b**2*gam**2 - b**4*gam**2 + 2.d0*b**4*gam**4 +
698  $ a1**4*gam**2*(-1.d0 + 2.d0*gam**2) +
699  $ 2.d0*a1**2*b**2*gam**2*(-1.d0 + 2.d0*gam**2) +
700  $ (a1**2 + b**2)**2*gam**2*(-1.d0 + gam**2)*dcos(2.d0*theta)) +
701  $ 2.d0*gam**2*(-1.d0 + gam**2)*
702  $ (2.d0 + 2.d0*a1 + a1**2*gam**2 + b**2*gam**2)**2*
703  $ dsin(theta)**2 - 2.d0*gam**4*v**2*dsin(2.d0*theta)**2 -
704  $ 4.d0*a1*gam**4*v**2*dsin(2.d0*theta)**2 -
705  $ 4.d0*a1**2*gam**4*v**2*dsin(2.d0*theta)**2 -
706  $ 2.d0*a1**3*gam**4*v**2*dsin(2.d0*theta)**2 -
707  $ a1**4*gam**4*v**2*dsin(2.d0*theta)**2 -
708  $ 2.d0*b**2*gam**4*v**2*dsin(2.d0*theta)**2 -
709  $ 2.d0*a1*b**2*gam**4*v**2*dsin(2.d0*theta)**2 -
710  $ 2.d0*a1**2*b**2*gam**4*v**2*dsin(2.d0*theta)**2 -
711  $ b**4*gam**4*v**2*dsin(2.d0*theta)**2 -
712  $ 2.d0*a1**2*gam**6*v**2*dsin(2.d0*theta)**2 -
713  $ 2.d0*a1**3*gam**6*v**2*dsin(2.d0*theta)**2 +
714  $ a1**4*gam**6*v**2*dsin(2.d0*theta)**2 -
715  $ 2.d0*b**2*gam**6*v**2*dsin(2.d0*theta)**2 -
716  $ 2.d0*a1*b**2*gam**6*v**2*dsin(2.d0*theta)**2 +
717  $ 2.d0*a1**2*b**2*gam**6*v**2*dsin(2.d0*theta)**2 +
718  $ b**4*gam**6*v**2*dsin(2.d0*theta)**2 -
719  $ a1**4*gam**8*v**2*dsin(2.d0*theta)**2 -
720  $ 2.d0*a1**2*b**2*gam**8*v**2*dsin(2.d0*theta)**2 -
721  $ b**4*gam**8*v**2*dsin(2.d0*theta)**2))/
722  $ (4.d0*gam**4*(-1.d0 + v**2*dcos(theta)**2)**2)
723 
724 
725  r(1,2) =
726  $ (b*e4*v*(-32.d0 + 32.d0*a1 - 48.d0*gam**2 - 192.d0*a1*gam**2 -
727  $ 52.d0*a1**2*gam**2 - 20.d0*b**2*gam**2 + 112.d0*a1*gam**4 +
728  $ 44.d0*a1**2*gam**4 + 12.d0*b**2*gam**4 + 28.d0*gam**2*v**2 +
729  $ 120.d0*a1*gam**2*v**2 + 27.d0*a1**2*gam**2*v**2 +
730  $ 3.d0*b**2*gam**2*v**2 - 172.d0*a1*gam**4*v**2 -
731  $ 65.d0*a1**2*gam**4*v**2 - 9.d0*b**2*gam**4*v**2 +
732  $ 72.d0*a1*gam**4*v**4 + 27.d0*a1**2*gam**4*v**4 +
733  $ 3.d0*b**2*gam**4*v**4 +
734  $ 4.d0*(8.d0 + b**2*gam**4*(-1.d0 - 2.d0*v**2 + v**4) +
735  $ a1**2*gam**2*
736  $ (-5.d0 + 7.d0*gam**2 + (9.d0 - 18.d0*gam**2)*v**2 +
737  $ 9.d0*gam**2*v**4) +
738  $ gam**2*(-4.d0 + 8.d0*v**2 + b**2*(3.d0 + v**2)) +
739  $ 4.d0*a1*(2.d0 + gam**2*(-6.d0 + 9.d0*v**2) +
740  $ gam**4*(5.d0 - 12.d0*v**2 + 6.d0*v**4)))*dcos(2.d0*theta) +
741  $ gam**2*v**2*(4.d0 + b**2*(1.d0 + gam**2*(1.d0 + v**2)) +
742  $ 4.d0*a1*(6.d0 + gam**2*(-5.d0 + 6.d0*v**2)) +
743  $ a1**2*(9.d0 + gam**2*(-7.d0 + 9.d0*v**2)))*dcos(4.d0*theta)))/
744  $ (16.d0*gam**2*(-1.d0 + v**2*dcos(theta)**2)**2)
745 
746 
747  r(1,3) =
748  $ -(e4*dcos(theta)*(4.d0 - 4.d0*gam**2 + 4.d0*b**2*gam**2 -
749  $ 4.d0*b**2*gam**4 + b**4*gam**4 - b**4*gam**6 - 4.d0*v**2 -
750  $ 4.d0*b**2*v**2 + 4.d0*gam**2*v**2 - 6.d0*b**2*gam**2*v**2 -
751  $ 2.d0*b**4*gam**2*v**2 + 4.d0*b**2*gam**4*v**2 +
752  $ b**4*gam**6*v**2 +
753  $ 4.d0*a1**3*gam**2*(1.d0 - gam**2 + (-2.d0 + gam**2)*v**2) +
754  $ a1**4*(gam**4 - gam**6 +
755  $ gam**2*(-2.d0 + gam**4)*v**2) +
756  $ 2.d0*a1**2*(-((-1.d0 + gam**2)*
757  $ (2.d0 + 2.d0*gam**2 + b**2*gam**4)) +
758  $ gam**2*(-3.d0 + 2.d0*gam**2 + b**2*(-2.d0 + gam**4))*v**2) +
759  $ 4.d0*a1*(2.d0 + b**2*gam**4*(-1.d0 + v**2) +
760  $ gam**2*(-2.d0 + v**2 + b**2*(1.d0 - 2.d0*v**2))) -
761  $ 2.d0*gam**2*v**2*
762  $ (2.d0*(-1.d0 + v**2) +
763  $ b**2*(-2.d0*(1.d0 + gam**2) + (1.d0 + 2.d0*gam**2)*v**2) +
764  $ 2.d0*a1**3*(-1.d0 + gam**2*(-1.d0 + v**2)) +
765  $ a1**4*(-1.d0 + gam**2 + gam**4*(-1.d0 + v**2)) +
766  $ b**4*(-1.d0 + gam**2 + gam**4*(-1.d0 + v**2)) +
767  $ 2.d0*a1*(-2.d0 + v**2 +
768  $ b**2*(-1.d0 + gam**2*(-1.d0 + v**2))) +
769  $ a1**2*(-4.d0 - 2.d0*gam**2 + v**2 + 2.d0*gam**2*v**2 +
770  $ 2.d0*b**2*(-1.d0 + gam**2 + gam**4*(-1.d0 + v**2))))*
771  $ dcos(theta)**2 +
772  $ (a1**2 + b**2)**2*gam**4*v**4*
773  $ (1.d0 + gam**2*(-1.d0 + v**2))*dcos(theta)**4)*dsin(theta))/
774  $ (2.d0*gam*(-1.d0 + v**2*dcos(theta)**2)**2)
775 
776 
777  r(1,4) = 0.d0
778 
779  r(2,1) = -r(1,2)
780 
781 
782  r(2,2) =
783  $ (e4*(-8.d0 + (8.d0 - 16.d0*b**2)*gam**2 -
784  $ 2.d0*(2.d0 + 4.d0*a1**3 + a1**4 - 6.d0*b**2 + b**4 +
785  $ 4.d0*a1*(2.d0 + b**2) + 2.d0*a1**2*(5.d0 + b**2))*gam**4 +
786  $ 2.d0*(a1**4 - 4.d0*b**2 + 2.d0*a1**2*b**2 + b**4)*gam**6 -
787  $ (a1**2 + b**2)**2*gam**8 +
788  $ gam**2*(-8.d0 + 2.d0*
789  $ (4.d0 + 4.d0*a1**3 + a1**4 - 6.d0*b**2 + b**4 +
790  $ 4.d0*a1*(2.d0 + b**2) + 2.d0*a1**2*(5.d0 + b**2))*gam**2 -
791  $ 4.d0*(a1**4 - 4.d0*b**2 + 2.d0*a1**2*b**2 + b**4)*gam**4 +
792  $ 3.d0*(a1**2 + b**2)**2*gam**6)*v**2*dcos(theta)**2 -
793  $ gam**4*(4.d0 + 8.d0*b**2*gam**2 +
794  $ a1**4*gam**2*(-2.d0 + 3.d0*gam**2) +
795  $ 2.d0*a1**2*b**2*gam**2*(-2.d0 + 3.d0*gam**2) +
796  $ b**4*gam**2*(-2.d0 + 3.d0*gam**2))*v**4*dcos(theta)**4 +
797  $ (a1**2 + b**2)**2*gam**8*v**6*dcos(theta)**6))/
798  $ (4.d0*gam**4*(-1.d0 + v**2*dcos(theta)**2)**2)
799 
800 
801  r(2,3) =
802  $ -(b*e4*v*(-16.d0 + 8.d0*gam**2 - 6.d0*b**2*gam**2 +
803  $ 2.d0*b**2*gam**4 -
804  $ 4.d0*gam**2*v**2 + 7.d0*b**2*gam**2*v**2 -
805  $ 3.d0*b**2*gam**4*v**2 + b**2*gam**4*v**4 -
806  $ a1**2*gam**2*(-10.d0 + 14.d0*gam**2 + v**2 -
807  $ 21.d0*gam**2*v**2 + 7.d0*gam**2*v**4) -
808  $ 4.d0*a1*(4.d0 + 6.d0*gam**2*(-2.d0 + v**2) +
809  $ 5.d0*gam**4*(2.d0 - 3.d0*v**2 + v**4)) -
810  $ gam**2*v**2*(4.d0 + b**2*(1.d0 + gam**2 - gam**2*v**2) +
811  $ 4.d0*a1*(6.d0 + 5.d0*gam**2*(-1.d0 + v**2)) +
812  $ a1**2*(9.d0 + 7.d0*gam**2*(-1.d0 + v**2)))*dcos(2.d0*theta))*
813  $ dsin(2.d0*theta))/(8.d0*gam*(-1.d0 + v**2*dcos(theta)**2)**2)
814 
815 
816  r(2,4) = 0.d0
817 
818  r(3,1) = r(1,3)
819 
820  r(3,2) = -r(2,3)
821 
822 
823  r(3,3) =
824  $ -(e4*(-8.d0 - 8.d0*(-3.d0 + 2.d0*b**2)*gam**2 -
825  $ 2.d0*(a1**2 + b**2)**2*gam**10*(-1.d0 + v**2) -
826  $ 2.d0*gam**4*(10.d0 + 4.d0*a1**3 + a1**4 + b**4 - 4.d0*v**2 +
827  $ 4.d0*a1*(2.d0 + b**2 + 2.d0*v**2) +
828  $ 2.d0*a1**2*(5.d0 + b**2 + 2.d0*v**2) +
829  $ b**2*(-22.d0 + 8.d0*v**2)) +
830  $ 2.d0*gam**6*(4.d0 + 3.d0*a1**4 + 3.d0*b**4 - 4.d0*v**2 +
831  $ 2.d0*a1**2*(10.d0 + 3.d0*b**2 - 8.d0*v**2) -
832  $ 4.d0*a1**3*(-2.d0 + v**2) + 4.d0*b**2*(-4.d0 + 3.d0*v**2) -
833  $ 4.d0*a1*(b**2*(-2.d0 + v**2) + 4.d0*(-1.d0 + v**2))) +
834  $ gam**8*(a1**4*(-5.d0 + 2.d0*v**2) +
835  $ 2.d0*a1**2*b**2*(-5.d0 + 2.d0*v**2) +
836  $ b**2*(-16.d0*(-1.d0 + v**2) + b**2*(-5.d0 + 2.d0*v**2))) +
837  $ gam**2*(-2.d0*gam**2*(-1.d0 + gam**2)*
838  $ (2.d0 + 2.d0*a1 + a1**2*gam**2 + b**2*gam**2)**2 +
839  $ (-8.d0 + 2.d0*(4.d0 + 4.d0*a1**3 + a1**4 -
840  $ 14.d0*b**2 + b**4 +
841  $ 4.d0*a1*(2.d0 + b**2) + 2.d0*a1**2*(5.d0 + b**2))*gam**2 -
842  $ 16.d0*(a1 + 3.d0*a1**3 + a1**4 + 3.d0*a1*b**2 +
843  $ b**2*(-1.d0 + b**2) + 2*a1**2*(2.d0 + b**2))*
844  $ gam**4 +
845  $ (16.d0*a1**3 + 11.d0*a1**4 + 16.d0*a1*b**2 +
846  $ b**2*(-16.d0 + 11.d0*b**2) + 2.d0*a1**2*(8.d0 + 11.d0*b**2))*
847  $ gam**6 - 2.d0*(a1**2 + b**2)**2*gam**8)*v**2 +
848  $ 4.d0*gam**4*
849  $ (2.d0 - 2.d0*a1**3*(-3.d0 + gam**2) +
850  $ b**2*(-2.d0 + 6.d0*gam**2) +
851  $ a1**4*(1.d0 - gam**2 + gam**4) +
852  $ b**4*(1.d0 - gam**2 + gam**4) +
853  $ a1*(8.d0 - 2.d0*b**2*(-3.d0 + gam**2)) +
854  $ 2.d0*a1**2*
855  $ (6.d0 - gam**2 + b**2*(1.d0 - gam**2 + gam**4)))*
856  $ v**4)*dcos(theta)**2 -
857  $ gam**4*v**2*(-4.d0*gam**2*
858  $ (2.d0 + 2.d0*a1**3*(1.d0 + gam**2) +
859  $ 2.d0*b**2*(1.d0 + gam**2) +
860  $ a1**4*(1.d0 - gam**2 + gam**4) +
861  $ b**4*(1.d0 - gam**2 + gam**4) +
862  $ 2.d0*a1*(2.d0 + b**2*(1.d0 + gam**2)) +
863  $ 2.d0*a1**2*
864  $ (2.d0 + gam**2 + b**2*(1.d0 - gam**2 + gam**4))) +
865  $ (4.d0 - 2.d0*(-4.d0 - 8.d0*a1 + a1**4 - 8.d0*b**2 + b**4 +
866  $ 2.d0*a1**2*(-2.d0 + b**2))*gam**2 +
867  $ (16.d0*a1**3 + 7.d0*a1**4 + 16.d0*a1*b**2 + 7.d0*b**4 +
868  $ 2.d0*a1**2*(8.d0 + 7.d0*b**2))*gam**4 +
869  $ 2.d0*(a1**2 + b**2)**2*gam**6)*v**2 +
870  $ 2.d0*gam**4*(-4.d0*a1**3 - 4.d0*a1*b**2 +
871  $ a1**4*(-1.d0 + gam**2) +
872  $ 2.d0*a1**2*(-2.d0 + b**2*(-1.d0 + gam**2)) +
873  $ b**2*(4.d0 + b**2*(-1.d0 + gam**2)))*v**4)*
874  $ dcos(theta)**4 +
875  $ (a1**2 + b**2)**2*gam**8*v**4*
876  $ (2.d0 - 2.d0*gam**2 +
877  $ (1.d0 + 2.d0*gam**2)*v**2)*dcos(theta)**6))/
878  $ (4.d0*gam**4*(-1.d0 + v**2*dcos(theta)**2)**2)
879 
880 
881  r(3,4) = 0.d0
882 
883  r(4,1) = 0.d0
884 
885  r(4,2) = 0.d0
886 
887  r(4,3) = 0.d0
888 
889 
890  r(4,4) =
891  $ -(e4*(8.d0 - 8.d0*gam**2 +
892  $ 2.d0*(-2.d0 + 4.d0*a1**3 + a1**4 + 2.d0*b**2 + b**4 +
893  $ 4.d0*a1*(-2.d0 + b**2) + 2.d0*a1**2*(-1.d0 + b**2))*gam**4 -
894  $ 2.d0*(8.d0*a1**3 + a1**4 + 8.d0*a1*b**2 +
895  $ 2.d0*a1**2*(4.d0 + b**2) +
896  $ b**2*(8.d0 + b**2))*gam**6 -
897  $ (a1**2 + b**2)**2*gam**8 +
898  $ gam**2*(8.d0 + 8.d0*a1**3*gam**2*(-1.d0 + 4.d0*gam**2) +
899  $ 4.d0*b**2*gam**2*(-1.d0 + 8.d0*gam**2) +
900  $ a1**4*gam**2*(-2.d0 + 4.d0*gam**2 + 3.d0*gam**4) +
901  $ b**4*gam**2*(-2.d0 + 4.d0*gam**2 + 3.d0*gam**4) +
902  $ 8.d0*a1*gam**2*(2.d0 + b**2*(-1.d0 + 4.d0*gam**2)) +
903  $ 2.d0*a1**2*gam**2*
904  $ (2.d0 + 16.d0*gam**2 + b**2*(-2.d0 + 4.d0*gam**2 +
905  $ 3.d0*gam**4)))*
906  $ v**2*dcos(theta)**2 -
907  $ gam**4*(-4.d0 + 16.d0*a1**3*gam**2 + 16.d0*b**2*gam**2 +
908  $ 16.d0*a1*b**2*gam**2 + a1**4*gam**2*(2.d0 + 3.d0*gam**2) +
909  $ b**4*gam**2*(2.d0 + 3.d0*gam**2) +
910  $ 2.d0*a1**2*gam**2*(8.d0 + b**2*(2.d0 + 3.d0*gam**2)))*v**4*
911  $ dcos(theta)**4 +
912  $ (a1**2 + b**2)**2*gam**8*v**6*dcos(theta)**6))/
913  $ (4.d0*gam**4*(-1.d0 + v**2*dcos(theta)**2)**2)
914 
915 ! to normalize as in tauspinner
916  ! FNORM=1.0
917  do k=1,4
918  do l=1,4
919  r(k,l)=r(k,l)/fnorm
920  enddo
921  enddo
922 
923 
924  return
925  end
926 
927 
928  Subroutine dipolqq(iqed,Energy,theta,channel,Amz0,Gamz0,sin2W0,alphaQED,ReA0,ImA0, ReB,ImB, ReX,ImX,ReY,ImY,GSWr,GSWi,R)
929 !
930 ! This routine performs adjustment of frame orientation between conventions used by Korchin and the one of TauSpinner
931 ! which is inherited from conventions of KORALB-KKMC: that consist of rotations
932 ! Also index shift is performed from xyzt (1234) to txyz [in c++ that will be seen as (0123)]
933 
934  implicit none
935  real*8 amz0,gamz0,sin2w0,alphaqed,gswr(7),gswi(7)
936  integer iqed,channel,k,j,k0,j0
937  real*8 energy,theta,rea,ima,rea0,ima0,reb,imb,rex,imx,rey,imy,xnor,thet, buf
938  real*8 r(1:4, 1:4),r0(1:4, 1:4),sign(4)
939 ! ReA0=0
940 ! ImA0=0
941 ! ReB=0
942 ! ImB=0
943 ! ReX=0
944 ! ImX=0
945 ! ReY=0
946 ! ImY=0
947  sign(1)=-1.0
948  sign(2)= 1.0
949  sign(3)=-1.0
950  sign(4)= 1.0
951  thet=acos(cos(theta))
952 
953  call dipolqqrijradcor (iqed,energy, thet,amz0,gamz0,sin2w0,alphaqed,
954  # gswr,gswi, rea0, ima0, reb, imb, rex, imx, rey, imy, r0, channel)
955 
956 
957 
958 ! rotation from q ~q --> tau- tau+ to ~q q --> tau+ tau- frames.
959 ! WARNING: we normalize R0, but x-section is set into R0(4,4).
960  xnor=r0(4,4)
961  do k=1,4
962  do j=1,4
963  r0(k,j)=r0(k,j)/xnor*sign(k)*sign(j)
964  enddo
965  enddo
966  r0(4,4)=xnor
967 
968  ! rotation in second index by angle pi/2 around z-axis
969  do k=1,4
970  do j=1,4
971  if(j.eq.1) then
972  buf=r0(k,2)
973  r0(k,2)=r0(k,j)
974  r0(k,1)=-buf
975  endif
976 
977  enddo
978  enddo
979 
980 ! rotation in first index by angle pi/2 around z-axis
981  do j=1,4
982  do k=1,4
983  if(k.eq.1) then
984  buf=r0(2,j)
985  r0(2,j)=r0(k,j)
986  r0(1,j)=-buf
987  endif
988 
989  enddo
990  enddo
991 
992 
993 ! shift order of entries from x,y,z,t to t,x,y,z (necessary for C++ conventions).
994  do k0=1,4
995  do j0=1,4
996  k=k0+1
997  j=j0+1
998  if(j.eq.5) j=1
999  if(k.eq.5) k=1
1000  r(k,j)=r0(k0,j0)
1001  enddo
1002  enddo
1003  ! if(R0(3,3)/R0(4,4).gt.1d0) then
1004  ! write(*,*) "z fortranu thet=", thet," energy= ",energy, " Rzz=", R0(3,3)/R0(4,4)," ", R0(3,3)," ",R0(4,4)
1005 
1006  ! write(*,*) "iqed,Energy, thet, channel ",iqed,Energy, thet, channel
1007  ! write(*,*) "Amz0,Gamz0,sin2W0,alphaQED ", Amz0,Gamz0,sin2W0,alphaQED
1008  ! write(*,*) "ReA0, ImA0, ReB, ImB ", ReA0, ImA0, ReB, ImB
1009  ! write(*,*) "ReX, ImX, ReY, ImY ", ReX, ImX, ReY, ImY
1010  ! write(*,*) "GSWr ", GSWr
1011  ! write(*,*) "GSWi ", GSWi
1012  ! endif
1013  ! IF(R(4,4)/R(1,1).gt.1.0) then
1014  ! write(*,*) "==============",R(4,4)/R(1,1)
1015  ! write(*,*) 'energuy,theta,channel= ',Energy,' ',cos(theta),' ',channel
1016  ! do k=1,4
1017  ! 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)
1018  ! enddo
1019  ! endif
1020 
1021 
1022  7 format(1x,i4, 4(2x,f12.6))
1023 
1024  end
1025 
1026  Subroutine dipolqqrijradcor (iqed,Energy,theta,Amz0,Gamz0,sin2W0,alphaQED,
1027  # GSWr,GSWi,ReA0,ImA0,ReB,ImB,ReX,ImX,ReY,ImY,R,channel)
1028 C version of dipole routine with electroweak form factors included,
1029 C If ported to different applications, beware of electroweak schemes, constant
1030 C and form-factors initialization
1031 C qed anomalous magnetic moment coded, but commented out
1032 C (comments) #####################################
1033 C 1) code is adopted to run with KKMC
1034 C 2) and EW for channel=1 (initial state e, final state tau or mu.
1035 C 3) internal activation of electroweak KKMC libs with IfGSW=KeyElw, call BornV_GetKeyElw(KeyElw)
1036 C 4) note distinct conventions, papers: e-Print: 2307.03526 versus Eur.Phys.J.C 79 (2019) 6, 480
1037 C that is about photon/Z normalization too
1038 C ########################################################
1039 
1040 c PROCESS fermion_i + fermionbar_i -> fermion_f- + fermionbar_f
1041 c
1042 c IMPROVED BORN APPROXIMATION (IBA) with radiative corrections
1043 c Integer parameter channel = 1 for LEPTONS, 2 -for UP QUARK, 3 - for DOWN QUARK
1044 c
1045 c Photon and Z-boson exchanges are included.
1046 c Calculates spin-correlation coefficients R(i,j) as functions of beam
1047 c energy: Energy = sqrt(s)/2 (in GeV), and scattaring angle theta.
1048 c
1049 c Anomalous magnetic formfactor A = ReA +i*ImA,
1050 c and electric dipole formfactor B = ReB + i*ImB are complex.
1051 c Also weak anomalous magnetic formfactor X = ReX + i*ImX and
1052 c weak electric dipole formfactor Y = ReY + i*ImY are complex.
1053 c
1054 c V is velocity of tau, gam is Lorentz factor, alpha is fine-structure constant,
1055 c G_F is Fermi constant of weak interaction.
1056 c RSM(i,j) - spin-correlation matrix in IBA without anomalous moments.
1057 c RDM(i,j) - spin-correlation matrix, linear in anomalous dipole moments A, B, X, Y.
1058 c R(i,j) - total spin-correlation matrix.
1059 c Order of coefficients R(ij): i = 1,2,3 correspond to S_x, S_y, S_z for tau-,
1060 c j = 1,2,3 correspond to S'_x, S'_y, S'_z for tau+,
1061 c i = j = 4 corresponds to 1 (no spin dependence).
1062 c
1063 c Functions below and their notation are from the paper:
1064 c E. Richter-Was, Z. Was Eur. Phys.J. C (2019) 79, 480
1065 c
1066 c indices i = 1, 2, 3 correspond to lepton, up quark, down quark, respectively
1067 c Gamma_vp - vacuum polarization factor (function of s)
1068 c Rho_11, Rho_12, Rho_13 =rho_{lepton,i}(s,t) for i= L,Up,Down
1069 c K_1,K_2,K_3 = K_i(s,t) for i= L, Up, Down
1070 c K_11,K_12,K_13 = K_{lepton,i}(s,t) for i = L,Up,Down
1071 
1072  Implicit none
1073 ! INCLUDE '../../basf2/KK2f/GPS.fi'
1074  real*8 m_sw2,m_mz,m_gammz,m_gmu,m_alfinv,m_swsq,m_pi
1075  ! INCLUDE '../../basf2/bornv/BornV.fi'
1076 cc not needed external function K_1, K_2, K_3, K_11, K_12, K_13
1077 cc not needed external function Rho_11,Rho_12,Rho_13,Gamma_vp
1078  complex*16 A,B,X,Y, Pg, Pz, Den, vi, vf
1079  complex*16 K1,K2,K3,K11,K12,K13,Gammavp,Rho11,Rho12,Rho13
1080  DOUBLE COMPLEX GSW(100)
1081  real*8 amz0,gamz0,sin2w0,alphaqed,gswr(7),gswi(7)
1082  DOUBLE PRECISION Svar,CosThetD
1083  DOUBLE COMPLEX RhoEW, VPgamma, CorEle, CorFin, CorEleFin, VVCef
1084  integer i, j, channel,iqed,KFf,IfGSW,KeyElw,IfPrint
1085  real*8 energy,theta,rea,ima,rea0,ima0,reb,imb,rex,imx,rey,imy,gam,v
1086  real*8 arqed, aiqed, ar1, ai1
1087  real*8 e, f, m, qi, qf, ai, af, sw2, sw, cw2,cw, s, t, fvivf
1088  real*8 r(1:4, 1:4), rsm(1:4, 1:4), rdm(1:4, 1:4)
1089  real*8 pi/3.141592653589793238d0/,alpha/7.2973525693d-3/
1090  real*8 mz/91.1876d0/,gz/2.4952d0/ ! Mass and width of Z in GeV
1091 cc real*8 v0/-0.03783d0/,a0/-0.50123d0/ ! vector and axial couplings of leptons
1092  real*8 mtau/1.77686d0/ ! mass of tau in GeV
1093 
1094  real*8 g_f/1.1663788d-5/ ! from PDG Fermi constant, GeV^{-2}
1095  real*8 g_mu/1.166389d-5/,mw/80.353d0/ ! G_mu and M_W from Eur.Phys.J. C (2019) 79, 480
1096 
1097 C ==================================================================
1098 C definition of form-factors as <in-line> functions. Complex form-factors
1099 C initialized from KKMC libraries will be used. Note that electroweak
1100 C scheme need to be reconsidered if function is used for other purposes.
1101  k3(s,t) = corele
1102  k1(s,t) = corele
1103  k13(s,t) = corelefin
1104  k2(s,t) = corele
1105  k12(s,t) = corelefin
1106  k11(s,t) = corelefin
1107 
1108  gammavp(s) = vpgamma
1109  rho11(s,t) = rhoew
1110  rho12(s,t) = rhoew
1111  rho13(s,t) = rhoew
1112 C ===================================================================
1113 
1114 
1115  kff=15 ! 415-400
1116  svar=4*energy**2
1117  costhetd=cos(theta)
1118 ! RhoEW = GSW(1)
1119 ! VPgamma = GSW(6)
1120 ! CorEle = GSW(2)
1121 ! CorFin = GSW(3)
1122 ! CorEleFin = GSW(4)
1123  ifgsw=1 ! switch to activate electroweak formfactor
1124 ! call BornV_GetKeyElw(KeyElw)
1125  ifgsw=keyelw
1126  ifgsw=1
1127  m_pi=pi
1128  m_gmu=g_mu
1129  ifprint=0 ! switch to activate prints
1130  If(ifgsw.eq.1) then
1131 ! CALL BornV_InterpoGSW(KFf,Svar,CosThetD)
1132 ! CALL BornV_GetGSW(GSW)
1133  rhoew = dcmplx(gswr(1),gswi(1))
1134  vpgamma = dcmplx(gswr(7),gswi(7))
1135  vpgamma = 1.d0 /(2.d0-dcmplx(gswr(7),gswi(7)))
1136  corele = dcmplx(gswr(2),gswi(2))
1137  corfin = dcmplx(gswr(3),gswi(3))
1138  corelefin = dcmplx(gswr(4),gswi(4))
1139  sw2 = sin2w0 ! m_Sw2 ! 0.22351946d0 ! sin(theta_w)^2 from Eur.Phys.J. C (2019) 79, 480
1140  m_swsq=sw2
1141  m_sw2=sw2
1142 ! alpha=1.d0/m_AlfInv
1143  alpha=alphaqed
1144  m_alfinv=1.d0/alpha
1145  mz=amz0 ! m_MZ
1146  m_mz=mz
1147  gz=m_gammz*svar/mz**2 ! running width
1148  m_gammz=gamz0
1149  gz=gamz0*svar/mz**2 ! running width
1150  if (ifprint.eq.1) then
1151  write(*,*) ' m_KeyElw= ', keyelw
1152  write(*,*) 'sw2,alpha,Mz,Gz,svar=', sw2,alpha,mz,gz,svar
1153  write(*,*) 'G_mu, m_gmu= ', g_mu, m_gmu
1154  write(*,*) 'RhoEW = ',rhoew!, GSWr(1)
1155  write(*,*) 'VPgamma = ',vpgamma!, GSWr(6)
1156  write(*,*) 'gsw(6) (7)= ', gswr(6), gswr(7)
1157  write(*,*) 'CorEle = ', corele!, GSWr(2)
1158  write(*,*) 'CorFin = ',corfin !, GSWr(3)
1159  write(*,*) 'CorEleFin = ',corelefin!, GSWr(5)
1160  write(*,*) 'CorEleFin -CorEle*CorFin = ',corelefin-corele*corfin
1161 ! stop
1162  endif
1163  else
1164  rhoew = 1.d0
1165  vpgamma = 1.d0
1166  corele = 1.d0
1167  corfin = 1.d0
1168  corelefin = 1.d0
1169  sw2 = 0.22351946d0 ! sin(theta_w)^2 from Eur.Phys.J. C (2019) 79, 480
1170  m_swsq=sw2
1171  endif
1172 
1173  v= dsqrt(1.d0 -(mtau/energy)**2) ! tau velocity
1174 
1175 c Contributions to ReA(s) and ImA(s) from QED in one loop
1176  if(iqed.eq.10) then ! Warning: temporarily this contribution is blocked
1177  arqed = -alpha*mtau**2 /(pi*v*4.d0*energy**2)*dlog((1.d0+v)/(1.d0-v))
1178  aiqed = alpha *mtau**2 /(v *4.d0 *energy**2)
1179  else
1180  arqed = 0.d0 !switch for dipole QED part
1181  aiqed = 0.d0
1182  endif
1183 
1184 
1185 
1186  rea = arqed + rea0 ! QED + new physics
1187  ima = aiqed + ima0 ! QED + new physics
1188 
1189 
1190  m= mtau
1191  gam= energy/mtau ! Lorentz factor of tau
1192  e= dsqrt(4.0d0 *pi *alpha) ! electric charge
1193 c Below coupling f=g_w/(2*cos(theta_w))=e/(2*cos(theta_w)*sin(theta_w)),
1194 c which is expressed via Fermi constant:
1195  f = mz *dsqrt(g_mu *dsqrt(2.d0))
1196 
1197 
1198  sw = dsqrt(sw2) ! sin(theta_w)
1199  cw2 = 1.d0 - sw2
1200  cw = dsqrt(cw2)
1201 
1202 c Mandelstam variables s and t (mass of tau is included in t)
1203  s = 4.d0 *energy**2
1204  t = m**2 -0.5d0 *s *(1.d0 -v*dcos(theta))
1205 
1206 c Denominator of the Z propagator with running width
1207  den = s - mz**2 + (0.d0, 1.d0) *gz *mz
1208 c Constant for effective photon like couplings correction coming from Z boson
1209  fvivf = 4.d0 *(f**2/e**2) *sw**4 ! WARNING: the *sw**4 factor come from vector couplings numerators Fvivf
1210 c Note: Fvivf can also be written as Fvivf= sw**2/cw**2, or through the constant G_mu:
1211 c Fvivf= sqrt(2) *G_mu *Mz**2 *sw**4 /(pi*alpha) This agrees with definition below.
1212 
1213 c FINAL FERMION IS TAU LEPTON (f=lepton) with couplings:
1214  qf= -1.d0
1215  vf= -0.03783d0 ! from PDG, used in previous code (v0 coupling)
1216  af= -0.50123d0 ! from PDG, used in previous code (a0 coupling)
1217  if(ifgsw.eq.1) then
1218  vf= -0.5d0 - 2.d0 *qf *sw2 * k1(s,t) ! is function for lepto
1219  af= -0.5d0
1220  endif
1221 
1222 c
1223 c COUPLINGS FOR INITIAL FERMIONS e u or d
1224  IF (channel.EQ.1) THEN ! LEPTONS
1225  qi= -1.d0
1226  vi= -0.03783d0 ! from PDG, used in the previous code
1227  ai= -0.50123d0 ! from PDG, used in the previous code
1228 c
1229 c Effective Z propagator (reduces to 1/Den without rad. corr.):
1230  pz = rho11(s,t) /den
1231  if(ifgsw.eq.1) then
1232  vi= -0.5d0 - 2.d0 *qi *sw2 *k1(s,t)
1233  ai= -0.5d0
1234  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'
1235  fvivf=
1236  $ m_gmu *m_mz**2 *m_alfinv /(dsqrt(2.d0)*8.d0*m_pi)
1237  $ *(m_sw2*(1.d0-m_sw2)) *16.d0
1238  fvivf=fvivf *m_sw2/(1.d0-m_sw2) ! why this factor? Because vi ai couplings in KKMC
1239  ! are divided by deno= 4 sqrt(m_swsq*(1d0-m_swsq)) in addition multiplied by 2 and we are not using it here
1240  ! Ve = (2*T3e -4*Qe*m_swsq)/deno and Ae = 2*T3e/deno.
1241  if (ifprint.eq.1) then
1242  write(*,*) 'with EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
1243  stop
1244  endif
1245  endif
1246 
1247 C Effective photon propagator with v_{if}-v_i*v_f correction to Z exchange
1248 C (NOTE: Pz reduces to 1/s if rad. corr. are off):
1249  pg = 1.d0/s *(gammavp(s) + fvivf *s/den*
1250  $ rho11(s,t)*(k11(s,t)-k1(s,t)*k1(s,t)))
1251 
1252 c
1253  ELSE IF (channel.EQ.2) THEN ! UP QUARK
1254  qi= +2.d0/3.d0
1255  vi= +0.266d0 ! from PDG, used in the previous code
1256  ai= +0.519d0 ! from PDG, used in the previous code
1257  if(ifgsw.eq.1) then
1258  vi= +0.5d0 - 2.d0 *qi *sw2 *k2(s,t)
1259  ai= +0.5d0
1260  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'
1261  fvivf=
1262  $ m_gmu *m_mz**2 *m_alfinv /(dsqrt(2.d0)*8.d0*m_pi)
1263  $ *(m_sw2*(1.d0-m_sw2)) *16.d0
1264  fvivf=fvivf *m_sw2/(1.d0-m_sw2) ! why this factor? Because vi ai couplings in KKMC
1265  ! are divided by deno= 4 sqrt(m_swsq*(1d0-m_swsq)) in addition multiplied by 2.
1266  ! Ve = (2*T3e -4*Qe*m_swsq)/deno and Ae = 2*T3e/deno. Note that Fvivf is used only for calculation of
1267  ! (vv_if-v_i*v_f)
1268  if (ifprint.eq.1) then
1269  write(*,*) 'with EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
1270  stop
1271  endif
1272  endif
1273 
1274 c
1275 c Effective Z propagator (reduces to 1/Den without rad. corr.):
1276  pz = rho12(s,t) /den
1277 c Effective photon propagator with v_{if}-v_i*v_f correction to Z exchange
1278 C (NOTE: Pz reduces to 1/s if rad. corr. are off):
1279  pg = 1.d0/s *(gammavp(s) + fvivf *s/den*
1280  $ rho12(s,t)*(k12(s,t)-k1(s,t)*k2(s,t)))
1281 c
1282  ELSE IF (channel.EQ.3) THEN ! DOWN QUARK
1283  qi= -1.d0/3.d0
1284  vi= -0.38d0 ! from PDG, used in the previous code
1285  ai= -0.527d0 ! from PDG, used in the previous code
1286  if(ifgsw.eq.1) then
1287  vi= -0.5d0 - 2.d0 *qi *sw2 *k3(s,t)
1288  ai= -0.5d0
1289  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'
1290  fvivf=
1291  $ m_gmu *m_mz**2 *m_alfinv /(dsqrt(2.d0)*8.d0*m_pi)
1292  $ *(m_sw2*(1.d0-m_sw2)) *16.d0
1293  fvivf=fvivf *m_sw2/(1.d0-m_sw2) ! why this factor? Because vi ai couplings in KKMC
1294  ! are divided by deno= 4 sqrt(m_swsq*(1d0-m_swsq)) in addition multiplied by 2.
1295  ! Ve = (2*T3e -4*Qe*m_swsq)/deno and Ae = 2*T3e/deno. Note that Fvivf is used only for calculation of
1296  ! (vv_if-v_i*v_f)
1297  if (ifprint.eq.1) then
1298  write(*,*) 'with EW-corr Fvivf=',fvivf,' Fvivf is for (v_if-v_i*v_f) implemented as add-up to photon'
1299  stop
1300  endif
1301  endif
1302 
1303 c
1304 c Effective Z propagator (reduces to 1/Den without rad. corr.):
1305  pz = rho13(s,t) /den
1306 c Effective photon propagator (reduces to 1/s without rad. corr.):
1307  pg = 1.d0/s *(gammavp(s) + fvivf *s/den*
1308  $ rho13(s,t) *(k13(s,t)-k1(s,t)*k3(s,t)))
1309 c
1310  ELSE ! OUTPUT WILL BE ZERO
1311  qi= 0.d0
1312  vi= (0.d0, 0.d0)
1313  ai= 0.d0
1314  pz= (0.d0, 0.d0)
1315  pg= (0.d0, 0.d0)
1316  END IF
1317 
1318 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1319 c Complex magnetic and electric form-factors
1320  a = dcmplx(rea, ima)
1321  b = dcmplx(reb, imb)
1322  x = dcmplx(rex, imx)
1323  y = dcmplx(rey, imy)
1324 
1325 c Contributions to A(s), X(s) from QED in one loop: NOT INCLUDED
1326 cc ArQED = -alpha *mtau**2 /(PI *V *4.d0 *E**2) *dlog((1.d0+V)/(1.d0-V))
1327 cc AiQED = alpha *mtau**2 /(V *4.d0 *E**2)
1328 cc Ar1 = ArQED + Ar ! QED + new physics
1329 cc Ai1 = AiQED + Ai ! QED + new physics
1330 
1331 
1332 
1333 c ----------------------------------------------------------------
1334 c CONTRIBUTIONS to R(i,j) IN STANDARD MODEL NO DIPOLE MOMENTS
1335 c ----------------------------------------------------------------
1336 c
1337 
1338  rsm(1,1)= 4.d0*gam**2*m**4*(e**2*(1.d0 + gam**2)*qf*qi*
1339  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(pg) +
1340  - f**2*dconjg(pz)*
1341  - (-(af**2*f**2*(-1.d0 + gam**2)*pz*
1342  - (ai**2 + vi*dconjg(vi))) +
1343  - (1.d0 + gam**2)*dconjg(vf)*
1344  - (ai**2*f**2*pz*vf +
1345  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))))*
1346  - dsin(theta)**2
1347 
1348  rsm(1,2)= (0.d0, -4.d0)*af*f**2*gam**4*m**4*v*
1349  - (e**2*pz*qf*qi*vi*dconjg(pg) +
1350  - dconjg(pz)*(-(ai**2*f**2*pz*vf) -
1351  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi) +
1352  - f**2*pz*dconjg(vf)*(ai**2 + vi*dconjg(vi))))*
1353  - dsin(theta)**2
1354 
1355  rsm(2,1)= rsm(1,2)
1356 
1357  rsm(2,2)= 4.d0*gam**2*(-1.d0 + gam**2)*m**4*
1358  - (-(e**2*qf*qi*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1359  - dconjg(pg)) +
1360  - dconjg(pz)*(af**2*f**4*pz*
1361  - (ai**2 + vi*dconjg(vi)) -
1362  - f**2*dconjg(vf)*
1363  - (ai**2*f**2*pz*vf +
1364  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))))*
1365  - dsin(theta)**2
1366 
1367  rsm(1,3)= 4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1368  - (af*ai*f**2*pz*v -
1369  - 2.d0*gam**2*(-1.d0 + v**2)*
1370  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dcos(theta)) +
1371  - f**2*dconjg(pz)*
1372  - (af*(ai*v*(e**2*pg*qf*qi +
1373  - f**2*pz*vf*(vi + dconjg(vi))) +
1374  - 2.d0*af*f**2*pz*(1.d0 + gam**2*(-1.d0 + v**2))*
1375  - (ai**2 + vi*dconjg(vi))*dcos(theta)) +
1376  - dconjg(vf)*
1377  - (ai*f**2*pz*(af*v*vi -
1378  - 2.d0*ai*gam**2*(-1.d0 + v**2)*vf*dcos(theta)) +
1379  - dconjg(vi)*
1380  - (af*ai*f**2*pz*v -
1381  - 2.d0*gam**2*(-1.d0 + v**2)*
1382  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dcos(theta))))
1383  - )*dsin(theta)
1384 
1385  rsm(3,1)= rsm(1,3)
1386 
1387  rsm(2,3)= (0.d0, -2.d0)*af*f**2*gam**3*m**4*v*
1388  - (e**2*pz*qf*qi*vi*dconjg(pg) +
1389  - dconjg(pz)*(-(ai**2*f**2*pz*vf) -
1390  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi) +
1391  - f**2*pz*dconjg(vf)*(ai**2 + vi*dconjg(vi))))*
1392  - dsin(2.d0*theta)
1393 
1394  rsm(3,2)= rsm(2,3)
1395 
1396  rsm(3,3)= 2.d0*gam**2*m**4*(-(e**4*pg*qf**2*qi**2*dconjg(pg)*
1397  - (1.d0 - 3.d0*gam**2 +
1398  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1399  - dcos(2.d0*theta))) -
1400  - e**2*f**2*qf*qi*(pz*dconjg(pg)*
1401  - (-4.d0*af*ai*gam**2*v*dcos(theta) +
1402  - vf*vi*(1.d0 - 3.d0*gam**2 +
1403  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1404  - dcos(2.d0*theta))) +
1405  - pg*dconjg(pz)*
1406  - (-4.d0*af*ai*gam**2*v*dcos(theta) +
1407  - dconjg(vf)*dconjg(vi)*
1408  - (1.d0 - 3.d0*gam**2 +
1409  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1410  - dcos(2.d0*theta)))) +
1411  - f**4*pz*dconjg(pz)*
1412  - (af*(4*ai*gam**2*v*vf*(vi + dconjg(vi))*
1413  - dcos(theta) +
1414  - af*(ai**2 + vi*dconjg(vi))*
1415  - (1.d0 + gam**2*(-1.d0 + 4.d0*v**2) +
1416  - (-1.d0 + 5.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1417  - dcos(2.d0*theta))) +
1418  - dconjg(vf)*
1419  - (dconjg(vi)*
1420  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1421  - vf*vi*(-1.d0 + 3.d0*gam**2 +
1422  - (1.d0 - 3.d0*gam**2 - 4.d0*gam**4*(-1.d0 + v**2))*
1423  - dcos(2.d0*theta))) +
1424  - ai*(4.d0*af*gam**2*v*vi*dcos(theta) -
1425  - ai*vf*(1.d0 - 3.d0*gam**2 +
1426  - (-1.d0 + 3.d0*gam**2 + 4.d0*gam**4*(-1.d0 + v**2))*
1427  - dcos(2.d0*theta))))))
1428 
1429  rsm(1,4)= -4.d0*f**2*gam**3*m**4*
1430  - (e**2*pz*qf*qi*dconjg(pg)*
1431  - (2.d0*ai*vf + af*v*vi*dcos(theta)) +
1432  - dconjg(pz)*(af*v*
1433  - (ai**2*f**2*pz*vf +
1434  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1435  - dcos(theta) +
1436  - dconjg(vf)*
1437  - (2.d0*ai*(e**2*pg*qf*qi +
1438  - f**2*pz*vf*(vi + dconjg(vi))) +
1439  - af*f**2*pz*v*(ai**2 + vi*dconjg(vi))*
1440  - dcos(theta))))*dsin(theta)
1441 
1442  rsm(4,1)= rsm(1,4)
1443 
1444  rsm(2,4)= (0.d0, 4.d0)*af*ai*f**2*gam**3*m**4*v*
1445  - (e**2*pz*qf*qi*dconjg(pg) +
1446  - dconjg(pz)*(-(e**2*pg*qf*qi) -
1447  - f**2*pz*vf*(vi + dconjg(vi)) +
1448  - f**2*pz*dconjg(vf)*(vi + dconjg(vi))))*
1449  - dsin(theta)
1450 
1451  rsm(4,2)= rsm(2,4)
1452 
1453  rsm(3,4)= -4.d0*f**2*gam**2*m**4*
1454  - (e**2*gam**2*pz*qf*qi*dconjg(pg)*
1455  - (2.d0*ai*vf*dcos(theta) + af*v*vi*(1.d0 + dcos(theta)**2)) +
1456  - (dconjg(pz)*(dconjg(vi)*
1457  - (4.d0*af**2*ai*f**2*(-1.d0 + gam**2)*pz*dcos(theta) +
1458  - 4.d0*ai*f**2*gam**2*pz*vf*dconjg(vf)*
1459  - dcos(theta) +
1460  - af*gam**2*v*
1461  - (e**2*pg*qf*qi +
1462  - f**2*pz*vi*(vf + dconjg(vf)))*
1463  - (3.d0 + dcos(2.d0*theta))) +
1464  - ai*(gam**2*dconjg(vf)*
1465  - (4.d0*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1466  - dcos(theta) +
1467  - af*ai*f**2*pz*v*(3.d0 + dcos(2.d0*theta))) +
1468  - af*f**2*pz*
1469  - (4.d0*af*(-1.d0 + gam**2)*vi*dcos(theta) +
1470  - ai*gam**2*v*vf*(3.d0 + dcos(2.d0*theta))))))/2.d0)
1471 
1472  rsm(4,3)= rsm(3,4)
1473 
1474 
1475  rsm(4,4)= 2.d0*gam**2*m**4*(e**4*pg*qf**2*qi**2*dconjg(pg)*
1476  - (1.d0 + 3.d0*gam**2 + (-1.d0 + gam**2)*dcos(2.d0*theta)) +
1477  - e**2*f**2*qf*qi*(pz*dconjg(pg)*
1478  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1479  - vf*vi*(1.d0 + 3.d0*gam**2 + (-1.d0 + gam**2)*dcos(2.d0*theta))
1480  - ) + pg*dconjg(pz)*
1481  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1482  - dconjg(vf)*dconjg(vi)*
1483  - (1.d0 + 3.d0*gam**2 + (-1.d0 + gam**2)*dcos(2.d0*theta)))) -
1484  - f**4*pz*dconjg(pz)*
1485  - (-(af*(4.d0*ai*gam**2*v*vf*(vi + dconjg(vi))*
1486  - dcos(theta) +
1487  - af*(-1.d0 + gam**2)*(ai**2 + vi*dconjg(vi))*
1488  - (3.d0 + dcos(2.d0*theta)))) -
1489  - dconjg(vf)*
1490  - (ai*(4.d0*af*gam**2*v*vi*dcos(theta) +
1491  - ai*vf*(1.d0 + 3.d0*gam**2 +
1492  - (-1.d0 + gam**2)*dcos(2.d0*theta))) +
1493  - dconjg(vi)*
1494  - (4.d0*af*ai*gam**2*v*dcos(theta) +
1495  - vf*vi*(1.d0 + 3.d0*gam**2 +
1496  - (-1.d0 + gam**2)*dcos(2.d0*theta))))))
1497 
1498 
1499 c -------------------------------------------------------------
1500 c CONTRIBUTION TO R(i,j), LINEAR IN DIPOLE MOMENTS A, B, X, Y
1501 c -------------------------------------------------------------
1502 
1503  rdm(1,1)= 8.d0*gam**4*m**4*(e**2*qf*qi*
1504  - (e**2*pg*qf*qi*(a + dconjg(a)) +
1505  - f**2*pz*vi*(x + vf*dconjg(a)))*dconjg(pg) +
1506  - f**2*dconjg(pz)*
1507  - (dconjg(vf)*(ai**2*f**2*pz*x +
1508  - (a*e**2*pg*qf*qi + f**2*pz*vi*x)*dconjg(vi)) +
1509  - (ai**2*f**2*pz*vf +
1510  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1511  - dconjg(x)))*dsin(theta)**2
1512 
1513  rdm(1,2)= 4.d0*gam**4*m**4*v*(e**2*qf*qi*
1514  - (e**2*pg*qf*qi*(b + dconjg(b)) +
1515  - f**2*pz*vi*(y - (0.d0, 1.d0)*af*dconjg(a) +
1516  - vf*dconjg(b)))*dconjg(pg) +
1517  - f**2*dconjg(pz)*
1518  - (dconjg(vf)*(ai**2*f**2*pz*y +
1519  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dconjg(vi)) +
1520  - (0.d0, 1.d0)*af*(dconjg(vi)*
1521  - (a*e**2*pg*qf*qi +
1522  - f**2*pz*vi*(x - dconjg(x))) +
1523  - ai**2*f**2*pz*(x - dconjg(x))) +
1524  - (ai**2*f**2*pz*vf +
1525  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1526  - dconjg(y)))*dsin(theta)**2
1527 
1528  rdm(2,1)= -4.d0*gam**4*m**4*v*(e**2*qf*qi*
1529  - (e**2*pg*qf*qi*(b + dconjg(b)) +
1530  - f**2*pz*vi*(y + (0,1)*af*dconjg(a) +
1531  - vf*dconjg(b)))*dconjg(pg) +
1532  - f**2*dconjg(pz)*
1533  - (dconjg(vf)*(ai**2*f**2*pz*y +
1534  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dconjg(vi)) -
1535  - (0.d0 ,1.d0)*af*(dconjg(vi)*
1536  - (a*e**2*pg*qf*qi +
1537  - f**2*pz*vi*(x - dconjg(x))) +
1538  - ai**2*f**2*pz*(x - dconjg(x))) +
1539  - (ai**2*f**2*pz*vf +
1540  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))*
1541  - dconjg(y)))*dsin(theta)**2
1542 
1543  rdm(2,2)= 0.d0
1544 
1545  rdm(1,3)= 4.d0*gam**5*m**4*(-(e**2*qf*qi*dconjg(pg)*
1546  - ((0.d0, -1.d0)*ai*f**2*pz*v*
1547  - (y - (0.d0, 1.d0)*af*dconjg(a) - vf*dconjg(b)) +
1548  - (e**2*pg*qf*qi*(-2 + v**2)*(a + dconjg(a)) +
1549  - f**2*pz*vi*
1550  - ((-2.d0 + v**2)*(x + vf*dconjg(a)) +
1551  - (0.d0, 1.d0)*af*v**2*dconjg(b)))*dcos(theta))) +
1552  - f**2*dconjg(pz)*
1553  - (dconjg(vf)*(ai*
1554  - ((0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) -
1555  - ai*f**2*pz*(-2.d0 + v**2)*x*dcos(theta)) +
1556  - dconjg(vi)*
1557  - ((0.d0, 1.d0)*ai*f**2*pz*v*y -
1558  - (-2.d0 + v**2)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1559  - dcos(theta))) +
1560  - dconjg(vi)*
1561  - ((0.d0,-1.d0)*ai*f**2*pz*v*vf*dconjg(y) -
1562  - (-2.d0 + v**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1563  - dconjg(x)*dcos(theta) +
1564  - af*v*(ai*f**2*pz*(x + dconjg(x)) +
1565  - (0.d0, 1.d0)*v*
1566  - (b*e**2*pg*qf*qi +
1567  - f**2*pz*vi*(y - dconjg(y)))*dcos(theta)))
1568  - + ai*((0.d0, -1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1569  - dconjg(y) +
1570  - ai*f**2*pz*(2.d0 - v**2)*vf*dconjg(x)*
1571  - dcos(theta) +
1572  - af*v*(a*e**2*pg*qf*qi +
1573  - f**2*pz*vi*(x + dconjg(x)) +
1574  - (0.d0, 1.d0)*ai*f**2*pz*v*(y - dconjg(y))*
1575  - dcos(theta)))))*dsin(theta)
1576 
1577  rdm(3,1)= 4.d0*gam**5*m**4*(e**2*qf*qi*dconjg(pg)*
1578  - (ai*f**2*pz*v*(af*dconjg(a) -
1579  - (0.d0, 1.d0)*(y - vf*dconjg(b))) -
1580  - (e**2*pg*qf*qi*(-2.d0 + v**2)*(a + dconjg(a)) +
1581  - f**2*pz*vi*
1582  - ((-2.d0 + v**2)*(x + vf*dconjg(a)) -
1583  - (0.d0, 1.d0)*af*v**2*dconjg(b)))*dcos(theta)) +
1584  - f**2*dconjg(pz)*
1585  - (dconjg(vf)*(-(ai*
1586  - ((0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) +
1587  - ai*f**2*pz*(-2.d0 + v**2)*x*dcos(theta))) +
1588  - dconjg(vi)*
1589  - ((0.d0, -1.d0)*ai*f**2*pz*v*y -
1590  - (-2.d0 + v**2)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1591  - dcos(theta))) +
1592  - dconjg(vi)*
1593  - ((0.d0, 1.d0)*ai*f**2*pz*v*vf*dconjg(y) -
1594  - (-2.d0 + v**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1595  - dconjg(x)*dcos(theta) +
1596  - af*v*(ai*f**2*pz*(x + dconjg(x)) -
1597  - (0.d0, 1.d0)*v*
1598  - (b*e**2*pg*qf*qi +
1599  - f**2*pz*vi*(y - dconjg(y)))*dcos(theta)))
1600  - + ai*((0.d0,1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1601  - dconjg(y) +
1602  - ai*f**2*pz*(2.d0 - v**2)*vf*dconjg(x)*
1603  - dcos(theta) +
1604  - af*v*(a*e**2*pg*qf*qi +
1605  - f**2*pz*vi*(x + dconjg(x)) -
1606  - (0.d0, 1.d0)*ai*f**2*pz*v*(y - dconjg(y))*
1607  - dcos(theta)))))*dsin(theta)
1608 
1609  rdm(2,3)= (0.d0, 4.d0)*gam**5*m**4*v*
1610  - (e**2*qf*qi*dconjg(pg)*
1611  - (ai*f**2*pz*v*(x - vf*dconjg(a) +
1612  - (0.d0, 1.d0)*af*dconjg(b)) +
1613  - (0.d0, 1.d0)*(e**2*pg*qf*qi*(b + dconjg(b)) +
1614  - f**2*pz*vi*
1615  - (y + (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b)))*
1616  - dcos(theta)) +
1617  - f**2*dconjg(pz)*
1618  - (dconjg(vf)*(ai*
1619  - (a*e**2*pg*qf*qi*v + f**2*pz*v*vi*x +
1620  - (0.d0, 1.d0)*ai*f**2*pz*y*dcos(theta)) +
1621  - dconjg(vi)*
1622  - (ai*f**2*pz*v*x +
1623  - (0.d0, 1.d0)*(b*e**2*pg*qf*qi + f**2*pz*vi*y)*
1624  - dcos(theta))) +
1625  - (0.d0, 1.d0)*(dconjg(vi)*
1626  - ((0.d0, 1.d0)*ai*f**2*pz*v*vf*dconjg(x) +
1627  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(y)*
1628  - dcos(theta) +
1629  - af*(ai*f**2*pz*v*(y + dconjg(y)) -
1630  - (0.d0, 1.d0)*
1631  - (a*e**2*pg*qf*qi +
1632  - f**2*pz*vi*(x - dconjg(x)))*
1633  - dcos(theta))) +
1634  - ai*((0.d0, 1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1635  - dconjg(x) +
1636  - ai*f**2*pz*vf*dconjg(y)*dcos(theta) +
1637  - af*(v*(b*e**2*pg*qf*qi +
1638  - f**2*pz*vi*(y + dconjg(y))) -
1639  - (0.d0, 1.d0)*ai*f**2*pz*(x - dconjg(x))*
1640  - dcos(theta))))))*dsin(theta)
1641 
1642  rdm(3,2)= 4.d0*gam**5*m**4*v*(e**2*qf*qi*dconjg(pg)*
1643  - (ai*f**2*pz*v*((0.d0, 1.d0)*(x - vf*dconjg(a)) +
1644  - af*dconjg(b)) +
1645  - (e**2*pg*qf*qi*(b + dconjg(b)) +
1646  - f**2*pz*vi*
1647  - (y - (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b)))*
1648  - dcos(theta)) +
1649  - f**2*dconjg(pz)*
1650  - (dconjg(vf)*(ai*
1651  - ((0.d0, 1.d0)*v*(a*e**2*pg*qf*qi + f**2*pz*vi*x) +
1652  - ai*f**2*pz*y*dcos(theta)) +
1653  - dconjg(vi)*
1654  - ((0.d0, 1.d0)*ai*f**2*pz*v*x +
1655  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dcos(theta))) +
1656  - dconjg(vi)*
1657  - ((0.d0, -1.d0)*ai*f**2*pz*v*vf*dconjg(x) +
1658  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(y)*
1659  - dcos(theta) +
1660  - af*(ai*f**2*pz*v*(y + dconjg(y)) +
1661  - (0.d0, 1.d0)*(a*e**2*pg*qf*qi +
1662  - f**2*pz*vi*(x - dconjg(x)))*dcos(theta)))
1663  - + ai*((0.d0, -1.d0)*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1664  - dconjg(x) +
1665  - ai*f**2*pz*vf*dconjg(y)*dcos(theta) +
1666  - af*(v*(b*e**2*pg*qf*qi +
1667  - f**2*pz*vi*(y + dconjg(y))) +
1668  - (0.d0, 1.d0)*ai*f**2*pz*(x - dconjg(x))*dcos(theta)
1669  - ))))*dsin(theta)
1670 
1671  rdm(3,3)= -8.d0*gam**6*m**4*(-1.d0 + v**2)*dcos(theta)*
1672  - (e**2*qf*qi*dconjg(pg)*
1673  - ((a*e**2*pg*qf*qi + f**2*pz*vi*x)*dcos(theta) +
1674  - dconjg(a)*(af*ai*f**2*pz*v +
1675  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dcos(theta))) +
1676  - f**2*dconjg(pz)*
1677  - (ai*(af*v*(a*e**2*pg*qf*qi +
1678  - f**2*pz*vi*(x + dconjg(x))) +
1679  - ai*f**2*pz*x*dconjg(vf)*dcos(theta) +
1680  - ai*f**2*pz*vf*dconjg(x)*dcos(theta)) +
1681  - dconjg(vi)*
1682  - (af*ai*f**2*pz*v*(x + dconjg(x)) +
1683  - ((a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1684  - dconjg(vf) +
1685  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(x))*
1686  - dcos(theta))))
1687 
1688  rdm(1,4)= -4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1689  - (ai*f**2*pz*((1 + gam**2)*(x + vf*dconjg(a)) -
1690  - (0.d0, 1.d0)*af*(-1.d0 + gam**2)*dconjg(b)) +
1691  - (0.d0, 1.d0)*gam**2*v*
1692  - (e**2*pg*qf*qi*(b - dconjg(b)) +
1693  - f**2*pz*vi*
1694  - (y - (0.d0, 1.d0)*af*dconjg(a) - vf*dconjg(b)))*
1695  - dcos(theta)) +
1696  - f**2*dconjg(pz)*
1697  - (dconjg(vf)*(ai*
1698  - ((1.d0 + gam**2)*
1699  - (a*e**2*pg*qf*qi + f**2*pz*vi*x) +
1700  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*y*dcos(theta)) +
1701  - dconjg(vi)*
1702  - (ai*f**2*(1.d0 + gam**2)*pz*x +
1703  - (0.d0, 1.d0)*gam**2*v*
1704  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dcos(theta)))
1705  - + ai*((1.d0 + gam**2)*
1706  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(x) -
1707  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*vf*dconjg(y)*
1708  - dcos(theta) +
1709  - af*((0.d0, 1.d0)*(-1.d0 + gam**2)*
1710  - (b*e**2*pg*qf*qi +
1711  - f**2*pz*vi*(y - dconjg(y))) +
1712  - ai*f**2*gam**2*pz*v*(x + dconjg(x))*
1713  - dcos(theta))) +
1714  - dconjg(vi)*
1715  - (ai*f**2*(1.d0 + gam**2)*pz*vf*dconjg(x) -
1716  - (0.d0, 1.d0)*gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1717  - dconjg(y)*dcos(theta) +
1718  - af*((0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*
1719  - (y - dconjg(y)) +
1720  - gam**2*v*
1721  - (a*e**2*pg*qf*qi +
1722  - f**2*pz*vi*(x + dconjg(x)))*dcos(theta)))
1723  - ))*dsin(theta)
1724 
1725  rdm(4,1)= -4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1726  - (ai*f**2*pz*((1.d0 + gam**2)*(x + vf*dconjg(a)) +
1727  - (0.d0, 1.d0)*af*(-1.d0 + gam**2)*dconjg(b)) -
1728  - (0.d0, 1.d0)*gam**2*v*
1729  - (e**2*pg*qf*qi*(b - dconjg(b)) +
1730  - f**2*pz*vi*
1731  - (y + (0.d0, 1.d0)*af*dconjg(a) - vf*dconjg(b)))*
1732  - dcos(theta)) +
1733  - f**2*dconjg(pz)*
1734  - (dconjg(vf)*(ai*
1735  - ((1.d0 + gam**2)*
1736  - (a*e**2*pg*qf*qi + f**2*pz*vi*x) -
1737  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*y*dcos(theta)) +
1738  - dconjg(vi)*
1739  - (ai*f**2*(1.d0 + gam**2)*pz*x -
1740  - (0.d0, 1.d0)*gam**2*v*
1741  - (b*e**2*pg*qf*qi + f**2*pz*vi*y)*dcos(theta)))
1742  - + ai*((1.d0 + gam**2)*
1743  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(x) +
1744  - (0.d0, 1.d0)*ai*f**2*gam**2*pz*v*vf*dconjg(y)*
1745  - dcos(theta) +
1746  - af*((0.d0, -1.d0)*(-1.d0 + gam**2)*
1747  - (b*e**2*pg*qf*qi +
1748  - f**2*pz*vi*(y - dconjg(y))) +
1749  - ai*f**2*gam**2*pz*v*(x + dconjg(x))*
1750  - dcos(theta))) +
1751  - dconjg(vi)*
1752  - (ai*f**2*(1.d0 + gam**2)*pz*vf*dconjg(x) +
1753  - (0.d0, 1.d0)*gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1754  - dconjg(y)*dcos(theta) +
1755  - af*((0.d0, -1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*
1756  - (y - dconjg(y)) +
1757  - gam**2*v*
1758  - (a*e**2*pg*qf*qi +
1759  - f**2*pz*vi*(x + dconjg(x)))*dcos(theta)))
1760  - ))*dsin(theta)
1761 
1762  rdm(2,4)= 4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1763  - (ai*f**2*gam**2*pz*v*
1764  - (y + (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b)) -
1765  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1766  - (e**2*pg*qf*qi*(a - dconjg(a)) +
1767  - f**2*pz*vi*
1768  - (x - vf*dconjg(a) + (0.d0, 1.d0)*af*dconjg(b)))*
1769  - dcos(theta)) +
1770  - f**2*dconjg(pz)*
1771  - (dconjg(vf)*(ai*
1772  - (gam**2*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) -
1773  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*x*dcos(theta)) +
1774  - dconjg(vi)*
1775  - (ai*f**2*gam**2*pz*v*y -
1776  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1777  - (a*e**2*pg*qf*qi + f**2*pz*vi*x)*dcos(theta)))
1778  - + ai*(gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1779  - dconjg(y) +
1780  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*vf*dconjg(x)*
1781  - dcos(theta) +
1782  - af*((0.d0, -1.d0)*gam**2*v*
1783  - (a*e**2*pg*qf*qi +
1784  - f**2*pz*vi*(x - dconjg(x))) +
1785  - ai*f**2*(-1.d0 + gam**2)*pz*(y + dconjg(y))*
1786  - dcos(theta))) +
1787  - dconjg(vi)*
1788  - ((0.d0, 1.d0)*((0.d0, -1.d0)*ai*f**2*gam**2*pz*v*vf*
1789  - dconjg(y) +
1790  - (-1.d0 + gam**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1791  - dconjg(x)*dcos(theta)) +
1792  - af*((0.d0, -1.d0)*ai*f**2*gam**2*pz*v*
1793  - (x - dconjg(x)) +
1794  - (-1.d0 + gam**2)*
1795  - (b*e**2*pg*qf*qi +
1796  - f**2*pz*vi*(y + dconjg(y)))*dcos(theta)))
1797  - ))*dsin(theta)
1798 
1799  rdm(4,2)= 4.d0*gam**3*m**4*(e**2*qf*qi*dconjg(pg)*
1800  - (-(ai*f**2*gam**2*pz*v*
1801  - (y - (0.d0, 1.d0)*af*dconjg(a) + vf*dconjg(b))) -
1802  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1803  - (e**2*pg*qf*qi*(a - dconjg(a)) +
1804  - f**2*pz*vi*
1805  - (x - vf*dconjg(a) - (0.d0, 1.d0)*af*dconjg(b)))*
1806  - dcos(theta)) -
1807  - f**2*dconjg(pz)*
1808  - (dconjg(vf)*(ai*
1809  - (gam**2*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y) +
1810  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*x*dcos(theta)) +
1811  - dconjg(vi)*
1812  - (ai*f**2*gam**2*pz*v*y +
1813  - (0.d0, 1.d0)*(-1.d0 + gam**2)*
1814  - (a*e**2*pg*qf*qi + f**2*pz*vi*x)*dcos(theta)))
1815  - + ai*(gam**2*v*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1816  - dconjg(y) -
1817  - (0.d0, 1.d0)*ai*f**2*(-1.d0 + gam**2)*pz*vf*dconjg(x)*
1818  - dcos(theta) +
1819  - af*((0.d0, 1.d0)*gam**2*v*
1820  - (a*e**2*pg*qf*qi +
1821  - f**2*pz*vi*(x - dconjg(x))) +
1822  - ai*f**2*(-1.d0 + gam**2)*pz*(y + dconjg(y))*
1823  - dcos(theta))) +
1824  - dconjg(vi)*
1825  - ((0.d0, -1.d0)*((0.d0, 1.d0)*ai*f**2*gam**2*pz*v*vf*
1826  - dconjg(y) +
1827  - (-1.d0 + gam**2)*(e**2*pg*qf*qi + f**2*pz*vf*vi)*
1828  - dconjg(x)*dcos(theta)) +
1829  - af*((0.d0, 1.d0)*ai*f**2*gam**2*pz*v*
1830  - (x - dconjg(x)) +
1831  - (-1.d0 + gam**2)*
1832  - (b*e**2*pg*qf*qi +
1833  - f**2*pz*vi*(y + dconjg(y)))*dcos(theta)))
1834  - ))*dsin(theta)
1835 
1836  rdm(3,4)= -2.d0*gam**4*m**4*(f**2*dconjg(pz)*
1837  - (-(af*v*(ai**2*f**2*pz*(x + dconjg(x)) +
1838  - dconjg(vi)*
1839  - (a*e**2*pg*qf*qi +
1840  - f**2*pz*vi*(x + dconjg(x))))*
1841  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1842  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1843  - (0.d0, 1.d0)*((0.d0, -4.d0)*ai*
1844  - (e**2*pg*qf*qi +
1845  - f**2*pz*vf*(vi + dconjg(vi)))*dconjg(x)*
1846  - dcos(theta) +
1847  - v*(ai**2*f**2*pz*vf +
1848  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))
1849  - *dconjg(y)*
1850  - (2.d0 + gam**2*(-1.d0 + v**2) +
1851  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1852  - dconjg(vf)*
1853  - ((0.d0, -1.d0)*ai*((0.d0, 4.d0)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1854  - dcos(theta) +
1855  - ai*f**2*pz*v*y*
1856  - (2.d0 + gam**2*(-1.d0 + v**2) +
1857  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1858  - dconjg(vi)*
1859  - (4.d0*ai*f**2*pz*x*dcos(theta) -
1860  - (0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y)*
1861  - (2.d0 + gam**2*(-1.d0 + v**2) +
1862  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))) +
1863  - e**2*qf*qi*dconjg(pg)*
1864  - (4.d0*ai*f**2*pz*(x + vf*dconjg(a))*dcos(theta) -
1865  - (0.d0, 1.d0)*v*(e**2*pg*qf*qi*(b - dconjg(b))*
1866  - (2.d0 + gam**2*(-1.d0 + v**2) +
1867  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1868  - f**2*pz*vi*
1869  - ((0.d0, -1.d0)*af*dconjg(a)*
1870  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1871  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1872  - (y - vf*dconjg(b))*
1873  - (2.d0 + gam**2*(-1.d0 + v**2) +
1874  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))))
1875 
1876  rdm(4,3)= -2.d0*gam**4*m**4*(f**2*dconjg(pz)*
1877  - (-(af*v*(ai**2*f**2*pz*(x + dconjg(x)) +
1878  - dconjg(vi)*
1879  - (a*e**2*pg*qf*qi +
1880  - f**2*pz*vi*(x + dconjg(x))))*
1881  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1882  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) -
1883  - (0.d0, 1.d0)*((0.d0, 4.d0)*ai*
1884  - (e**2*pg*qf*qi +
1885  - f**2*pz*vf*(vi + dconjg(vi)))*dconjg(x)*
1886  - dcos(theta) +
1887  - v*(ai**2*f**2*pz*vf +
1888  - (e**2*pg*qf*qi + f**2*pz*vf*vi)*dconjg(vi))
1889  - *dconjg(y)*
1890  - (2.d0 + gam**2*(-1.d0 + v**2) +
1891  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1892  - dconjg(vf)*
1893  - ((0.d0, 1.d0)*ai*((0.d0, -4.d0)*(a*e**2*pg*qf*qi + f**2*pz*vi*x)*
1894  - dcos(theta) +
1895  - ai*f**2*pz*v*y*
1896  - (2.d0 + gam**2*(-1.d0 + v**2) +
1897  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))) +
1898  - dconjg(vi)*
1899  - (4.d0*ai*f**2*pz*x*dcos(theta) +
1900  - (0.d0, 1.d0)*v*(b*e**2*pg*qf*qi + f**2*pz*vi*y)*
1901  - (2.d0 + gam**2*(-1.d0 + v**2) +
1902  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))) +
1903  - e**2*qf*qi*dconjg(pg)*
1904  - (4.d0*ai*f**2*pz*(x + vf*dconjg(a))*dcos(theta) +
1905  - (0.d0, 1.d0)*v*(e**2*pg*qf*qi*(b - dconjg(b))*
1906  - (2.d0 + gam**2*(-1.d0 + v**2) +
1907  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1908  - f**2*pz*vi*
1909  - ((0.d0, 1.d0)*af*dconjg(a)*
1910  - (-2.d0 + gam**2*(-1.d0 + v**2) +
1911  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta)) +
1912  - (y - vf*dconjg(b))*
1913  - (2.d0 + gam**2*(-1.d0 + v**2) +
1914  - gam**2*(-1.d0 + v**2)*dcos(2.d0*theta))))))
1915 
1916  rdm(4,4)= 8.d0*gam**4*m**4*(e**4*pg*qf**2*qi**2*(a + dconjg(a))*
1917  - dconjg(pg) +
1918  - e**2*f**2*qf*qi*(pz*vi*x*dconjg(pg) +
1919  - a*pg*dconjg(pz)*dconjg(vf)*dconjg(vi) +
1920  - pg*dconjg(pz)*dconjg(vi)*dconjg(x) +
1921  - a*af*ai*pg*v*dconjg(pz)*dcos(theta) +
1922  - pz*dconjg(a)*dconjg(pg)*
1923  - (vf*vi + af*ai*v*dcos(theta))) +
1924  - f**4*pz*dconjg(pz)*
1925  - (ai*(ai*x*dconjg(vf) + ai*vf*dconjg(x) +
1926  - af*v*vi*x*dcos(theta) +
1927  - af*v*vi*dconjg(x)*dcos(theta)) +
1928  - dconjg(vi)*
1929  - (vi*(x*dconjg(vf) + vf*dconjg(x)) +
1930  - af*ai*v*(x + dconjg(x))*dcos(theta))))
1931 
1932 c-------------------------------------------------------------
1933 c TOTAL SPIN-CORRELATION MATRIX: STANDARD MODEL + DIPOLE MOMENTS
1934 c-------------------------------------------------------------
1935  do i=1,4
1936  do j=1,4
1937  r(i,j)= rsm(i,j) + rdm(i,j)
1938  enddo
1939  enddo
1940 
1941  return
1942  end
Generated by doxygen 1.9.1