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