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)
4c
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
14c
15 pi = 3.141592653589793238d0
16 alpha = 1.d0/137.0359895d0
17c
18 gam1 = 0.38d0
19 gam2 = 0.38d0
20 fpi = 0.0933d0
21c 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
39c
40c the omega coupling changed
41c
42 sgo = 1.55d0/sqrt(2.d0)
43CC sgo = 1.4d0/sqrt(2.d0)
44c
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
53c
54 a1m2 = a1m**2
55 con = a1g*a1m/gfun8(a1m2)
56c
57 amrop = 1.7d0
58 gamrop = 0.26d0
59 sig = -0.1d0
60 amrop_2 = amrop**2
61 amropg = amrop*gamrop
62c
63 beta = -0.145d0
64 rho1m_t = 1.37d0
65 rho1g_t = 0.51d0
66c
67 rho1m2_t = rho1m_t**2
68 rho1mg_t = rho1m_t*rho1g_t
69
70 return
71 end
72c*************************************************************************
73 complex*16 function anom_bwg(q1_2,q2_2)
74 implicit real*8 (a-h,o-z)
75c
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
81c
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
87c*************************************************************************
88 complex*16 function bwga1(q1_2)
89 implicit real*8 (a-h,o-z)
90c
91 common /cbwga1/ a1m2,con
92c
93 ggm = gfun8(q1_2)*con
94 bwga1 = dcmplx(a1m2,0.d0)/dcmplx(a1m2-q1_2,-ggm)
95c
96 return
97 end
98c*************************************************************************
99 real*8 function gfun8(q1_2)
100 implicit real*8 (a-h,o-z)
101c
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
104c
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
111c
112 return
113 end
114c*************************************************************************
115 complex*16 function bwgrho(q1_2)
116 implicit real*8 (a-h,o-z)
117c
118 complex*16 cbw,cbw1,cbw2,cbwo
119c
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
123c
124 c2 = 4.d0*pim**2/q1_2
125 if(c2.lt.1.d0)then
126c
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
141c
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)
148c
149 return
150 end
151c*************************************************************************
152 complex*16 function bwgrho_t(q1_2)
153 implicit real*8 (a-h,o-z)
154c
155 complex*16 cbw,cbw1,cbw2,cbwo
156c
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
162c
163 c2 = 4.d0*pim**2/q1_2
164c
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
177c
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)
182c
183 return
184 end
185c ************************************************************************
186 complex*16 function bwgf0(q1_2)
187 implicit real*8 (a-h,o-z)
188c
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
191c
192 f0m2 = f0m**2
193 f0mg = f0m*f0g
194 bwgf0 = dcmplx(f0m2,-f0mg)/dcmplx(f0m2-q1_2,-f0mg)
195c
196 return
197 end
198c***************************************************************************
199c*************************************************************************
200c the file contains all currents in 4pi mode
201c the basic building block is rho(0) -> pi+ pi- 2pi0 mode
202c other modes: rho(0) -> 2pi+ 2pi-
203c rho(-) -> 3pi0 pi-
204c rho(-) -> 2pi- pi+ pi0
205c
206c*************************************************************************
207c this is a code of hadronic current rho(0) -> 2pi+ 2pi-
208c
209c q1,q4 : pi+'s four momenta
210c q2,q3 : pi-'s four momenta
211c
212 subroutine had1(qq2,q1,q2,q3,q4,hadr)
213 implicit real*8 (a-h,o-z)
214c
215 complex*16 hadr(4),hadr1(4),hadr2(4),hadr3(4),hadr4(4)
216 dimension q1(4),q2(4),q3(4),q4(4)
217c
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)
222c
223 do i=1,4
224 hadr(i) = hadr1(i)+hadr2(i)+hadr3(i)+hadr4(i)
225 enddo
226c
227 return
228 end
229c*************************************************************************
230c this is a code of hadronic current rho(-) -> 3pi0 pi-
231c
232c q1,q2,q3 : pi0's four momenta
233c q4 : pi-'s four momentum
234c
235c
236 subroutine had3(qq2,q1,q2,q3,q4,hadr)
237 implicit real*8 (a-h,o-z)
238c
239 complex*16 hadr(4),hadr1(4),hadr2(4),hadr3(4)
240 dimension q1(4),q2(4),q3(4),q4(4)
241c
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)
245c
246 do i=1,4
247 hadr(i) = (hadr1(i)+hadr2(i)+hadr3(i))*sqrt(2.d0)
248 enddo
249c
250 return
251 end
252c*************************************************************************
253c this is a code of hadronic current rho(-) -> 2pi- pi+ pi0
254c
255c q1,q2 : pi-'s four momenta
256c q3 : pi0 four momentum
257c q4 : pi+ four momentum
258c
259c
260 subroutine had4(qq2,q1,q2,q3,q4,hadr)
261c
262 implicit real*8 (a-h,o-z)
263c
264 complex*16 hadr(4),hadr1(4),hadr2(4)
265 dimension q1(4),q2(4),q3(4),q4(4)
266c
267 call had2(qq2,q3,q1,q2,q4,hadr1)
268 call had2(qq2,q3,q2,q1,q4,hadr2)
269c
270 do i=1,4
271 hadr(i) = (hadr1(i)+hadr2(i))*sqrt(2.d0)
272 enddo
273c
274 return
275 end
276c*************************************************************************
277c*************************************************************************
278c this is a code of hadronic current rho(0) -> pi+ pi- 2pi0
279c
280c the basic building block for other currents
281c
282c q1,q2 : pi0's four momenta
283c q3 : pi- four momentum
284c q4 : pi+ four momentum
285c
286c the current was obtained in h1_t_f0.f(log)
287c
288 subroutine had2(qq2,q1,q2,q3,q4,hadr)
289 implicit real*8 (a-h,o-z)
290c
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)
297c
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
300c
301c the dot products:
302c
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
343c
344 c0 = bwgrho(qq2)*coupl
345c c0 = coupl
346c
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)
351c
352 c5 = bwga1(qmq3_2)
353 c6 = bwga1(qmq4_2)
354c
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
359c
360 ss(3,4,1,2) = bwgrho(q3p4_2)*bwgf0(q1p2_2)*gam2
361c
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 )
366c
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 )
371c
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)
377c
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
387c
388 do i=1,4
389 cfac(i) = cfac(i)*c0
390 enddo
391c
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
396c
397c from here Omega current
398c
399 fac3 = sgo * 1475.98d0*12.924d0 * 0.266d0 *rhom**2
400c
401c the dot products:
402c
403 do i=1,4
404 q134(i) = q1p3(i)+q4(i)
405 q234(i) = q2p4(i)+q3(i)
406 enddo
407c
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
422c
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)
429c
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
434c
435 return
436 end
437c*************************************************************************
438c this is a code of hadronic current rho(0) -> pi+ pi- 2pi0
439c
440c the basic building block for other currents: omega part
441c
442c q1,q2 : pi0's four momenta
443c q3 : pi- four momentum
444c q4 : pi+ four momentum
445c
446c the current was obtained in h1_t_f0.f(log)
447c
448 subroutine had2_om(qq2,q1,q2,q3,q4,hadr)
449 implicit real*8 (a-h,o-z)
450c
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)
457c
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
460c
461c the dot products:
462c
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
503c
504c
505c from here Omega current
506c
507 fac3 = sgo * 1475.98d0*12.924d0 * 0.266d0 *rhom**2
508c
509c the dot products:
510c
511 do i=1,4
512 q134(i) = q1p3(i)+q4(i)
513 q234(i) = q2p4(i)+q3(i)
514 enddo
515c
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
530c
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)
537c
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
542c
543 return
544 end
545c************************************************************************
546
547
548
549