C++ Interface to Tauola
tauola-F/formf.F
1  FUNCTION formom(XMAA,XMOM)
2 C ==================================================================
3 C formfactorfor pi-pi0 gamma final state
4 C R. Decker, Z. Phys C36 (1987) 487.
5 C ==================================================================
6  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
7  * ,ampiz,ampi,amro,gamro,ama1,gama1
8  * ,amk,amkz,amkst,gamkst
9 C
10  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
11  * ,ampiz,ampi,amro,gamro,ama1,gama1
12  * ,amk,amkz,amkst,gamkst
13  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
14  real*4 gfermi,gv,ga,ccabib,scabib,gamel
15  COMMON /testa1/ keya1
16  COMPLEX BWIGN,FORMOM
17  DATA icont /1/
18 * THIS INLINE FUNCT. CALCULATES THE SCALAR PART OF THE PROPAGATOR
19  bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
20 * HADRON CURRENT
21  fro =0.266*amro**2
22  elpha=- 0.1
23  amrop = 1.7
24  gamrop= 0.26
25  amom =0.782
26  gamom =0.0085
27  aromeg= 1.0
28  gcoup=12.924
29  gcoup=gcoup*aromeg
30  fqed =sqrt(4.0*3.1415926535/137.03604)
31  formom=fqed*fro**2/sqrt(2.0)*gcoup**2*bwign(xmom,amom,gamom)
32  $ *(bwign(xmaa,amro,gamro)+elpha*bwign(xmaa,amrop,gamrop))
33  $ *(bwign( 0.0,amro,gamro)+elpha*bwign( 0.0,amrop,gamrop))
34  END
35 #if defined (CePeCe)
36 #elif defined(ALEPH)
37 #else
38 C=======================================================================
39  COMPLEX FUNCTION fk1ab(XMSQ,INDX)
40 C ==================================================================
41 C complex form-factor for a1+a1prime. AJW 1/98
42 C ==================================================================
43 
44  COMPLEX F1,F2,AMPA,AMPB
45  INTEGER IFIRST,INDX
46  DATA ifirst/0/
47 
48  IF (ifirst.EQ.0) THEN
49  ifirst = 1
50  xm1 = pkorb(1,19)
51  xg1 = pkorb(2,19)
52  xm2 = pkorb(1,20)
53  xg2 = pkorb(2,20)
54 
55  xm1sq = xm1*xm1
56  gf1 = gfun(xm1sq)
57  gg1 = xm1*xg1/gf1
58  xm2sq = xm2*xm2
59  gf2 = gfun(xm2sq)
60  gg2 = xm2*xg2/gf2
61  END IF
62 
63  IF (indx.EQ.1) THEN
64  ampa = cmplx(pkorb(3,81),0.)
65  ampb = cmplx(pkorb(3,82),0.)
66  ELSE IF (indx.EQ.2) THEN
67  ampa = cmplx(pkorb(3,83),0.)
68  ampb = cmplx(pkorb(3,84),0.)
69  ELSEIF (indx.EQ.3) THEN
70  ampa = cmplx(pkorb(3,85),0.)
71  ampb = cmplx(pkorb(3,86),0.)
72  ELSEIF (indx.EQ.4) THEN
73  ampa = cmplx(pkorb(3,87),0.)
74  ampb = cmplx(pkorb(3,88),0.)
75  END IF
76 
77  gf = gfun(xmsq)
78  fg1 = gg1*gf
79  fg2 = gg2*gf
80  f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
81  f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
82  fk1ab = ampa*f1+ampb*f2
83 
84  RETURN
85  END
86 #endif
87  FUNCTION form1(MNUM,QQ,S1,SDWA)
88 C ==================================================================
89 C formfactorfor F1 for 3 scalar final state
90 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
91 C H. Georgi, Weak interactions and modern particle theory,
92 C The Benjamin/Cummings Pub. Co., Inc. 1984.
93 C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
94 C and erratum !!!!!!
95 C ==================================================================
96 C
97  COMPLEX FORM1,WIGNER,WIGFOR,FPIKM,BWIGM
98  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
99  * ,ampiz,ampi,amro,gamro,ama1,gama1
100  * ,amk,amkz,amkst,gamkst
101 C
102  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
103  * ,ampiz,ampi,amro,gamro,ama1,gama1
104  * ,amk,amkz,amkst,gamkst
105 #if defined (CLEO)
106  COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
107  COMPLEX FA1A1P,FK1AB,F3PI
108 C
109  IF (mnum.EQ.0) THEN
110 C ------------ 3 pi hadronic state (a1)
111 C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
112 C FORMRO = F3PI(1,QQ,S1,SDWA)
113 C FORMA1 = FA1A1P(QQ)
114 C FORM1 = FORMA1*FORMRO
115  form1 = f3pi(1,qq,s1,sdwa)
116 
117  ELSEIF (mnum.EQ.1) THEN
118 C ------------ K- pi- K+ (K*0 K-)
119  formks = bwigm(s1,amkst,gamkst,ampi,amk)
120  forma1 = fa1a1p(qq)
121  form1 = forma1*formks
122 
123  ELSEIF (mnum.EQ.2) THEN
124 C ------------ K0 pi- K0B (K*- K0)
125  formks = bwigm(s1,amkst,gamkst,ampi,amk)
126  forma1 = fa1a1p(qq)
127  form1 = forma1*formks
128 
129  ELSEIF (mnum.EQ.3) THEN
130 C ------------ K- pi0 K0 (K*0 K-)
131  formks = bwigm(s1,amkst,gamkst,ampi,amk)
132  forma1 = fa1a1p(qq)
133  form1 = forma1*formks
134 
135  ELSEIF (mnum.EQ.4) THEN
136 C ------------ pi0 pi0 K- (K*-pi0)
137  formks = bwigm(s1,amkst,gamkst,ampi,amk)
138  formk1 = fk1ab(qq,3)
139  form1 = formk1*formks
140 
141  ELSEIF (mnum.EQ.5) THEN
142 C ------------ K- pi- pi+ (rho0 K-)
143  formk1 = fk1ab(qq,4)
144  formro = fpikm(sqrt(s1),ampi,ampi)
145  form1 = formk1*formro
146 
147  ELSEIF (mnum.EQ.6) THEN
148 C ------------ pi- K0B pi0 (pi- K*0B)
149  formk1 = fk1ab(qq,1)
150  formks = bwigm(s1,amkst,gamkst,amk,ampi)
151  form1 = formk1*formks
152 
153  ELSEIF (mnum.EQ.7) THEN
154 C -------------- eta pi- pi0 final state
155  form1=0.0
156  ENDIF
157 #else
158  wigner(a,b,c)= cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
159  IF (mnum.EQ.0) THEN
160 C ------------ 3 pi hadronic state (a1)
161  gamax=gama1*gfun(qq)/gfun(ama1**2)
162  form1=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
163  ELSEIF (mnum.EQ.1) THEN
164 C ------------ K- pi- K+
165 #if defined (ALEPH)
166  form1=bwigm(s1,amkst,gamkst,ampi,amkz)
167 #else
168  form1=bwigm(s1,amkst,gamkst,ampi,amk)
169 #endif
170  gamax=gama1*gfun(qq)/gfun(ama1**2)
171  form1=ama1**2*wigner(qq,ama1,gamax)*form1
172  ELSEIF (mnum.EQ.2) THEN
173 C ------------ K0 pi- K0B
174 #if defined (ALEPH)
175  form1=bwigm(s1,amkst,gamkst,ampi,amkz)
176 #else
177  form1=bwigm(s1,amkst,gamkst,ampi,amk)
178 #endif
179  gamax=gama1*gfun(qq)/gfun(ama1**2)
180  form1=ama1**2*wigner(qq,ama1,gamax)*form1
181  ELSEIF (mnum.EQ.3) THEN
182 C ------------ K- K0 pi0
183  form1=0.0
184  gamax=gama1*gfun(qq)/gfun(ama1**2)
185  form1=ama1**2*wigner(qq,ama1,gamax)*form1
186  ELSEIF (mnum.EQ.4) THEN
187 C ------------ pi0 pi0 K-
188  xm2=1.402
189  gam2=0.174
190 #if defined (ALEPH)
191  form1=bwigm(s1,amkst,gamkst,amk,ampiz)
192 #else
193  form1=bwigm(s1,amkst,gamkst,amk,ampi)
194 #endif
195  form1=wigfor(qq,xm2,gam2)*form1
196  ELSEIF (mnum.EQ.5) THEN
197 C ------------ K- pi- pi+
198  xm2=1.402
199  gam2=0.174
200  form1=wigfor(qq,xm2,gam2)*fpikm(sqrt(s1),ampi,ampi)
201  ELSEIF (mnum.EQ.6) THEN
202  form1=0.0
203  ELSEIF (mnum.EQ.7) THEN
204 C -------------- eta pi- pi0 final state
205  form1=0.0
206  ENDIF
207 #endif
208  END
209  FUNCTION form2(MNUM,QQ,S1,SDWA)
210 C ==================================================================
211 C formfactorfor F2 for 3 scalar final state
212 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
213 C H. Georgi, Weak interactions and modern particle theory,
214 C The Benjamin/Cummings Pub. Co., Inc. 1984.
215 C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
216 C and erratum !!!!!!
217 C ==================================================================
218 C
219  COMPLEX FORM2,WIGNER,WIGFOR,FPIKM,BWIGM
220  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
221  * ,ampiz,ampi,amro,gamro,ama1,gama1
222  * ,amk,amkz,amkst,gamkst
223 C
224  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
225  * ,ampiz,ampi,amro,gamro,ama1,gama1
226  * ,amk,amkz,amkst,gamkst
227 #if defined (CLEO)
228  COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
229  COMPLEX FA1A1P,FK1AB,F3PI
230 
231  IF (mnum.EQ.0) THEN
232 C ------------ 3 pi hadronic state (a1)
233 C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
234 C FORMRO = F3PI(2,QQ,S1,SDWA)
235 C FORMA1 = FA1A1P(QQ)
236 C FORM2 = FORMA1*FORMRO
237  form2 = f3pi(2,qq,s1,sdwa)
238 
239  ELSEIF (mnum.EQ.1) THEN
240 C ------------ K- pi- K+ (rho0 pi-)
241  formro = fpikm(sqrt(s1),amk,amk)
242  forma1 = fa1a1p(qq)
243  form2 = forma1*formro
244 
245  ELSEIF (mnum.EQ.2) THEN
246 C ------------ K0 pi- K0B (rho0 pi-)
247  formro = fpikm(sqrt(s1),amk,amk)
248  forma1 = fa1a1p(qq)
249  form2 = forma1*formro
250 
251  ELSEIF (mnum.EQ.3) THEN
252 C ------------ K- pi0 K0 (rho- pi0)
253  formro = fpikm(sqrt(s1),amk,amk)
254  forma1 = fa1a1p(qq)
255  form2 = forma1*formro
256 
257  ELSEIF (mnum.EQ.4) THEN
258 C ------------ pi0 pi0 K- (K*-pi0)
259  formks = bwigm(s1,amkst,gamkst,ampi,amk)
260  formk1 = fk1ab(qq,3)
261  form2 = formk1*formks
262 
263  ELSEIF (mnum.EQ.5) THEN
264 C ------------ K- pi- pi+ (K*0B pi-)
265  formks = bwigm(s1,amkst,gamkst,ampi,amk)
266  formk1 = fk1ab(qq,1)
267  form2 = formk1*formks
268 C
269  ELSEIF (mnum.EQ.6) THEN
270 C ------------ pi- K0B pi0 (rho- K0B)
271  formro = fpikm(sqrt(s1),ampi,ampi)
272  formk1 = fk1ab(qq,2)
273  form2 = formk1*formro
274 C
275  ELSEIF (mnum.EQ.7) THEN
276 C -------------- eta pi- pi0 final state
277  form2=0.0
278  ENDIF
279 C
280 #else
281  wigner(a,b,c)= cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
282  IF (mnum.EQ.0) THEN
283 C ------------ 3 pi hadronic state (a1)
284  gamax=gama1*gfun(qq)/gfun(ama1**2)
285  form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
286  ELSEIF (mnum.EQ.1) THEN
287 C ------------ K- pi- K+
288  gamax=gama1*gfun(qq)/gfun(ama1**2)
289  form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
290  ELSEIF (mnum.EQ.2) THEN
291 C ------------ K0 pi- K0B
292  gamax=gama1*gfun(qq)/gfun(ama1**2)
293  form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
294  ELSEIF (mnum.EQ.3) THEN
295 C ------------ K- K0 pi0
296  gamax=gama1*gfun(qq)/gfun(ama1**2)
297  form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
298  ELSEIF (mnum.EQ.4) THEN
299 C ------------ pi0 pi0 K-
300  xm2=1.402
301  gam2=0.174
302 #if defined (ALEPH)
303  form2=bwigm(s1,amkst,gamkst,amk,ampiz)
304 #else
305  form2=bwigm(s1,amkst,gamkst,amk,ampi)
306 #endif
307  form2=wigfor(qq,xm2,gam2)*form2
308  ELSEIF (mnum.EQ.5) THEN
309 C ------------ K- pi- pi+
310  xm2=1.402
311  gam2=0.174
312 #if defined (ALEPH)
313  form2=bwigm(s1,amkst,gamkst,amk,ampiz)
314 #else
315  form2=bwigm(s1,amkst,gamkst,amk,ampi)
316 #endif
317  form2=wigfor(qq,xm2,gam2)*form2
318 C
319  ELSEIF (mnum.EQ.6) THEN
320  xm2=1.402
321  gam2=0.174
322  form2=wigfor(qq,xm2,gam2)*fpikm(sqrt(s1),ampi,ampi)
323 C
324  ELSEIF (mnum.EQ.7) THEN
325 C -------------- eta pi- pi0 final state
326  form2=0.0
327  ENDIF
328 C
329 #endif
330  END
331  COMPLEX FUNCTION bwigm(S,M,G,XM1,XM2)
332 C **********************************************************
333 C P-WAVE BREIT-WIGNER FOR RHO
334 C **********************************************************
335  REAL S,M,G,XM1,XM2
336  REAL PI,QS,QM,W,GS
337  DATA init /0/
338 C ------------ PARAMETERS --------------------
339  IF (init.EQ.0) THEN
340  init=1
341  pi=3.141592654
342 C ------- BREIT-WIGNER -----------------------
343  ENDIF
344  IF (s.GT.(xm1+xm2)**2) THEN
345  qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
346  qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
347  w=sqrt(s)
348  gs=g*(m/w)**2*(qs/qm)**3
349  ELSE
350  gs=0.0
351  ENDIF
352  bwigm=m**2/cmplx(m**2-s,-sqrt(s)*gs)
353  RETURN
354  END
355  COMPLEX FUNCTION fpikm(W,XM1,XM2)
356 C **********************************************************
357 C PION FORM FACTOR
358 C **********************************************************
359  COMPLEX BWIGM
360  REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
361  EXTERNAL bwig
362  DATA init /0/
363 C
364 C ------------ PARAMETERS --------------------
365  IF (init.EQ.0 ) THEN
366  init=1
367  pi=3.141592654
368  pim=.140
369  rom=0.773
370  rog=0.145
371  rom1=1.370
372  rog1=0.510
373  beta1=-0.145
374  ENDIF
375 C -----------------------------------------------
376  s=w**2
377  fpikm=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
378  & /(1+beta1)
379  RETURN
380  END
381  COMPLEX FUNCTION fpikmd(W,XM1,XM2)
382 C **********************************************************
383 C PION FORM FACTOR
384 C **********************************************************
385  COMPLEX BWIGM
386  REAL ROM,ROG,ROM1,ROG1,PI,PIM,S,W
387  EXTERNAL bwig
388  DATA init /0/
389 C
390 C ------------ PARAMETERS --------------------
391  IF (init.EQ.0 ) THEN
392  init=1
393  pi=3.141592654
394  pim=.140
395  rom=0.773
396  rog=0.145
397  rom1=1.500
398  rog1=0.220
399  rom2=1.750
400  rog2=0.120
401  beta=6.5
402  delta=-26.0
403  ENDIF
404 C -----------------------------------------------
405  s=w**2
406  fpikmd=(delta*bwigm(s,rom,rog,xm1,xm2)
407  $ +beta*bwigm(s,rom1,rog1,xm1,xm2)
408  $ + bwigm(s,rom2,rog2,xm1,xm2))
409  & /(1+beta+delta)
410  RETURN
411  END
412 
413  FUNCTION form3(MNUM,QQ,S1,SDWA)
414 C ==================================================================
415 C formfactorfor F3 for 3 scalar final state
416 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
417 C H. Georgi, Weak interactions and modern particle theory,
418 C The Benjamin/Cummings Pub. Co., Inc. 1984.
419 C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
420 C and erratum !!!!!!
421 C ==================================================================
422 C
423  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
424  * ,ampiz,ampi,amro,gamro,ama1,gama1
425  * ,amk,amkz,amkst,gamkst
426 C
427  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
428  * ,ampiz,ampi,amro,gamro,ama1,gama1
429  * ,amk,amkz,amkst,gamkst
430 #if defined (CLEO)
431  COMPLEX FORM3,BWIGM
432  COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
433  COMPLEX FA1A1P,FK1AB,F3PI
434 C
435  IF (mnum.EQ.0) THEN
436 C ------------ 3 pi hadronic state (a1)
437 C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
438 C FORMRO = F3PI(3,QQ,S1,SDWA)
439 C FORMA1 = FA1A1P(QQ)
440 C FORM3 = FORMA1*FORMRO
441  form3 = f3pi(3,qq,s1,sdwa)
442 
443  ELSEIF (mnum.EQ.3) THEN
444 C ------------ K- pi0 K0 (K*- K0)
445  formks = bwigm(s1,amkst,gamkst,ampiz,amk)
446  forma1 = fa1a1p(qq)
447  form3 = forma1*formks
448 
449  ELSEIF (mnum.EQ.6) THEN
450 C ------------ pi- K0B pi0 (K*- pi0)
451  formks = bwigm(s1,amkst,gamkst,amk,ampi)
452  formk1 = fk1ab(qq,3)
453  form3 = formk1*formks
454 
455  ELSE
456  form3=cmplx(0.,0.)
457  ENDIF
458 #else
459  COMPLEX FORM3
460  IF (mnum.EQ.6) THEN
461  form3=cmplx(0.0)
462  ELSE
463  form3=cmplx(0.0)
464  ENDIF
465  form3=0
466 
467 #endif
468  END
469  FUNCTION form4(MNUM,QQ,S1,S2,S3)
470 C ==================================================================
471 C formfactorfor F4 for 3 scalar final state
472 C R. Decker, in preparation
473 C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
474 C and erratum !!!!!!
475 C ==================================================================
476 C
477  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
478  * ,ampiz,ampi,amro,gamro,ama1,gama1
479  * ,amk,amkz,amkst,gamkst
480 C
481  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
482  * ,ampiz,ampi,amro,gamro,ama1,gama1
483  * ,amk,amkz,amkst,gamkst
484  COMPLEX FORM4,WIGNER,FPIKM
485  real*4 m
486 #if defined (CLEO)
487 #else
488  wigner(a,b,c)=cmplx(1.0,0.0) /cmplx(a-b**2,b*c)
489  IF (mnum.EQ.0) THEN
490 C ------------ 3 pi hadronic state (a1)
491  g1=5.8
492  g2=6.08
493  fpip=0.02
494  ampip=1.3
495  gampip=0.3
496  s=qq
497  g=gampip
498  xm1=ampiz
499  xm2=amro
500  m =ampip
501  IF (s.GT.(xm1+xm2)**2) THEN
502  qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
503  qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
504  w=sqrt(s)
505  gs=g*(m/w)**2*(qs/qm)**5
506  ELSE
507  gs=0.0
508  ENDIF
509  gamx=gs*w/m
510  form4=g1*g2*fpip/amro**4/ampip**2
511  $ *ampip**2*wigner(qq,ampip,gamx)
512  $ *( s1*(s2-s3)*fpikm(sqrt(s1),ampiz,ampiz)
513  $ +s2*(s1-s3)*fpikm(sqrt(s2),ampiz,ampiz) )
514  ELSEIF (mnum.EQ.1) THEN
515 C ------------ K- pi- K+
516  g1=5.8
517  g2=6.08
518  fpip=0.02
519  ampip=1.3
520  gampip=0.3
521  s=qq
522  g=gampip
523  xm1=ampiz
524  xm2=amro
525  m =ampip
526  IF (s.GT.(xm1+xm2)**2) THEN
527  qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
528  qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
529  w=sqrt(s)
530  gs=g*(m/w)**2*(qs/qm)**5
531  ELSE
532  gs=0.0
533  ENDIF
534  gamx=gs*w/m
535  form4=g1*g2*fpip/amro**4/ampip**2
536  $ *ampip**2*wigner(qq,ampip,gamx)
537  $ *( s1*(s2-s3)*fpikm(sqrt(s1),ampiz,ampiz)
538  $ +s2*(s1-s3)*fpikm(sqrt(s2),ampiz,ampiz) )
539  ELSE
540  form4=cmplx(0.0,0.0)
541  ENDIF
542 #endif
543 #if defined (ALEPH)
544 C ---- this formfactor is switched off .. .
545 cam FORM4=CMPLX(0.0,0.0)
546 #else
547 C ---- this formfactor is switched off .. .
548  form4=cmplx(0.0,0.0)
549 #endif
550  END
551  FUNCTION form5(MNUM,QQ,S1,S2)
552 C ==================================================================
553 C formfactorfor F5 for 3 scalar final state
554 C G. Kramer, W. Palmer, S. Pinsky, Phys. Rev. D30 (1984) 89.
555 C G. Kramer, W. Palmer Z. Phys. C25 (1984) 195.
556 C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
557 C and erratum !!!!!!
558 C ==================================================================
559 C
560  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
561  * ,ampiz,ampi,amro,gamro,ama1,gama1
562  * ,amk,amkz,amkst,gamkst
563 C
564  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
565  * ,ampiz,ampi,amro,gamro,ama1,gama1
566  * ,amk,amkz,amkst,gamkst
567  COMPLEX FORM5,WIGNER,FPIKM,FPIKMD,BWIGM
568 #if defined (CePeCe)
569 C to satisfy compiler (since it isn-t actually used here)
570  wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
571 #elif defined (ALEPH)
572  wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
573 #endif
574  IF (mnum.EQ.0) THEN
575 C ------------ 3 pi hadronic state (a1)
576  form5=0.0
577  ELSEIF (mnum.EQ.1) THEN
578 C ------------ K- pi- K+
579  elpha=-0.2
580  form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
581  $ *( fpikm(sqrt(s2),ampi,ampi)
582  $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
583  ELSEIF (mnum.EQ.2) THEN
584 C ------------ K0 pi- K0B
585  elpha=-0.2
586  form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
587  $ *( fpikm(sqrt(s2),ampi,ampi)
588  $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
589  ELSEIF (mnum.EQ.3) THEN
590 C ------------ K- K0 pi0
591  form5=0.0
592  ELSEIF (mnum.EQ.4) THEN
593 C ------------ pi0 pi0 K-
594  form5=0.0
595  ELSEIF (mnum.EQ.5) THEN
596 C ------------ K- pi- pi+
597  elpha=-0.2
598  form5=bwigm(qq,amkst,gamkst,ampi,amk)/(1+elpha)
599  $ *( fpikm(sqrt(s1),ampi,ampi)
600  $ +elpha*bwigm(s2,amkst,gamkst,ampi,amk))
601  ELSEIF (mnum.EQ.6) THEN
602 C ------------ pi- K0B pi0
603  elpha=-0.2
604  form5=bwigm(qq,amkst,gamkst,ampi,amkz)/(1+elpha)
605  $ *( fpikm(sqrt(s2),ampi,ampi)
606  $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
607  ELSEIF (mnum.EQ.7) THEN
608 C -------------- eta pi- pi0 final state
609  form5=fpikmd(sqrt(qq),ampi,ampi)*fpikm(sqrt(s1),ampi,ampi)
610  ENDIF
611 C
612  END
613 #if defined (ALEPH)
614  SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
615 #else
616  SUBROUTINE currx(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
617 #endif
618 C ==================================================================
619 C hadronic current for 4 pi final state
620 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
621 C R. Decker Z. Phys C36 (1987) 487.
622 C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
623 C ==================================================================
624 
625  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
626  * ,ampiz,ampi,amro,gamro,ama1,gama1
627  * ,amk,amkz,amkst,gamkst
628 C
629  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
630  * ,ampiz,ampi,amro,gamro,ama1,gama1
631  * ,amk,amkz,amkst,gamkst
632  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
633  real*4 gfermi,gv,ga,ccabib,scabib,gamel
634 #if defined (ALEPH)
635 #else
636 C ARBITRARY FIXING OF THE FOUR PI X-SECTION NORMALIZATION
637  COMMON /arbit/ arflat,aromeg
638 #endif
639  REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
640 #if defined (ALEPH)
641 cam COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
642  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,WIGFOR
643 #else
644  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
645 #endif
646  COMPLEX BWIGN
647  REAL PA(4),PB(4)
648  REAL AA(4,4),PP(4,4)
649  DATA pi /3.141592653589793238462643/
650  DATA fpi /93.3e-3/
651  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
652 C
653 C --- masses and constants
654 #if defined (ALEPH)
655 cam rho-prim taken as in Dolinsky et al (PhysLett B174 (1986) 453)
656 cam (best fit to Argus data)
657 #endif
658  g1=12.924
659  g2=1475.98
660  g =g1*g2
661 #if defined (ALEPH)
662 cam ELPHA=-.1
663 cam AMROP=1.7
664 cam GAMROP=0.26
665  elpha= .02
666  amrop=1.250
667  gamrop=0.125
668 #else
669  elpha=-.1
670  amrop=1.7
671  gamrop=0.26
672 #endif
673  amom=.782
674  gamom=0.0085
675 #if defined (ALEPH)
676 cam ARFLAT=1.0
677 cam AROMEG=1.0
678  arflat=1.3
679  aromeg=2.0
680 #else
681  arflat=1.0
682  aromeg=1.0
683 #endif
684 C
685  fro=0.266*amro**2
686  coef1=2.0*sqrt(3.0)/fpi**2*arflat
687  coef2=fro*g*aromeg
688 C --- initialization of four vectors
689  DO 7 k=1,4
690  DO 8 l=1,4
691  8 aa(k,l)=0.0
692  hadcur(k)=cmplx(0.0)
693  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
694  pp(1,k)=pim1(k)
695  pp(2,k)=pim2(k)
696  pp(3,k)=pim3(k)
697  7 pp(4,k)=pim4(k)
698 C
699  IF (mnum.EQ.1) THEN
700 C ===================================================================
701 C pi- pi- p0 pi+ case ====
702 C ===================================================================
703  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
704 C --- loop over thre contribution of the non-omega current
705  DO 201 k=1,3
706  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
707  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
708 C -- definition of AA matrix
709 C -- cronecker delta
710  DO 202 i=1,4
711  DO 203 j=1,4
712  203 aa(i,j)=0.0
713  202 aa(i,i)=1.0
714 C ... and the rest ...
715  DO 204 l=1,3
716  IF (l.NE.k) THEN
717  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
718  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
719  DO 205 i=1,4
720  DO 205 j=1,4
721  sig= 1.0
722  IF(j.NE.4) sig=-sig
723  aa(i,j)=aa(i,j)
724  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
725  205 CONTINUE
726  ENDIF
727  204 CONTINUE
728 C --- lets add something to HADCURR
729 #if defined (ALEPH)
730 cam FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
731 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
732  form1=wigfor(sk,amro,gamro)
733 #else
734  form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
735 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
736 CCCCCCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
737 #endif
738 C
739  fix=1.0
740  IF (k.EQ.3) fix=-2.0
741  DO 206 i=1,4
742  DO 206 j=1,4
743  hadcur(i)=
744  $ hadcur(i)+cmplx(fix*coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
745  206 CONTINUE
746 C --- end of the non omega current (3 possibilities)
747  201 CONTINUE
748 C
749 C
750 C --- there are two possibilities for omega current
751 C --- PA PB are corresponding first and second pi-s
752  DO 301 kk=1,2
753  DO 302 i=1,4
754  pa(i)=pp(kk,i)
755  pb(i)=pp(3-kk,i)
756  302 CONTINUE
757 C --- lorentz invariants
758  qqa=0.0
759  ss23=0.0
760  ss24=0.0
761  ss34=0.0
762  qp1p2=0.0
763  qp1p3=0.0
764  qp1p4=0.0
765  p1p2 =0.0
766  p1p3 =0.0
767  p1p4 =0.0
768  DO 303 k=1,4
769  sign=-1.0
770  IF (k.EQ.4) sign= 1.0
771  qqa=qqa+sign*(paa(k)-pa(k))**2
772  ss23=ss23+sign*(pb(k) +pim3(k))**2
773  ss24=ss24+sign*(pb(k) +pim4(k))**2
774  ss34=ss34+sign*(pim3(k)+pim4(k))**2
775  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
776  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
777  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
778  p1p2=p1p2+sign*pa(k)*pb(k)
779  p1p3=p1p3+sign*pa(k)*pim3(k)
780  p1p4=p1p4+sign*pa(k)*pim4(k)
781  303 CONTINUE
782 C
783  form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
784 C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
785 C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
786  form3=bwign(qqa,amom,gamom)
787 C
788  DO 304 k=1,4
789  hadcur(k)=hadcur(k)+form2*form3*(
790  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
791  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
792  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
793  304 CONTINUE
794  301 CONTINUE
795 C
796  ELSE
797 C ===================================================================
798 C pi0 pi0 p0 pi- case ====
799 C ===================================================================
800  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
801  DO 101 k=1,3
802 C --- loop over thre contribution of the non-omega current
803  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
804  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
805 C -- definition of AA matrix
806 C -- cronecker delta
807  DO 102 i=1,4
808  DO 103 j=1,4
809  103 aa(i,j)=0.0
810  102 aa(i,i)=1.0
811 C
812 C ... and the rest ...
813  DO 104 l=1,3
814  IF (l.NE.k) THEN
815  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
816  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
817  DO 105 i=1,4
818  DO 105 j=1,4
819  sig=1.0
820  IF(j.NE.4) sig=-sig
821  aa(i,j)=aa(i,j)
822  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
823  105 CONTINUE
824  ENDIF
825  104 CONTINUE
826 C --- lets add something to HADCURR
827 #if defined (ALEPH)
828 cam FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
829 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
830  form1=wigfor(sk,amro,gamro)
831 #else
832  form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
833 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
834 CCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
835 #endif
836  DO 106 i=1,4
837  DO 106 j=1,4
838  hadcur(i)=
839  $ hadcur(i)+cmplx(coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
840  106 CONTINUE
841 C --- end of the non omega current (3 possibilities)
842  101 CONTINUE
843  ENDIF
844  END
845  FUNCTION wigfor(S,XM,XGAM)
846  COMPLEX WIGFOR,WIGNOR
847  wignor=cmplx(-xm**2,xm*xgam)
848  wigfor=wignor/cmplx(s-xm**2,xm*xgam)
849  END
850 #if defined (ALEPH)
851 #else
852  SUBROUTINE curinf
853 C HERE the form factors of M. Finkemeier et al. start
854 C it ends with the string: M. Finkemeier et al. END
855  COMMON /inout/ inut, iout
856  WRITE (unit = iout,fmt = 99)
857  WRITE (unit = iout,fmt = 98)
858 c print *, 'here is curinf'
859  99 FORMAT(
860  . /, ' *************************************************** ',
861  . /, ' YOU ARE USING THE 4 PION DECAY MODE FORM FACTORS ',
862  . /, ' WHICH HAVE BEEN DESCRIBED IN:',
863  . /, ' R. DECKER, M. FINKEMEIER, P. HEILIGER AND H.H. JONSSON',
864  . /, ' "TAU DECAYS INTO FOUR PIONS" ',
865  . /, ' UNIVERSITAET KARLSRUHE PREPRINT TTP 94-13 (1994);',
866  . /, ' LNF-94/066(IR); HEP-PH/9410260 ',
867  . /, ' ',
868  . /, ' PLEASE NOTE THAT THIS ROUTINE IS USING PARAMETERS',
869  . /, ' RELATED TO THE 3 PION DECAY MODE (A1 MODE), SUCH AS',
870  . /, ' THE A1 MASS AND WIDTH (TAKEN FROM THE COMMON /PARMAS/)',
871  . /, ' AND THE 2 PION VECTOR RESONANCE FORM FACTOR (BY USING',
872  . /, ' THE ROUTINE FPIKM)' ,
873  . /, ' THUS IF YOU DECIDE TO CHANGE ANY OF THESE, YOU WILL' ,
874  . /, ' HAVE TO REFIT THE 4 PION PARAMETERS IN THE COMMON' )
875  98 FORMAT(
876  . ' BLOCK /TAU4PI/, OR YOU MIGHT GET A BAD DISCRIPTION' ,
877  . /, ' OF TAU -> 4 PIONS' ,
878  . /, ' for these formfactors set in routine CHOICE for',
879  . /, .eq.' mnum102 -- AMRX=1.42 and GAMRX=.21',
880  . /, .eq.' mnum101 -- AMRX=1.3 and GAMRX=.46 PROB1,PROB2=0.2',
881  . /, ' to optimize phase space parametrization',
882  . /, ' *************************************************** ',
883  . /, ' coded by M. Finkemeier and P. Heiliger, 29. sept. 1994',
884  . /, ' incorporated to TAUOLA by Z. Was 17. jan. 1995',
885 c . /, ' fitted on (day/month/year) by ... ',
886 c . /, ' to .... data ',
887  . /, ' changed by: Z. Was on 17.01.95',
888  . /, ' changes by: M. Finkemeier on 30.01.95' )
889  END
890 C
891  SUBROUTINE curini
892  COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
893  . rom2,rog2,beta2
894  real*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
895  . rom2,rog2,beta2
896  gomega = 1.4
897  gamma1 = 0.38
898  gamma2 = 0.38
899  rom1 = 1.35
900  rog1 = 0.3
901  beta1 = 0.08
902  rom2 = 1.70
903  rog2 = 0.235
904  beta2 = -0.0075
905  END
906  COMPLEX FUNCTION bwiga1(QA)
907 C ================================================================
908 C breit-wigner enhancement of a1
909 C ================================================================
910  COMPLEX WIGNER
911  COMMON / parmas/ amtau,amnuta,amel,amnue,ammu,amnumu,
912  % AMPIZ,ampi,amro,gamro,ama1,gama1,
913  % AMK,amkz,amkst,gamkst
914 
915 C
916  real*4 amtau,amnuta,amel,amnue,ammu,amnumu,
917  % AMPIZ,ampi,amro,gamro,ama1,gama1,
918  % AMK,amkz,amkst,gamkst
919 
920  wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
921  gamax=gama1*gfun(qa)/gfun(ama1**2)
922  bwiga1=-ama1**2*wigner(qa,ama1,gamax)
923  RETURN
924  END
925  COMPLEX FUNCTION bwigeps(QEPS)
926 C =============================================================
927 C breit-wigner enhancement of epsilon
928 C =============================================================
929  REAL QEPS,MEPS,GEPS
930  meps=1.300
931  geps=.600
932  bwigeps=cmplx(meps**2,-meps*geps)/
933  % CMPLX(meps**2-qeps,-meps*geps)
934  RETURN
935  END
936  COMPLEX FUNCTION frho4(W,XM1,XM2)
937 C ===========================================================
938 C rho-type resonance factor with higher radials, to be used
939 C by CURR for the four pion mode
940 C ===========================================================
941  COMPLEX BWIGM
942  COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
943  . rom2,rog2,beta2
944  real*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
945  . rom2,rog2,beta2
946  REAL ROM,ROG,PI,PIM,S,W
947  EXTERNAL bwig
948  DATA init /0/
949 C
950 C ------------ PARAMETERS --------------------
951  IF (init.EQ.0 ) THEN
952  init=1
953  pi=3.141592654
954  pim=.140
955  rom=0.773
956  rog=0.145
957  ENDIF
958 C -----------------------------------------------
959  s=w**2
960 c print *,'rom2,rog2 =',rom2,rog2
961  frho4=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2)
962  & +beta2*bwigm(s,rom2,rog2,xm1,xm2))
963  & /(1+beta1+beta2)
964  RETURN
965  END
966  SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
967 C ==================================================================
968 C Hadronic current for 4 pi final state, according to:
969 C R. Decker, M. Finkemeier, P. Heiliger, H.H.Jonsson, TTP94-13
970 C
971 C See also:
972 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
973 C R. Decker Z. Phys C36 (1987) 487.
974 C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
975 C ==================================================================
976 
977  COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
978  . rom2,rog2,beta2
979  real*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
980  . rom2,rog2,beta2
981  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
982  * ,ampiz,ampi,amro,gamro,ama1,gama1
983  * ,amk,amkz,amkst,gamkst
984 C
985  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
986  * ,ampiz,ampi,amro,gamro,ama1,gama1
987  * ,amk,amkz,amkst,gamkst
988  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
989  real*4 gfermi,gv,ga,ccabib,scabib,gamel
990  REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
991  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
992  COMPLEX BWIGN,FRHO4
993  COMPLEX BWIGEPS,BWIGA1
994  COMPLEX HCOMP1(4),HCOMP2(4),HCOMP3(4),HCOMP4(4)
995 
996  COMPLEX T243,T213,T143,T123,T341,T342
997  COMPLEX T124,T134,T214,T234,T314,T324
998  COMPLEX S2413,S2314,S1423,S1324,S34
999  COMPLEX S2431,S3421
1000  COMPLEX BRACK1,BRACK2,BRACK3,BRACK4A,BRACK4B,BRACK4
1001 
1002  REAL QMP1,QMP2,QMP3,QMP4
1003  REAL PS43,PS41,PS42,PS34,PS14,PS13,PS24,PS23
1004  REAL PS21,PS31
1005 
1006  REAL PD243,PD241,PD213,PD143,PD142
1007  REAL PD123,PD341,PD342,PD413,PD423
1008  REAL PD124,PD134,PD214,PD234,PD314,PD324
1009  REAL QP1,QP2,QP3,QP4
1010 
1011  REAL PA(4),PB(4)
1012  REAL AA(4,4),PP(4,4)
1013  DATA pi /3.141592653589793238462643/
1014  DATA fpi /93.3e-3/
1015  DATA init /0/
1016  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
1017 C
1018  IF (init.EQ.0) THEN
1019  CALL curini
1020  CALL curinf
1021  init = 1
1022  ENDIF
1023 C
1024 C --- MASSES AND CONSTANTS
1025  g1=12.924
1026  g2=1475.98 * gomega
1027  g =g1*g2
1028  elpha=-.1
1029  amrop=1.7
1030  gamrop=0.26
1031  amom=.782
1032  gamom=0.0085
1033  arflat=1.0
1034  aromeg=1.0
1035 C
1036  fro=0.266*amro**2
1037  coef1=2.0*sqrt(3.0)/fpi**2*arflat
1038  coef2=fro*g*aromeg
1039 C --- INITIALIZATION OF FOUR VECTORS
1040  DO 7 k=1,4
1041  DO 8 l=1,4
1042  8 aa(k,l)=0.0
1043  hadcur(k)=cmplx(0.0)
1044  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
1045  pp(1,k)=pim1(k)
1046  pp(2,k)=pim2(k)
1047  pp(3,k)=pim3(k)
1048  7 pp(4,k)=pim4(k)
1049 C
1050  IF (mnum.EQ.1) THEN
1051 C ===================================================================
1052 C PI- PI- P0 PI+ CASE ====
1053 C ===================================================================
1054  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1055 
1056 C FIRST DEFINITION OF SCALAR PRODUCTS OF MOMENTUM VECTORS
1057 
1058 C DEFINE (Q-PI)**2 AS QPI:
1059 
1060  qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1061  % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1062 
1063  qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1064  % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1065 
1066  qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1067  % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1068 
1069  qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1070  % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1071 
1072 
1073 C DEFINE (PI+PK)**2 AS PSIK:
1074 
1075  ps43=(pim4(4)+pim3(4))**2-(pim4(3)+pim3(3))**2
1076  % -(pim4(2)+pim3(2))**2-(pim4(1)+pim3(1))**2
1077 
1078  ps41=(pim4(4)+pim1(4))**2-(pim4(3)+pim1(3))**2
1079  % -(pim4(2)+pim1(2))**2-(pim4(1)+pim1(1))**2
1080 
1081  ps42=(pim4(4)+pim2(4))**2-(pim4(3)+pim2(3))**2
1082  % -(pim4(2)+pim2(2))**2-(pim4(1)+pim2(1))**2
1083 
1084  ps34=ps43
1085 
1086  ps14=ps41
1087 
1088  ps13=(pim1(4)+pim3(4))**2-(pim1(3)+pim3(3))**2
1089  % -(pim1(2)+pim3(2))**2-(pim1(1)+pim3(1))**2
1090 
1091  ps24=ps42
1092 
1093  ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1094  % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1095 
1096  pd243=pim2(4)*(pim4(4)-pim3(4))-pim2(3)*(pim4(3)-pim3(3))
1097  % -pim2(2)*(pim4(2)-pim3(2))-pim2(1)*(pim4(1)-pim3(1))
1098 
1099  pd241=pim2(4)*(pim4(4)-pim1(4))-pim2(3)*(pim4(3)-pim1(3))
1100  % -pim2(2)*(pim4(2)-pim1(2))-pim2(1)*(pim4(1)-pim1(1))
1101 
1102  pd213=pim2(4)*(pim1(4)-pim3(4))-pim2(3)*(pim1(3)-pim3(3))
1103  % -pim2(2)*(pim1(2)-pim3(2))-pim2(1)*(pim1(1)-pim3(1))
1104 
1105  pd143=pim1(4)*(pim4(4)-pim3(4))-pim1(3)*(pim4(3)-pim3(3))
1106  % -pim1(2)*(pim4(2)-pim3(2))-pim1(1)*(pim4(1)-pim3(1))
1107 
1108  pd142=pim1(4)*(pim4(4)-pim2(4))-pim1(3)*(pim4(3)-pim2(3))
1109  % -pim1(2)*(pim4(2)-pim2(2))-pim1(1)*(pim4(1)-pim2(1))
1110 
1111  pd123=pim1(4)*(pim2(4)-pim3(4))-pim1(3)*(pim2(3)-pim3(3))
1112  % -pim1(2)*(pim2(2)-pim3(2))-pim1(1)*(pim2(1)-pim3(1))
1113 
1114  pd341=pim3(4)*(pim4(4)-pim1(4))-pim3(3)*(pim4(3)-pim1(3))
1115  % -pim3(2)*(pim4(2)-pim1(2))-pim3(1)*(pim4(1)-pim1(1))
1116 
1117  pd342=pim3(4)*(pim4(4)-pim2(4))-pim3(3)*(pim4(3)-pim2(3))
1118  % -pim3(2)*(pim4(2)-pim2(2))-pim3(1)*(pim4(1)-pim2(1))
1119 
1120  pd413=pim4(4)*(pim1(4)-pim3(4))-pim4(3)*(pim1(3)-pim3(3))
1121  % -pim4(2)*(pim1(2)-pim3(2))-pim4(1)*(pim1(1)-pim3(1))
1122 
1123  pd423=pim4(4)*(pim2(4)-pim3(4))-pim4(3)*(pim2(3)-pim3(3))
1124  % -pim4(2)*(pim2(2)-pim3(2))-pim4(1)*(pim2(1)-pim3(1))
1125 
1126 C DEFINE Q*PI = QPI:
1127 
1128  qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1129  % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1130  % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1131  % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1132 
1133  qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1134  % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1135  % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1136  % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1137 
1138  qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1139  % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1140  % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1141  % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1142 
1143  qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1144  % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1145  % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1146  % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1147 
1148 
1149 
1150 C DEFINE T(PI;PJ,PK)= TIJK:
1151 
1152  t243=bwiga1(qmp2)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
1153  t213=bwiga1(qmp2)*fpikm(sqrt(ps13),ampi,ampi)*gamma1
1154  t143=bwiga1(qmp1)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
1155  t123=bwiga1(qmp1)*fpikm(sqrt(ps23),ampi,ampi)*gamma1
1156  t341=bwiga1(qmp3)*fpikm(sqrt(ps41),ampi,ampi)*gamma1
1157  t342=bwiga1(qmp3)*fpikm(sqrt(ps42),ampi,ampi)*gamma1
1158 
1159 C DEFINE S(I,J;K,L)= SIJKL:
1160 
1161  s2413=frho4(sqrt(ps24),ampi,ampi)*gamma2
1162  s2314=frho4(sqrt(ps23),ampi,ampi)*bwigeps(ps14)*gamma2
1163  s1423=frho4(sqrt(ps14),ampi,ampi)*gamma2
1164  s1324=frho4(sqrt(ps13),ampi,ampi)*bwigeps(ps24)*gamma2
1165  s34=frho4(sqrt(ps34),ampi,ampi)*gamma2
1166 
1167 C DEFINITION OF AMPLITUDE, FIRST THE [] BRACKETS:
1168 
1169  brack1=2.*t143+2.*t243+t123+t213
1170  % +t341*(pd241/qmp3-1.)+t342*(pd142/qmp3-1.)
1171  % +3./4.*(s1423+s2413-s2314-s1324)-3.*s34
1172 
1173  brack2=2.*t143*pd243/qmp1+3.*t213
1174  % +t123*(2.*pd423/qmp1+1.)+t341*(pd241/qmp3+3.)
1175  % +t342*(pd142/qmp3+1.)
1176  % -3./4.*(s2314+3.*s1324+3.*s1423+s2413)
1177 
1178  brack3=2.*t243*pd143/qmp2+3.*t123
1179  % +t213*(2.*pd413/qmp2+1.)+t341*(pd241/qmp3+1.)
1180  % +t342*(pd142/qmp3+3.)
1181  % -3./4.*(3.*s2314+s1324+s1423+3.*s2413)
1182 
1183  brack4a=2.*t143*(pd243/qq*(qp1/qmp1+1.)+pd143/qq)
1184  % +2.*t243*(pd143/qq*(qp2/qmp2+1.)+pd243/qq)
1185  % +t123+t213
1186  % +2.*t123*(pd423/qq*(qp1/qmp1+1.)+pd123/qq)
1187  % +2.*t213*(pd413/qq*(qp2/qmp2+1.)+pd213/qq)
1188  % +t341*(pd241/qmp3+1.-2.*pd241/qq*(qp3/qmp3+1.)
1189  % -2.*pd341/qq)
1190  % +t342*(pd142/qmp3+1.-2.*pd142/qq*(qp3/qmp3+1.)
1191  % -2.*pd342/qq)
1192 
1193  brack4b=-3./4.*(s2314*(2.*(qp2-qp3)/qq+1.)
1194  % +s1324*(2.*(qp1-qp3)/qq+1.)
1195  % +s1423*(2.*(qp1-qp4)/qq+1.)
1196  % +s2413*(2.*(qp2-qp4)/qq+1.)
1197  % +4.*s34*(qp4-qp3)/qq)
1198 
1199  brack4=brack4a+brack4b
1200 
1201  DO 208 k=1,4
1202 
1203  hcomp1(k)=(pim3(k)-pim4(k))*brack1
1204  hcomp2(k)=pim1(k)*brack2
1205  hcomp3(k)=pim2(k)*brack3
1206  hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1207 
1208  208 CONTINUE
1209 
1210  DO 209 i=1,4
1211 
1212  hadcur(i)=hcomp1(i)-hcomp2(i)-hcomp3(i)+hcomp4(i)
1213  hadcur(i)=-coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1214 
1215  209 CONTINUE
1216 
1217 
1218 C --- END OF THE NON OMEGA CURRENT (3 POSSIBILITIES)
1219  201 CONTINUE
1220 C
1221 C
1222 C --- THERE ARE TWO POSSIBILITIES FOR OMEGA CURRENT
1223 C --- PA PB ARE CORRESPONDING FIRST AND SECOND PI-S
1224  DO 301 kk=1,2
1225  DO 302 i=1,4
1226  pa(i)=pp(kk,i)
1227  pb(i)=pp(3-kk,i)
1228  302 CONTINUE
1229 C --- LORENTZ INVARIANTS
1230  qqa=0.0
1231  ss23=0.0
1232  ss24=0.0
1233  ss34=0.0
1234  qp1p2=0.0
1235  qp1p3=0.0
1236  qp1p4=0.0
1237  p1p2 =0.0
1238  p1p3 =0.0
1239  p1p4 =0.0
1240  DO 303 k=1,4
1241  sign=-1.0
1242  IF (k.EQ.4) sign= 1.0
1243  qqa=qqa+sign*(paa(k)-pa(k))**2
1244  ss23=ss23+sign*(pb(k) +pim3(k))**2
1245  ss24=ss24+sign*(pb(k) +pim4(k))**2
1246  ss34=ss34+sign*(pim3(k)+pim4(k))**2
1247  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
1248  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
1249  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
1250  p1p2=p1p2+sign*pa(k)*pb(k)
1251  p1p3=p1p3+sign*pa(k)*pim3(k)
1252  p1p4=p1p4+sign*pa(k)*pim4(k)
1253  303 CONTINUE
1254 C
1255  form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
1256 C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
1257 C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
1258  form3=bwign(qqa,amom,gamom)
1259 C
1260  DO 304 k=1,4
1261  hadcur(k)=hadcur(k)+form2*form3*(
1262  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
1263  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
1264  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
1265  304 CONTINUE
1266  301 CONTINUE
1267 C
1268  ELSE
1269 C ===================================================================
1270 C PI0 PI0 P0 PI- CASE ====
1271 C ===================================================================
1272  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1273 
1274 
1275 C FIRST DEFINITION OF SCALAR PRODUCTS OF MOMENTUM VECTORS
1276 
1277 C DEFINE (Q-PI)**2 AS QPI:
1278 
1279  qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1280  % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1281 
1282  qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1283  % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1284 
1285  qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1286  % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1287 
1288  qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1289  % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1290 
1291 
1292 C DEFINE (PI+PK)**2 AS PSIK:
1293 
1294  ps14=(pim1(4)+pim4(4))**2-(pim1(3)+pim4(3))**2
1295  % -(pim1(2)+pim4(2))**2-(pim1(1)+pim4(1))**2
1296 
1297  ps21=(pim2(4)+pim1(4))**2-(pim2(3)+pim1(3))**2
1298  % -(pim2(2)+pim1(2))**2-(pim2(1)+pim1(1))**2
1299 
1300  ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1301  % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1302 
1303  ps24=(pim2(4)+pim4(4))**2-(pim2(3)+pim4(3))**2
1304  % -(pim2(2)+pim4(2))**2-(pim2(1)+pim4(1))**2
1305 
1306  ps31=(pim3(4)+pim1(4))**2-(pim3(3)+pim1(3))**2
1307  % -(pim3(2)+pim1(2))**2-(pim3(1)+pim1(1))**2
1308 
1309  ps34=(pim3(4)+pim4(4))**2-(pim3(3)+pim4(3))**2
1310  % -(pim3(2)+pim4(2))**2-(pim3(1)+pim4(1))**2
1311 
1312 
1313 
1314  pd324=pim3(4)*(pim2(4)-pim4(4))-pim3(3)*(pim2(3)-pim4(3))
1315  % -pim3(2)*(pim2(2)-pim4(2))-pim3(1)*(pim2(1)-pim4(1))
1316 
1317  pd314=pim3(4)*(pim1(4)-pim4(4))-pim3(3)*(pim1(3)-pim4(3))
1318  % -pim3(2)*(pim1(2)-pim4(2))-pim3(1)*(pim1(1)-pim4(1))
1319 
1320  pd234=pim2(4)*(pim3(4)-pim4(4))-pim2(3)*(pim3(3)-pim4(3))
1321  % -pim2(2)*(pim3(2)-pim4(2))-pim2(1)*(pim3(1)-pim4(1))
1322 
1323  pd214=pim2(4)*(pim1(4)-pim4(4))-pim2(3)*(pim1(3)-pim4(3))
1324  % -pim2(2)*(pim1(2)-pim4(2))-pim2(1)*(pim1(1)-pim4(1))
1325 
1326  pd134=pim1(4)*(pim3(4)-pim4(4))-pim1(3)*(pim3(3)-pim4(3))
1327  % -pim1(2)*(pim3(2)-pim4(2))-pim1(1)*(pim3(1)-pim4(1))
1328 
1329  pd124=pim1(4)*(pim2(4)-pim4(4))-pim1(3)*(pim2(3)-pim4(3))
1330  % -pim1(2)*(pim2(2)-pim4(2))-pim1(1)*(pim2(1)-pim4(1))
1331 
1332 C DEFINE Q*PI = QPI:
1333 
1334  qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1335  % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1336  % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1337  % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1338 
1339  qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1340  % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1341  % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1342  % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1343 
1344  qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1345  % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1346  % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1347  % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1348 
1349  qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1350  % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1351  % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1352  % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1353 
1354 
1355 C DEFINE T(PI;PJ,PK)= TIJK:
1356 
1357  t324=bwiga1(qmp3)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1358  t314=bwiga1(qmp3)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1359  t234=bwiga1(qmp2)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1360  t214=bwiga1(qmp2)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1361  t134=bwiga1(qmp1)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1362  t124=bwiga1(qmp1)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1363 
1364 C DEFINE S(I,J;K,L)= SIJKL:
1365 
1366  s1423=frho4(sqrt(ps14),ampi,ampi)*bwigeps(ps23)*gamma2
1367  s2431=frho4(sqrt(ps24),ampi,ampi)*bwigeps(ps31)*gamma2
1368  s3421=frho4(sqrt(ps34),ampi,ampi)*bwigeps(ps21)*gamma2
1369 
1370 
1371 C DEFINITION OF AMPLITUDE, FIRST THE [] BRACKETS:
1372 
1373  brack1=t234+t324+2.*t314+t134+2.*t214+t124
1374  % +t134*pd234/qmp1+t124*pd324/qmp1
1375  % -3./2.*(s3421+s2431+2.*s1423)
1376 
1377 
1378  brack2=t234*(1.+2.*pd134/qmp2)+3.*t324+3.*t124
1379  % +t134*(1.-pd234/qmp1)+2.*t214*pd314/qmp2
1380  % -t124*pd324/qmp1
1381  % -3./2.*(s3421+3.*s2431)
1382 
1383  brack3=t324*(1.+2.*pd124/qmp3)+3.*t234+3.*t134
1384  % +t124*(1.-pd324/qmp1)+2.*t314*pd214/qmp3
1385  % -t134*pd234/qmp1
1386  % -3./2.*(3.*s3421+s2431)
1387 
1388  brack4a=2.*t234*(1./2.+pd134/qq*(qp2/qmp2+1.)+pd234/qq)
1389  % +2.*t324*(1./2.+pd124/qq*(qp3/qmp3+1.)+pd324/qq)
1390  % +2.*t134*(1./2.+pd234/qq*(qp1/qmp1+1.)
1391  % -1./2.*pd234/qmp1+pd134/qq)
1392  % +2.*t124*(1./2.+pd324/qq*(qp1/qmp1+1.)
1393  % -1./2.*pd324/qmp1+pd124/qq)
1394  % +2.*t214*(pd314/qq*(qp2/qmp2+1.)+pd214/qq)
1395  % +2.*t314*(pd214/qq*(qp3/qmp3+1.)+pd314/qq)
1396 
1397  brack4b=-3./2.*(s3421*(2.*(qp3-qp4)/qq+1.)
1398  % +s2431*(2.*(qp2-qp4)/qq+1.)
1399  % +s1423*2.*(qp1-qp4)/qq)
1400 
1401 
1402  brack4=brack4a+brack4b
1403 
1404  DO 308 k=1,4
1405 
1406  hcomp1(k)=(pim1(k)-pim4(k))*brack1
1407  hcomp2(k)=pim2(k)*brack2
1408  hcomp3(k)=pim3(k)*brack3
1409  hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1410 
1411  308 CONTINUE
1412 
1413  DO 309 i=1,4
1414 
1415  hadcur(i)=hcomp1(i)+hcomp2(i)+hcomp3(i)-hcomp4(i)
1416  hadcur(i)=coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1417 
1418  309 CONTINUE
1419 
1420  101 CONTINUE
1421  ENDIF
1422 C M. Finkemeier et al. END
1423  END
1424 #endif