/* 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