C++InterfacetoTauola
f3pi.f
1 /* copyright(c) 1991-2021 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)
45 C.......................................................................
46 C.
47 C. F3PI - 1 version of a1 form factor, used in TAUOLA
48 C.
49 C. Inputs : None
50 C. :
51 C. Outputs : None
52 C.
53 C. COMMON : None
54 C.
55 C. Calls :
56 C. Called : by FORM1-FORM3 in $C_CVSSRC/korb/koralb/formf.F
57 C. Author : Alan Weinstein 2/98
58 C.
59 C. Detailed description
60 C. First determine whether we are doing pi-2pi0 or 3pi.
61 C. Then implement full form-factor from fit:
62 C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
63 C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
64 C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
65 C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
66 C. All the parameters in this routine are hard-coded!!
67 C.
68 C.......................................................................
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 *
102 C. Hard-code the fit parameters:
103 .EQ. IF (IFIRST0) THEN
104  IFIRST = 1
105 C 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 
119 C 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 
134 C Initialize to 0:
135  F3PI = CMPLX(0.,0.)
136 
137 C. First determine whether we are doing pi-2pi0 or 3pi.
138 C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
139 C since KORALB doesnt bother to remember!!
140  R = PKORB(4,11)
141 .EQ. IF (R0.) THEN
142 C it is 2pi0pi-
143  IDK = 1
144  M1 = MPIZ
145  M2 = MPIZ
146  M3 = MPIC
147  ELSE
148 C 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 
158 C. Then implement full form-factor from fit:
159 C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
160 C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
161 C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
162 C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
163 
164 C Note that for FORM1, the arguments are S1, S2;
165 C for FORM2, the arguments are S2, S1;
166 C for FORM3, the arguments are S3, S1.
167 C Here, we implement FORM1 and FORM2 at the same time,
168 C 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
177 C it is 2pi0pi-
178 C 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 
183 C 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 
197 C F3PI = FPIKM(SQRT(S1),M2,M3)
198 .EQ. ELSEIF (IDK2) THEN
199 C it is 3pi
200 C 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 
206 C 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 
221 C 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
231 C it is 2pi0pi-
232 C 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 
237 C 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 
249 C F3PI = CMPLX(0.,0.)
250 .EQ. ELSEIF (IDK2) THEN
251 C it is 3pi
252 C 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 
260 C 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 
279 C F3PI = CMPLX(0.,0.)
280  END IF
281  END IF
282 
283 C Add overall a1/a1prime:
284  FORMA1 = FA1A1P(QQ)
285  F3PI = F3PI*FORMA1
286 
287  RETURN
288  END
289 C **********************************************************
290  COMPLEX FUNCTION BWIGML(S,M,G,M1,M2,L)
291 C **********************************************************
292 C L-WAVE BREIT-WIGNER
293 C **********************************************************
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
314 C=======================================================================
315  COMPLEX FUNCTION FA1A1P(XMSQ)
316 C ==================================================================
317 C complex form-factor for a1+a1prime. AJW 1/98
318 C ==================================================================
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 
329 C 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.)
335 C 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
352 C=======================================================================
353  FUNCTION WGA1(QQ)
354 
355 C mass-dependent M*Gamma of a1 through its decays to
356 C. [(rho-pi S-wave) + (rho-pi D-wave) +
357 C. (f2 pi D-wave) + (f0pi S-wave)]
358 C. 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
364 C-----------------------------------------------------------------------
365 C
366 .NE. IF (IFIRST987) THEN
367  IFIRST = 987
368 C
369 C 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
374 C coupling constants squared:
375  C3PI = 0.2384D0**2
376  CKST = 4.7621D0**2*C3PI
377  END IF
378 
379 C-----------------------------------------------------------------------
380 C Parameterization of numerical integral of total width of a1 to 3pi.
381 C From M. Schmidtler, CBX-97-64-Update.
382  S = DBLE(QQ)
383  WG3PIC = WGA1C(S)
384  WG3PIN = WGA1N(S)
385 
386 C 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
394 C=======================================================================
395  DOUBLE PRECISION FUNCTION WGA1C(S)
396 C
397 C parameterization of m*Gamma(m^2) for pi-2pi0 system
398 C
399  DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
400 C
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)
404 C
405  PARAMETER (STH = 0.1753D0)
406 C---------------------------------------------------------------------
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
419 C=======================================================================
420  DOUBLE PRECISION FUNCTION WGA1N(S)
421 C
422 C parameterization of m*Gamma(m^2) for pi-pi+pi- system
423 C
424  DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
425 C
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)
429 C
430  PARAMETER (STH = 0.1676D0)
431 C---------------------------------------------------------------------
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
Generated on Tue Mar 12 2024 22:57:50 for C++InterfacetoTauola by   doxygen 1.8.7