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