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