1 /* copyright(c) 1991-2012 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public
6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public
15 license along with the gnu c library;
if not, see
16 <http://www.gnu.org/licenses/>. */
19 /* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27 /* we
do support the iec 559 math functionality,
real and complex. */
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
32 /* we
do not support c11 <threads.h>. */
34 FUNCTION formom(XMAA,XMOM)
35 c ==================================================================
36 c formfactorfor pi-pi0 gamma final state
37 c r. decker, z. phys c36(1987) 487.
38 c ==================================================================
39 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
40 * ,ampiz,ampi,amro,gamro,ama1,gama1
41 * ,amk,amkz,amkst,gamkst
43 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
44 * ,ampiz,ampi,amro,gamro,ama1,gama1
45 * ,amk,amkz,amkst,gamkst
46 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
47 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
51 * this inline funct. calculates the scalar part of the propagator
52 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
63 fqed =sqrt(4.0*3.1415926535/137.03604)
64 formom=fqed*fro**2/sqrt(2.0)*gcoup**2*bwign(xmom,amom,gamom)
65 $ *(bwign(xmaa,amro,gamro)+elpha*bwign(xmaa,amrop,gamrop))
66 $ *(bwign( 0.0,amro,gamro)+elpha*bwign( 0.0,amrop,gamrop))
68 c=======================================================================
69 COMPLEX FUNCTION fk1ab(XMSQ,INDX)
70 c ==================================================================
71 c
complex form-factor for a1+a1prime. ajw 1/98
72 c ==================================================================
74 COMPLEX f1,f2,ampa,ampb
94 ampa = cmplx(pkorb(3,81),0.)
95 ampb = cmplx(pkorb(3,82),0.)
96 ELSE IF (indx.EQ.2)
THEN
97 ampa = cmplx(pkorb(3,83),0.)
98 ampb = cmplx(pkorb(3,84),0.)
99 ELSEIF (indx.EQ.3)
THEN
100 ampa = cmplx(pkorb(3,85),0.)
101 ampb = cmplx(pkorb(3,86),0.)
102 ELSEIF (indx.EQ.4)
THEN
103 ampa = cmplx(pkorb(3,87),0.)
104 ampb = cmplx(pkorb(3,88),0.)
110 f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
111 f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
112 fk1ab = ampa*f1+ampb*f2
116 FUNCTION form1(MNUM,QQ,S1,SDWA)
117 c ==================================================================
118 c formfactorfor f1 for 3 scalar final state
119 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
120 c h. georgi, weak interactions and modern particle theory,
121 c the benjamin/cummings pub. co., inc. 1984.
122 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
124 c ==================================================================
126 COMPLEX form1,wigner,wigfor,fpikm,bwigm
127 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
128 * ,ampiz,ampi,amro,gamro,ama1,gama1
129 * ,amk,amkz,amkst,gamkst
131 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
132 * ,ampiz,ampi,amro,gamro,ama1,gama1
133 * ,amk,amkz,amkst,gamkst
134 COMPLEX forma1,formk1,formro,formks
135 COMPLEX fa1a1p,fk1ab,f3pi
138 c ------------ 3 pi hadronic state(a1)
139 c formro = fpikm(sqrt(s1),ampi,ampi)
140 c formro = f3pi(1,qq,s1,sdwa)
141 c forma1 = fa1a1p(qq)
142 c form1 = forma1*formro
143 form1 = f3pi(1,qq,s1,sdwa)
145 ELSEIF (mnum.EQ.1)
THEN
146 c ------------ k- pi- k+ (k*0 k-)
147 formks = bwigm(s1,amkst,gamkst,ampi,amk)
149 form1 = forma1*formks
151 ELSEIF (mnum.EQ.2)
THEN
152 c ------------ k0 pi- k0b(k*- k0)
153 formks = bwigm(s1,amkst,gamkst,ampi,amk)
155 form1 = forma1*formks
157 ELSEIF (mnum.EQ.3)
THEN
158 c ------------ k- pi0 k0(k*0 k-)
159 formks = bwigm(s1,amkst,gamkst,ampi,amk)
161 form1 = forma1*formks
163 ELSEIF (mnum.EQ.4)
THEN
164 c ------------ pi0 pi0 k- (k*-pi0)
165 formks = bwigm(s1,amkst,gamkst,ampi,amk)
167 form1 = formk1*formks
169 ELSEIF (mnum.EQ.5)
THEN
170 c ------------ k- pi- pi+ (rho0 k-)
172 formro = fpikm(sqrt(s1),ampi,ampi)
173 form1 = formk1*formro
175 ELSEIF (mnum.EQ.6)
THEN
176 c ------------ pi- k0b pi0(pi- k*0b)
178 formks = bwigm(s1,amkst,gamkst,amk,ampi)
179 form1 = formk1*formks
181 ELSEIF (mnum.EQ.7)
THEN
182 c -------------- eta pi- pi0 final state
186 FUNCTION form2(MNUM,QQ,S1,SDWA)
187 c ==================================================================
188 c formfactorfor f2 for 3 scalar final state
189 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
190 c h. georgi, weak interactions and modern particle theory,
191 c the benjamin/cummings pub. co., inc. 1984.
192 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
194 c ==================================================================
196 COMPLEX form2,wigner,wigfor,fpikm,bwigm
197 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
198 * ,ampiz,ampi,amro,gamro,ama1,gama1
199 * ,amk,amkz,amkst,gamkst
201 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
202 * ,ampiz,ampi,amro,gamro,ama1,gama1
203 * ,amk,amkz,amkst,gamkst
204 COMPLEX forma1,formk1,formro,formks
205 COMPLEX fa1a1p,fk1ab,f3pi
208 c ------------ 3 pi hadronic state(a1)
209 c formro = fpikm(sqrt(s1),ampi,ampi)
210 c formro = f3pi(2,qq,s1,sdwa)
211 c forma1 = fa1a1p(qq)
212 c form2 = forma1*formro
213 form2 = f3pi(2,qq,s1,sdwa)
215 ELSEIF (mnum.EQ.1)
THEN
216 c ------------ k- pi- k+ (rho0 pi-)
217 formro = fpikm(sqrt(s1),amk,amk)
219 form2 = forma1*formro
221 ELSEIF (mnum.EQ.2)
THEN
222 c ------------ k0 pi- k0b(rho0 pi-)
223 formro = fpikm(sqrt(s1),amk,amk)
225 form2 = forma1*formro
227 ELSEIF (mnum.EQ.3)
THEN
228 c ------------ k- pi0 k0(rho- pi0)
229 formro = fpikm(sqrt(s1),amk,amk)
231 form2 = forma1*formro
233 ELSEIF (mnum.EQ.4)
THEN
234 c ------------ pi0 pi0 k- (k*-pi0)
235 formks = bwigm(s1,amkst,gamkst,ampi,amk)
237 form2 = formk1*formks
239 ELSEIF (mnum.EQ.5)
THEN
240 c ------------ k- pi- pi+ (k*0b pi-)
241 formks = bwigm(s1,amkst,gamkst,ampi,amk)
243 form2 = formk1*formks
245 ELSEIF (mnum.EQ.6)
THEN
246 c ------------ pi- k0b pi0(rho- k0b)
247 formro = fpikm(sqrt(s1),ampi,ampi)
249 form2 = formk1*formro
251 ELSEIF (mnum.EQ.7)
THEN
252 c -------------- eta pi- pi0 final state
257 COMPLEX FUNCTION bwigm(S,M,G,XM1,XM2)
258 c **********************************************************
259 c p-wave breit-wigner for rho
260 c **********************************************************
264 c ------------ parameters --------------------
268 c ------- breit-wigner -----------------------
270 IF (s.GT.(xm1+xm2)**2)
THEN
271 qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
272 qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
274 gs=g*(m/w)**2*(qs/qm)**3
278 bwigm=m**2/cmplx(m**2-s,-sqrt(s)*gs)
281 COMPLEX FUNCTION fpikm(W,XM1,XM2)
282 c **********************************************************
284 c **********************************************************
286 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
290 c ------------ parameters --------------------
301 c -----------------------------------------------
303 fpikm=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
307 COMPLEX FUNCTION fpikmd(W,XM1,XM2)
308 c **********************************************************
310 c **********************************************************
312 REAL rom,rog,rom1,rog1,pi,pim,s,w
316 c ------------ parameters --------------------
330 c -----------------------------------------------
332 fpikmd=(delta*bwigm(s,rom,rog,xm1,xm2)
333 $ +beta*bwigm(s,rom1,rog1,xm1,xm2)
334 $ + bwigm(s,rom2,rog2,xm1,xm2))
339 FUNCTION form3(MNUM,QQ,S1,SDWA)
340 c ==================================================================
341 c formfactorfor f3 for 3 scalar final state
342 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
343 c h. georgi, weak interactions and modern particle theory,
344 c the benjamin/cummings pub. co., inc. 1984.
345 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
347 c ==================================================================
349 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
350 * ,ampiz,ampi,amro,gamro,ama1,gama1
351 * ,amk,amkz,amkst,gamkst
353 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
354 * ,ampiz,ampi,amro,gamro,ama1,gama1
355 * ,amk,amkz,amkst,gamkst
357 COMPLEX forma1,formk1,formro,formks
358 COMPLEX fa1a1p,fk1ab,f3pi
361 c ------------ 3 pi hadronic state(a1)
362 c formro = fpikm(sqrt(s1),ampi,ampi)
363 c formro = f3pi(3,qq,s1,sdwa)
364 c forma1 = fa1a1p(qq)
365 c form3 = forma1*formro
366 form3 = f3pi(3,qq,s1,sdwa)
368 ELSEIF (mnum.EQ.3)
THEN
369 c ------------ k- pi0 k0(k*- k0)
370 formks = bwigm(s1,amkst,gamkst,ampiz,amk)
372 form3 = forma1*formks
374 ELSEIF (mnum.EQ.6)
THEN
375 c ------------ pi- k0b pi0(k*- pi0)
376 formks = bwigm(s1,amkst,gamkst,amk,ampi)
378 form3 = formk1*formks
384 FUNCTION form4(MNUM,QQ,S1,S2,S3)
385 c ==================================================================
386 c formfactorfor f4 for 3 scalar final state
387 c r. decker, in preparation
388 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
390 c ==================================================================
392 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
393 * ,ampiz,ampi,amro,gamro,ama1,gama1
394 * ,amk,amkz,amkst,gamkst
396 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
397 * ,ampiz,ampi,amro,gamro,ama1,gama1
398 * ,amk,amkz,amkst,gamkst
399 COMPLEX form4,wigner,fpikm
401 c ---- this formfactor is switched off .. .
404 FUNCTION form5(MNUM,QQ,S1,S2)
405 c ==================================================================
406 c formfactorfor f5 for 3 scalar final state
407 c g. kramer, w. palmer, s. pinsky, phys. rev. d30(1984) 89.
408 c g. kramer, w. palmer z. phys. c25(1984) 195.
409 c r. decker, e. mirkes, r. sauer, z. was karlsruhe preprint ttp92-25
411 c ==================================================================
413 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
414 * ,ampiz,ampi,amro,gamro,ama1,gama1
415 * ,amk,amkz,amkst,gamkst
417 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
418 * ,ampiz,ampi,amro,gamro,ama1,gama1
419 * ,amk,amkz,amkst,gamkst
420 COMPLEX form5,wigner,fpikm,fpikmd,bwigm
422 c ------------ 3 pi hadronic state(a1)
424 ELSEIF (mnum.EQ.1)
THEN
425 c ------------ k- pi- k+
427 form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
428 $ *( fpikm(sqrt(s2),ampi,ampi)
429 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
430 ELSEIF (mnum.EQ.2)
THEN
431 c ------------ k0 pi- k0b
433 form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
434 $ *( fpikm(sqrt(s2),ampi,ampi)
435 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
436 ELSEIF (mnum.EQ.3)
THEN
437 c ------------ k- k0 pi0
439 ELSEIF (mnum.EQ.4)
THEN
440 c ------------ pi0 pi0 k-
442 ELSEIF (mnum.EQ.5)
THEN
443 c ------------ k- pi- pi+
445 form5=bwigm(qq,amkst,gamkst,ampi,amk)/(1+elpha)
446 $ *( fpikm(sqrt(s1),ampi,ampi)
447 $ +elpha*bwigm(s2,amkst,gamkst,ampi,amk))
448 ELSEIF (mnum.EQ.6)
THEN
449 c ------------ pi- k0b pi0
451 form5=bwigm(qq,amkst,gamkst,ampi,amkz)/(1+elpha)
452 $ *( fpikm(sqrt(s2),ampi,ampi)
453 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
454 ELSEIF (mnum.EQ.7)
THEN
455 c -------------- eta pi- pi0 final state
456 form5=fpikmd(sqrt(qq),ampi,ampi)*fpikm(sqrt(s1),ampi,ampi)
460 SUBROUTINE currx(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
461 c ==================================================================
462 c hadronic current for 4 pi final state
463 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
464 c r. decker z. phys c36(1987) 487.
465 c m. gell-mann, d. sharp, w. wagner phys. rev. lett 8 (1962) 261.
466 c ==================================================================
468 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
469 * ,ampiz,ampi,amro,gamro,ama1,gama1
470 * ,amk,amkz,amkst,gamkst
472 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
473 * ,ampiz,ampi,amro,gamro,ama1,gama1
474 * ,amk,amkz,amkst,gamkst
475 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
476 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
477 c arbitrary fixing of the four pi x-section normalization
478 COMMON /arbit/ arflat,aromeg
479 REAL pim1(4),pim2(4),pim3(4),pim4(4),paa(4)
480 COMPLEX hadcur(4),form1,form2,form3,fpikm
484 DATA pi /3.141592653589793238462643/
486 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
488 c --- masses and constants
501 coef1=2.0*sqrt(3.0)/fpi**2*arflat
503 c --- initialization of four vectors
508 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
515 c ===================================================================
516 c pi- pi- p0 pi+
case ====
517 c ===================================================================
518 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
519 c --- loop over thre contribution of the non-omega current
521 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
522 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
523 c -- definition of aa matrix
529 c ... and the rest ...
532 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
533 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
539 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
543 c --- lets add something to hadcurr
544 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
545 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikmd(sqrt(qq),ampi,ampi)
546 ccccccccccccccccc form1=wigfor(sk,amro,gamro) (tests)
553 $ hadcur(i)+cmplx(fix*coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
555 c ---
end of the non omega current(3 possibilities)
559 c --- there are two possibilities for omega current
560 c --- pa pb are corresponding first and second pi-s
566 c --- lorentz invariants
579 IF (k.EQ.4) sign= 1.0
580 qqa=qqa+sign*(paa(k)-pa(k))**2
581 ss23=ss23+sign*(pb(k) +pim3(k))**2
582 ss24=ss24+sign*(pb(k) +pim4(k))**2
583 ss34=ss34+sign*(pim3(k)+pim4(k))**2
584 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
585 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
586 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
587 p1p2=p1p2+sign*pa(k)*pb(k)
588 p1p3=p1p3+sign*pa(k)*pim3(k)
589 p1p4=p1p4+sign*pa(k)*pim4(k)
592 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
593 c form3=bwign(qqa,amom,gamom)*(bwign(ss23,amro,gamro)+
594 c $ bwign(ss24,amro,gamro)+bwign(ss34,amro,gamro))
595 form3=bwign(qqa,amom,gamom)
598 hadcur(k)=hadcur(k)+form2*form3*(
599 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
600 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
601 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
606 c ===================================================================
607 c pi0 pi0 p0 pi-
case ====
608 c ===================================================================
609 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
611 c --- loop over thre contribution of the non-omega current
612 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
613 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
614 c -- definition of aa matrix
621 c ... and the rest ...
624 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
625 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
631 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
635 c --- lets add something to hadcurr
636 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
637 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikmd(sqrt(qq),ampi,ampi)
638 ccccccccccccc form1=wigfor(sk,amro,gamro) (tests)
642 $ hadcur(i)+cmplx(coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
644 c ---
end of the non omega current(3 possibilities)
648 FUNCTION wigfor(S,XM,XGAM)
649 COMPLEX wigfor,wignor
650 wignor=cmplx(-xm**2,xm*xgam)
651 wigfor=wignor/cmplx(s-xm**2,xm*xgam)
654 c here the form factors of m. finkemeier et al. start
655 c it ends with the string: m. finkemeier et al.
END
656 COMMON /inout/ inut, iout
657 WRITE (unit = iout,fmt = 99)
658 WRITE (unit = iout,fmt = 98)
659 c print *,
'here is curinf'
661 . /,
' *************************************************** ',
662 . /,
' YOU ARE USING THE 4 PION DECAY MODE FORM FACTORS ',
663 . /,
' WHICH HAVE BEEN DESCRIBED IN:',
664 . /,
' R. DECKER, M. FINKEMEIER, P. HEILIGER AND H.H. JONSSON',
665 . /,
' "TAU DECAYS INTO FOUR PIONS" ',
666 . /,
' UNIVERSITAET KARLSRUHE PREPRINT TTP 94-13 (1994);',
667 . /,
' LNF-94/066(IR); HEP-PH/9410260 ',
669 . /,
' PLEASE NOTE THAT THIS ROUTINE IS USING PARAMETERS',
670 . /,
' RELATED TO THE 3 PION DECAY MODE (A1 MODE), SUCH AS',
671 . /,
' THE A1 MASS AND WIDTH (TAKEN FROM THE COMMON /PARMAS/)',
672 . /,
' AND THE 2 PION VECTOR RESONANCE FORM FACTOR (BY USING',
673 . /,
' THE ROUTINE FPIKM)' ,
674 . /,
' THUS IF YOU DECIDE TO CHANGE ANY OF THESE, YOU WILL' ,
675 . /,
' HAVE TO REFIT THE 4 PION PARAMETERS IN THE COMMON' )
677 .
' BLOCK /TAU4PI/, OR YOU MIGHT GET A BAD DISCRIPTION' ,
678 . /,
' OF TAU -> 4 PIONS' ,
679 . /,
' for these formfactors set in routine CHOICE for',
680 . /, .eq.
' mnum102 -- AMRX=1.42 and GAMRX=.21',
681 . /, .eq.
' mnum101 -- AMRX=1.3 and GAMRX=.46 PROB1,PROB2=0.2',
682 . /,
' to optimize phase space parametrization',
683 . /,
' *************************************************** ',
684 . /,
' coded by M. Finkemeier and P. Heiliger, 29. sept. 1994',
685 . /,
' incorporated to TAUOLA by Z. Was 17. jan. 1995',
686 c . /,
' fitted on (day/month/year) by ... ',
687 c . /,
' to .... data ',
688 . /,
' changed by: Z. Was on 17.01.95',
689 . /,
' changes by: M. Finkemeier on 30.01.95' )
693 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
695 REAL*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
707 COMPLEX FUNCTION bwiga1(QA)
708 c ================================================================
709 c breit-wigner enhancement of a1
710 c ================================================================
712 COMMON / parmas/ amtau,amnuta,amel,amnue,ammu,amnumu,
713 % AMPIZ,ampi,amro,gamro,ama1,gama1,
714 % AMK,amkz,amkst,gamkst
717 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu,
718 % AMPIZ,ampi,amro,gamro,ama1,gama1,
719 % AMK,amkz,amkst,gamkst
721 wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
722 gamax=gama1*gfun(qa)/gfun(ama1**2)
723 bwiga1=-ama1**2*wigner(qa,ama1,gamax)
726 COMPLEX FUNCTION bwigeps(QEPS)
727 c =============================================================
728 c breit-wigner enhancement of epsilon
729 c =============================================================
733 bwigeps=cmplx(meps**2,-meps*geps)/
734 % CMPLX(meps**2-qeps,-meps*geps)
737 COMPLEX FUNCTION frho4(W,XM1,XM2)
738 c ===========================================================
739 c rho-
type resonance factor with higher radials, to be used
740 c by curr for the four pion mode
741 c ===========================================================
743 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
745 REAL*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
747 REAL rom,rog,pi,pim,s,w
751 c ------------ parameters --------------------
759 c -----------------------------------------------
761 c print *,
'rom2,rog2 =',rom2,rog2
762 frho4=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2)
763 & +beta2*bwigm(s,rom2,rog2,xm1,xm2))
767 SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
768 c ==================================================================
769 c hadronic current for 4 pi final state, according to:
770 c r. decker, m. finkemeier, p. heiliger, h.h.jonsson, ttp94-13
773 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
774 c r. decker z. phys c36(1987) 487.
775 c m. gell-mann, d. sharp, w. wagner phys. rev. lett 8 (1962) 261.
776 c ==================================================================
778 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
780 REAL*4 gomega,gamma1,gamma2,rom1,rog1,beta1,
782 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
783 * ,ampiz,ampi,amro,gamro,ama1,gama1
784 * ,amk,amkz,amkst,gamkst
786 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
787 * ,ampiz,ampi,amro,gamro,ama1,gama1
788 * ,amk,amkz,amkst,gamkst
789 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
790 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
791 REAL pim1(4),pim2(4),pim3(4),pim4(4),paa(4)
792 COMPLEX hadcur(4),form1,form2,form3,fpikm
794 COMPLEX bwigeps,bwiga1
795 COMPLEX hcomp1(4),hcomp2(4),hcomp3(4),hcomp4(4)
797 COMPLEX t243,t213,t143,t123,t341,t342
798 COMPLEX t124,t134,t214,t234,t314,t324
799 COMPLEX s2413,s2314,s1423,s1324,s34
801 COMPLEX brack1,brack2,brack3,brack4a,brack4b,brack4
803 REAL qmp1,qmp2,qmp3,qmp4
804 REAL ps43,ps41,ps42,ps34,ps14,ps13,ps24,ps23
807 REAL pd243,pd241,pd213,pd143,pd142
808 REAL pd123,pd341,pd342,pd413,pd423
809 REAL pd124,pd134,pd214,pd234,pd314,pd324
814 DATA pi /3.141592653589793238462643/
817 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
825 c --- masses and constants
838 coef1=2.0*sqrt(3.0)/fpi**2*arflat
840 c --- initialization of four vectors
845 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
852 c ===================================================================
853 c pi- pi- p0 pi+
CASE ====
854 c ===================================================================
855 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
857 c first definition of scalar products of momentum vectors
859 c define(q-pi)**2 as qpi:
861 qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
862 % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
864 qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
865 % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
867 qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
868 % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
870 qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
871 % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
874 c define(pi+pk)**2 as psik:
876 ps43=(pim4(4)+pim3(4))**2-(pim4(3)+pim3(3))**2
877 % -(pim4(2)+pim3(2))**2-(pim4(1)+pim3(1))**2
879 ps41=(pim4(4)+pim1(4))**2-(pim4(3)+pim1(3))**2
880 % -(pim4(2)+pim1(2))**2-(pim4(1)+pim1(1))**2
882 ps42=(pim4(4)+pim2(4))**2-(pim4(3)+pim2(3))**2
883 % -(pim4(2)+pim2(2))**2-(pim4(1)+pim2(1))**2
889 ps13=(pim1(4)+pim3(4))**2-(pim1(3)+pim3(3))**2
890 % -(pim1(2)+pim3(2))**2-(pim1(1)+pim3(1))**2
894 ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
895 % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
897 pd243=pim2(4)*(pim4(4)-pim3(4))-pim2(3)*(pim4(3)-pim3(3))
898 % -pim2(2)*(pim4(2)-pim3(2))-pim2(1)*(pim4(1)-pim3(1))
900 pd241=pim2(4)*(pim4(4)-pim1(4))-pim2(3)*(pim4(3)-pim1(3))
901 % -pim2(2)*(pim4(2)-pim1(2))-pim2(1)*(pim4(1)-pim1(1))
903 pd213=pim2(4)*(pim1(4)-pim3(4))-pim2(3)*(pim1(3)-pim3(3))
904 % -pim2(2)*(pim1(2)-pim3(2))-pim2(1)*(pim1(1)-pim3(1))
906 pd143=pim1(4)*(pim4(4)-pim3(4))-pim1(3)*(pim4(3)-pim3(3))
907 % -pim1(2)*(pim4(2)-pim3(2))-pim1(1)*(pim4(1)-pim3(1))
909 pd142=pim1(4)*(pim4(4)-pim2(4))-pim1(3)*(pim4(3)-pim2(3))
910 % -pim1(2)*(pim4(2)-pim2(2))-pim1(1)*(pim4(1)-pim2(1))
912 pd123=pim1(4)*(pim2(4)-pim3(4))-pim1(3)*(pim2(3)-pim3(3))
913 % -pim1(2)*(pim2(2)-pim3(2))-pim1(1)*(pim2(1)-pim3(1))
915 pd341=pim3(4)*(pim4(4)-pim1(4))-pim3(3)*(pim4(3)-pim1(3))
916 % -pim3(2)*(pim4(2)-pim1(2))-pim3(1)*(pim4(1)-pim1(1))
918 pd342=pim3(4)*(pim4(4)-pim2(4))-pim3(3)*(pim4(3)-pim2(3))
919 % -pim3(2)*(pim4(2)-pim2(2))-pim3(1)*(pim4(1)-pim2(1))
921 pd413=pim4(4)*(pim1(4)-pim3(4))-pim4(3)*(pim1(3)-pim3(3))
922 % -pim4(2)*(pim1(2)-pim3(2))-pim4(1)*(pim1(1)-pim3(1))
924 pd423=pim4(4)*(pim2(4)-pim3(4))-pim4(3)*(pim2(3)-pim3(3))
925 % -pim4(2)*(pim2(2)-pim3(2))-pim4(1)*(pim2(1)-pim3(1))
929 qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
930 % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
931 % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
932 % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
934 qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
935 % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
936 % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
937 % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
939 qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
940 % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
941 % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
942 % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
944 qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
945 % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
946 % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
947 % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
951 c define t(pi;pj,pk)= tijk:
953 t243=bwiga1(qmp2)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
954 t213=bwiga1(qmp2)*fpikm(sqrt(ps13),ampi,ampi)*gamma1
955 t143=bwiga1(qmp1)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
956 t123=bwiga1(qmp1)*fpikm(sqrt(ps23),ampi,ampi)*gamma1
957 t341=bwiga1(qmp3)*fpikm(sqrt(ps41),ampi,ampi)*gamma1
958 t342=bwiga1(qmp3)*fpikm(sqrt(ps42),ampi,ampi)*gamma1
960 c define s(i,j;k,l)= sijkl:
962 s2413=frho4(sqrt(ps24),ampi,ampi)*gamma2
963 s2314=frho4(sqrt(ps23),ampi,ampi)*bwigeps(ps14)*gamma2
964 s1423=frho4(sqrt(ps14),ampi,ampi)*gamma2
965 s1324=frho4(sqrt(ps13),ampi,ampi)*bwigeps(ps24)*gamma2
966 s34=frho4(sqrt(ps34),ampi,ampi)*gamma2
968 c definition of amplitude, first the [] brackets:
970 brack1=2.*t143+2.*t243+t123+t213
971 % +t341*(pd241/qmp3-1.)+t342*(pd142/qmp3-1.)
972 % +3./4.*(s1423+s2413-s2314-s1324)-3.*s34
974 brack2=2.*t143*pd243/qmp1+3.*t213
975 % +t123*(2.*pd423/qmp1+1.)+t341*(pd241/qmp3+3.)
976 % +t342*(pd142/qmp3+1.)
977 % -3./4.*(s2314+3.*s1324+3.*s1423+s2413)
979 brack3=2.*t243*pd143/qmp2+3.*t123
980 % +t213*(2.*pd413/qmp2+1.)+t341*(pd241/qmp3+1.)
981 % +t342*(pd142/qmp3+3.)
982 % -3./4.*(3.*s2314+s1324+s1423+3.*s2413)
984 brack4a=2.*t143*(pd243/qq*(qp1/qmp1+1.)+pd143/qq)
985 % +2.*t243*(pd143/qq*(qp2/qmp2+1.)+pd243/qq)
987 % +2.*t123*(pd423/qq*(qp1/qmp1+1.)+pd123/qq)
988 % +2.*t213*(pd413/qq*(qp2/qmp2+1.)+pd213/qq)
989 % +t341*(pd241/qmp3+1.-2.*pd241/qq*(qp3/qmp3+1.)
991 % +t342*(pd142/qmp3+1.-2.*pd142/qq*(qp3/qmp3+1.)
994 brack4b=-3./4.*(s2314*(2.*(qp2-qp3)/qq+1.)
995 % +s1324*(2.*(qp1-qp3)/qq+1.)
996 % +s1423*(2.*(qp1-qp4)/qq+1.)
997 % +s2413*(2.*(qp2-qp4)/qq+1.)
998 % +4.*s34*(qp4-qp3)/qq)
1000 brack4=brack4a+brack4b
1004 hcomp1(k)=(pim3(k)-pim4(k))*brack1
1005 hcomp2(k)=pim1(k)*brack2
1006 hcomp3(k)=pim2(k)*brack3
1007 hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1013 hadcur(i)=hcomp1(i)-hcomp2(i)-hcomp3(i)+hcomp4(i)
1014 hadcur(i)=-coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1019 c ---
END of the non omega current(3 possibilities)
1023 c --- there are two possibilities for omega current
1024 c --- pa pb are corresponding first and second pi-s
1030 c --- lorentz invariants
1043 IF (k.EQ.4) sign= 1.0
1044 qqa=qqa+sign*(paa(k)-pa(k))**2
1045 ss23=ss23+sign*(pb(k) +pim3(k))**2
1046 ss24=ss24+sign*(pb(k) +pim4(k))**2
1047 ss34=ss34+sign*(pim3(k)+pim4(k))**2
1048 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
1049 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
1050 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
1051 p1p2=p1p2+sign*pa(k)*pb(k)
1052 p1p3=p1p3+sign*pa(k)*pim3(k)
1053 p1p4=p1p4+sign*pa(k)*pim4(k)
1056 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
1057 c form3=bwign(qqa,amom,gamom)*(bwign(ss23,amro,gamro)+
1058 c $ bwign(ss24,amro,gamro)+bwign(ss34,amro,gamro))
1059 form3=bwign(qqa,amom,gamom)
1062 hadcur(k)=hadcur(k)+form2*form3*(
1063 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
1064 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
1065 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
1070 c ===================================================================
1071 c pi0 pi0 p0 pi-
CASE ====
1072 c ===================================================================
1073 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1076 c first definition of scalar products of momentum vectors
1078 c define(q-pi)**2 as qpi:
1080 qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1081 % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1083 qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1084 % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1086 qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1087 % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1089 qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1090 % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1093 c define(pi+pk)**2 as psik:
1095 ps14=(pim1(4)+pim4(4))**2-(pim1(3)+pim4(3))**2
1096 % -(pim1(2)+pim4(2))**2-(pim1(1)+pim4(1))**2
1098 ps21=(pim2(4)+pim1(4))**2-(pim2(3)+pim1(3))**2
1099 % -(pim2(2)+pim1(2))**2-(pim2(1)+pim1(1))**2
1101 ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1102 % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1104 ps24=(pim2(4)+pim4(4))**2-(pim2(3)+pim4(3))**2
1105 % -(pim2(2)+pim4(2))**2-(pim2(1)+pim4(1))**2
1107 ps31=(pim3(4)+pim1(4))**2-(pim3(3)+pim1(3))**2
1108 % -(pim3(2)+pim1(2))**2-(pim3(1)+pim1(1))**2
1110 ps34=(pim3(4)+pim4(4))**2-(pim3(3)+pim4(3))**2
1111 % -(pim3(2)+pim4(2))**2-(pim3(1)+pim4(1))**2
1115 pd324=pim3(4)*(pim2(4)-pim4(4))-pim3(3)*(pim2(3)-pim4(3))
1116 % -pim3(2)*(pim2(2)-pim4(2))-pim3(1)*(pim2(1)-pim4(1))
1118 pd314=pim3(4)*(pim1(4)-pim4(4))-pim3(3)*(pim1(3)-pim4(3))
1119 % -pim3(2)*(pim1(2)-pim4(2))-pim3(1)*(pim1(1)-pim4(1))
1121 pd234=pim2(4)*(pim3(4)-pim4(4))-pim2(3)*(pim3(3)-pim4(3))
1122 % -pim2(2)*(pim3(2)-pim4(2))-pim2(1)*(pim3(1)-pim4(1))
1124 pd214=pim2(4)*(pim1(4)-pim4(4))-pim2(3)*(pim1(3)-pim4(3))
1125 % -pim2(2)*(pim1(2)-pim4(2))-pim2(1)*(pim1(1)-pim4(1))
1127 pd134=pim1(4)*(pim3(4)-pim4(4))-pim1(3)*(pim3(3)-pim4(3))
1128 % -pim1(2)*(pim3(2)-pim4(2))-pim1(1)*(pim3(1)-pim4(1))
1130 pd124=pim1(4)*(pim2(4)-pim4(4))-pim1(3)*(pim2(3)-pim4(3))
1131 % -pim1(2)*(pim2(2)-pim4(2))-pim1(1)*(pim2(1)-pim4(1))
1133 c define q*pi = qpi:
1135 qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1136 % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1137 % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1138 % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1140 qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1141 % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1142 % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1143 % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1145 qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1146 % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1147 % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1148 % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1150 qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1151 % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1152 % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1153 % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1156 c define t(pi;pj,pk)= tijk:
1158 t324=bwiga1(qmp3)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1159 t314=bwiga1(qmp3)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1160 t234=bwiga1(qmp2)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1161 t214=bwiga1(qmp2)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1162 t134=bwiga1(qmp1)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1163 t124=bwiga1(qmp1)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1165 c define s(i,j;k,l)= sijkl:
1167 s1423=frho4(sqrt(ps14),ampi,ampi)*bwigeps(ps23)*gamma2
1168 s2431=frho4(sqrt(ps24),ampi,ampi)*bwigeps(ps31)*gamma2
1169 s3421=frho4(sqrt(ps34),ampi,ampi)*bwigeps(ps21)*gamma2
1172 c definition of amplitude, first the [] brackets:
1174 brack1=t234+t324+2.*t314+t134+2.*t214+t124
1175 % +t134*pd234/qmp1+t124*pd324/qmp1
1176 % -3./2.*(s3421+s2431+2.*s1423)
1179 brack2=t234*(1.+2.*pd134/qmp2)+3.*t324+3.*t124
1180 % +t134*(1.-pd234/qmp1)+2.*t214*pd314/qmp2
1182 % -3./2.*(s3421+3.*s2431)
1184 brack3=t324*(1.+2.*pd124/qmp3)+3.*t234+3.*t134
1185 % +t124*(1.-pd324/qmp1)+2.*t314*pd214/qmp3
1187 % -3./2.*(3.*s3421+s2431)
1189 brack4a=2.*t234*(1./2.+pd134/qq*(qp2/qmp2+1.)+pd234/qq)
1190 % +2.*t324*(1./2.+pd124/qq*(qp3/qmp3+1.)+pd324/qq)
1191 % +2.*t134*(1./2.+pd234/qq*(qp1/qmp1+1.)
1192 % -1./2.*pd234/qmp1+pd134/qq)
1193 % +2.*t124*(1./2.+pd324/qq*(qp1/qmp1+1.)
1194 % -1./2.*pd324/qmp1+pd124/qq)
1195 % +2.*t214*(pd314/qq*(qp2/qmp2+1.)+pd214/qq)
1196 % +2.*t314*(pd214/qq*(qp3/qmp3+1.)+pd314/qq)
1198 brack4b=-3./2.*(s3421*(2.*(qp3-qp4)/qq+1.)
1199 % +s2431*(2.*(qp2-qp4)/qq+1.)
1200 % +s1423*2.*(qp1-qp4)/qq)
1203 brack4=brack4a+brack4b
1207 hcomp1(k)=(pim1(k)-pim4(k))*brack1
1208 hcomp2(k)=pim2(k)*brack2
1209 hcomp3(k)=pim3(k)*brack3
1210 hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1216 hadcur(i)=hcomp1(i)+hcomp2(i)+hcomp3(i)-hcomp4(i)
1217 hadcur(i)=coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1223 c m. finkemeier et al.
END