C++ Interface to Tauola
tauola/curr_cleo.f
1/* copyright(c) 1991-2023 free software foundation, inc.
2 this file is part of the gnu c library.
3
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.
8
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.
13
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 <https://www.gnu.org/licenses/>. */
17
18
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. */
26
27/* glibc's intent is to support the IEC 559 math functionality, real
28 and complex. If the GCC (4.9 and later) predefined macros
29 specifying compiler intent are available, use them to determine
30 whether the overall intent is to support these features; otherwise,
31 presume an older compiler has intent to support these features and
32 define these macros by default. */
33
34
35
36/* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is
37 synchronized with ISO/IEC 10646:2017, fifth edition, plus
38 the following additions from Amendment 1 to the fifth edition:
39 - 56 emoji characters
40 - 285 hentaigana
41 - 3 additional Zanabazar Square characters */
42
43*AJW 1 version of CURR from KORALB.
44 SUBROUTINE CURR_CLEO(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
45C ==================================================================
46C AJW, 11/97 - based on original CURR from TAUOLA:
47C hadronic current for 4 pi final state
48C R. Fisher, J. Wess and F. Wagner Z. Phys C3 (1980) 313
49C R. Decker Z. Phys C36 (1987) 487.
50C M. Gell-Mann, D. Sharp, W. Wagner Phys. Rev. Lett 8 (1962) 261.
51C BUT, rewritten to be more general and less "theoretical",
52C using parameters tuned by Vasia and DSC.
53C ==================================================================
54
55 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
56 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
57 * ,AMK,AMKZ,AMKST,GAMKST
58C
59 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
60 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
61 * ,AMK,AMKZ,AMKST,GAMKST
62C
63 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4)
64 COMPLEX HADCUR(4)
65
66 INTEGER K,L,MNUM,K1,K2,IRO,I,J,KK
67 REAL PA(4),PB(4),PAA(4)
68 REAL AA(4,4),PP(4,4)
69 REAL A,XM,XG,G1,G2,G,AMRO2,GAMRO2,AMRO3,GAMRO3,AMOM,GAMOM
70 REAL FRO,COEF1,FPI,COEF2,QQ,SK,DENOM,SIG,QQA,SS23,SS24,SS34,QP1P2
71 REAL QP1P3,QP1P4,P1P2,P1P3,P1P4,SIGN
72 REAL PKORB,AMPA
73 COMPLEX ALF0,ALF1,ALF2,ALF3
74 COMPLEX LAM0,LAM1,LAM2,LAM3
75 COMPLEX BET1,BET2,BET3
76 COMPLEX FORM1,FORM2,FORM3,FORM4,FORM2PI
77 COMPLEX BWIGM,WIGFOR,FPIKM,FPIKMD
78 COMPLEX AMPL(7),AMPR
79 COMPLEX BWIGN
80C
81 BWIGN(A,XM,XG)=1.0/CMPLX(A-XM**2,XM*XG)
82C*******************************************************************************
83C
84C --- masses and constants
85.NE. IF (G112.924) THEN
86 G1=12.924
87 G2=1475.98
88 FPI=93.3E-3
89 G =G1*G2
90 FRO=0.266*AMRO**2
91 COEF1=2.0*SQRT(3.0)/FPI**2
92 COEF2=FRO*G ! overall constant for the omega current
93 COEF2= COEF2*0.56 ! factor 0.56 reduces contribution of omega from 68% to 40 %
94
95C masses and widths for for rho-prim and rho-bis:
96 AMRO2 = 1.465
97 GAMRO2= 0.310
98 AMRO3=1.700
99 GAMRO3=0.235
100C
101 AMOM = PKORB(1,14)
102 GAMOM = PKORB(2,14)
103 AMRO2 = PKORB(1,21)
104 GAMRO2= PKORB(2,21)
105 AMRO3 = PKORB(1,22)
106 GAMRO3= PKORB(2,22)
107C
108C Amplitudes for (pi-pi-pi0pi+) -> PS, rho0, rho-, rho+, omega.
109 AMPL(1) = CMPLX(PKORB(3,31)*COEF1,0.)
110 AMPL(2) = CMPLX(PKORB(3,32)*COEF1,0.)*CEXP(CMPLX(0.,PKORB(3,42)))
111 AMPL(3) = CMPLX(PKORB(3,33)*COEF1,0.)*CEXP(CMPLX(0.,PKORB(3,43)))
112 AMPL(4) = CMPLX(PKORB(3,34)*COEF1,0.)*CEXP(CMPLX(0.,PKORB(3,44)))
113 AMPL(5) = CMPLX(PKORB(3,35)*COEF2,0.)*CEXP(CMPLX(0.,PKORB(3,45)))
114C Amplitudes for (pi0pi0pi0pi-) -> PS, rho-.
115 AMPL(6) = CMPLX(PKORB(3,36)*COEF1)
116 AMPL(7) = CMPLX(PKORB(3,37)*COEF1)
117C
118C rho' contributions to rho' -> pi-omega:
119 ALF0 = CMPLX(PKORB(3,51),0.0)
120 ALF1 = CMPLX(PKORB(3,52)*AMRO**2,0.0)
121 ALF2 = CMPLX(PKORB(3,53)*AMRO2**2,0.0)
122 ALF3 = CMPLX(PKORB(3,54)*AMRO3**2,0.0)
123C rho' contribtions to rho' -> rhopipi:
124 LAM0 = CMPLX(PKORB(3,55),0.0)
125 LAM1 = CMPLX(PKORB(3,56)*AMRO**2,0.0)
126 LAM2 = CMPLX(PKORB(3,57)*AMRO2**2,0.0)
127 LAM3 = CMPLX(PKORB(3,58)*AMRO3**2,0.0)
128C rho contributions to rhopipi, rho -> 2pi:
129 BET1 = CMPLX(PKORB(3,59)*AMRO**2,0.0)
130 BET2 = CMPLX(PKORB(3,60)*AMRO2**2,0.0)
131 BET3 = CMPLX(PKORB(3,61)*AMRO3**2,0.0)
132C
133 END IF
134C**************************************************
135C
136C --- initialization of four vectors
137 DO 7 K=1,4
138 DO 8 L=1,4
139 8 AA(K,L)=0.0
140 HADCUR(K)=CMPLX(0.0)
141 PAA(K)=PIM1(K)+PIM2(K)+PIM3(K)+PIM4(K)
142 PP(1,K)=PIM1(K)
143 PP(2,K)=PIM2(K)
144 PP(3,K)=PIM3(K)
145 7 PP(4,K)=PIM4(K)
146C
147.EQ. IF (MNUM1) THEN
148C ===================================================================
149C pi- pi- p0 pi+ case ====
150C ===================================================================
151 QQ=PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2
152
153C Add M(4pi)-dependence to rhopipi channels:
154 FORM4= LAM0+LAM1*BWIGN(QQ,AMRO,GAMRO)
155 * +LAM2*BWIGN(QQ,AMRO2,GAMRO2)
156 * +LAM3*BWIGN(QQ,AMRO3,GAMRO3)
157
158C --- loop over five contributions of the rho-pi-pi
159 DO 201 K1=1,3
160 DO 201 K2=3,4
161C
162.EQ. IF (K2K1) THEN
163 GOTO 201
164.EQ. ELSEIF (K23) THEN
165C rho-
166 AMPR = AMPL(3)
167 AMPA = AMPIZ
168.EQ. ELSEIF (K13) THEN
169C rho+
170 AMPR = AMPL(4)
171 AMPA = AMPIZ
172 ELSE
173C rho0
174 AMPR = AMPL(2)
175 AMPA = AMPI
176 END IF
177C
178 SK=(PP(K1,4)+PP(K2,4))**2-(PP(K1,3)+PP(K2,3))**2
179 $ -(PP(K1,2)+PP(K2,2))**2-(PP(K1,1)+PP(K2,1))**2
180
181C -- definition of AA matrix
182C -- cronecker delta
183 DO 202 I=1,4
184 DO 203 J=1,4
185 203 AA(I,J)=0.0
186 202 AA(I,I)=1.0
187C ... and the rest ...
188 DO 204 L=1,4
189.NE..AND..NE. IF (LK1LK2) THEN
190 DENOM=(PAA(4)-PP(L,4))**2-(PAA(3)-PP(L,3))**2
191 $ -(PAA(2)-PP(L,2))**2-(PAA(1)-PP(L,1))**2
192 DO 205 I=1,4
193 DO 205 J=1,4
194 SIG= 1.0
195.NE. IF(J4) SIG=-SIG
196 AA(I,J)=AA(I,J)
197 $ -SIG*(PAA(I)-2.0*PP(L,I))*(PAA(J)-PP(L,J))/DENOM
198 205 CONTINUE
199 ENDIF
200 204 CONTINUE
201C
202C --- lets add something to HADCURR
203C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
204C FORM1= AMPL(1)+AMPR*FPIKM(SQRT(SK),AMPI,AMPI)
205
206 FORM2PI= BET1*BWIGM(SK,AMRO,GAMRO,AMPA,AMPI)
207 1 +BET2*BWIGM(SK,AMRO2,GAMRO2,AMPA,AMPI)
208 2 +BET3*BWIGM(SK,AMRO3,GAMRO3,AMPA,AMPI)
209 FORM1= AMPL(1)+AMPR*FORM2PI
210C
211 DO 206 I=1,4
212 DO 206 J=1,4
213 HADCUR(I)=HADCUR(I)+FORM1*FORM4*AA(I,J)*(PP(K1,J)-PP(K2,J))
214 206 CONTINUE
215C --- end of the rho-pi-pi current (5 possibilities)
216 201 CONTINUE
217C
218C ===================================================================
219C Now modify the coefficient for the omega-pi current: =
220C ===================================================================
221.EQ. IF (AMPL(5)CMPLX(0.,0.)) GOTO 311
222
223C Overall rho+rhoprime for the 4pi system:
224C FORM2=AMPL(5)*(BWIGN(QQ,AMRO,GAMRO)+ELPHA*BWIGN(QQ,AMROP,GAMROP))
225C Modified M(4pi)-dependence:
226 FORM2=AMPL(5)*(ALF0+ALF1*BWIGN(QQ,AMRO,GAMRO)
227 * +ALF2*BWIGN(QQ,AMRO2,GAMRO2)
228 * +ALF3*BWIGN(QQ,AMRO3,GAMRO3))
229C
230C --- there are two possibilities for omega current
231C --- PA PB are corresponding first and second pi-s
232 DO 301 KK=1,2
233 DO 302 I=1,4
234 PA(I)=PP(KK,I)
235 PB(I)=PP(3-KK,I)
236 302 CONTINUE
237C --- lorentz invariants
238 QQA=0.0
239 SS23=0.0
240 SS24=0.0
241 SS34=0.0
242 QP1P2=0.0
243 QP1P3=0.0
244 QP1P4=0.0
245 P1P2 =0.0
246 P1P3 =0.0
247 P1P4 =0.0
248 DO 303 K=1,4
249 SIGN=-1.0
250.EQ. IF (K4) SIGN= 1.0
251 QQA=QQA+SIGN*(PAA(K)-PA(K))**2
252 SS23=SS23+SIGN*(PB(K) +PIM3(K))**2
253 SS24=SS24+SIGN*(PB(K) +PIM4(K))**2
254 SS34=SS34+SIGN*(PIM3(K)+PIM4(K))**2
255 QP1P2=QP1P2+SIGN*(PAA(K)-PA(K))*PB(K)
256 QP1P3=QP1P3+SIGN*(PAA(K)-PA(K))*PIM3(K)
257 QP1P4=QP1P4+SIGN*(PAA(K)-PA(K))*PIM4(K)
258 P1P2=P1P2+SIGN*PA(K)*PB(K)
259 P1P3=P1P3+SIGN*PA(K)*PIM3(K)
260 P1P4=P1P4+SIGN*PA(K)*PIM4(K)
261 303 CONTINUE
262C
263C omega -> rho pi for the 3pi system:
264C FORM3=BWIGN(QQA,AMOM,GAMOM)*(BWIGN(SS23,AMRO,GAMRO)+
265C $ BWIGN(SS24,AMRO,GAMRO)+BWIGN(SS34,AMRO,GAMRO))
266C No omega -> rho pi; just straight omega:
267 FORM3=BWIGN(QQA,AMOM,GAMOM)
268C
269 DO 304 K=1,4
270 HADCUR(K)=HADCUR(K)+FORM2*FORM3*(
271 $ PB (K)*(QP1P3*P1P4-QP1P4*P1P3)
272 $ +PIM3(K)*(QP1P4*P1P2-QP1P2*P1P4)
273 $ +PIM4(K)*(QP1P2*P1P3-QP1P3*P1P2) )
274 304 CONTINUE
275 301 CONTINUE
276 311 CONTINUE
277C
278 ELSE
279C ===================================================================
280C pi0 pi0 p0 pi- case ====
281C ===================================================================
282 QQ=PAA(4)**2-PAA(3)**2-PAA(2)**2-PAA(1)**2
283
284C --- loop over three contribution of the non-omega current
285 DO 101 K=1,3
286 SK=(PP(K,4)+PIM4(4))**2-(PP(K,3)+PIM4(3))**2
287 $ -(PP(K,2)+PIM4(2))**2-(PP(K,1)+PIM4(1))**2
288
289C -- definition of AA matrix
290C -- cronecker delta
291 DO 102 I=1,4
292 DO 103 J=1,4
293 103 AA(I,J)=0.0
294 102 AA(I,I)=1.0
295C
296C ... and the rest ...
297 DO 104 L=1,3
298.NE. IF (LK) THEN
299 DENOM=(PAA(4)-PP(L,4))**2-(PAA(3)-PP(L,3))**2
300 $ -(PAA(2)-PP(L,2))**2-(PAA(1)-PP(L,1))**2
301 DO 105 I=1,4
302 DO 105 J=1,4
303 SIG=1.0
304.NE. IF(J4) SIG=-SIG
305 AA(I,J)=AA(I,J)
306 $ -SIG*(PAA(I)-2.0*PP(L,I))*(PAA(J)-PP(L,J))/DENOM
307 105 CONTINUE
308 ENDIF
309 104 CONTINUE
310
311C --- lets add something to HADCURR
312C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKMD(SQRT(QQ),AMPI,AMPI)
313CCCCCCCCCCCCC FORM1=WIGFOR(SK,AMRO,GAMRO) (tests)
314C FORM1= FPIKM(SQRT(SK),AMPI,AMPI) *FPIKM(SQRT(QQ),AMPI,AMPI)
315 FORM1 = AMPL(6)+AMPL(7)*FPIKM(SQRT(SK),AMPI,AMPI)
316
317 DO 106 I=1,4
318 DO 106 J=1,4
319 HADCUR(I)=HADCUR(I)+FORM1*AA(I,J)*(PP(K,J)-PP(4,J))
320 106 CONTINUE
321C --- end of the non omega current (3 possibilities)
322 101 CONTINUE
323
324 ENDIF
325 END
326
327
328