C++InterfacetoTauola
F/prod/formf.f
1 /* copyright(c) 1991-2012 free software foundation, inc.
2  this file is part of the gnu c library.
3 
4  the gnu c library is free software; you can redistribute it and/or
5  modify it under the terms of the gnu lesser general Public
6  license as published by the free software foundation; either
7  version 2.1 of the license, or(at your option) any later version.
8 
9  the gnu c library is distributed in the hope that it will be useful,
10  but without any warranty; without even the implied warranty of
11  merchantability or fitness for a particular purpose. see the gnu
12  lesser general Public license for more details.
13 
14  you should have received a copy of the gnu lesser general Public
15  license along with the gnu c library; if not, see
16  <http://www.gnu.org/licenses/>. */
17 
18 
19 /* this header is separate from features.h so that the compiler can
20  include it implicitly at the start of every compilation. it must
21  not itself include <features.h> or any other header that includes
22  <features.h> because the implicit include comes before any feature
23  test macros that may be defined in a source file before it first
24  explicitly includes a system header. gcc knows the name of this
25  header in order to preinclude it. */
26 
27 /* we do support the iec 559 math functionality, real and complex. */
28 
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
30  unicode 6.0. */
31 
32 /* we do not support c11 <threads.h>. */
33 
34  FUNCTION formom(XMAA,XMOM)
35 c ==================================================================
36 c formfactorfor pi-pi0 gamma final state
37 c r. decker, z. phys c36(1987) 487.
38 c ==================================================================
39  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
40  * ,ampiz,ampi,amro,gamro,ama1,gama1
41  * ,amk,amkz,amkst,gamkst
42 c
43  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
44  * ,ampiz,ampi,amro,gamro,ama1,gama1
45  * ,amk,amkz,amkst,gamkst
46  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
47  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
48  COMMON /testa1/ keya1
49  COMPLEX bwign,formom
50  DATA icont /1/
51 * this inline funct. calculates the scalar part of the propagator
52  bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
53 * hadron current
54  fro =0.266*amro**2
55  elpha=- 0.1
56  amrop = 1.7
57  gamrop= 0.26
58  amom =0.782
59  gamom =0.0085
60  aromeg= 1.0
61  gcoup=12.924
62  gcoup=gcoup*aromeg
63  fqed =sqrt(4.0*3.1415926535/137.03604)
64  formom=fqed*fro**2/sqrt(2.0)*gcoup**2*bwign(xmom,amom,gamom)
65  $ *(bwign(xmaa,amro,gamro)+elpha*bwign(xmaa,amrop,gamrop))
66  $ *(bwign( 0.0,amro,gamro)+elpha*bwign( 0.0,amrop,gamrop))
67  END
68 c=======================================================================
69  COMPLEX FUNCTION fk1ab(XMSQ,INDX)
70 c ==================================================================
71 c complex form-factor for a1+a1prime. ajw 1/98
72 c ==================================================================
73 
74  COMPLEX f1,f2,ampa,ampb
75  INTEGER ifirst,indx
76  DATA ifirst/0/
77 
78  IF (ifirst.EQ.0) THEN
79  ifirst = 1
80  xm1 = pkorb(1,19)
81  xg1 = pkorb(2,19)
82  xm2 = pkorb(1,20)
83  xg2 = pkorb(2,20)
84 
85  xm1sq = xm1*xm1
86  gf1 = gfun(xm1sq)
87  gg1 = xm1*xg1/gf1
88  xm2sq = xm2*xm2
89  gf2 = gfun(xm2sq)
90  gg2 = xm2*xg2/gf2
91  END IF
92 
93  IF (indx.EQ.1) THEN
94  ampa = cmplx(pkorb(3,81),0.)
95  ampb = cmplx(pkorb(3,82),0.)
96  ELSE IF (indx.EQ.2) THEN
97  ampa = cmplx(pkorb(3,83),0.)
98  ampb = cmplx(pkorb(3,84),0.)
99  ELSEIF (indx.EQ.3) THEN
100  ampa = cmplx(pkorb(3,85),0.)
101  ampb = cmplx(pkorb(3,86),0.)
102  ELSEIF (indx.EQ.4) THEN
103  ampa = cmplx(pkorb(3,87),0.)
104  ampb = cmplx(pkorb(3,88),0.)
105  END IF
106 
107  gf = gfun(xmsq)
108  fg1 = gg1*gf
109  fg2 = gg2*gf
110  f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
111  f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
112  fk1ab = ampa*f1+ampb*f2
113 
114  RETURN
115  END
116  FUNCTION form1(MNUM,QQ,S1,SDWA)
117 c ==================================================================
118 c formfactorfor f1 for 3 scalar final state
119 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
120 c h. georgi, weak interactions and modern particle theory,
121 c the benjamin/cummings pub. co., inc. 1984.
122 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
123 c and erratum !!!!!!
124 c ==================================================================
125 c
126  COMPLEX form1,wigner,wigfor,fpikm,bwigm
127  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
128  * ,ampiz,ampi,amro,gamro,ama1,gama1
129  * ,amk,amkz,amkst,gamkst
130 c
131  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
132  * ,ampiz,ampi,amro,gamro,ama1,gama1
133  * ,amk,amkz,amkst,gamkst
134  COMPLEX forma1,formk1,formro,formks
135  COMPLEX fa1a1p,fk1ab,f3pi
136 c
137  IF (mnum.EQ.0) THEN
138 c ------------ 3 pi hadronic state(a1)
139 c formro = fpikm(sqrt(s1),ampi,ampi)
140 c formro = f3pi(1,qq,s1,sdwa)
141 c forma1 = fa1a1p(qq)
142 c form1 = forma1*formro
143  form1 = f3pi(1,qq,s1,sdwa)
144 
145  ELSEIF (mnum.EQ.1) THEN
146 c ------------ k- pi- k+ (k*0 k-)
147  formks = bwigm(s1,amkst,gamkst,ampi,amk)
148  forma1 = fa1a1p(qq)
149  form1 = forma1*formks
150 
151  ELSEIF (mnum.EQ.2) THEN
152 c ------------ k0 pi- k0b(k*- k0)
153  formks = bwigm(s1,amkst,gamkst,ampi,amk)
154  forma1 = fa1a1p(qq)
155  form1 = forma1*formks
156 
157  ELSEIF (mnum.EQ.3) THEN
158 c ------------ k- pi0 k0(k*0 k-)
159  formks = bwigm(s1,amkst,gamkst,ampi,amk)
160  forma1 = fa1a1p(qq)
161  form1 = forma1*formks
162 
163  ELSEIF (mnum.EQ.4) THEN
164 c ------------ pi0 pi0 k- (k*-pi0)
165  formks = bwigm(s1,amkst,gamkst,ampi,amk)
166  formk1 = fk1ab(qq,3)
167  form1 = formk1*formks
168 
169  ELSEIF (mnum.EQ.5) THEN
170 c ------------ k- pi- pi+ (rho0 k-)
171  formk1 = fk1ab(qq,4)
172  formro = fpikm(sqrt(s1),ampi,ampi)
173  form1 = formk1*formro
174 
175  ELSEIF (mnum.EQ.6) THEN
176 c ------------ pi- k0b pi0(pi- k*0b)
177  formk1 = fk1ab(qq,1)
178  formks = bwigm(s1,amkst,gamkst,amk,ampi)
179  form1 = formk1*formks
180 
181  ELSEIF (mnum.EQ.7) THEN
182 c -------------- eta pi- pi0 final state
183  form1=0.0
184  ENDIF
185  END
186  FUNCTION form2(MNUM,QQ,S1,SDWA)
187 c ==================================================================
188 c formfactorfor f2 for 3 scalar final state
189 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
190 c h. georgi, weak interactions and modern particle theory,
191 c the benjamin/cummings pub. co., inc. 1984.
192 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
193 c and erratum !!!!!!
194 c ==================================================================
195 c
196  COMPLEX form2,wigner,wigfor,fpikm,bwigm
197  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
198  * ,ampiz,ampi,amro,gamro,ama1,gama1
199  * ,amk,amkz,amkst,gamkst
200 c
201  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
202  * ,ampiz,ampi,amro,gamro,ama1,gama1
203  * ,amk,amkz,amkst,gamkst
204  COMPLEX forma1,formk1,formro,formks
205  COMPLEX fa1a1p,fk1ab,f3pi
206 
207  IF (mnum.EQ.0) THEN
208 c ------------ 3 pi hadronic state(a1)
209 c formro = fpikm(sqrt(s1),ampi,ampi)
210 c formro = f3pi(2,qq,s1,sdwa)
211 c forma1 = fa1a1p(qq)
212 c form2 = forma1*formro
213  form2 = f3pi(2,qq,s1,sdwa)
214 
215  ELSEIF (mnum.EQ.1) THEN
216 c ------------ k- pi- k+ (rho0 pi-)
217  formro = fpikm(sqrt(s1),amk,amk)
218  forma1 = fa1a1p(qq)
219  form2 = forma1*formro
220 
221  ELSEIF (mnum.EQ.2) THEN
222 c ------------ k0 pi- k0b(rho0 pi-)
223  formro = fpikm(sqrt(s1),amk,amk)
224  forma1 = fa1a1p(qq)
225  form2 = forma1*formro
226 
227  ELSEIF (mnum.EQ.3) THEN
228 c ------------ k- pi0 k0(rho- pi0)
229  formro = fpikm(sqrt(s1),amk,amk)
230  forma1 = fa1a1p(qq)
231  form2 = forma1*formro
232 
233  ELSEIF (mnum.EQ.4) THEN
234 c ------------ pi0 pi0 k- (k*-pi0)
235  formks = bwigm(s1,amkst,gamkst,ampi,amk)
236  formk1 = fk1ab(qq,3)
237  form2 = formk1*formks
238 
239  ELSEIF (mnum.EQ.5) THEN
240 c ------------ k- pi- pi+ (k*0b pi-)
241  formks = bwigm(s1,amkst,gamkst,ampi,amk)
242  formk1 = fk1ab(qq,1)
243  form2 = formk1*formks
244 c
245  ELSEIF (mnum.EQ.6) THEN
246 c ------------ pi- k0b pi0(rho- k0b)
247  formro = fpikm(sqrt(s1),ampi,ampi)
248  formk1 = fk1ab(qq,2)
249  form2 = formk1*formro
250 c
251  ELSEIF (mnum.EQ.7) THEN
252 c -------------- eta pi- pi0 final state
253  form2=0.0
254  ENDIF
255 c
256  END
257  COMPLEX FUNCTION bwigm(S,M,G,XM1,XM2)
258 c **********************************************************
259 c p-wave breit-wigner for rho
260 c **********************************************************
261  REAL s,m,g,xm1,xm2
262  REAL pi,qs,qm,w,gs
263  DATA init /0/
264 c ------------ parameters --------------------
265  IF (init.EQ.0) THEN
266  init=1
267  pi=3.141592654
268 c ------- breit-wigner -----------------------
269  ENDIF
270  IF (s.GT.(xm1+xm2)**2) THEN
271  qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
272  qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
273  w=sqrt(s)
274  gs=g*(m/w)**2*(qs/qm)**3
275  ELSE
276  gs=0.0
277  ENDIF
278  bwigm=m**2/cmplx(m**2-s,-sqrt(s)*gs)
279  RETURN
280  END
281  COMPLEX FUNCTION fpikm(W,XM1,XM2)
282 c **********************************************************
283 c pion form factor
284 c **********************************************************
285  COMPLEX bwigm
286  REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
287  EXTERNAL bwig
288  DATA init /0/
289 c
290 c ------------ parameters --------------------
291  IF (init.EQ.0 ) THEN
292  init=1
293  pi=3.141592654
294  pim=.140
295  rom=0.773
296  rog=0.145
297  rom1=1.370
298  rog1=0.510
299  beta1=-0.145
300  ENDIF
301 c -----------------------------------------------
302  s=w**2
303  fpikm=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
304  & /(1+beta1)
305  RETURN
306  END
307  COMPLEX FUNCTION fpikmd(W,XM1,XM2)
308 c **********************************************************
309 c pion form factor
310 c **********************************************************
311  COMPLEX bwigm
312  REAL rom,rog,rom1,rog1,pi,pim,s,w
313  EXTERNAL bwig
314  DATA init /0/
315 c
316 c ------------ parameters --------------------
317  IF (init.EQ.0 ) THEN
318  init=1
319  pi=3.141592654
320  pim=.140
321  rom=0.773
322  rog=0.145
323  rom1=1.500
324  rog1=0.220
325  rom2=1.750
326  rog2=0.120
327  beta=6.5
328  delta=-26.0
329  ENDIF
330 c -----------------------------------------------
331  s=w**2
332  fpikmd=(delta*bwigm(s,rom,rog,xm1,xm2)
333  $ +beta*bwigm(s,rom1,rog1,xm1,xm2)
334  $ + bwigm(s,rom2,rog2,xm1,xm2))
335  & /(1+beta+delta)
336  RETURN
337  END
338 
339  FUNCTION form3(MNUM,QQ,S1,SDWA)
340 c ==================================================================
341 c formfactorfor f3 for 3 scalar final state
342 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
343 c h. georgi, weak interactions and modern particle theory,
344 c the benjamin/cummings pub. co., inc. 1984.
345 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
346 c and erratum !!!!!!
347 c ==================================================================
348 c
349  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
350  * ,ampiz,ampi,amro,gamro,ama1,gama1
351  * ,amk,amkz,amkst,gamkst
352 c
353  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
354  * ,ampiz,ampi,amro,gamro,ama1,gama1
355  * ,amk,amkz,amkst,gamkst
356  COMPLEX form3,bwigm
357  COMPLEX forma1,formk1,formro,formks
358  COMPLEX fa1a1p,fk1ab,f3pi
359 c
360  IF (mnum.EQ.0) THEN
361 c ------------ 3 pi hadronic state(a1)
362 c formro = fpikm(sqrt(s1),ampi,ampi)
363 c formro = f3pi(3,qq,s1,sdwa)
364 c forma1 = fa1a1p(qq)
365 c form3 = forma1*formro
366  form3 = f3pi(3,qq,s1,sdwa)
367 
368  ELSEIF (mnum.EQ.3) THEN
369 c ------------ k- pi0 k0(k*- k0)
370  formks = bwigm(s1,amkst,gamkst,ampiz,amk)
371  forma1 = fa1a1p(qq)
372  form3 = forma1*formks
373 
374  ELSEIF (mnum.EQ.6) THEN
375 c ------------ pi- k0b pi0(k*- pi0)
376  formks = bwigm(s1,amkst,gamkst,amk,ampi)
377  formk1 = fk1ab(qq,3)
378  form3 = formk1*formks
379 
380  ELSE
381  form3=cmplx(0.,0.)
382  ENDIF
383  END
384  FUNCTION form4(MNUM,QQ,S1,S2,S3)
385 c ==================================================================
386 c formfactorfor f4 for 3 scalar final state
387 c r. decker, in preparation
388 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
389 c and erratum !!!!!!
390 c ==================================================================
391 c
392  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
393  * ,ampiz,ampi,amro,gamro,ama1,gama1
394  * ,amk,amkz,amkst,gamkst
395 c
396  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
397  * ,ampiz,ampi,amro,gamro,ama1,gama1
398  * ,amk,amkz,amkst,gamkst
399  COMPLEX form4,wigner,fpikm
400  REAL*4 m
401 c ---- this formfactor is switched off .. .
402  form4=cmplx(0.0,0.0)
403  END
404  FUNCTION form5(MNUM,QQ,S1,S2)
405 c ==================================================================
406 c formfactorfor f5 for 3 scalar final state
407 c g. kramer, w. palmer, s. pinsky, phys. rev. d30(1984) 89.
408 c g. kramer, w. palmer z. phys. c25(1984) 195.
409 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
410 c and erratum !!!!!!
411 c ==================================================================
412 c
413  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
414  * ,ampiz,ampi,amro,gamro,ama1,gama1
415  * ,amk,amkz,amkst,gamkst
416 c
417  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
418  * ,ampiz,ampi,amro,gamro,ama1,gama1
419  * ,amk,amkz,amkst,gamkst
420  COMPLEX form5,wigner,fpikm,fpikmd,bwigm
421  IF (mnum.EQ.0) THEN
422 c ------------ 3 pi hadronic state(a1)
423  form5=0.0
424  ELSEIF (mnum.EQ.1) THEN
425 c ------------ k- pi- k+
426  elpha=-0.2
427  form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
428  $ *( fpikm(sqrt(s2),ampi,ampi)
429  $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
430  ELSEIF (mnum.EQ.2) THEN
431 c ------------ k0 pi- k0b
432  elpha=-0.2
433  form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
434  $ *( fpikm(sqrt(s2),ampi,ampi)
435  $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
436  ELSEIF (mnum.EQ.3) THEN
437 c ------------ k- k0 pi0
438  form5=0.0
439  ELSEIF (mnum.EQ.4) THEN
440 c ------------ pi0 pi0 k-
441  form5=0.0
442  ELSEIF (mnum.EQ.5) THEN
443 c ------------ k- pi- pi+
444  elpha=-0.2
445  form5=bwigm(qq,amkst,gamkst,ampi,amk)/(1+elpha)
446  $ *( fpikm(sqrt(s1),ampi,ampi)
447  $ +elpha*bwigm(s2,amkst,gamkst,ampi,amk))
448  ELSEIF (mnum.EQ.6) THEN
449 c ------------ pi- k0b pi0
450  elpha=-0.2
451  form5=bwigm(qq,amkst,gamkst,ampi,amkz)/(1+elpha)
452  $ *( fpikm(sqrt(s2),ampi,ampi)
453  $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
454  ELSEIF (mnum.EQ.7) THEN
455 c -------------- eta pi- pi0 final state
456  form5=fpikmd(sqrt(qq),ampi,ampi)*fpikm(sqrt(s1),ampi,ampi)
457  ENDIF
458 c
459  END
460  SUBROUTINE currx(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
461 c ==================================================================
462 c hadronic current for 4 pi final state
463 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
464 c r. decker z. phys c36(1987) 487.
465 c m. gell-mann, d. sharp, w. wagner phys. rev. lett 8 (1962) 261.
466 c ==================================================================
467 
468  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
469  * ,ampiz,ampi,amro,gamro,ama1,gama1
470  * ,amk,amkz,amkst,gamkst
471 c
472  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
473  * ,ampiz,ampi,amro,gamro,ama1,gama1
474  * ,amk,amkz,amkst,gamkst
475  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
476  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
477 c arbitrary fixing of the four pi x-section normalization
478  COMMON /arbit/ arflat,aromeg
479  REAL pim1(4),pim2(4),pim3(4),pim4(4),paa(4)
480  COMPLEX hadcur(4),form1,form2,form3,fpikm
481  COMPLEX bwign
482  REAL pa(4),pb(4)
483  REAL aa(4,4),pp(4,4)
484  DATA pi /3.141592653589793238462643/
485  DATA fpi /93.3e-3/
486  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
487 c
488 c --- masses and constants
489  g1=12.924
490  g2=1475.98
491  g =g1*g2
492  elpha=-.1
493  amrop=1.7
494  gamrop=0.26
495  amom=.782
496  gamom=0.0085
497  arflat=1.0
498  aromeg=1.0
499 c
500  fro=0.266*amro**2
501  coef1=2.0*sqrt(3.0)/fpi**2*arflat
502  coef2=fro*g*aromeg
503 c --- initialization of four vectors
504  DO 7 k=1,4
505  DO 8 l=1,4
506  8 aa(k,l)=0.0
507  hadcur(k)=cmplx(0.0)
508  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
509  pp(1,k)=pim1(k)
510  pp(2,k)=pim2(k)
511  pp(3,k)=pim3(k)
512  7 pp(4,k)=pim4(k)
513 c
514  IF (mnum.EQ.1) THEN
515 c ===================================================================
516 c pi- pi- p0 pi+ case ====
517 c ===================================================================
518  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
519 c --- loop over thre contribution of the non-omega current
520  DO 201 k=1,3
521  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
522  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
523 c -- definition of aa matrix
524 c -- cronecker delta
525  DO 202 i=1,4
526  DO 203 j=1,4
527  203 aa(i,j)=0.0
528  202 aa(i,i)=1.0
529 c ... and the rest ...
530  DO 204 l=1,3
531  IF (l.NE.k) THEN
532  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
533  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
534  DO 205 i=1,4
535  DO 205 j=1,4
536  sig= 1.0
537  IF(j.NE.4) sig=-sig
538  aa(i,j)=aa(i,j)
539  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
540  205 CONTINUE
541  ENDIF
542  204 CONTINUE
543 c --- lets add something to hadcurr
544  form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
545 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikmd(sqrt(qq),ampi,ampi)
546 ccccccccccccccccc form1=wigfor(sk,amro,gamro) (tests)
547 c
548  fix=1.0
549  IF (k.EQ.3) fix=-2.0
550  DO 206 i=1,4
551  DO 206 j=1,4
552  hadcur(i)=
553  $ hadcur(i)+cmplx(fix*coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
554  206 CONTINUE
555 c --- end of the non omega current(3 possibilities)
556  201 CONTINUE
557 c
558 c
559 c --- there are two possibilities for omega current
560 c --- pa pb are corresponding first and second pi-s
561  DO 301 kk=1,2
562  DO 302 i=1,4
563  pa(i)=pp(kk,i)
564  pb(i)=pp(3-kk,i)
565  302 CONTINUE
566 c --- lorentz invariants
567  qqa=0.0
568  ss23=0.0
569  ss24=0.0
570  ss34=0.0
571  qp1p2=0.0
572  qp1p3=0.0
573  qp1p4=0.0
574  p1p2 =0.0
575  p1p3 =0.0
576  p1p4 =0.0
577  DO 303 k=1,4
578  sign=-1.0
579  IF (k.EQ.4) sign= 1.0
580  qqa=qqa+sign*(paa(k)-pa(k))**2
581  ss23=ss23+sign*(pb(k) +pim3(k))**2
582  ss24=ss24+sign*(pb(k) +pim4(k))**2
583  ss34=ss34+sign*(pim3(k)+pim4(k))**2
584  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
585  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
586  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
587  p1p2=p1p2+sign*pa(k)*pb(k)
588  p1p3=p1p3+sign*pa(k)*pim3(k)
589  p1p4=p1p4+sign*pa(k)*pim4(k)
590  303 CONTINUE
591 c
592  form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
593 c form3=bwign(qqa,amom,gamom)*(bwign(ss23,amro,gamro)+
594 c $ bwign(ss24,amro,gamro)+bwign(ss34,amro,gamro))
595  form3=bwign(qqa,amom,gamom)
596 c
597  DO 304 k=1,4
598  hadcur(k)=hadcur(k)+form2*form3*(
599  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
600  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
601  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
602  304 CONTINUE
603  301 CONTINUE
604 c
605  ELSE
606 c ===================================================================
607 c pi0 pi0 p0 pi- case ====
608 c ===================================================================
609  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
610  DO 101 k=1,3
611 c --- loop over thre contribution of the non-omega current
612  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
613  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
614 c -- definition of aa matrix
615 c -- cronecker delta
616  DO 102 i=1,4
617  DO 103 j=1,4
618  103 aa(i,j)=0.0
619  102 aa(i,i)=1.0
620 c
621 c ... and the rest ...
622  DO 104 l=1,3
623  IF (l.NE.k) THEN
624  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
625  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
626  DO 105 i=1,4
627  DO 105 j=1,4
628  sig=1.0
629  IF(j.NE.4) sig=-sig
630  aa(i,j)=aa(i,j)
631  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
632  105 CONTINUE
633  ENDIF
634  104 CONTINUE
635 c --- lets add something to hadcurr
636  form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
637 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikmd(sqrt(qq),ampi,ampi)
638 ccccccccccccc form1=wigfor(sk,amro,gamro) (tests)
639  DO 106 i=1,4
640  DO 106 j=1,4
641  hadcur(i)=
642  $ hadcur(i)+cmplx(coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
643  106 CONTINUE
644 c --- end of the non omega current(3 possibilities)
645  101 CONTINUE
646  ENDIF
647  END
648  FUNCTION wigfor(S,XM,XGAM)
649  COMPLEX wigfor,wignor
650  wignor=cmplx(-xm**2,xm*xgam)
651  wigfor=wignor/cmplx(s-xm**2,xm*xgam)
652  END
653  SUBROUTINE curinf
654 c here the form factors of m. finkemeier et al. start
655 c it ends with the string: m. finkemeier et al. END
656  COMMON /inout/ inut, iout
657  WRITE (unit = iout,fmt = 99)
658  WRITE (unit = iout,fmt = 98)
659 c print *, 'here is curinf'
660  99 FORMAT(
661  . /, ' *************************************************** ',
662  . /, ' YOU ARE USING THE 4 PION DECAY MODE FORM FACTORS ',
663  . /, ' WHICH HAVE BEEN DESCRIBED IN:',
664  . /, ' R. DECKER, M. FINKEMEIER, P. HEILIGER AND H.H. JONSSON',
665  . /, ' "TAU DECAYS INTO FOUR PIONS" ',
666  . /, ' UNIVERSITAET KARLSRUHE PREPRINT TTP 94-13 (1994);',
667  . /, ' LNF-94/066(IR); HEP-PH/9410260 ',
668  . /, ' ',
669  . /, ' PLEASE NOTE THAT THIS ROUTINE IS USING PARAMETERS',
670  . /, ' RELATED TO THE 3 PION DECAY MODE (A1 MODE), SUCH AS',
671  . /, ' THE A1 MASS AND WIDTH (TAKEN FROM THE COMMON /PARMAS/)',
672  . /, ' AND THE 2 PION VECTOR RESONANCE FORM FACTOR (BY USING',
673  . /, ' THE ROUTINE FPIKM)' ,
674  . /, ' THUS IF YOU DECIDE TO CHANGE ANY OF THESE, YOU WILL' ,
675  . /, ' HAVE TO REFIT THE 4 PION PARAMETERS IN THE COMMON' )
676  98 FORMAT(
677  . ' BLOCK /TAU4PI/, OR YOU MIGHT GET A BAD DISCRIPTION' ,
678  . /, ' OF TAU -> 4 PIONS' ,
679  . /, ' for these formfactors set in routine CHOICE for',
680  . /, .eq.' mnum102 -- AMRX=1.42 and GAMRX=.21',
681  . /, .eq.' mnum101 -- AMRX=1.3 and GAMRX=.46 PROB1,PROB2=0.2',
682  . /, ' to optimize phase space parametrization',
683  . /, ' *************************************************** ',
684  . /, ' coded by M. Finkemeier and P. Heiliger, 29. sept. 1994',
685  . /, ' incorporated to TAUOLA by Z. Was 17. jan. 1995',
686 c . /, ' fitted on (day/month/year) by ... ',
687 c . /, ' to .... data ',
688  . /, ' changed by: Z. Was on 17.01.95',
689  . /, ' changes by: M. Finkemeier on 30.01.95' )
690  END
691 c
692  SUBROUTINE curini
693  COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
694  . rom2,rog2,beta2
695  REAL*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
696  . rom2,rog2,beta2
697  gomega = 1.4
698  gamma1 = 0.38
699  gamma2 = 0.38
700  rom1 = 1.35
701  rog1 = 0.3
702  beta1 = 0.08
703  rom2 = 1.70
704  rog2 = 0.235
705  beta2 = -0.0075
706  END
707  COMPLEX FUNCTION bwiga1(QA)
708 c ================================================================
709 c breit-wigner enhancement of a1
710 c ================================================================
711  COMPLEX wigner
712  COMMON / parmas/ amtau,amnuta,amel,amnue,ammu,amnumu,
713  % AMPIZ,ampi,amro,gamro,ama1,gama1,
714  % AMK,amkz,amkst,gamkst
715 
716 c
717  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu,
718  % AMPIZ,ampi,amro,gamro,ama1,gama1,
719  % AMK,amkz,amkst,gamkst
720 
721  wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
722  gamax=gama1*gfun(qa)/gfun(ama1**2)
723  bwiga1=-ama1**2*wigner(qa,ama1,gamax)
724  RETURN
725  END
726  COMPLEX FUNCTION bwigeps(QEPS)
727 c =============================================================
728 c breit-wigner enhancement of epsilon
729 c =============================================================
730  REAL qeps,meps,geps
731  meps=1.300
732  geps=.600
733  bwigeps=cmplx(meps**2,-meps*geps)/
734  % CMPLX(meps**2-qeps,-meps*geps)
735  RETURN
736  END
737  COMPLEX FUNCTION frho4(W,XM1,XM2)
738 c ===========================================================
739 c rho-type resonance factor with higher radials, to be used
740 c by curr for the four pion mode
741 c ===========================================================
742  COMPLEX bwigm
743  COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
744  . rom2,rog2,beta2
745  REAL*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
746  . rom2,rog2,beta2
747  REAL rom,rog,pi,pim,s,w
748  EXTERNAL bwig
749  DATA init /0/
750 c
751 c ------------ parameters --------------------
752  IF (init.EQ.0 ) THEN
753  init=1
754  pi=3.141592654
755  pim=.140
756  rom=0.773
757  rog=0.145
758  ENDIF
759 c -----------------------------------------------
760  s=w**2
761 c print *,'rom2,rog2 =',rom2,rog2
762  frho4=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2)
763  & +beta2*bwigm(s,rom2,rog2,xm1,xm2))
764  & /(1+beta1+beta2)
765  RETURN
766  END
767  SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
768 c ==================================================================
769 c hadronic current for 4 pi final state, according to:
770 c r. decker, m. finkemeier, p. heiliger, h.h.jonsson, ttp94-13
771 c
772 c see also:
773 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
774 c r. decker z. phys c36(1987) 487.
775 c m. gell-mann, d. sharp, w. wagner phys. rev. lett 8 (1962) 261.
776 c ==================================================================
777 
778  COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
779  . rom2,rog2,beta2
780  REAL*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
781  . rom2,rog2,beta2
782  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
783  * ,ampiz,ampi,amro,gamro,ama1,gama1
784  * ,amk,amkz,amkst,gamkst
785 c
786  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
787  * ,ampiz,ampi,amro,gamro,ama1,gama1
788  * ,amk,amkz,amkst,gamkst
789  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
790  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
791  REAL pim1(4),pim2(4),pim3(4),pim4(4),paa(4)
792  COMPLEX hadcur(4),form1,form2,form3,fpikm
793  COMPLEX bwign,frho4
794  COMPLEX bwigeps,bwiga1
795  COMPLEX hcomp1(4),hcomp2(4),hcomp3(4),hcomp4(4)
796 
797  COMPLEX t243,t213,t143,t123,t341,t342
798  COMPLEX t124,t134,t214,t234,t314,t324
799  COMPLEX s2413,s2314,s1423,s1324,s34
800  COMPLEX s2431,s3421
801  COMPLEX brack1,brack2,brack3,brack4a,brack4b,brack4
802 
803  REAL qmp1,qmp2,qmp3,qmp4
804  REAL ps43,ps41,ps42,ps34,ps14,ps13,ps24,ps23
805  REAL ps21,ps31
806 
807  REAL pd243,pd241,pd213,pd143,pd142
808  REAL pd123,pd341,pd342,pd413,pd423
809  REAL pd124,pd134,pd214,pd234,pd314,pd324
810  REAL qp1,qp2,qp3,qp4
811 
812  REAL pa(4),pb(4)
813  REAL aa(4,4),pp(4,4)
814  DATA pi /3.141592653589793238462643/
815  DATA fpi /93.3e-3/
816  DATA init /0/
817  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
818 c
819  IF (init.EQ.0) THEN
820  CALL curini
821  CALL curinf
822  init = 1
823  ENDIF
824 c
825 c --- masses and constants
826  g1=12.924
827  g2=1475.98 * gomega
828  g =g1*g2
829  elpha=-.1
830  amrop=1.7
831  gamrop=0.26
832  amom=.782
833  gamom=0.0085
834  arflat=1.0
835  aromeg=1.0
836 c
837  fro=0.266*amro**2
838  coef1=2.0*sqrt(3.0)/fpi**2*arflat
839  coef2=fro*g*aromeg
840 c --- initialization of four vectors
841  DO 7 k=1,4
842  DO 8 l=1,4
843  8 aa(k,l)=0.0
844  hadcur(k)=cmplx(0.0)
845  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
846  pp(1,k)=pim1(k)
847  pp(2,k)=pim2(k)
848  pp(3,k)=pim3(k)
849  7 pp(4,k)=pim4(k)
850 c
851  IF (mnum.EQ.1) THEN
852 c ===================================================================
853 c pi- pi- p0 pi+ CASE ====
854 c ===================================================================
855  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
856 
857 c first definition of scalar products of momentum vectors
858 
859 c define(q-pi)**2 as qpi:
860 
861  qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
862  % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
863 
864  qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
865  % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
866 
867  qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
868  % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
869 
870  qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
871  % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
872 
873 
874 c define(pi+pk)**2 as psik:
875 
876  ps43=(pim4(4)+pim3(4))**2-(pim4(3)+pim3(3))**2
877  % -(pim4(2)+pim3(2))**2-(pim4(1)+pim3(1))**2
878 
879  ps41=(pim4(4)+pim1(4))**2-(pim4(3)+pim1(3))**2
880  % -(pim4(2)+pim1(2))**2-(pim4(1)+pim1(1))**2
881 
882  ps42=(pim4(4)+pim2(4))**2-(pim4(3)+pim2(3))**2
883  % -(pim4(2)+pim2(2))**2-(pim4(1)+pim2(1))**2
884 
885  ps34=ps43
886 
887  ps14=ps41
888 
889  ps13=(pim1(4)+pim3(4))**2-(pim1(3)+pim3(3))**2
890  % -(pim1(2)+pim3(2))**2-(pim1(1)+pim3(1))**2
891 
892  ps24=ps42
893 
894  ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
895  % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
896 
897  pd243=pim2(4)*(pim4(4)-pim3(4))-pim2(3)*(pim4(3)-pim3(3))
898  % -pim2(2)*(pim4(2)-pim3(2))-pim2(1)*(pim4(1)-pim3(1))
899 
900  pd241=pim2(4)*(pim4(4)-pim1(4))-pim2(3)*(pim4(3)-pim1(3))
901  % -pim2(2)*(pim4(2)-pim1(2))-pim2(1)*(pim4(1)-pim1(1))
902 
903  pd213=pim2(4)*(pim1(4)-pim3(4))-pim2(3)*(pim1(3)-pim3(3))
904  % -pim2(2)*(pim1(2)-pim3(2))-pim2(1)*(pim1(1)-pim3(1))
905 
906  pd143=pim1(4)*(pim4(4)-pim3(4))-pim1(3)*(pim4(3)-pim3(3))
907  % -pim1(2)*(pim4(2)-pim3(2))-pim1(1)*(pim4(1)-pim3(1))
908 
909  pd142=pim1(4)*(pim4(4)-pim2(4))-pim1(3)*(pim4(3)-pim2(3))
910  % -pim1(2)*(pim4(2)-pim2(2))-pim1(1)*(pim4(1)-pim2(1))
911 
912  pd123=pim1(4)*(pim2(4)-pim3(4))-pim1(3)*(pim2(3)-pim3(3))
913  % -pim1(2)*(pim2(2)-pim3(2))-pim1(1)*(pim2(1)-pim3(1))
914 
915  pd341=pim3(4)*(pim4(4)-pim1(4))-pim3(3)*(pim4(3)-pim1(3))
916  % -pim3(2)*(pim4(2)-pim1(2))-pim3(1)*(pim4(1)-pim1(1))
917 
918  pd342=pim3(4)*(pim4(4)-pim2(4))-pim3(3)*(pim4(3)-pim2(3))
919  % -pim3(2)*(pim4(2)-pim2(2))-pim3(1)*(pim4(1)-pim2(1))
920 
921  pd413=pim4(4)*(pim1(4)-pim3(4))-pim4(3)*(pim1(3)-pim3(3))
922  % -pim4(2)*(pim1(2)-pim3(2))-pim4(1)*(pim1(1)-pim3(1))
923 
924  pd423=pim4(4)*(pim2(4)-pim3(4))-pim4(3)*(pim2(3)-pim3(3))
925  % -pim4(2)*(pim2(2)-pim3(2))-pim4(1)*(pim2(1)-pim3(1))
926 
927 c define q*pi = qpi:
928 
929  qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
930  % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
931  % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
932  % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
933 
934  qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
935  % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
936  % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
937  % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
938 
939  qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
940  % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
941  % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
942  % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
943 
944  qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
945  % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
946  % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
947  % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
948 
949 
950 
951 c define t(pi;pj,pk)= tijk:
952 
953  t243=bwiga1(qmp2)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
954  t213=bwiga1(qmp2)*fpikm(sqrt(ps13),ampi,ampi)*gamma1
955  t143=bwiga1(qmp1)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
956  t123=bwiga1(qmp1)*fpikm(sqrt(ps23),ampi,ampi)*gamma1
957  t341=bwiga1(qmp3)*fpikm(sqrt(ps41),ampi,ampi)*gamma1
958  t342=bwiga1(qmp3)*fpikm(sqrt(ps42),ampi,ampi)*gamma1
959 
960 c define s(i,j;k,l)= sijkl:
961 
962  s2413=frho4(sqrt(ps24),ampi,ampi)*gamma2
963  s2314=frho4(sqrt(ps23),ampi,ampi)*bwigeps(ps14)*gamma2
964  s1423=frho4(sqrt(ps14),ampi,ampi)*gamma2
965  s1324=frho4(sqrt(ps13),ampi,ampi)*bwigeps(ps24)*gamma2
966  s34=frho4(sqrt(ps34),ampi,ampi)*gamma2
967 
968 c definition of amplitude, first the [] brackets:
969 
970  brack1=2.*t143+2.*t243+t123+t213
971  % +t341*(pd241/qmp3-1.)+t342*(pd142/qmp3-1.)
972  % +3./4.*(s1423+s2413-s2314-s1324)-3.*s34
973 
974  brack2=2.*t143*pd243/qmp1+3.*t213
975  % +t123*(2.*pd423/qmp1+1.)+t341*(pd241/qmp3+3.)
976  % +t342*(pd142/qmp3+1.)
977  % -3./4.*(s2314+3.*s1324+3.*s1423+s2413)
978 
979  brack3=2.*t243*pd143/qmp2+3.*t123
980  % +t213*(2.*pd413/qmp2+1.)+t341*(pd241/qmp3+1.)
981  % +t342*(pd142/qmp3+3.)
982  % -3./4.*(3.*s2314+s1324+s1423+3.*s2413)
983 
984  brack4a=2.*t143*(pd243/qq*(qp1/qmp1+1.)+pd143/qq)
985  % +2.*t243*(pd143/qq*(qp2/qmp2+1.)+pd243/qq)
986  % +t123+t213
987  % +2.*t123*(pd423/qq*(qp1/qmp1+1.)+pd123/qq)
988  % +2.*t213*(pd413/qq*(qp2/qmp2+1.)+pd213/qq)
989  % +t341*(pd241/qmp3+1.-2.*pd241/qq*(qp3/qmp3+1.)
990  % -2.*pd341/qq)
991  % +t342*(pd142/qmp3+1.-2.*pd142/qq*(qp3/qmp3+1.)
992  % -2.*pd342/qq)
993 
994  brack4b=-3./4.*(s2314*(2.*(qp2-qp3)/qq+1.)
995  % +s1324*(2.*(qp1-qp3)/qq+1.)
996  % +s1423*(2.*(qp1-qp4)/qq+1.)
997  % +s2413*(2.*(qp2-qp4)/qq+1.)
998  % +4.*s34*(qp4-qp3)/qq)
999 
1000  brack4=brack4a+brack4b
1001 
1002  DO 208 k=1,4
1003 
1004  hcomp1(k)=(pim3(k)-pim4(k))*brack1
1005  hcomp2(k)=pim1(k)*brack2
1006  hcomp3(k)=pim2(k)*brack3
1007  hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1008 
1009  208 CONTINUE
1010 
1011  DO 209 i=1,4
1012 
1013  hadcur(i)=hcomp1(i)-hcomp2(i)-hcomp3(i)+hcomp4(i)
1014  hadcur(i)=-coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1015 
1016  209 CONTINUE
1017 
1018 
1019 c --- END of the non omega current(3 possibilities)
1020  201 CONTINUE
1021 c
1022 c
1023 c --- there are two possibilities for omega current
1024 c --- pa pb are corresponding first and second pi-s
1025  DO 301 kk=1,2
1026  DO 302 i=1,4
1027  pa(i)=pp(kk,i)
1028  pb(i)=pp(3-kk,i)
1029  302 CONTINUE
1030 c --- lorentz invariants
1031  qqa=0.0
1032  ss23=0.0
1033  ss24=0.0
1034  ss34=0.0
1035  qp1p2=0.0
1036  qp1p3=0.0
1037  qp1p4=0.0
1038  p1p2 =0.0
1039  p1p3 =0.0
1040  p1p4 =0.0
1041  DO 303 k=1,4
1042  sign=-1.0
1043  IF (k.EQ.4) sign= 1.0
1044  qqa=qqa+sign*(paa(k)-pa(k))**2
1045  ss23=ss23+sign*(pb(k) +pim3(k))**2
1046  ss24=ss24+sign*(pb(k) +pim4(k))**2
1047  ss34=ss34+sign*(pim3(k)+pim4(k))**2
1048  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
1049  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
1050  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
1051  p1p2=p1p2+sign*pa(k)*pb(k)
1052  p1p3=p1p3+sign*pa(k)*pim3(k)
1053  p1p4=p1p4+sign*pa(k)*pim4(k)
1054  303 CONTINUE
1055 c
1056  form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
1057 c form3=bwign(qqa,amom,gamom)*(bwign(ss23,amro,gamro)+
1058 c $ bwign(ss24,amro,gamro)+bwign(ss34,amro,gamro))
1059  form3=bwign(qqa,amom,gamom)
1060 c
1061  DO 304 k=1,4
1062  hadcur(k)=hadcur(k)+form2*form3*(
1063  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
1064  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
1065  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
1066  304 CONTINUE
1067  301 CONTINUE
1068 c
1069  ELSE
1070 c ===================================================================
1071 c pi0 pi0 p0 pi- CASE ====
1072 c ===================================================================
1073  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1074 
1075 
1076 c first definition of scalar products of momentum vectors
1077 
1078 c define(q-pi)**2 as qpi:
1079 
1080  qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1081  % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1082 
1083  qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1084  % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1085 
1086  qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1087  % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1088 
1089  qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1090  % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1091 
1092 
1093 c define(pi+pk)**2 as psik:
1094 
1095  ps14=(pim1(4)+pim4(4))**2-(pim1(3)+pim4(3))**2
1096  % -(pim1(2)+pim4(2))**2-(pim1(1)+pim4(1))**2
1097 
1098  ps21=(pim2(4)+pim1(4))**2-(pim2(3)+pim1(3))**2
1099  % -(pim2(2)+pim1(2))**2-(pim2(1)+pim1(1))**2
1100 
1101  ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1102  % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1103 
1104  ps24=(pim2(4)+pim4(4))**2-(pim2(3)+pim4(3))**2
1105  % -(pim2(2)+pim4(2))**2-(pim2(1)+pim4(1))**2
1106 
1107  ps31=(pim3(4)+pim1(4))**2-(pim3(3)+pim1(3))**2
1108  % -(pim3(2)+pim1(2))**2-(pim3(1)+pim1(1))**2
1109 
1110  ps34=(pim3(4)+pim4(4))**2-(pim3(3)+pim4(3))**2
1111  % -(pim3(2)+pim4(2))**2-(pim3(1)+pim4(1))**2
1112 
1113 
1114 
1115  pd324=pim3(4)*(pim2(4)-pim4(4))-pim3(3)*(pim2(3)-pim4(3))
1116  % -pim3(2)*(pim2(2)-pim4(2))-pim3(1)*(pim2(1)-pim4(1))
1117 
1118  pd314=pim3(4)*(pim1(4)-pim4(4))-pim3(3)*(pim1(3)-pim4(3))
1119  % -pim3(2)*(pim1(2)-pim4(2))-pim3(1)*(pim1(1)-pim4(1))
1120 
1121  pd234=pim2(4)*(pim3(4)-pim4(4))-pim2(3)*(pim3(3)-pim4(3))
1122  % -pim2(2)*(pim3(2)-pim4(2))-pim2(1)*(pim3(1)-pim4(1))
1123 
1124  pd214=pim2(4)*(pim1(4)-pim4(4))-pim2(3)*(pim1(3)-pim4(3))
1125  % -pim2(2)*(pim1(2)-pim4(2))-pim2(1)*(pim1(1)-pim4(1))
1126 
1127  pd134=pim1(4)*(pim3(4)-pim4(4))-pim1(3)*(pim3(3)-pim4(3))
1128  % -pim1(2)*(pim3(2)-pim4(2))-pim1(1)*(pim3(1)-pim4(1))
1129 
1130  pd124=pim1(4)*(pim2(4)-pim4(4))-pim1(3)*(pim2(3)-pim4(3))
1131  % -pim1(2)*(pim2(2)-pim4(2))-pim1(1)*(pim2(1)-pim4(1))
1132 
1133 c define q*pi = qpi:
1134 
1135  qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1136  % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1137  % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1138  % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1139 
1140  qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1141  % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1142  % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1143  % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1144 
1145  qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1146  % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1147  % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1148  % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1149 
1150  qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1151  % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1152  % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1153  % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1154 
1155 
1156 c define t(pi;pj,pk)= tijk:
1157 
1158  t324=bwiga1(qmp3)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1159  t314=bwiga1(qmp3)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1160  t234=bwiga1(qmp2)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1161  t214=bwiga1(qmp2)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1162  t134=bwiga1(qmp1)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1163  t124=bwiga1(qmp1)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1164 
1165 c define s(i,j;k,l)= sijkl:
1166 
1167  s1423=frho4(sqrt(ps14),ampi,ampi)*bwigeps(ps23)*gamma2
1168  s2431=frho4(sqrt(ps24),ampi,ampi)*bwigeps(ps31)*gamma2
1169  s3421=frho4(sqrt(ps34),ampi,ampi)*bwigeps(ps21)*gamma2
1170 
1171 
1172 c definition of amplitude, first the [] brackets:
1173 
1174  brack1=t234+t324+2.*t314+t134+2.*t214+t124
1175  % +t134*pd234/qmp1+t124*pd324/qmp1
1176  % -3./2.*(s3421+s2431+2.*s1423)
1177 
1178 
1179  brack2=t234*(1.+2.*pd134/qmp2)+3.*t324+3.*t124
1180  % +t134*(1.-pd234/qmp1)+2.*t214*pd314/qmp2
1181  % -t124*pd324/qmp1
1182  % -3./2.*(s3421+3.*s2431)
1183 
1184  brack3=t324*(1.+2.*pd124/qmp3)+3.*t234+3.*t134
1185  % +t124*(1.-pd324/qmp1)+2.*t314*pd214/qmp3
1186  % -t134*pd234/qmp1
1187  % -3./2.*(3.*s3421+s2431)
1188 
1189  brack4a=2.*t234*(1./2.+pd134/qq*(qp2/qmp2+1.)+pd234/qq)
1190  % +2.*t324*(1./2.+pd124/qq*(qp3/qmp3+1.)+pd324/qq)
1191  % +2.*t134*(1./2.+pd234/qq*(qp1/qmp1+1.)
1192  % -1./2.*pd234/qmp1+pd134/qq)
1193  % +2.*t124*(1./2.+pd324/qq*(qp1/qmp1+1.)
1194  % -1./2.*pd324/qmp1+pd124/qq)
1195  % +2.*t214*(pd314/qq*(qp2/qmp2+1.)+pd214/qq)
1196  % +2.*t314*(pd214/qq*(qp3/qmp3+1.)+pd314/qq)
1197 
1198  brack4b=-3./2.*(s3421*(2.*(qp3-qp4)/qq+1.)
1199  % +s2431*(2.*(qp2-qp4)/qq+1.)
1200  % +s1423*2.*(qp1-qp4)/qq)
1201 
1202 
1203  brack4=brack4a+brack4b
1204 
1205  DO 308 k=1,4
1206 
1207  hcomp1(k)=(pim1(k)-pim4(k))*brack1
1208  hcomp2(k)=pim2(k)*brack2
1209  hcomp3(k)=pim3(k)*brack3
1210  hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1211 
1212  308 CONTINUE
1213 
1214  DO 309 i=1,4
1215 
1216  hadcur(i)=hcomp1(i)+hcomp2(i)+hcomp3(i)-hcomp4(i)
1217  hadcur(i)=coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1218 
1219  309 CONTINUE
1220 
1221  101 CONTINUE
1222  ENDIF
1223 c m. finkemeier et al. END
1224  END