C++ Interface to Tauola
tauola-BBB/curr_cleo.F
1
2 subroutine testresu
3C this routine calculates Gamma_X/Gamma_e for tau decay into five
4C masless pions via narrow a_2 and later omega +rho resonances
5C (also narrow). Karlsruhe 16 Feb 2005. Test worked down to 1 % level
6C to go beyond, rho presampler must be implemented in phase space.
7C otherwise one can not get rid off the tails of distribution tails
8C 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
12C
13 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
14 * ,ampiz,ampi,amro,gamro,ama1,gama1
15 * ,amk,amkz,amkst,gamkst
16C
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
66C --
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
83C --
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)
153C ==================================================================
154C ZBW, 02/2005 - prototype current for 5 pione, several options.
155C ==================================================================
156
157 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
158 * ,ampiz,ampi,amro,gamro,ama1,gama1
159 * ,amk,amkz,amkst,gamkst
160C
161 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
162 * ,ampiz,ampi,amro,gamro,ama1,gama1
163 * ,amk,amkz,amkst,gamkst
164C
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
182C
183 DATA pi /3.141592653589793238462643/
184 bwign(a,xm,xg)=xm**2/cmplx(a-xm**2,xm*xg)
185C*******************************************************************************
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
191C
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
328C*AJW 1 version of CURR from KORALB.
329 SUBROUTINE curr_cleo(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
330C ==================================================================
331C AJW, 11/97 - based on original CURR from TAUOLA:
332C hadronic current for 4 pi final state
333C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
334C R. Decker Z. Phys C36 (1987) 487.
335C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
336C BUT, rewritten to be more general and less "theoretical",
337C using parameters tuned by Vasia and DSC.
338C ==================================================================
339
340 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
341 * ,ampiz,ampi,amro,gamro,ama1,gama1
342 * ,amk,amkz,amkst,gamkst
343C
344 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
345 * ,ampiz,ampi,amro,gamro,ama1,gama1
346 * ,amk,amkz,amkst,gamkst
347C
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
365C
366 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
367C*******************************************************************************
368C
369C --- 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
380C 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
385C
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)
392C
393C 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)))
399C Amplitudes for (pi0pi0pi0pi-) -> PS, rho-.
400 ampl(6) = cmplx(pkorb(3,36)*coef1)
401 ampl(7) = cmplx(pkorb(3,37)*coef1)
402C
403C 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)
408C 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)
413C 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)
417C
418 END IF
419C**************************************************
420C
421C --- 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)
431C
432! IF (mnum.gt.2) write(*,*) 'curr cleo mnum=',mnum
433! IF (MNUM.gt.2) goto 389
434 IF (mnum.EQ.1) THEN
435C ===================================================================
436C pi- pi- p0 pi+ case ====
437C ===================================================================
438 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
439
440C 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
445C --- loop over five contributions of the rho-pi-pi
446 DO 201 k1=1,3
447 DO 201 k2=3,4
448C
449 IF (k2.EQ.k1) THEN
450 GOTO 201
451 ELSEIF (k2.EQ.3) THEN
452C rho-
453 ampr = ampl(3)
454 ampa = ampiz
455 ELSEIF (k1.EQ.3) THEN
456C rho+
457 ampr = ampl(4)
458 ampa = ampiz
459 ELSE
460C rho0
461 ampr = ampl(2)
462 ampa = ampi
463 END IF
464C
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
468C -- definition of AA matrix
469C -- 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
474C ... 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
488C
489C --- lets add something to HADCURR
490C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
491C 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
497C
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
502C --- end of the rho-pi-pi current (5 possibilities)
503 201 CONTINUE
504C
505C ===================================================================
506C Now modify the coefficient for the omega-pi current: =
507C ===================================================================
508 IF (ampl(5).EQ.cmplx(0.,0.)) GOTO 311
509
510C Overall rho+rhoprime for the 4pi system:
511C FORM2=AMPL(5)*(BWIGN(QQ,AMRO,GAMRO)+ELPHA*BWIGN(QQ,AMROP,GAMROP))
512C 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))
516C
517C --- there are two possibilities for omega current
518C --- 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
524C --- 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
549C
550C omega -> rho pi for the 3pi system:
551C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
552C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
553C No omega -> rho pi; just straight omega:
554 form3=bwign(qqa,amom,gamom)
555C
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
564C
565 ELSE
566C ===================================================================
567C pi0 pi0 p0 pi- case ====
568C ===================================================================
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
572C --- 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
577C -- definition of AA matrix
578C -- 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
583C
584C ... 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
599C --- lets add something to HADCURR
600C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
601CCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
602C 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
609C --- end of the non omega current (3 possibilities)
610 101 CONTINUE
611
612 ENDIF
613 END
614
615
616