C++ Interface to Tauola
tauola-F/prod/f3pi.f
1/* copyright(c) 1991-2024 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 a1 form factor
44 COMPLEX FUNCTION F3PI(IFORM,QQ,SA,SB)
45C.......................................................................
46C.
47C. F3PI - 1 version of a1 form factor, used in TAUOLA
48C.
49C. Inputs : None
50C. :
51C. Outputs : None
52C.
53C. COMMON : None
54C.
55C. Calls :
56C. Called : by FORM1-FORM3 in $C_CVSSRC/korb/koralb/formf.F
57C. Author : Alan Weinstein 2/98
58C.
59C. Detailed description
60C. First determine whether we are doing pi-2pi0 or 3pi.
61C. Then implement full form-factor from fit:
62C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
63C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
64C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
65C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
66C. All the parameters in this routine are hard-coded!!
67C.
68C.......................................................................
69* -------------------- Argument declarations ---------------
70
71 INTEGER IFORM
72 REAL QQ,SA,SB
73* -------------------- EXTERNAL declarations ---------------
74*
75 REAL PKORB
76 COMPLEX BWIGML
77* -------------------- SEQUENCE declarations ---------------
78*
79* -------------------- Local declarations ---------------
80*
81 CHARACTER*(*) CRNAME
82 PARAMETER( CRNAME = 'f3pi' )
83*
84 INTEGER IFIRST,IDK
85 REAL MRO,GRO,MRP,GRP,MF2,GF2,MF0,GF0,MSG,GSG
86 REAL M1,M2,M3,M1SQ,M2SQ,M3SQ,MPIZ,MPIC
87 REAL S1,S2,S3,R,PI
88 REAL F134,F150,F15A,F15B,F167
89 REAL F34A,F34B,F35,F35A,F35B,F36A,F36B
90 COMPLEX BT1,BT2,BT3,BT4,BT5,BT6,BT7
91 COMPLEX FRO1,FRO2,FRP1,FRP2
92 COMPLEX FF21,FF22,FF23,FSG1,FSG2,FSG3,FF01,FF02,FF03
93 COMPLEX FA1A1P,FORMA1
94
95* -------------------- SAVE declarations ---------------
96*
97* -------------------- DATA initializations ---------------
98*
99 DATA IFIRST/0/
100* ----------------- Executable code starts here ------------
101*
102C. Hard-code the fit parameters:
103.EQ. IF (IFIRST0) THEN
104 IFIRST = 1
105C rho, rhoprime, f2(1275), f0(1186), sigma(made up!)
106 MRO = 0.7743
107 GRO = 0.1491
108 MRP = 1.370
109 GRP = 0.386
110 MF2 = 1.275
111 GF2 = 0.185
112 MF0 = 1.186
113 GF0 = 0.350
114 MSG = 0.860
115 GSG = 0.880
116 MPIZ = PKORB(1,7)
117 MPIC = PKORB(1,8)
118
119C Fit coefficients for each of the contributions:
120 PI = 3.14159
121 BT1 = CMPLX(1.,0.)
122 BT2 = CMPLX(0.12,0.)*CEXP(CMPLX(0., 0.99*PI))
123 BT3 = CMPLX(0.37,0.)*CEXP(CMPLX(0.,-0.15*PI))
124 BT4 = CMPLX(0.87,0.)*CEXP(CMPLX(0., 0.53*PI))
125 BT5 = CMPLX(0.71,0.)*CEXP(CMPLX(0., 0.56*PI))
126 BT6 = CMPLX(2.10,0.)*CEXP(CMPLX(0., 0.23*PI))
127 BT7 = CMPLX(0.77,0.)*CEXP(CMPLX(0.,-0.54*PI))
128
129 PRINT *,' in f3pi: add(rho-pi s-wave) + (rhop-pi s-wave) +'
130 PRINT *,' (rho-pi d-wave) + (rhop-pi d-wave) +'
131 PRINT *,' (f2 pi d-wave) + (sigmapi s-wave) + (f0pi s-wave)'
132 END IF
133
134C Initialize to 0:
135 F3PI = CMPLX(0.,0.)
136
137C. First determine whether we are doing pi-2pi0 or 3pi.
138C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
139C since KORALB doesnt bother to remember!!
140 R = PKORB(4,11)
141.EQ. IF (R0.) THEN
142C it is 2pi0pi-
143 IDK = 1
144 M1 = MPIZ
145 M2 = MPIZ
146 M3 = MPIC
147 ELSE
148C it is 3pi
149 IDK = 2
150 M1 = MPIC
151 M2 = MPIC
152 M3 = MPIC
153 END IF
154 M1SQ = M1*M1
155 M2SQ = M2*M2
156 M3SQ = M3*M3
157
158C. Then implement full form-factor from fit:
159C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
160C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
161C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
162C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
163
164C Note that for FORM1, the arguments are S1, S2;
165C for FORM2, the arguments are S2, S1;
166C for FORM3, the arguments are S3, S1.
167C Here, we implement FORM1 and FORM2 at the same time,
168C so the above switch is just what we need!
169
170.EQ..OR..EQ. IF (IFORM1IFORM2) THEN
171 S1 = SA
172 S2 = SB
173 S3 = QQ-SA-SB+M1SQ+M2SQ+M3SQ
174.LE..OR..LE. IF (S30.S20.) RETURN
175
176.EQ. IF (IDK1) THEN
177C it is 2pi0pi-
178C Lorentz invariants for all the contributions:
179 F134 = -(1./3.)*((S3-M3SQ)-(S1-M1SQ))
180 F150 = (1./18.)*(QQ-M3SQ+S3)*(2.*M1SQ+2.*M2SQ-S3)/S3
181 F167 = (2./3.)
182
183C Breit Wigners for all the contributions:
184 FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
185 FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
186 FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
187 FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
188 FF23 = BWIGML(S3,MF2,GF2,M1,M2,2)
189 FSG3 = BWIGML(S3,MSG,GSG,M1,M2,0)
190 FF03 = BWIGML(S3,MF0,GF0,M1,M2,0)
191
192 F3PI = BT1*FRO1+BT2*FRP1+
193 1 BT3*CMPLX(F134,0.)*FRO2+BT4*CMPLX(F134,0.)*FRP2+
194 1 BT5*CMPLX(F150,0.)*FF23+
195 1 BT6*CMPLX(F167,0.)*FSG3+BT7*CMPLX(F167,0.)*FF03
196
197C F3PI = FPIKM(SQRT(S1),M2,M3)
198.EQ. ELSEIF (IDK2) THEN
199C it is 3pi
200C Lorentz invariants for all the contributions:
201 F134 = -(1./3.)*((S3-M3SQ)-(S1-M1SQ))
202 F15A = -(1./2.)*((S2-M2SQ)-(S3-M3SQ))
203 F15B = -(1./18.)*(QQ-M2SQ+S2)*(2.*M1SQ+2.*M3SQ-S2)/S2
204 F167 = -(2./3.)
205
206C Breit Wigners for all the contributions:
207 FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
208 FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
209 FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
210 FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
211 FF21 = BWIGML(S1,MF2,GF2,M2,M3,2)
212 FF22 = BWIGML(S2,MF2,GF2,M3,M1,2)
213 FSG2 = BWIGML(S2,MSG,GSG,M3,M1,0)
214 FF02 = BWIGML(S2,MF0,GF0,M3,M1,0)
215
216 F3PI = BT1*FRO1+BT2*FRP1+
217 1 BT3*CMPLX(F134,0.)*FRO2+BT4*CMPLX(F134,0.)*FRP2
218 1 -BT5*CMPLX(F15A,0.)*FF21-BT5*CMPLX(F15B,0.)*FF22
219 1 -BT6*CMPLX(F167,0.)*FSG2-BT7*CMPLX(F167,0.)*FF02
220
221C F3PI = FPIKM(SQRT(S1),M2,M3)
222 END IF
223
224.EQ. ELSE IF (IFORM3) THEN
225 S3 = SA
226 S1 = SB
227 S2 = QQ-SA-SB+M1SQ+M2SQ+M3SQ
228.LE..OR..LE. IF (S10.S20.) RETURN
229
230.EQ. IF (IDK1) THEN
231C it is 2pi0pi-
232C Lorentz invariants for all the contributions:
233 F34A = (1./3.)*((S2-M2SQ)-(S3-M3SQ))
234 F34B = (1./3.)*((S3-M3SQ)-(S1-M1SQ))
235 F35 =-(1./2.)*((S1-M1SQ)-(S2-M2SQ))
236
237C Breit Wigners for all the contributions:
238 FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
239 FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
240 FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
241 FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
242 FF23 = BWIGML(S3,MF2,GF2,M1,M2,2)
243
244 F3PI =
245 1 BT3*(CMPLX(F34A,0.)*FRO1+CMPLX(F34B,0.)*FRO2)+
246 1 BT4*(CMPLX(F34A,0.)*FRP1+CMPLX(F34B,0.)*FRP2)+
247 1 BT5*CMPLX(F35,0.)*FF23
248
249C F3PI = CMPLX(0.,0.)
250.EQ. ELSEIF (IDK2) THEN
251C it is 3pi
252C Lorentz invariants for all the contributions:
253 F34A = (1./3.)*((S2-M2SQ)-(S3-M3SQ))
254 F34B = (1./3.)*((S3-M3SQ)-(S1-M1SQ))
255 F35A = -(1./18.)*(QQ-M1SQ+S1)*(2.*M2SQ+2.*M3SQ-S1)/S1
256 F35B = (1./18.)*(QQ-M2SQ+S2)*(2.*M3SQ+2.*M1SQ-S2)/S2
257 F36A = -(2./3.)
258 F36B = (2./3.)
259
260C Breit Wigners for all the contributions:
261 FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
262 FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
263 FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
264 FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
265 FF21 = BWIGML(S1,MF2,GF2,M2,M3,2)
266 FF22 = BWIGML(S2,MF2,GF2,M3,M1,2)
267 FSG1 = BWIGML(S1,MSG,GSG,M2,M3,0)
268 FSG2 = BWIGML(S2,MSG,GSG,M3,M1,0)
269 FF01 = BWIGML(S1,MF0,GF0,M2,M3,0)
270 FF02 = BWIGML(S2,MF0,GF0,M3,M1,0)
271
272 F3PI =
273 1 BT3*(CMPLX(F34A,0.)*FRO1+CMPLX(F34B,0.)*FRO2)+
274 1 BT4*(CMPLX(F34A,0.)*FRP1+CMPLX(F34B,0.)*FRP2)
275 1 -BT5*(CMPLX(F35A,0.)*FF21+CMPLX(F35B,0.)*FF22)
276 1 -BT6*(CMPLX(F36A,0.)*FSG1+CMPLX(F36B,0.)*FSG2)
277 1 -BT7*(CMPLX(F36A,0.)*FF01+CMPLX(F36B,0.)*FF02)
278
279C F3PI = CMPLX(0.,0.)
280 END IF
281 END IF
282
283C Add overall a1/a1prime:
284 FORMA1 = FA1A1P(QQ)
285 F3PI = F3PI*FORMA1
286
287 RETURN
288 END
289C **********************************************************
290 COMPLEX FUNCTION BWIGML(S,M,G,M1,M2,L)
291C **********************************************************
292C L-WAVE BREIT-WIGNER
293C **********************************************************
294 REAL S,M,G,M1,M2
295 INTEGER L,IPOW
296 REAL MSQ,W,WGS,MP,MM,QS,QM
297
298 MP = (M1+M2)**2
299 MM = (M1-M2)**2
300 MSQ = M*M
301 W = SQRT(S)
302 WGS = 0.0
303.GT. IF (W(M1+M2)) THEN
304 QS=SQRT(ABS((S -MP)*(S -MM)))/W
305 QM=SQRT(ABS((MSQ -MP)*(MSQ -MM)))/M
306 IPOW = 2*L+1
307 WGS=G*(MSQ/W)*(QS/QM)**IPOW
308 ENDIF
309
310 BWIGML=CMPLX(MSQ,0.)/CMPLX(MSQ-S,-WGS)
311
312 RETURN
313 END
314C=======================================================================
315 COMPLEX FUNCTION FA1A1P(XMSQ)
316C ==================================================================
317C complex form-factor for a1+a1prime. AJW 1/98
318C ==================================================================
319
320 REAL XMSQ
321 REAL PKORB,WGA1
322 REAL XM1,XG1,XM2,XG2,XM1SQ,XM2SQ,GG1,GG2,GF,FG1,FG2
323 COMPLEX BET,F1,F2
324 INTEGER IFIRST/0/
325
326.EQ. IF (IFIRST0) THEN
327 IFIRST = 1
328
329C The user may choose masses and widths that differ from nominal:
330 XM1 = PKORB(1,10)
331 XG1 = PKORB(2,10)
332 XM2 = PKORB(1,17)
333 XG2 = PKORB(2,17)
334 BET = CMPLX(PKORB(3,17),0.)
335C scale factors relative to nominal:
336 GG1 = XM1*XG1/(1.3281*0.806)
337 GG2 = XM2*XG2/(1.3281*0.806)
338
339 XM1SQ = XM1*XM1
340 XM2SQ = XM2*XM2
341 END IF
342
343 GF = WGA1(XMSQ)
344 FG1 = GG1*GF
345 FG2 = GG2*GF
346 F1 = CMPLX(-XM1SQ,0.0)/CMPLX(XMSQ-XM1SQ,FG1)
347 F2 = CMPLX(-XM2SQ,0.0)/CMPLX(XMSQ-XM2SQ,FG2)
348 FA1A1P = F1+BET*F2
349
350 RETURN
351 END
352C=======================================================================
353 FUNCTION WGA1(QQ)
354
355C mass-dependent M*Gamma of a1 through its decays to
356C. [(rho-pi S-wave) + (rho-pi D-wave) +
357C. (f2 pi D-wave) + (f0pi S-wave)]
358C. AND simple K*K S-wave
359
360 REAL QQ,WGA1
361 DOUBLE PRECISION MKST,MK,MK1SQ,MK2SQ,C3PI,CKST
362 DOUBLE PRECISION S,WGA1C,WGA1N,WG3PIC,WG3PIN,GKST
363 INTEGER IFIRST
364C-----------------------------------------------------------------------
365C
366.NE. IF (IFIRST987) THEN
367 IFIRST = 987
368C
369C Contribution to M*Gamma(m(3pi)^2) from S-wave K*K:
370 MKST = 0.894D0
371 MK = 0.496D0
372 MK1SQ = (MKST+MK)**2
373 MK2SQ = (MKST-MK)**2
374C coupling constants squared:
375 C3PI = 0.2384D0**2
376 CKST = 4.7621D0**2*C3PI
377 END IF
378
379C-----------------------------------------------------------------------
380C Parameterization of numerical integral of total width of a1 to 3pi.
381C From M. Schmidtler, CBX-97-64-Update.
382 S = DBLE(QQ)
383 WG3PIC = WGA1C(S)
384 WG3PIN = WGA1N(S)
385
386C Contribution to M*Gamma(m(3pi)^2) from S-wave K*K, if above threshold
387 GKST = 0.D0
388.GT. IF (SMK1SQ) GKST = SQRT((S-MK1SQ)*(S-MK2SQ))/(2.*S)
389
390 WGA1 = SNGL(C3PI*(WG3PIC+WG3PIN)+CKST*GKST)
391
392 RETURN
393 END
394C=======================================================================
395 DOUBLE PRECISION FUNCTION WGA1C(S)
396C
397C parameterization of m*Gamma(m^2) for pi-2pi0 system
398C
399 DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
400C
401 PARAMETER(Q0 = 5.80900D0,Q1 = -3.00980D0,Q2 = 4.57920D0,
402 1 P0 = -13.91400D0,P1 = 27.67900D0,P2 = -13.39300D0,
403 2 P3 = 3.19240D0,P4 = -0.10487D0)
404C
405 PARAMETER (STH = 0.1753D0)
406C---------------------------------------------------------------------
407
408.LT. IF(SSTH) THEN
409 G1_IM = 0.D0
410.GT..AND..LT. ELSEIF((SSTH)(S0.823D0)) THEN
411 G1_IM = Q0*(S-STH)**3*(1. + Q1*(S-STH) + Q2*(S-STH)**2)
412 ELSE
413 G1_IM = P0 + P1*S + P2*S**2+ P3*S**3 + P4*S**4
414 ENDIF
415
416 WGA1C = G1_IM
417 RETURN
418 END
419C=======================================================================
420 DOUBLE PRECISION FUNCTION WGA1N(S)
421C
422C parameterization of m*Gamma(m^2) for pi-pi+pi- system
423C
424 DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
425C
426 PARAMETER(Q0 = 6.28450D0,Q1 = -2.95950D0,Q2 = 4.33550D0,
427 1 P0 = -15.41100D0,P1 = 32.08800D0,P2 = -17.66600D0,
428 2 P3 = 4.93550D0,P4 = -0.37498D0)
429C
430 PARAMETER (STH = 0.1676D0)
431C---------------------------------------------------------------------
432
433.LT. IF(SSTH) THEN
434 G1_IM = 0.D0
435.GT..AND..LT. ELSEIF((SSTH)(S0.823D0)) THEN
436 G1_IM = Q0*(S-STH)**3*(1. + Q1*(S-STH) + Q2*(S-STH)**2)
437 ELSE
438 G1_IM = P0 + P1*S + P2*S**2+ P3*S**3 + P4*S**4
439 ENDIF
440
441 WGA1N = G1_IM
442 RETURN
443 END