/* Copyright (C) 1991-2012 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
. */
/* This header is separate from features.h so that the compiler can
include it implicitly at the start of every compilation. It must
not itself include or any other header that includes
because the implicit include comes before any feature
test macros that may be defined in a source file before it first
explicitly includes a system header. GCC knows the name of this
header in order to preinclude it. */
/* We do support the IEC 559 math functionality, real and complex. */
/* wchar_t uses ISO/IEC 10646 (2nd ed., published 2011-03-15) /
Unicode 6.0. */
/* We do not support C11 . */
*AJW 1 version of a1 form factor
COMPLEX FUNCTION F3PI(IFORM,QQ,SA,SB)
C.......................................................................
C.
C. F3PI - 1 version of a1 form factor, used in TAUOLA
C.
C. Inputs : None
C. :
C. Outputs : None
C.
C. COMMON : None
C.
C. Calls :
C. Called : by FORM1-FORM3 in $C_CVSSRC/korb/koralb/formf.F
C. Author : Alan Weinstein 2/98
C.
C. Detailed description
C. First determine whether we are doing pi-2pi0 or 3pi.
C. Then implement full form-factor from fit:
C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
C. All the parameters in this routine are hard-coded!!
C.
C.......................................................................
* -------------------- Argument declarations ---------------
INTEGER IFORM
REAL QQ,SA,SB
* -------------------- EXTERNAL declarations ---------------
*
REAL PKORB
COMPLEX BWIGML
* -------------------- SEQUENCE declarations ---------------
*
* -------------------- Local declarations ---------------
*
CHARACTER*(*) CRNAME
PARAMETER( CRNAME = 'F3PI' )
*
INTEGER IFIRST,IDK
REAL MRO,GRO,MRP,GRP,MF2,GF2,MF0,GF0,MSG,GSG
REAL M1,M2,M3,M1SQ,M2SQ,M3SQ,MPIZ,MPIC
REAL S1,S2,S3,R,PI
REAL F134,F150,F15A,F15B,F167
REAL F34A,F34B,F35,F35A,F35B,F36A,F36B
COMPLEX BT1,BT2,BT3,BT4,BT5,BT6,BT7
COMPLEX FRO1,FRO2,FRP1,FRP2
COMPLEX FF21,FF22,FF23,FSG1,FSG2,FSG3,FF01,FF02,FF03
COMPLEX FA1A1P,FORMA1
* -------------------- SAVE declarations ---------------
*
* -------------------- DATA initializations ---------------
*
DATA IFIRST/0/
* ----------------- Executable code starts here ------------
*
C. Hard-code the fit parameters:
IF (IFIRST.EQ.0) THEN
IFIRST = 1
C rho, rhoprime, f2(1275), f0(1186), sigma(made up!)
MRO = 0.7743
GRO = 0.1491
MRP = 1.370
GRP = 0.386
MF2 = 1.275
GF2 = 0.185
MF0 = 1.186
GF0 = 0.350
MSG = 0.860
GSG = 0.880
MPIZ = PKORB(1,7)
MPIC = PKORB(1,8)
C Fit coefficients for each of the contributions:
PI = 3.14159
BT1 = CMPLX(1.,0.)
BT2 = CMPLX(0.12,0.)*CEXP(CMPLX(0., 0.99*PI))
BT3 = CMPLX(0.37,0.)*CEXP(CMPLX(0.,-0.15*PI))
BT4 = CMPLX(0.87,0.)*CEXP(CMPLX(0., 0.53*PI))
BT5 = CMPLX(0.71,0.)*CEXP(CMPLX(0., 0.56*PI))
BT6 = CMPLX(2.10,0.)*CEXP(CMPLX(0., 0.23*PI))
BT7 = CMPLX(0.77,0.)*CEXP(CMPLX(0.,-0.54*PI))
PRINT *,' In F3pi: add (rho-pi S-wave) + (rhop-pi S-wave) +'
PRINT *,' (rho-pi D-wave) + (rhop-pi D-wave) +'
PRINT *,' (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)'
END IF
C Initialize to 0:
F3PI = CMPLX(0.,0.)
C. First determine whether we are doing pi-2pi0 or 3pi.
C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
C since KORALB doesnt bother to remember!!
R = PKORB(4,11)
IF (R.EQ.0.) THEN
C it is 2pi0pi-
IDK = 1
M1 = MPIZ
M2 = MPIZ
M3 = MPIC
ELSE
C it is 3pi
IDK = 2
M1 = MPIC
M2 = MPIC
M3 = MPIC
END IF
M1SQ = M1*M1
M2SQ = M2*M2
M3SQ = M3*M3
C. Then implement full form-factor from fit:
C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
C Note that for FORM1, the arguments are S1, S2;
C for FORM2, the arguments are S2, S1;
C for FORM3, the arguments are S3, S1.
C Here, we implement FORM1 and FORM2 at the same time,
C so the above switch is just what we need!
IF (IFORM.EQ.1.OR.IFORM.EQ.2) THEN
S1 = SA
S2 = SB
S3 = QQ-SA-SB+M1SQ+M2SQ+M3SQ
IF (S3.LE.0..OR.S2.LE.0.) RETURN
IF (IDK.EQ.1) THEN
C it is 2pi0pi-
C Lorentz invariants for all the contributions:
F134 = -(1./3.)*((S3-M3SQ)-(S1-M1SQ))
F150 = (1./18.)*(QQ-M3SQ+S3)*(2.*M1SQ+2.*M2SQ-S3)/S3
F167 = (2./3.)
C Breit Wigners for all the contributions:
FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
FF23 = BWIGML(S3,MF2,GF2,M1,M2,2)
FSG3 = BWIGML(S3,MSG,GSG,M1,M2,0)
FF03 = BWIGML(S3,MF0,GF0,M1,M2,0)
F3PI = BT1*FRO1+BT2*FRP1+
1 BT3*CMPLX(F134,0.)*FRO2+BT4*CMPLX(F134,0.)*FRP2+
1 BT5*CMPLX(F150,0.)*FF23+
1 BT6*CMPLX(F167,0.)*FSG3+BT7*CMPLX(F167,0.)*FF03
C F3PI = FPIKM(SQRT(S1),M2,M3)
ELSEIF (IDK.EQ.2) THEN
C it is 3pi
C Lorentz invariants for all the contributions:
F134 = -(1./3.)*((S3-M3SQ)-(S1-M1SQ))
F15A = -(1./2.)*((S2-M2SQ)-(S3-M3SQ))
F15B = -(1./18.)*(QQ-M2SQ+S2)*(2.*M1SQ+2.*M3SQ-S2)/S2
F167 = -(2./3.)
C Breit Wigners for all the contributions:
FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
FF21 = BWIGML(S1,MF2,GF2,M2,M3,2)
FF22 = BWIGML(S2,MF2,GF2,M3,M1,2)
FSG2 = BWIGML(S2,MSG,GSG,M3,M1,0)
FF02 = BWIGML(S2,MF0,GF0,M3,M1,0)
F3PI = BT1*FRO1+BT2*FRP1+
1 BT3*CMPLX(F134,0.)*FRO2+BT4*CMPLX(F134,0.)*FRP2
1 -BT5*CMPLX(F15A,0.)*FF21-BT5*CMPLX(F15B,0.)*FF22
1 -BT6*CMPLX(F167,0.)*FSG2-BT7*CMPLX(F167,0.)*FF02
C F3PI = FPIKM(SQRT(S1),M2,M3)
END IF
ELSE IF (IFORM.EQ.3) THEN
S3 = SA
S1 = SB
S2 = QQ-SA-SB+M1SQ+M2SQ+M3SQ
IF (S1.LE.0..OR.S2.LE.0.) RETURN
IF (IDK.EQ.1) THEN
C it is 2pi0pi-
C Lorentz invariants for all the contributions:
F34A = (1./3.)*((S2-M2SQ)-(S3-M3SQ))
F34B = (1./3.)*((S3-M3SQ)-(S1-M1SQ))
F35 =-(1./2.)*((S1-M1SQ)-(S2-M2SQ))
C Breit Wigners for all the contributions:
FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
FF23 = BWIGML(S3,MF2,GF2,M1,M2,2)
F3PI =
1 BT3*(CMPLX(F34A,0.)*FRO1+CMPLX(F34B,0.)*FRO2)+
1 BT4*(CMPLX(F34A,0.)*FRP1+CMPLX(F34B,0.)*FRP2)+
1 BT5*CMPLX(F35,0.)*FF23
C F3PI = CMPLX(0.,0.)
ELSEIF (IDK.EQ.2) THEN
C it is 3pi
C Lorentz invariants for all the contributions:
F34A = (1./3.)*((S2-M2SQ)-(S3-M3SQ))
F34B = (1./3.)*((S3-M3SQ)-(S1-M1SQ))
F35A = -(1./18.)*(QQ-M1SQ+S1)*(2.*M2SQ+2.*M3SQ-S1)/S1
F35B = (1./18.)*(QQ-M2SQ+S2)*(2.*M3SQ+2.*M1SQ-S2)/S2
F36A = -(2./3.)
F36B = (2./3.)
C Breit Wigners for all the contributions:
FRO1 = BWIGML(S1,MRO,GRO,M2,M3,1)
FRP1 = BWIGML(S1,MRP,GRP,M2,M3,1)
FRO2 = BWIGML(S2,MRO,GRO,M3,M1,1)
FRP2 = BWIGML(S2,MRP,GRP,M3,M1,1)
FF21 = BWIGML(S1,MF2,GF2,M2,M3,2)
FF22 = BWIGML(S2,MF2,GF2,M3,M1,2)
FSG1 = BWIGML(S1,MSG,GSG,M2,M3,0)
FSG2 = BWIGML(S2,MSG,GSG,M3,M1,0)
FF01 = BWIGML(S1,MF0,GF0,M2,M3,0)
FF02 = BWIGML(S2,MF0,GF0,M3,M1,0)
F3PI =
1 BT3*(CMPLX(F34A,0.)*FRO1+CMPLX(F34B,0.)*FRO2)+
1 BT4*(CMPLX(F34A,0.)*FRP1+CMPLX(F34B,0.)*FRP2)
1 -BT5*(CMPLX(F35A,0.)*FF21+CMPLX(F35B,0.)*FF22)
1 -BT6*(CMPLX(F36A,0.)*FSG1+CMPLX(F36B,0.)*FSG2)
1 -BT7*(CMPLX(F36A,0.)*FF01+CMPLX(F36B,0.)*FF02)
C F3PI = CMPLX(0.,0.)
END IF
END IF
C Add overall a1/a1prime:
FORMA1 = FA1A1P(QQ)
F3PI = F3PI*FORMA1
RETURN
END
C **********************************************************
COMPLEX FUNCTION BWIGML(S,M,G,M1,M2,L)
C **********************************************************
C L-WAVE BREIT-WIGNER
C **********************************************************
REAL S,M,G,M1,M2
INTEGER L,IPOW
REAL MSQ,W,WGS,MP,MM,QS,QM
MP = (M1+M2)**2
MM = (M1-M2)**2
MSQ = M*M
W = SQRT(S)
WGS = 0.0
IF (W.GT.(M1+M2)) THEN
QS=SQRT(ABS((S -MP)*(S -MM)))/W
QM=SQRT(ABS((MSQ -MP)*(MSQ -MM)))/M
IPOW = 2*L+1
WGS=G*(MSQ/W)*(QS/QM)**IPOW
ENDIF
BWIGML=CMPLX(MSQ,0.)/CMPLX(MSQ-S,-WGS)
RETURN
END
C=======================================================================
COMPLEX FUNCTION FA1A1P(XMSQ)
C ==================================================================
C complex form-factor for a1+a1prime. AJW 1/98
C ==================================================================
REAL XMSQ
REAL PKORB,WGA1
REAL XM1,XG1,XM2,XG2,XM1SQ,XM2SQ,GG1,GG2,GF,FG1,FG2
COMPLEX BET,F1,F2
INTEGER IFIRST/0/
IF (IFIRST.EQ.0) THEN
IFIRST = 1
C The user may choose masses and widths that differ from nominal:
XM1 = PKORB(1,10)
XG1 = PKORB(2,10)
XM2 = PKORB(1,17)
XG2 = PKORB(2,17)
BET = CMPLX(PKORB(3,17),0.)
C scale factors relative to nominal:
GG1 = XM1*XG1/(1.3281*0.806)
GG2 = XM2*XG2/(1.3281*0.806)
XM1SQ = XM1*XM1
XM2SQ = XM2*XM2
END IF
GF = WGA1(XMSQ)
FG1 = GG1*GF
FG2 = GG2*GF
F1 = CMPLX(-XM1SQ,0.0)/CMPLX(XMSQ-XM1SQ,FG1)
F2 = CMPLX(-XM2SQ,0.0)/CMPLX(XMSQ-XM2SQ,FG2)
FA1A1P = F1+BET*F2
RETURN
END
C=======================================================================
FUNCTION WGA1(QQ)
C mass-dependent M*Gamma of a1 through its decays to
C. [(rho-pi S-wave) + (rho-pi D-wave) +
C. (f2 pi D-wave) + (f0pi S-wave)]
C. AND simple K*K S-wave
REAL QQ,WGA1
DOUBLE PRECISION MKST,MK,MK1SQ,MK2SQ,C3PI,CKST
DOUBLE PRECISION S,WGA1C,WGA1N,WG3PIC,WG3PIN,GKST
INTEGER IFIRST
C-----------------------------------------------------------------------
C
IF (IFIRST.NE.987) THEN
IFIRST = 987
C
C Contribution to M*Gamma(m(3pi)^2) from S-wave K*K:
MKST = 0.894D0
MK = 0.496D0
MK1SQ = (MKST+MK)**2
MK2SQ = (MKST-MK)**2
C coupling constants squared:
C3PI = 0.2384D0**2
CKST = 4.7621D0**2*C3PI
END IF
C-----------------------------------------------------------------------
C Parameterization of numerical integral of total width of a1 to 3pi.
C From M. Schmidtler, CBX-97-64-Update.
S = DBLE(QQ)
WG3PIC = WGA1C(S)
WG3PIN = WGA1N(S)
C Contribution to M*Gamma(m(3pi)^2) from S-wave K*K, if above threshold
GKST = 0.D0
IF (S.GT.MK1SQ) GKST = SQRT((S-MK1SQ)*(S-MK2SQ))/(2.*S)
WGA1 = SNGL(C3PI*(WG3PIC+WG3PIN)+CKST*GKST)
RETURN
END
C=======================================================================
DOUBLE PRECISION FUNCTION WGA1C(S)
C
C parameterization of m*Gamma(m^2) for pi-2pi0 system
C
DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
C
PARAMETER(Q0 = 5.80900D0,Q1 = -3.00980D0,Q2 = 4.57920D0,
1 P0 = -13.91400D0,P1 = 27.67900D0,P2 = -13.39300D0,
2 P3 = 3.19240D0,P4 = -0.10487D0)
C
PARAMETER (STH = 0.1753D0)
C---------------------------------------------------------------------
IF(S.LT.STH) THEN
G1_IM = 0.D0
ELSEIF((S.GT.STH).AND.(S.LT.0.823D0)) THEN
G1_IM = Q0*(S-STH)**3*(1. + Q1*(S-STH) + Q2*(S-STH)**2)
ELSE
G1_IM = P0 + P1*S + P2*S**2+ P3*S**3 + P4*S**4
ENDIF
WGA1C = G1_IM
RETURN
END
C=======================================================================
DOUBLE PRECISION FUNCTION WGA1N(S)
C
C parameterization of m*Gamma(m^2) for pi-pi+pi- system
C
DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
C
PARAMETER(Q0 = 6.28450D0,Q1 = -2.95950D0,Q2 = 4.33550D0,
1 P0 = -15.41100D0,P1 = 32.08800D0,P2 = -17.66600D0,
2 P3 = 4.93550D0,P4 = -0.37498D0)
C
PARAMETER (STH = 0.1676D0)
C---------------------------------------------------------------------
IF(S.LT.STH) THEN
G1_IM = 0.D0
ELSEIF((S.GT.STH).AND.(S.LT.0.823D0)) THEN
G1_IM = Q0*(S-STH)**3*(1. + Q1*(S-STH) + Q2*(S-STH)**2)
ELSE
G1_IM = P0 + P1*S + P2*S**2+ P3*S**3 + P4*S**4
ENDIF
WGA1N = G1_IM
RETURN
END