C++ Interface to Tauola
tauola-BBB/prod/curr_cleo.f
1 
2  subroutine testresu
3 C this routine calculates Gamma_X/Gamma_e for tau decay into five
4 C masless pions via narrow a_2 and later omega +rho resonances
5 C (also narrow). Karlsruhe 16 Feb 2005. Test worked down to 1 % level
6 C to go beyond, rho presampler must be implemented in phase space.
7 C otherwise one can not get rid off the tails of distribution tails
8 C which remain with sizable effect,
9  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
10  * ,ampiz,ampi,amro,gamro,ama1,gama1
11  * ,amk,amkz,amkst,gamkst
12 C
13  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
14  * ,ampiz,ampi,amro,gamro,ama1,gama1
15  * ,amk,amkz,amkst,gamkst
16 C
17  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
18  real*4 gfermi,gv,ga,ccabib,scabib,gamel
19 
20  xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
21  DATA pi /3.141592653589793238462643/
22  amom=.782
23  gamom= 0.0085
24  ama2=1.260 ! for the test to work ama2> AMOM+ AMRO
25  gama2=0.001 !0.400
26 
27 
28  carb=3000
29  gropp=6
30  fomega=0.07
31  wyni= 1.0d0/4.0d0/amtau**2*(1-ama2**2/amtau**2)**2
32  wyni=wyni*(1+2*ama2**2/amtau**2)*carb**2*ama2**4/ama2/gama2*pi
33  dwpynd=xlam(ama2**2,amro**2,amom**2)/ama2
34  eff=dwpynd/ama2*(0.5*dwpynd**2*(1d0/amro**2+1d0/amom**2)+6)
35  wyni=wyni*eff
36  gam3pi= 1d0/3.0d0/128d0/(2*pi)**3*amom**7*
37  $ (fomega*gropp/amro**2)**2/120d0
38  gam2pi= gropp**2/48d0/pi*amro
39  wyni=wyni* gam3pi/gamom
40  wyni=wyni* gam2pi/gamro
41  wyni=wyni*ccabib**2
42  write(*,*) 'testresu=',wyni
43  end
44 
45  SUBROUTINE curr5(MNUM,PIM1,PIM2,PIM3,PIM4,PIM5,HADCUR)
46  REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PIM5(4)
47  COMPLEX HADCUR(4), HADCU(4)
48  IF (mnum.EQ.24) THEN ! --+00
49  CALL curr5x(mnum,pim1,pim2,pim3,pim4,pim5,hadcu)
50  DO k=1,4
51  hadcur(k)=hadcu(k)
52  ENDDO
53  ELSEIF (mnum.EQ.26) THEN ! +--00
54  CALL curr5x(mnum,pim2,pim3,pim1,pim4,pim5,hadcu)
55  DO k=1,4
56  hadcur(k)=hadcu(k)
57  ENDDO
58  CALL curr5x(mnum,pim4,pim5,pim2,pim3,pim1,hadcu)
59  DO k=1,4
60  hadcur(k)=hadcur(k)+hadcu(k)
61  ENDDO
62  CALL curr5x(mnum,pim4,pim5,pim3,pim2,pim1,hadcu)
63  DO k=1,4
64  hadcur(k)=hadcur(k)+hadcu(k)
65  ENDDO
66 C --
67  CALL curr5x(24,pim2,pim3,pim1,pim4,pim5,hadcu)
68  DO k=1,4
69  hadcur(k)=hadcu(k) +hadcur(k)
70  ENDDO
71  CALL curr5x(24,pim2,pim3,pim1,pim5,pim4,hadcu)
72  DO k=1,4
73  hadcur(k)=hadcur(k)+hadcu(k)
74  ENDDO
75  CALL curr5x(24,pim3,pim2,pim1,pim4,pim5,hadcu)
76  DO k=1,4
77  hadcur(k)=hadcur(k)+hadcu(k)
78  ENDDO
79  CALL curr5x(24,pim3,pim2,pim1,pim5,pim4,hadcu)
80  DO k=1,4
81  hadcur(k)=hadcur(k)+hadcu(k)
82  ENDDO
83 C --
84  DO k=1,4
85  hadcur(k)=hadcur(k)*sqrt(0.25) ! statistical factor
86  ENDDO
87  ELSEIF (mnum.EQ.27) THEN ! -0000
88  CALL curr5x(mnum,pim2,pim3,pim1,pim4,pim5,hadcu)
89  DO k=1,4
90  hadcur(k)=hadcu(k)
91  ENDDO
92  CALL curr5x(mnum,pim2,pim4,pim1,pim3,pim5,hadcu)
93  DO k=1,4
94  hadcur(k)=hadcur(k)+hadcu(k)
95  ENDDO
96  CALL curr5x(mnum,pim2,pim5,pim1,pim3,pim4,hadcu)
97  DO k=1,4
98  hadcur(k)=hadcur(k)+hadcu(k)
99  ENDDO
100  CALL curr5x(mnum,pim4,pim3,pim1,pim2,pim5,hadcu)
101  DO k=1,4
102  hadcur(k)=hadcur(k)+hadcu(k)
103  ENDDO
104  CALL curr5x(mnum,pim5,pim3,pim1,pim4,pim2,hadcu)
105  DO k=1,4
106  hadcur(k)=hadcur(k)+hadcu(k)
107  ENDDO
108  CALL curr5x(mnum,pim5,pim4,pim1,pim3,pim2,hadcu)
109  DO k=1,4
110  hadcur(k)=hadcur(k)+hadcu(k)
111  ENDDO
112  DO k=1,4
113  hadcur(k)=hadcur(k)*sqrt(1.0/24.0) ! statistical factor
114  ENDDO
115  ELSEIF (mnum.EQ.28) THEN ! -++--
116  CALL curr5x(mnum,pim4,pim5,pim2,pim3,pim1,hadcu)
117  DO k=1,4
118  hadcur(k)=hadcu(k)
119  ENDDO
120  CALL curr5x(mnum,pim1,pim5,pim2,pim3,pim4,hadcu)
121  DO k=1,4
122  hadcur(k)=hadcur(k)+hadcu(k)
123  ENDDO
124  CALL curr5x(mnum,pim1,pim4,pim2,pim3,pim5,hadcu)
125  DO k=1,4
126  hadcur(k)=hadcur(k)+hadcu(k)
127  ENDDO
128  CALL curr5x(mnum,pim4,pim5,pim3,pim2,pim1,hadcu)
129  DO k=1,4
130  hadcur(k)=hadcur(k)+hadcu(k)
131  ENDDO
132  CALL curr5x(mnum,pim1,pim5,pim3,pim2,pim4,hadcu)
133  DO k=1,4
134  hadcur(k)=hadcur(k)+hadcu(k)
135  ENDDO
136  CALL curr5x(mnum,pim1,pim4,pim3,pim2,pim5,hadcu)
137  DO k=1,4
138  hadcur(k)=hadcur(k)+hadcu(k)
139  ENDDO
140  DO k=1,4
141  hadcur(k)=hadcur(k)*sqrt(1.0/12.0) ! statistical factor
142  ENDDO
143 
144  ELSE
145  CALL curr5x(mnum,pim1,pim2,pim3,pim4,pim5,hadcu)
146  DO k=1,4
147  hadcur(k)=hadcu(k)
148  ENDDO
149  ENDIF
150  END
151 *AJW 1 version of CURR from KORALB.
152  SUBROUTINE curr5x(MNUM,PIM1,PIM2,PIM3,PIM4,PIM5,HADCUR)
153 C ==================================================================
154 C ZBW, 02/2005 - prototype current for 5 pione, several options.
155 C ==================================================================
156 
157  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
158  * ,ampiz,ampi,amro,gamro,ama1,gama1
159  * ,amk,amkz,amkst,gamkst
160 C
161  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
162  * ,ampiz,ampi,amro,gamro,ama1,gama1
163  * ,amk,amkz,amkst,gamkst
164 C
165  REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PIM5(4)
166  COMPLEX HADCUR(4)
167 
168  INTEGER K,L,MNUM,K1,K2,IRO,I,J,KK
169  REAL PA(4),PB(4),PAA(4),PC(4),PD(4)
170  REAL AA(4,4),PP(4,4)
171  REAL A,XM,XG,G1,G2,G,AMRO2,GAMRO2,AMRO3,GAMRO3,AMOM,GAMOM
172  REAL FRO,COEF1,FPI,COEF2,QQ,SK,DENOM,SIG,QQA,SS23,SS24,SS34,QP1P2
173  REAL QP1P3,QP1P4,P1P2,P1P3,P1P4,SIGN
174  REAL PKORB,AMPA
175  COMPLEX ALF0,ALF1,ALF2,ALF3
176  COMPLEX LAM0,LAM1,LAM2,LAM3
177  COMPLEX BET1,BET2,BET3
178  COMPLEX FORM1,FORM2,FORM3,FORM4,FORM2PI
179  COMPLEX BWIGM,WIGFOR,FPIKM,FPIKMD
180  COMPLEX AMPL(7),AMPR
181  COMPLEX BWIGN
182 C
183  DATA pi /3.141592653589793238462643/
184  bwign(a,xm,xg)=xm**2/cmplx(a-xm**2,xm*xg)
185 C*******************************************************************************
186  sa1=(pim1(4)+pim2(4)+pim3(4)+pim4(4)+pim5(4))**2
187  $ -(pim1(3)+pim2(3)+pim3(3)+pim4(3)+pim5(3))**2
188  $ -(pim1(2)+pim2(2)+pim3(2)+pim4(2)+pim5(2))**2
189  $ -(pim1(1)+pim2(1)+pim3(1)+pim4(1)+pim5(1))**2
190 
191 C
192  IF (mnum.EQ.24) THEN ! simple, semi realistic current for
193  ! pi+pi+pi-pi0pi0, saturated with a2 --> rho omega
194 
195  somega=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
196  $ -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
197  sp= (pim2(4)+pim4(4))**2-(pim2(3)+pim4(3))**2
198  $ -(pim2(2)+pim4(2))**2-(pim2(1)+pim4(1))**2
199  sm= (pim3(4)+pim4(4))**2-(pim3(3)+pim4(3))**2
200  $ -(pim3(2)+pim4(2))**2-(pim3(1)+pim4(1))**2
201  s0= (pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
202  $ -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
203 
204 
205 
206  srho=(pim1(4)+pim5(4))**2-(pim1(3)+pim5(3))**2
207  $ -(pim1(2)+pim5(2))**2-(pim1(1)+pim5(1))**2
208 
209  DO k=1,4
210  hadcur(k)=cmplx(0.0)
211  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)+pim5(k)
212  pa(k) =pim2(k)+pim3(k)+pim4(k)
213  pb(k)=pim1(k)-pim5(k)
214  ENDDO
215  CALL levici(pc,pim2,pim3,pim4)
216  CALL levici(pd,pb,pc,paa)
217 ! write(*,*) sqrt(somega),
218 
219  amom=.782
220  gamom= 0.0085
221  ama2=1.260
222  gama2=0.400
223 
224 
225  carb=3000
226  gropp=6
227  fomega=0.07
228  coef1=carb/amro**2/amom**2* gropp* (fomega*gropp/amro**2)
229  DO k=1,4
230  hadcur(k)=coef1*pd(k)
231  hadcur(k)=hadcur(k)*bwign(somega,amom,gamom)
232  $ *bwign(srho,amro,gamro)*bwign(sa1,ama2,gama2)
233 ! $ *(BWIGN(SP,AMRO,GAMRO)+BWIGN(SM,AMRO,GAMRO)+BWIGN(S0,AMRO,GAMRO))
234 ! write(*,*) sqrt(somega)-amom,amom
235  ENDDO
236  ELSEIF (mnum.EQ.26.OR.mnum.EQ.27.OR.mnum.EQ.28) THEN ! simple, semi realistic current for
237  ! pi-pi-pi+pi0pi0, saturated with a2 --> f0 a2
238 
239  DO k=1,4
240  hadcur(k)=cmplx(0.0)
241  ENDDO
242 
243  sa2=(pim2(4)+pim3(4)+pim1(4))**2-(pim2(3)+pim3(3)+pim1(3))**2
244  $ -(pim2(2)+pim3(2)+pim1(2))**2-(pim2(1)+pim3(1)+pim1(1))**2
245  s2= (pim1(4)+pim3(4))**2-(pim1(3)+pim3(3))**2
246  $ -(pim1(2)+pim3(2))**2-(pim1(1)+pim3(1))**2
247  s1= (pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
248  $ -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
249  s2x13=pim2(4)*(pim1(4)-pim3(4))-pim2(3)*(pim1(3)-pim3(3))
250  $ -pim2(2)*(pim1(2)-pim3(2))-pim2(1)*(pim1(1)-pim3(1))
251  s1x23=pim1(4)*(pim2(4)-pim3(4))-pim1(3)*(pim2(3)-pim3(3))
252  $ -pim1(2)*(pim2(2)-pim3(2))-pim1(1)*(pim2(1)-pim3(1))
253 
254 
255 
256  sf=(pim4(4)+pim5(4))**2-(pim4(3)+pim5(3))**2
257  $ -(pim4(2)+pim5(2))**2-(pim4(1)+pim5(1))**2
258 
259  DO k=1,4
260  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)+pim5(k)
261  pa(k) =pim1(k)+pim2(k)+pim3(k)
262  pb(k)=pa(k)*s2x13/sa2-(pim1(k)-pim3(k))
263  pc(k)=pa(k)*s1x23/sa2-(pim2(k)-pim3(k))
264  ENDDO
265  papb=paa(4)*pb(4)-paa(3)*pb(3)-paa(2)*pb(2)-paa(1)*pb(1)
266  papc=paa(4)*pc(4)-paa(3)*pc(3)-paa(2)*pc(2)-paa(1)*pc(1)
267  DO k=1,4
268  hadcur(k)=hadcur(k)+(paa(k)*papb/sa1-pb(k))*bwign(s2,amro,gamro)
269  hadcur(k)=hadcur(k)+(paa(k)*papc/sa1-pc(k))*bwign(s1,amro,gamro)
270  ENDDO
271 ! write(*,*) sqrt(somega),
272 
273  amf2=.800
274  gamf2=.600
275  ama2=1.260
276  gama2=.400
277 
278 
279  carb=4.
280  faaf=4.
281  fpp=5.
282  grorop=6.
283  gropp=6.
284  coef1=carb/ama2**4/amf2**2/amro**2*faaf*fpp* grorop* gropp
285  DO k=1,4
286  hadcur(k)=coef1*hadcur(k)
287  hadcur(k)=hadcur(k)*bwign(sf,amf2,gamf2)
288  $ *bwign(sa2,ama2,gama2)*bwign(sa1,ama2,gama2)
289  ENDDO
290 
291  ELSE !MNUM! not realistic current, for tests only !
292 
293  ! write(*,*) coef1
294  fpi=93.3e-3
295  coef1=2*2.0*sqrt(3.0)/fpi**3
296  coef1= 1d0/amtau**3 *(4*3*2*1)* (4*pi)**3* sqrt(20.0d0) ! this
297  ! normalization gives Gamma/Gamma_e=ccabib**2 in masless limit
298  ! Benchmark current !1
299  ! write(*,*) coef1
300  ! stop
301  DO k=1,4
302  hadcur(k)=cmplx(0.0)
303  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)+pim5(k)
304  hadcur(k)=coef1*paa(k)
305  ENDDO
306  ENDIF !MNUM!
307  END
308  SUBROUTINE levici(P,A,B,C)
309  REAL P(4),A(4),B(4),C(4)
310 
311  p(1)=a(2)*b(3)*c(4)+a(3)*b(4)*c(2)+a(4)*b(2)*c(3)
312  $ -a(2)*b(4)*c(3)-a(4)*b(3)*c(2)-a(3)*b(2)*c(4)
313 
314  p(2)=a(1)*b(4)*c(3)+a(3)*b(1)*c(4)+a(4)*b(3)*c(1)
315  $ -a(1)*b(3)*c(4)-a(4)*b(1)*c(3)-a(3)*b(4)*c(1)
316 
317  p(3)=a(1)*b(2)*c(4)+a(4)*b(1)*c(2)+a(2)*b(4)*c(1)
318  $ -a(1)*b(4)*c(2)-a(2)*b(1)*c(4)-a(4)*b(2)*c(1)
319 
320  p(4)=a(1)*b(3)*c(2)+a(2)*b(1)*c(3)+a(3)*b(2)*c(1)
321  $ -a(1)*b(2)*c(3)-a(3)*b(1)*c(2)-a(2)*b(3)*c(1)
322 
323  p(1)=-p(1)
324  p(2)=-p(2)
325  p(3)=-p(3)
326  END
327 
328 C*AJW 1 version of CURR from KORALB.
329  SUBROUTINE curr_cleo(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
330 C ==================================================================
331 C AJW, 11/97 - based on original CURR from TAUOLA:
332 C hadronic current for 4 pi final state
333 C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
334 C R. Decker Z. Phys C36 (1987) 487.
335 C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
336 C BUT, rewritten to be more general and less "theoretical",
337 C using parameters tuned by Vasia and DSC.
338 C ==================================================================
339 
340  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
341  * ,ampiz,ampi,amro,gamro,ama1,gama1
342  * ,amk,amkz,amkst,gamkst
343 C
344  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
345  * ,ampiz,ampi,amro,gamro,ama1,gama1
346  * ,amk,amkz,amkst,gamkst
347 C
348  REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4)
349  COMPLEX HADCUR(4)
350 
351  INTEGER K,L,MNUM,K1,K2,IRO,I,J,KK
352  REAL PA(4),PB(4),PAA(4)
353  REAL AA(4,4),PP(4,4)
354  REAL A,XM,XG,G1,G2,G,AMRO2,GAMRO2,AMRO3,GAMRO3,AMOM,GAMOM
355  REAL FRO,COEF1,FPI,COEF2,QQ,SK,DENOM,SIG,QQA,SS23,SS24,SS34,QP1P2
356  REAL QP1P3,QP1P4,P1P2,P1P3,P1P4,SIGN
357  REAL PKORB,AMPA
358  COMPLEX ALF0,ALF1,ALF2,ALF3
359  COMPLEX LAM0,LAM1,LAM2,LAM3
360  COMPLEX BET1,BET2,BET3
361  COMPLEX FORM1,FORM2,FORM3,FORM4,FORM2PI
362  COMPLEX BWIGM,WIGFOR,FPIKM,FPIKMD
363  COMPLEX AMPL(7),AMPR
364  COMPLEX BWIGN
365 C
366  bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
367 C*******************************************************************************
368 C
369 C --- masses and constants
370  IF (g1.NE.12.924) THEN
371  g1=12.924
372  g2=1475.98
373  fpi=93.3e-3
374  g =g1*g2
375  fro=0.266*amro**2
376  coef1=2.0*sqrt(3.0)/fpi**2
377  coef2=fro*g ! overall constant for the omega current
378  coef2= coef2*0.56 ! factor 0.56 reduces contribution of omega from 68% to 40 %
379 
380 C masses and widths for for rho-prim and rho-bis:
381  amro2 = 1.465
382  gamro2= 0.310
383  amro3=1.700
384  gamro3=0.235
385 C
386  amom = pkorb(1,14)
387  gamom = pkorb(2,14)
388  amro2 = pkorb(1,21)
389  gamro2= pkorb(2,21)
390  amro3 = pkorb(1,22)
391  gamro3= pkorb(2,22)
392 C
393 C Amplitudes for (pi-pi-pi0pi+) -> PS, rho0, rho-, rho+, omega.
394  ampl(1) = cmplx(pkorb(3,31)*coef1,0.)
395  ampl(2) = cmplx(pkorb(3,32)*coef1,0.)*cexp(cmplx(0.,pkorb(3,42)))
396  ampl(3) = cmplx(pkorb(3,33)*coef1,0.)*cexp(cmplx(0.,pkorb(3,43)))
397  ampl(4) = cmplx(pkorb(3,34)*coef1,0.)*cexp(cmplx(0.,pkorb(3,44)))
398  ampl(5) = cmplx(pkorb(3,35)*coef2,0.)*cexp(cmplx(0.,pkorb(3,45)))
399 C Amplitudes for (pi0pi0pi0pi-) -> PS, rho-.
400  ampl(6) = cmplx(pkorb(3,36)*coef1)
401  ampl(7) = cmplx(pkorb(3,37)*coef1)
402 C
403 C rho' contributions to rho' -> pi-omega:
404  alf0 = cmplx(pkorb(3,51),0.0)
405  alf1 = cmplx(pkorb(3,52)*amro**2,0.0)
406  alf2 = cmplx(pkorb(3,53)*amro2**2,0.0)
407  alf3 = cmplx(pkorb(3,54)*amro3**2,0.0)
408 C rho' contribtions to rho' -> rhopipi:
409  lam0 = cmplx(pkorb(3,55),0.0)
410  lam1 = cmplx(pkorb(3,56)*amro**2,0.0)
411  lam2 = cmplx(pkorb(3,57)*amro2**2,0.0)
412  lam3 = cmplx(pkorb(3,58)*amro3**2,0.0)
413 C rho contributions to rhopipi, rho -> 2pi:
414  bet1 = cmplx(pkorb(3,59)*amro**2,0.0)
415  bet2 = cmplx(pkorb(3,60)*amro2**2,0.0)
416  bet3 = cmplx(pkorb(3,61)*amro3**2,0.0)
417 C
418  END IF
419 C**************************************************
420 C
421 C --- initialization of four vectors
422  DO 7 k=1,4
423  DO 8 l=1,4
424  8 aa(k,l)=0.0
425  hadcur(k)=cmplx(0.0)
426  paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
427  pp(1,k)=pim1(k)
428  pp(2,k)=pim2(k)
429  pp(3,k)=pim3(k)
430  7 pp(4,k)=pim4(k)
431 C
432 ! IF (mnum.gt.2) write(*,*) 'curr cleo mnum=',mnum
433 ! IF (MNUM.gt.2) goto 389
434  IF (mnum.EQ.1) THEN
435 C ===================================================================
436 C pi- pi- p0 pi+ case ====
437 C ===================================================================
438  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
439 
440 C Add M(4pi)-dependence to rhopipi channels:
441  form4= lam0+lam1*bwign(qq,amro,gamro)
442  * +lam2*bwign(qq,amro2,gamro2)
443  * +lam3*bwign(qq,amro3,gamro3)
444 
445 C --- loop over five contributions of the rho-pi-pi
446  DO 201 k1=1,3
447  DO 201 k2=3,4
448 C
449  IF (k2.EQ.k1) THEN
450  GOTO 201
451  ELSEIF (k2.EQ.3) THEN
452 C rho-
453  ampr = ampl(3)
454  ampa = ampiz
455  ELSEIF (k1.EQ.3) THEN
456 C rho+
457  ampr = ampl(4)
458  ampa = ampiz
459  ELSE
460 C rho0
461  ampr = ampl(2)
462  ampa = ampi
463  END IF
464 C
465  sk=(pp(k1,4)+pp(k2,4))**2-(pp(k1,3)+pp(k2,3))**2
466  $ -(pp(k1,2)+pp(k2,2))**2-(pp(k1,1)+pp(k2,1))**2
467 
468 C -- definition of AA matrix
469 C -- cronecker delta
470  DO 202 i=1,4
471  DO 203 j=1,4
472  203 aa(i,j)=0.0
473  202 aa(i,i)=1.0
474 C ... and the rest ...
475  DO 204 l=1,4
476  IF (l.NE.k1.AND.l.NE.k2) THEN
477  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
478  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
479  DO 205 i=1,4
480  DO 205 j=1,4
481  sig= 1.0
482  IF(j.NE.4) sig=-sig
483  aa(i,j)=aa(i,j)
484  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
485  205 CONTINUE
486  ENDIF
487  204 CONTINUE
488 C
489 C --- lets add something to HADCURR
490 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
491 C FORM1= AMPL(1)+AMPR*FPIKM(SQRT(SK),AMPI,AMPI)
492 
493  form2pi= bet1*bwigm(sk,amro,gamro,ampa,ampi)
494  1 +bet2*bwigm(sk,amro2,gamro2,ampa,ampi)
495  2 +bet3*bwigm(sk,amro3,gamro3,ampa,ampi)
496  form1= ampl(1)+ampr*form2pi
497 C
498  DO 206 i=1,4
499  DO 206 j=1,4
500  hadcur(i)=hadcur(i)+form1*form4*aa(i,j)*(pp(k1,j)-pp(k2,j))
501  206 CONTINUE
502 C --- end of the rho-pi-pi current (5 possibilities)
503  201 CONTINUE
504 C
505 C ===================================================================
506 C Now modify the coefficient for the omega-pi current: =
507 C ===================================================================
508  IF (ampl(5).EQ.cmplx(0.,0.)) GOTO 311
509 
510 C Overall rho+rhoprime for the 4pi system:
511 C FORM2=AMPL(5)*(BWIGN(QQ,AMRO,GAMRO)+ELPHA*BWIGN(QQ,AMROP,GAMROP))
512 C Modified M(4pi)-dependence:
513  form2=ampl(5)*(alf0+alf1*bwign(qq,amro,gamro)
514  * +alf2*bwign(qq,amro2,gamro2)
515  * +alf3*bwign(qq,amro3,gamro3))
516 C
517 C --- there are two possibilities for omega current
518 C --- PA PB are corresponding first and second pi-s
519  DO 301 kk=1,2
520  DO 302 i=1,4
521  pa(i)=pp(kk,i)
522  pb(i)=pp(3-kk,i)
523  302 CONTINUE
524 C --- lorentz invariants
525  qqa=0.0
526  ss23=0.0
527  ss24=0.0
528  ss34=0.0
529  qp1p2=0.0
530  qp1p3=0.0
531  qp1p4=0.0
532  p1p2 =0.0
533  p1p3 =0.0
534  p1p4 =0.0
535  DO 303 k=1,4
536  sign=-1.0
537  IF (k.EQ.4) sign= 1.0
538  qqa=qqa+sign*(paa(k)-pa(k))**2
539  ss23=ss23+sign*(pb(k) +pim3(k))**2
540  ss24=ss24+sign*(pb(k) +pim4(k))**2
541  ss34=ss34+sign*(pim3(k)+pim4(k))**2
542  qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
543  qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
544  qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
545  p1p2=p1p2+sign*pa(k)*pb(k)
546  p1p3=p1p3+sign*pa(k)*pim3(k)
547  p1p4=p1p4+sign*pa(k)*pim4(k)
548  303 CONTINUE
549 C
550 C omega -> rho pi for the 3pi system:
551 C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
552 C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
553 C No omega -> rho pi; just straight omega:
554  form3=bwign(qqa,amom,gamom)
555 C
556  DO 304 k=1,4
557  hadcur(k)=hadcur(k)+form2*form3*(
558  $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
559  $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
560  $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
561  304 CONTINUE
562  301 CONTINUE
563  311 CONTINUE
564 C
565  ELSE
566 C ===================================================================
567 C pi0 pi0 p0 pi- case ====
568 C ===================================================================
569 ! 389 continue ! temporary solution for `new' 4-pi modes as well.
570  qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
571 
572 C --- loop over three contribution of the non-omega current
573  DO 101 k=1,3
574  sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
575  $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
576 
577 C -- definition of AA matrix
578 C -- cronecker delta
579  DO 102 i=1,4
580  DO 103 j=1,4
581  103 aa(i,j)=0.0
582  102 aa(i,i)=1.0
583 C
584 C ... and the rest ...
585  DO 104 l=1,3
586  IF (l.NE.k) THEN
587  denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
588  $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
589  DO 105 i=1,4
590  DO 105 j=1,4
591  sig=1.0
592  IF(j.NE.4) sig=-sig
593  aa(i,j)=aa(i,j)
594  $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
595  105 CONTINUE
596  ENDIF
597  104 CONTINUE
598 
599 C --- lets add something to HADCURR
600 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
601 CCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
602 C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
603  form1 = ampl(6)+ampl(7)*fpikm(sqrt(sk),ampi,ampi)
604 
605  DO 106 i=1,4
606  DO 106 j=1,4
607  hadcur(i)=hadcur(i)+form1*aa(i,j)*(pp(k,j)-pp(4,j))
608  106 CONTINUE
609 C --- end of the non omega current (3 possibilities)
610  101 CONTINUE
611 
612  ENDIF
613  END
614 
615 
616