C++ Interface to Tauola
curr_extracted.f
1  subroutine had1_init
2  implicit real*8 (a-h,o-z)
3  dimension p1(4),p2(4)
4 c
5  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
6  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
7  common /input/ su,su2,qq2,p1,p2,ngen,iseed,mode,iww,nhit
8  common /param/ pi,alpha,f_max
9  common /cbwgrho/ rhom2,rho1m2,rho2m2,omm2,rhomg,rho1mg,rho2mg
10  1 ,ommg
11  common /cbwga1/ a1m2,con
12  common /anom/amrop,gamrop,sig,amrop_2,amropg
13  common /cbwgrho_t/ rho1m2_t,rho1mg_t,beta
14 c
15  pi = 3.141592653589793238d0
16  alpha = 1.d0/137.0359895d0
17 c
18  gam1 = 0.38d0
19  gam2 = 0.38d0
20  fpi = 0.0933d0
21 c coupl = 2.d0*sqrt(3.d0)/fpi**2
22  coupl = sqrt(6.d0)/fpi**2 ! normalization change /sqrt(2)
23  a1m = 1.251d0
24  a1g = 0.599d0
25  rhom = 0.773d0
26  rhog = 0.145d0
27  rho1m = 1.35d0
28  rho1g = 0.3d0
29  rho2m = 1.7d0
30  rho2g = 0.235d0
31  omm = 0.782d0
32  omg = 0.0085d0
33  aa = 0.d0 ! to compare with Finkemeier (no omega)
34  bb1 = 0.08d0
35  bb2 = -0.0075d0
36  f0m = 1.3d0
37  f0g = 0.6d0
38  pim = 0.14d0
39 c
40 c the omega coupling changed
41 c
42  sgo = 1.55d0/sqrt(2.d0)
43 CC sgo = 1.4d0/sqrt(2.d0)
44 c
45  rhom2 = rhom**2
46  rho1m2 = rho1m**2
47  rho2m2 = rho2m**2
48  omm2 = omm**2
49  rhomg = rhom*rhog
50  rho1mg = rho1m*rho1g
51  rho2mg = rho2m*rho2g
52  ommg = omm*omg
53 c
54  a1m2 = a1m**2
55  con = a1g*a1m/gfun8(a1m2)
56 c
57  amrop = 1.7d0
58  gamrop = 0.26d0
59  sig = -0.1d0
60  amrop_2 = amrop**2
61  amropg = amrop*gamrop
62 c
63  beta = -0.145d0
64  rho1m_t = 1.37d0
65  rho1g_t = 0.51d0
66 c
67  rho1m2_t = rho1m_t**2
68  rho1mg_t = rho1m_t*rho1g_t
69 
70  return
71  end
72 c*************************************************************************
73  complex*16 function anom_bwg(q1_2,q2_2)
74  implicit real*8 (a-h,o-z)
75 c
76  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
77  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
78  common /cbwgrho/ rhom2,rho1m2,rho2m2,omm2,rhomg,rho1mg,rho2mg
79  1 ,ommg
80  common /anom/amrop,gamrop,sig,amrop_2,amropg
81 c
82  anom_bwg = (dcmplx(1.d0,0.d0)/dcmplx(rhom2-q1_2,-rhomg)
83  1 + dcmplx(sig,0.d0)/dcmplx(amrop_2-q1_2,-amropg) )
84  2 * dcmplx(1.d0,0.d0)/dcmplx(omm2-q2_2,-ommg)
85  return
86  end
87 c*************************************************************************
88  complex*16 function bwga1(q1_2)
89  implicit real*8 (a-h,o-z)
90 c
91  common /cbwga1/ a1m2,con
92 c
93  ggm = gfun8(q1_2)*con
94  bwga1 = dcmplx(a1m2,0.d0)/dcmplx(a1m2-q1_2,-ggm)
95 c
96  return
97  end
98 c*************************************************************************
99  real*8 function gfun8(q1_2)
100  implicit real*8 (a-h,o-z)
101 c
102  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
103  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
104 c
105  if(q1_2.gt.((rhom+pim)**2))then
106  gfun8 = q1_2*1.623d0 + 10.38d0 - 9.32d0/q1_2 + 0.65d0/q1_2**2
107  else
108  c1 = q1_2 - 9.d0*pim**2
109  gfun8 = 4.1d0 *c1**3 *(1.d0 - 3.3d0*c1 + 5.8d0*c1**2)
110  endif
111 c
112  return
113  end
114 c*************************************************************************
115  complex*16 function bwgrho(q1_2)
116  implicit real*8 (a-h,o-z)
117 c
118  complex*16 cbw,cbw1,cbw2,cbwo
119 c
120  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
121  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
122  common /cbwgrho/ rhom2,rho1m2,rho2m2,omm2,rhomg,rho1mg,rho2mg,ommg
123 c
124  c2 = 4.d0*pim**2/q1_2
125  if(c2.lt.1.d0)then
126 c
127  c1 = rhom2/q1_2
128  gamrho = rhomg*sqrt(c1*((1.d0-c2)/(c1-c2))**3)
129  c1 = rho1m2/q1_2
130  gamrho1 = rho1mg*sqrt(c1*((1.d0-c2)/(c1-c2))**3)
131  c1 = rho2m2/q1_2
132  gamrho2 = rho2mg*sqrt(c1*((1.d0-c2)/(c1-c2))**3)
133  c1 = omm2/q1_2
134  gamom = ommg*sqrt(c1*((1.d0-c2)/(c1-c2))**3)
135  else
136  gamrho =0.d0
137  gamrho1=0.d0
138  gamrho2=0.d0
139  gamom =0.d0
140  endif
141 c
142  cbw = dcmplx(rhom2,0.d0)/dcmplx(rhom2-q1_2,-gamrho)
143  cbw1 = dcmplx(rho1m2,0.d0)/dcmplx(rho1m2-q1_2,-gamrho1)
144  cbw2 = dcmplx(rho2m2,0.d0)/dcmplx(rho2m2-q1_2,-gamrho2)
145  cbwo = dcmplx(omm2,0.d0)/dcmplx(omm2-q1_2,-gamom)
146  bwgrho = ( cbw *(1.d0+aa*cbwo)/(1.d0+aa)
147  1 + bb1*cbw1+bb2*cbw2)/(1.d0+bb1+bb2)
148 c
149  return
150  end
151 c*************************************************************************
152  complex*16 function bwgrho_t(q1_2)
153  implicit real*8 (a-h,o-z)
154 c
155  complex*16 cbw,cbw1,cbw2,cbwo
156 c
157  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
158  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
159  common /cbwgrho/ rhom2,rho1m2,rho2m2,omm2,rhomg,rho1mg,rho2mg
160  1 ,ommg
161  common /cbwgrho_t/ rho1m2_t,rho1mg_t,beta
162 c
163  c2 = 4.d0*pim**2/q1_2
164 c
165  c1 = rhom2/q1_2
166  if(c2.gt.1.d0)then
167  gamrho = 0.d0
168  else
169  gamrho = rhomg*sqrt(c1*((1.d0-c2)/(c1-c2))**3)
170  endif
171  c1 = rho1m2_t/q1_2
172  if(c2.gt.1.d0)then
173  gamrho1 =0
174  else
175  gamrho1 = rho1mg_t*sqrt(c1*((1.d0-c2)/(c1-c2))**3)
176  endif
177 c
178  cbw = dcmplx(rhom2,0.d0)/dcmplx(rhom2-q1_2,-gamrho)
179  cbw1 = dcmplx(rho1m2,0.d0)/dcmplx(rho1m2-q1_2,-gamrho1)
180 
181  bwgrho_t = (cbw+beta*cbw1)/(1.d0+beta)
182 c
183  return
184  end
185 c ************************************************************************
186  complex*16 function bwgf0(q1_2)
187  implicit real*8 (a-h,o-z)
188 c
189  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
190  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
191 c
192  f0m2 = f0m**2
193  f0mg = f0m*f0g
194  bwgf0 = dcmplx(f0m2,-f0mg)/dcmplx(f0m2-q1_2,-f0mg)
195 c
196  return
197  end
198 c***************************************************************************
199 c*************************************************************************
200 c the file contains all currents in 4pi mode
201 c the basic building block is rho(0) -> pi+ pi- 2pi0 mode
202 c other modes: rho(0) -> 2pi+ 2pi-
203 c rho(-) -> 3pi0 pi-
204 c rho(-) -> 2pi- pi+ pi0
205 c
206 c*************************************************************************
207 c this is a code of hadronic current rho(0) -> 2pi+ 2pi-
208 c
209 c q1,q4 : pi+'s four momenta
210 c q2,q3 : pi-'s four momenta
211 c
212  subroutine had1(qq2,q1,q2,q3,q4,hadr)
213  implicit real*8 (a-h,o-z)
214 c
215  complex*16 hadr(4),hadr1(4),hadr2(4),hadr3(4),hadr4(4)
216  dimension q1(4),q2(4),q3(4),q4(4)
217 c
218  call had2(qq2,q1,q2,q3,q4,hadr1)
219  call had2(qq2,q4,q2,q3,q1,hadr2)
220  call had2(qq2,q1,q3,q2,q4,hadr3)
221  call had2(qq2,q4,q3,q2,q1,hadr4)
222 c
223  do i=1,4
224  hadr(i) = hadr1(i)+hadr2(i)+hadr3(i)+hadr4(i)
225  enddo
226 c
227  return
228  end
229 c*************************************************************************
230 c this is a code of hadronic current rho(-) -> 3pi0 pi-
231 c
232 c q1,q2,q3 : pi0's four momenta
233 c q4 : pi-'s four momentum
234 c
235 c
236  subroutine had3(qq2,q1,q2,q3,q4,hadr)
237  implicit real*8 (a-h,o-z)
238 c
239  complex*16 hadr(4),hadr1(4),hadr2(4),hadr3(4)
240  dimension q1(4),q2(4),q3(4),q4(4)
241 c
242  call had2(qq2,q1,q2,q3,q4,hadr1)
243  call had2(qq2,q1,q3,q2,q4,hadr2)
244  call had2(qq2,q3,q2,q1,q4,hadr3)
245 c
246  do i=1,4
247  hadr(i) = (hadr1(i)+hadr2(i)+hadr3(i))*sqrt(2.d0)
248  enddo
249 c
250  return
251  end
252 c*************************************************************************
253 c this is a code of hadronic current rho(-) -> 2pi- pi+ pi0
254 c
255 c q1,q2 : pi-'s four momenta
256 c q3 : pi0 four momentum
257 c q4 : pi+ four momentum
258 c
259 c
260  subroutine had4(qq2,q1,q2,q3,q4,hadr)
261 c
262  implicit real*8 (a-h,o-z)
263 c
264  complex*16 hadr(4),hadr1(4),hadr2(4)
265  dimension q1(4),q2(4),q3(4),q4(4)
266 c
267  call had2(qq2,q3,q1,q2,q4,hadr1)
268  call had2(qq2,q3,q2,q1,q4,hadr2)
269 c
270  do i=1,4
271  hadr(i) = (hadr1(i)+hadr2(i))*sqrt(2.d0)
272  enddo
273 c
274  return
275  end
276 c*************************************************************************
277 c*************************************************************************
278 c this is a code of hadronic current rho(0) -> pi+ pi- 2pi0
279 c
280 c the basic building block for other currents
281 c
282 c q1,q2 : pi0's four momenta
283 c q3 : pi- four momentum
284 c q4 : pi+ four momentum
285 c
286 c the current was obtained in h1_t_f0.f(log)
287 c
288  subroutine had2(qq2,q1,q2,q3,q4,hadr)
289  implicit real*8 (a-h,o-z)
290 c
291  complex*16 hadr(4),cfac(4),tt(4,4,4),ss(4,4,4,4)
292  complex*16 bwga1,bwgrho,bwgrho_t,bwgf0,c0,c5,c6
293  complex*16 c1_t,c2_t,c3_t,c4_t,anom_bwg
294  dimension q1(4),q2(4),q3(4),q4(4),q2m4(4),q3m1(4),q4m1(4),q3m2(4)
295  dimension q123(4),q124(4),qq(4),q3m4(4),q134(4),q234(4)
296  dimension q2p4(4),q1p3(4),q2p3(4),q1p4(4),q1p2(4),q3p4(4)
297 c
298  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
299  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
300 c
301 c the dot products:
302 c
303  do i=1,4
304  q2m4(i) = q2(i)-q4(i)
305  q3m1(i) = q3(i)-q1(i)
306  q3m4(i) = q3(i)-q4(i)
307  q4m1(i) = q4(i)-q1(i)
308  q3m2(i) = q3(i)-q2(i)
309  q2p4(i) = q2(i)+q4(i)
310  q1p3(i) = q1(i)+q3(i)
311  q1p2(i) = q1(i)+q2(i)
312  q2p3(i) = q2(i)+q3(i)
313  q1p4(i) = q1(i)+q4(i)
314  q3p4(i) = q3(i)+q4(i)
315  q123(i) = q2p3(i)+q1(i)
316  q124(i) = q2p4(i)+q1(i)
317  qq(i) = q123(i) + q4(i)
318  enddo
319  q1_2m4 = q1(1)*q2m4(1)-q1(2)*q2m4(2)-q1(3)*q2m4(3)-q1(4)*q2m4(4)
320  q1_3m2 = q1(1)*q3m2(1)-q1(2)*q3m2(2)-q1(3)*q3m2(3)-q1(4)*q3m2(4)
321  q3_2m4 = q3(1)*q2m4(1)-q3(2)*q2m4(2)-q3(3)*q2m4(3)-q3(4)*q2m4(4)
322  q2_3m1 = q2(1)*q3m1(1)-q2(2)*q3m1(2)-q2(3)*q3m1(3)-q2(4)*q3m1(4)
323  q2_4m1 = q2(1)*q4m1(1)-q2(2)*q4m1(2)-q2(3)*q4m1(3)-q2(4)*q4m1(4)
324  q3_4m1 = q3(1)*q4m1(1)-q3(2)*q4m1(2)-q3(3)*q4m1(3)-q3(4)*q4m1(4)
325  q4_3m1 = q4(1)*q3m1(1)-q4(2)*q3m1(2)-q4(3)*q3m1(3)-q4(4)*q3m1(4)
326  q4_3m2 = q4(1)*q3m2(1)-q4(2)*q3m2(2)-q4(3)*q3m2(3)-q4(4)*q3m2(4)
327  qmq3_2 = q124(1)**2 -q124(2)**2 -q124(3)**2 -q124(4)**2
328  qmq4_2 = q123(1)**2 -q123(2)**2 -q123(3)**2 -q123(4)**2
329  q_q3 = qq(1)*q3(1)-qq(2)*q3(2)-qq(3)*q3(3)-qq(4)*q3(4)
330  q_q4 = qq(1)*q4(1)-qq(2)*q4(2)-qq(3)*q4(3)-qq(4)*q4(4)
331  q2p4_2 = q2p4(1)**2 - q2p4(2)**2 - q2p4(3)**2 - q2p4(4)**2
332  q3p4_2 = q3p4(1)**2 - q3p4(2)**2 - q3p4(3)**2 - q3p4(4)**2
333  q1p3_2 = q1p3(1)**2 - q1p3(2)**2 - q1p3(3)**2 - q1p3(4)**2
334  q1p2_2 = q1p2(1)**2 - q1p2(2)**2 - q1p2(3)**2 - q1p2(4)**2
335  q2p3_2 = q2p3(1)**2 - q2p3(2)**2 - q2p3(3)**2 - q2p3(4)**2
336  q1p4_2 = q1p4(1)**2 - q1p4(2)**2 - q1p4(3)**2 - q1p4(4)**2
337  q1p2_3m4 = q1p2(1)*q3m4(1)
338  1 -q1p2(2)*q3m4(2)-q1p2(3)*q3m4(3)-q1p2(4)*q3m4(4)
339  q1p3_2m4 = q1_2m4 + q3_2m4
340  q1p4_3m2 = q1_3m2 + q4_3m2
341  q2p4_3m1 = q2_3m1 + q4_3m1
342  q2p3_4m1 = q2_4m1 + q3_4m1
343 c
344  c0 = bwgrho(qq2)*coupl
345 c c0 = coupl
346 c
347  c1_t = bwgrho_t(q2p4_2)
348  c2_t = bwgrho_t(q1p3_2)
349  c3_t = bwgrho_t(q2p3_2)
350  c4_t = bwgrho_t(q1p4_2)
351 c
352  c5 = bwga1(qmq3_2)
353  c6 = bwga1(qmq4_2)
354 c
355  tt(1,2,4) = c5*c1_t*gam1
356  tt(2,1,4) = c5*c4_t*gam1
357  tt(2,3,1) = c6*c2_t*gam1
358  tt(1,2,3) = c6*c3_t*gam1
359 c
360  ss(3,4,1,2) = bwgrho(q3p4_2)*bwgf0(q1p2_2)*gam2
361 c
362  cfac(1) = tt(1,2,3) * (-1.d0 - q1_3m2/qmq4_2 )
363  1 + tt(1,2,4) * ( 1.d0 - q1_2m4/qmq3_2 )
364  2 + tt(2,1,4) * ( 3.d0 + q2_4m1/qmq3_2 )
365  3 + tt(2,3,1) * (-3.d0 - q2_3m1/qmq4_2 )
366 c
367  cfac(2) = tt(1,2,3) * (-3.d0 - q1_3m2/qmq4_2 )
368  1 + tt(1,2,4) * ( 3.d0 - q1_2m4/qmq3_2 )
369  2 + tt(2,1,4) * ( 1.d0 + q2_4m1/qmq3_2 )
370  3 + tt(2,3,1) * (-1.d0 - q2_3m1/qmq4_2 )
371 c
372  cfac(3) = tt(1,2,3) * ( 1.d0 - q1_3m2/qmq4_2 )
373  1 + tt(1,2,4) * ( 1.d0 + q1_2m4/qmq3_2 )
374  2 + tt(2,1,4) * ( 1.d0 - q2_4m1/qmq3_2 )
375  3 + tt(2,3,1) * ( 1.d0 - q2_3m1/qmq4_2 )
376  4 -3.d0*ss(3,4,1,2)
377 c
378  cfac(4) = tt(1,2,3)
379  1 *(1.d0 -2.d0/qq2*(q_q4*q1_3m2/qmq4_2 +q1p4_3m2) +q1_3m2/qmq4_2 )
380  2 + tt(1,2,4)
381  3 *(-1.d0-2.d0/qq2*(q1_2m4/qmq3_2*q_q3 +q1p3_2m4) +q1_2m4/qmq3_2 )
382  4 + tt(2,1,4)
383  5 *(-1.d0+2.d0/qq2*(q_q3*q2_4m1/qmq3_2 +q2p3_4m1) -q2_4m1/qmq3_2 )
384  6 + tt(2,3,1)
385  7 *(1.d0 -2.d0/qq2*(q2_3m1/qmq4_2*q_q4 +q2p4_3m1) +q2_3m1/qmq4_2 )
386  8 +3.d0*ss(3,4,1,2)/qq2*q1p2_3m4
387 c
388  do i=1,4
389  cfac(i) = cfac(i)*c0
390  enddo
391 c
392  do i=1,4
393  hadr(i) = q1(i) *cfac(1) + q2(i)*cfac(2)
394  1 + q3m4(i)*cfac(3) + qq(i)*cfac(4)
395  enddo
396 c
397 c from here Omega current
398 c
399  fac3 = sgo * 1475.98d0*12.924d0 * 0.266d0 *rhom**2
400 c
401 c the dot products:
402 c
403  do i=1,4
404  q134(i) = q1p3(i)+q4(i)
405  q234(i) = q2p4(i)+q3(i)
406  enddo
407 c
408  q1_134 = q1(1)*q134(1)-q1(2)*q134(2)-q1(3)*q134(3)-q1(4)*q134(4)
409  q3_134 = q3(1)*q134(1)-q3(2)*q134(2)-q3(3)*q134(3)-q3(4)*q134(4)
410  q4_134 = q4(1)*q134(1)-q4(2)*q134(2)-q4(3)*q134(3)-q4(4)*q134(4)
411  q2_234 = q2(1)*q234(1)-q2(2)*q234(2)-q2(3)*q234(3)-q2(4)*q234(4)
412  q3_234 = q3(1)*q234(1)-q3(2)*q234(2)-q3(3)*q234(3)-q3(4)*q234(4)
413  q4_234 = q4(1)*q234(1)-q4(2)*q234(2)-q4(3)*q234(3)-q4(4)*q234(4)
414  q12 = q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) - q1(4)*q2(4)
415  q13 = q1(1)*q3(1) - q1(2)*q3(2) - q1(3)*q3(3) - q1(4)*q3(4)
416  q14 = q1(1)*q4(1) - q1(2)*q4(2) - q1(3)*q4(3) - q1(4)*q4(4)
417  q23 = q2(1)*q3(1) - q2(2)*q3(2) - q2(3)*q3(3) - q2(4)*q3(4)
418  q24 = q2(1)*q4(1) - q2(2)*q4(2) - q2(3)*q4(3) - q2(4)*q4(4)
419  q34 = q3(1)*q4(1) - q3(2)*q4(2) - q3(3)*q4(3) - q3(4)*q4(4)
420  q234_2 = q234(1)**2-q234(2)**2-q234(3)**2-q234(4)**2
421  q134_2 = q134(1)**2-q134(2)**2-q134(3)**2-q134(4)**2
422 c
423  cfac(1) = anom_bwg(qq2,q134_2) *(q3_134*q24 -q4_134*q23)
424  cfac(2) = anom_bwg(qq2,q234_2) *(q3_234*q14 -q4_234*q13)
425  cfac(3) = anom_bwg(qq2,q134_2) *(q4_134*q12 -q1_134*q24)
426  1 + anom_bwg(qq2,q234_2) *(q4_234*q12 -q2_234*q14)
427  cfac(4) = anom_bwg(qq2,q134_2) *(q1_134*q23 -q3_134*q12)
428  1 + anom_bwg(qq2,q234_2) *(q2_234*q13 -q3_234*q12)
429 c
430  do i =1,4
431  hadr(i) = hadr(i) + fac3* (q1(i)*cfac(1) + q2(i)*cfac(2)
432  1 + q3(i)*cfac(3) + q4(i)*cfac(4) )
433  enddo
434 c
435  return
436  end
437 c*************************************************************************
438 c this is a code of hadronic current rho(0) -> pi+ pi- 2pi0
439 c
440 c the basic building block for other currents: omega part
441 c
442 c q1,q2 : pi0's four momenta
443 c q3 : pi- four momentum
444 c q4 : pi+ four momentum
445 c
446 c the current was obtained in h1_t_f0.f(log)
447 c
448  subroutine had2_om(qq2,q1,q2,q3,q4,hadr)
449  implicit real*8 (a-h,o-z)
450 c
451  complex*16 hadr(4),cfac(4),tt(4,4,4),ss(4,4,4,4)
452  complex*16 bwga1,bwgrho,bwgrho_t,bwgf0,c0,c5,c6
453  complex*16 c1_t,c2_t,c3_t,c4_t,anom_bwg
454  dimension q1(4),q2(4),q3(4),q4(4),q2m4(4),q3m1(4),q4m1(4),q3m2(4)
455  dimension q123(4),q124(4),qq(4),q3m4(4),q134(4),q234(4)
456  dimension q2p4(4),q1p3(4),q2p3(4),q1p4(4),q1p2(4),q3p4(4)
457 c
458  common /had_par/ gam1,gam2,coupl,a1m,a1g,rhom,rhog,rho1m,rho1g
459  1 ,rho2m,rho2g,omm,omg,aa,bb1,bb2,f0m,f0g,pim,sgo
460 c
461 c the dot products:
462 c
463  do i=1,4
464  q2m4(i) = q2(i)-q4(i)
465  q3m1(i) = q3(i)-q1(i)
466  q3m4(i) = q3(i)-q4(i)
467  q4m1(i) = q4(i)-q1(i)
468  q3m2(i) = q3(i)-q2(i)
469  q2p4(i) = q2(i)+q4(i)
470  q1p3(i) = q1(i)+q3(i)
471  q1p2(i) = q1(i)+q2(i)
472  q2p3(i) = q2(i)+q3(i)
473  q1p4(i) = q1(i)+q4(i)
474  q3p4(i) = q3(i)+q4(i)
475  q123(i) = q2p3(i)+q1(i)
476  q124(i) = q2p4(i)+q1(i)
477  qq(i) = q123(i) + q4(i)
478  enddo
479  q1_2m4 = q1(1)*q2m4(1)-q1(2)*q2m4(2)-q1(3)*q2m4(3)-q1(4)*q2m4(4)
480  q1_3m2 = q1(1)*q3m2(1)-q1(2)*q3m2(2)-q1(3)*q3m2(3)-q1(4)*q3m2(4)
481  q3_2m4 = q3(1)*q2m4(1)-q3(2)*q2m4(2)-q3(3)*q2m4(3)-q3(4)*q2m4(4)
482  q2_3m1 = q2(1)*q3m1(1)-q2(2)*q3m1(2)-q2(3)*q3m1(3)-q2(4)*q3m1(4)
483  q2_4m1 = q2(1)*q4m1(1)-q2(2)*q4m1(2)-q2(3)*q4m1(3)-q2(4)*q4m1(4)
484  q3_4m1 = q3(1)*q4m1(1)-q3(2)*q4m1(2)-q3(3)*q4m1(3)-q3(4)*q4m1(4)
485  q4_3m1 = q4(1)*q3m1(1)-q4(2)*q3m1(2)-q4(3)*q3m1(3)-q4(4)*q3m1(4)
486  q4_3m2 = q4(1)*q3m2(1)-q4(2)*q3m2(2)-q4(3)*q3m2(3)-q4(4)*q3m2(4)
487  qmq3_2 = q124(1)**2 -q124(2)**2 -q124(3)**2 -q124(4)**2
488  qmq4_2 = q123(1)**2 -q123(2)**2 -q123(3)**2 -q123(4)**2
489  q_q3 = qq(1)*q3(1)-qq(2)*q3(2)-qq(3)*q3(3)-qq(4)*q3(4)
490  q_q4 = qq(1)*q4(1)-qq(2)*q4(2)-qq(3)*q4(3)-qq(4)*q4(4)
491  q2p4_2 = q2p4(1)**2 - q2p4(2)**2 - q2p4(3)**2 - q2p4(4)**2
492  q3p4_2 = q3p4(1)**2 - q3p4(2)**2 - q3p4(3)**2 - q3p4(4)**2
493  q1p3_2 = q1p3(1)**2 - q1p3(2)**2 - q1p3(3)**2 - q1p3(4)**2
494  q1p2_2 = q1p2(1)**2 - q1p2(2)**2 - q1p2(3)**2 - q1p2(4)**2
495  q2p3_2 = q2p3(1)**2 - q2p3(2)**2 - q2p3(3)**2 - q2p3(4)**2
496  q1p4_2 = q1p4(1)**2 - q1p4(2)**2 - q1p4(3)**2 - q1p4(4)**2
497  q1p2_3m4 = q1p2(1)*q3m4(1)
498  1 -q1p2(2)*q3m4(2)-q1p2(3)*q3m4(3)-q1p2(4)*q3m4(4)
499  q1p3_2m4 = q1_2m4 + q3_2m4
500  q1p4_3m2 = q1_3m2 + q4_3m2
501  q2p4_3m1 = q2_3m1 + q4_3m1
502  q2p3_4m1 = q2_4m1 + q3_4m1
503 c
504 c
505 c from here Omega current
506 c
507  fac3 = sgo * 1475.98d0*12.924d0 * 0.266d0 *rhom**2
508 c
509 c the dot products:
510 c
511  do i=1,4
512  q134(i) = q1p3(i)+q4(i)
513  q234(i) = q2p4(i)+q3(i)
514  enddo
515 c
516  q1_134 = q1(1)*q134(1)-q1(2)*q134(2)-q1(3)*q134(3)-q1(4)*q134(4)
517  q3_134 = q3(1)*q134(1)-q3(2)*q134(2)-q3(3)*q134(3)-q3(4)*q134(4)
518  q4_134 = q4(1)*q134(1)-q4(2)*q134(2)-q4(3)*q134(3)-q4(4)*q134(4)
519  q2_234 = q2(1)*q234(1)-q2(2)*q234(2)-q2(3)*q234(3)-q2(4)*q234(4)
520  q3_234 = q3(1)*q234(1)-q3(2)*q234(2)-q3(3)*q234(3)-q3(4)*q234(4)
521  q4_234 = q4(1)*q234(1)-q4(2)*q234(2)-q4(3)*q234(3)-q4(4)*q234(4)
522  q12 = q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) - q1(4)*q2(4)
523  q13 = q1(1)*q3(1) - q1(2)*q3(2) - q1(3)*q3(3) - q1(4)*q3(4)
524  q14 = q1(1)*q4(1) - q1(2)*q4(2) - q1(3)*q4(3) - q1(4)*q4(4)
525  q23 = q2(1)*q3(1) - q2(2)*q3(2) - q2(3)*q3(3) - q2(4)*q3(4)
526  q24 = q2(1)*q4(1) - q2(2)*q4(2) - q2(3)*q4(3) - q2(4)*q4(4)
527  q34 = q3(1)*q4(1) - q3(2)*q4(2) - q3(3)*q4(3) - q3(4)*q4(4)
528  q234_2 = q234(1)**2-q234(2)**2-q234(3)**2-q234(4)**2
529  q134_2 = q134(1)**2-q134(2)**2-q134(3)**2-q134(4)**2
530 c
531  cfac(1) = anom_bwg(qq2,q134_2) *(q3_134*q24 -q4_134*q23)
532  cfac(2) = anom_bwg(qq2,q234_2) *(q3_234*q14 -q4_234*q13)
533  cfac(3) = anom_bwg(qq2,q134_2) *(q4_134*q12 -q1_134*q24)
534  1 + anom_bwg(qq2,q234_2) *(q4_234*q12 -q2_234*q14)
535  cfac(4) = anom_bwg(qq2,q134_2) *(q1_134*q23 -q3_134*q12)
536  1 + anom_bwg(qq2,q234_2) *(q2_234*q13 -q3_234*q12)
537 c
538  do i =1,4
539  hadr(i) = fac3* (q1(i)*cfac(1) + q2(i)*cfac(2)
540  1 + q3(i)*cfac(3) + q4(i)*cfac(4) )
541  enddo
542 c
543  return
544  end
545 c************************************************************************
546 
547 
548 
549