C++ Interface to Tauola
tauola-modified/formf.f
1
2
3 FUNCTION formom(XMAA,XMOM)
4C ==================================================================
5C formfactorfor pi-pi0 gamma final state
6C R. Decker, Z. Phys C36 (1987) 487.
7C ==================================================================
8 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
9 * ,ampiz,ampi,amro,gamro,ama1,gama1
10 * ,amk,amkz,amkst,gamkst
11C
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
40C=======================================================================
41 COMPLEX FUNCTION fk1ab(XMSQ,INDX)
42C ==================================================================
43C complex form-factor for a1+a1prime. AJW 1/98
44C ==================================================================
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)
90C ==================================================================
91C formfactorfor F1 for 3 scalar final state
92C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
93C H. Georgi, Weak interactions and modern particle theory,
94C The Benjamin/Cummings Pub. Co., Inc. 1984.
95C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
96C and erratum !!!!!!
97C ==================================================================
98C
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
103C
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
112C
113 IF (mnum.EQ.0) THEN
114C ------------ 3 pi hadronic state (a1)
115C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
116C FORMRO = F3PI(1,QQ,S1,SDWA)
117C FORMA1 = FA1A1P(QQ)
118C 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
126C ------------ 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
132C ------------ 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
138C ------------ 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
144C ------------ 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
150C ------------ 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
156C ------------ 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
162C -------------- eta pi- pi0 final state
163 form1=0.0
164 ENDIF
165 END
166 FUNCTION form2(MNUM,QQ,S1,SDWA)
167C ==================================================================
168C formfactorfor F2 for 3 scalar final state
169C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
170C H. Georgi, Weak interactions and modern particle theory,
171C The Benjamin/Cummings Pub. Co., Inc. 1984.
172C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
173C and erratum !!!!!!
174C ==================================================================
175C
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
180C
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
191C ------------ 3 pi hadronic state (a1)
192C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
193C FORMRO = F3PI(2,QQ,S1,SDWA)
194C FORMA1 = FA1A1P(QQ)
195C FORM2 = FORMA1*FORMRO
196C 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
203C ------------ 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
209C ------------ 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
215C ------------ 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
221C ------------ 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
227C ------------ K- pi- pi+ (K*0B pi-)
228 formks = bwigm(s1,amkst,gamkst,ampi,amk)
229 formk1 = fk1ab(qq,1)
230 form2 = formk1*formks
231C
232 ELSEIF (mnum.EQ.6) THEN
233C ------------ pi- K0B pi0 (rho- K0B)
234 formro = fpikm(sqrt(s1),ampi,ampi)
235 formk1 = fk1ab(qq,2)
236 form2 = formk1*formro
237C
238 ELSEIF (mnum.EQ.7) THEN
239C -------------- eta pi- pi0 final state
240 form2=0.0
241 ENDIF
242C
243 END
244 COMPLEX FUNCTION bwigm(S,M,G,XM1,XM2)
245C **********************************************************
246C P-WAVE BREIT-WIGNER FOR RHO
247C **********************************************************
248 REAL S,M,G,XM1,XM2
249 REAL PI,QS,QM,W,GS
250 DATA init /0/
251C ------------ PARAMETERS --------------------
252 IF (init.EQ.0) THEN
253 init=1
254 pi=3.141592654
255C ------- 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)
269C **********************************************************
270C PION FORM FACTOR
271C **********************************************************
272 COMPLEX BWIGM
273 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
274 EXTERNAL bwig
275 DATA init /0/
276C
277C ------------ 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
288C -----------------------------------------------
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)
295C **********************************************************
296C PION FORM FACTOR
297C **********************************************************
298 COMPLEX BWIGM
299 REAL ROM,ROG,ROM1,ROG1,PI,PIM,S,W
300 EXTERNAL bwig
301 DATA init /0/
302C
303C ------------ 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
317C -----------------------------------------------
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)
327C ==================================================================
328C formfactorfor F3 for 3 scalar final state
329C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
330C H. Georgi, Weak interactions and modern particle theory,
331C The Benjamin/Cummings Pub. Co., Inc. 1984.
332C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
333C and erratum !!!!!!
334C ==================================================================
335C
336 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
337 * ,ampiz,ampi,amro,gamro,ama1,gama1
338 * ,amk,amkz,amkst,gamkst
339C
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
348C
349 IF (mnum.EQ.0) THEN
350C ------------ 3 pi hadronic state (a1)
351C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
352C FORMRO = F3PI(3,QQ,S1,SDWA)
353C FORMA1 = FA1A1P(QQ)
354C 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
362C ------------ 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
368C ------------ 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)
378C ==================================================================
379C formfactorfor F4 for 3 scalar final state
380C R. Decker, in preparation
381C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
382C and erratum !!!!!!
383C ==================================================================
384C
385 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
386 * ,ampiz,ampi,amro,gamro,ama1,gama1
387 * ,amk,amkz,amkst,gamkst
388C
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
397C ---- this formfactor is switched off for cleo version
398 form4=cmplx(0.0,0.0)
399 IF (mnum.EQ.0) THEN
400C ------------ 3 pi hadronic state (a1)
401C FORMRO = FPIKM(SQRT(S1),AMPI,AMPI)
402C FORMRO = F3PI(3,QQ,S1,SDWA)
403C FORMA1 = FA1A1P(QQ)
404C FORM3 = FORMA1*FORMRO
405C 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)
412C ==================================================================
413C formfactorfor F5 for 3 scalar final state
414C G. Kramer, W. Palmer, S. Pinsky, Phys. Rev. D30 (1984) 89.
415C G. Kramer, W. Palmer Z. Phys. C25 (1984) 195.
416C R. Decker, E. Mirkes, R. Sauer, Z. Was Karlsruhe preprint TTP92-25
417C and erratum !!!!!!
418C ==================================================================
419C
420 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
421 * ,ampiz,ampi,amro,gamro,ama1,gama1
422 * ,amk,amkz,amkst,gamkst
423C
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
429C ------------ 3 pi hadronic state (a1)
430 form5=0.0
431 ELSEIF (mnum.EQ.1) THEN
432C ------------ 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
438C ------------ 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
444C ------------ K- K0 pi0
445 form5=0.0
446 ELSEIF (mnum.EQ.4) THEN
447C ------------ pi0 pi0 K-
448 form5=0.0
449 ELSEIF (mnum.EQ.5) THEN
450C ------------ 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
456C ------------ 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
462C -------------- eta pi- pi0 final state
463 form5=fpikmd(sqrt(qq),ampi,ampi)*fpikm(sqrt(s1),ampi,ampi)
464 ENDIF
465C
466 END
467 SUBROUTINE currx(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
468C ==================================================================
469C hadronic current for 4 pi final state
470C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
471C R. Decker Z. Phys C36 (1987) 487.
472C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
473C ==================================================================
474
475 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
476 * ,ampiz,ampi,amro,gamro,ama1,gama1
477 * ,amk,amkz,amkst,gamkst
478C
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
484C 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)
494C
495C --- 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
506C
507 fro=0.266*amro**2
508 coef1=2.0*sqrt(3.0)/fpi**2*arflat
509 coef2=fro*g*aromeg
510C --- 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)
520C
521 IF (mnum.EQ.1) THEN
522C ===================================================================
523C pi- pi- p0 pi+ case ====
524C ===================================================================
525 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
526C --- 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
530C -- definition of AA matrix
531C -- 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
536C ... 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
550C --- lets add something to HADCURR
551 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
552C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
553CCCCCCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
554C
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
562C --- end of the non omega current (3 possibilities)
563 201 CONTINUE
564C
565C
566C --- there are two possibilities for omega current
567C --- 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
573C --- 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
598C
599 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
600C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
601C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
602 form3=bwign(qqa,amom,gamom)
603C
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
611C
612 ELSE
613C ===================================================================
614C pi0 pi0 p0 pi- case ====
615C ===================================================================
616 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
617 DO 101 k=1,3
618C --- 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
621C -- definition of AA matrix
622C -- 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
627C
628C ... 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
642C --- lets add something to HADCURR
643 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
644C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
645CCCCCCCCCCCCC 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
651C --- 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
661C HERE the form factors of M. Finkemeier et al. start
662C 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)
666c 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',
693c . /, ' fitted on (day/month/year) by ... ',
694c . /, ' to .... data ',
695 . /, ' changed by: Z. Was on 17.01.95',
696 . /, ' changes by: M. Finkemeier on 30.01.95' )
697 END
698C
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)
715C ================================================================
716C breit-wigner enhancement of a1
717C ================================================================
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
723C
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)
734C =============================================================
735C breit-wigner enhancement of epsilon
736C =============================================================
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)
745C ===========================================================
746C rho-type resonance factor with higher radials, to be used
747C by CURR for the four pion mode
748C ===========================================================
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/
757C
758C ------------ 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
766C -----------------------------------------------
767 s=w**2
768c 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)
775C ==================================================================
776C Hadronic current for 4 pi final state, according to:
777C R. Decker, M. Finkemeier, P. Heiliger, H.H.Jonsson, TTP94-13
778C
779C See also:
780C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
781C R. Decker Z. Phys C36 (1987) 487.
782C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
783C ==================================================================
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
792C
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)
825C
826 IF (init.EQ.0) THEN
827 CALL curini
828 CALL curinf
829 init = 1
830 ENDIF
831C
832C --- 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
843C
844 fro=0.266*amro**2
845 coef1=2.0*sqrt(3.0)/fpi**2*arflat
846 coef2=fro*g*aromeg
847C --- 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)
857C
858 IF (mnum.EQ.1) THEN
859C ===================================================================
860C PI- PI- P0 PI+ CASE ====
861C ===================================================================
862 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
863
864C FIRST DEFINITION OF SCALAR PRODUCTS OF MOMENTUM VECTORS
865
866C 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
881C 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
934C 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
958C 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
967C 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
975C 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
1026C --- END OF THE NON OMEGA CURRENT (3 POSSIBILITIES)
1027 201 CONTINUE
1028C
1029C
1030C --- THERE ARE TWO POSSIBILITIES FOR OMEGA CURRENT
1031C --- 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
1037C --- 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
1062C
1063 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
1064C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
1065C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
1066 form3=bwign(qqa,amom,gamom)
1067C
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
1075C
1076 ELSE
1077C ===================================================================
1078C PI0 PI0 P0 PI- CASE ====
1079C ===================================================================
1080 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1081
1082
1083C FIRST DEFINITION OF SCALAR PRODUCTS OF MOMENTUM VECTORS
1084
1085C 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
1100C 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
1140C 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
1163C 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
1172C 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
1179C 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
1230C M. Finkemeier et al. END
1231 END