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