From 7cd6c77260b06bbf9f01f2c39521452a99656637 Mon Sep 17 00:00:00 2001 From: Alexander Austregesilo Date: Mon, 21 Oct 2013 11:02:59 +0200 Subject: [PATCH] added: AMP_M (from ROOTPWA) AMP_K1 (my own) AMP_Kred (from M. Pennington) simplified Flatte --- bw.cc | 383 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- bw.h | 3 + 2 files changed, 375 insertions(+), 11 deletions(-) diff --git a/bw.cc b/bw.cc index 99c6553..a4ad76b 100644 --- a/bw.cc +++ b/bw.cc @@ -75,24 +75,385 @@ breakupMomentum2(double s, double m1, double m2) return .25*(s - pow(m1+m2, 2))*(s - pow(m1-m2,2)) / s; } - +// complex version with analytic continuation below threshold complex -Flatte(double s, double m1, double m2, double m0, double g1) +breakupMomentumComplex(double s, double m1, double m2) { + double q2 = breakupMomentum2(s, m1, m2); + double q = sqrt(fabs(q2)); + if (q2 < 0) + return complex(0, q); + else + return complex(q, 0); +} + +complex +AMP_M(double s, double m1, double m2, double alpha1, double alpha2) +{ + // parameter + //matrix > T(2, 2); + vector a(2, TMatrixD(2,2)); + vector c(5, TMatrixD(2,2)); + vector pol(3, TMatrixD(1,2)); + TMatrixD sP(1,2); + int vesSheet(1); + int kachaev(1); + + const double f[2] = {0.1968, -0.0154}; // AMP Table 1, M solution: f_1^1 and f_2^1 + + a[0](0, 0) = 0.1131; // AMP Table 1, M solution: f_2^2 // schmarrn, steht eigentlich fuer a11 + a[0](0, 1) = 0.0150; // AMP Table 1, M solution: f_1^3 // a12 + a[0](1, 0) = 0.0150; // AMP Table 1, M solution: f_1^3 // a12 + a[0](1, 1) = -0.3216; // AMP Table 1, M solution: f_2^3 // a22 + a[1](0, 0) = f[0] * f[0]; + a[1](0, 1) = f[0] * f[1]; + a[1](1, 0) = f[1] * f[0]; + a[1](1, 1) = f[1] * f[1]; + + c[0](0, 0) = 0.0337; // AMP Table 1, M solution: c_11^0 + c[1](0, 0) = -0.3185; // AMP Table 1, M solution: c_11^1 + c[2](0, 0) = -0.0942; // AMP Table 1, M solution: c_11^2 + c[3](0, 0) = -0.5927; // AMP Table 1, M solution: c_11^3 + c[4](0, 0) = 0.1957; // AMP Table 1, M solution: c_11^4 + c[0](0, 1) = c[0](1, 0) = -0.2826; // AMP Table 1, M solution: c_12^0 + c[1](0, 1) = c[1](1, 0) = 0.0918; // AMP Table 1, M solution: c_12^1 + c[2](0, 1) = c[2](1, 0) = 0.1669; // AMP Table 1, M solution: c_12^2 + c[3](0, 1) = c[3](1, 0) = -0.2082; // AMP Table 1, M solution: c_12^3 + c[4](0, 1) = c[4](1, 0) = -0.1386; // AMP Table 1, M solution: c_12^4 + c[0](1, 1) = 0.3010; // AMP Table 1, M solution: c_22^0 + c[1](1, 1) = -0.5140; // AMP Table 1, M solution: c_22^1 + c[2](1, 1) = 0.1176; // AMP Table 1, M solution: c_22^2 + c[3](1, 1) = 0.5204; // AMP Table 1, M solution: c_22^3 + c[4](1, 1) = -0.3977; // AMP Table 1, M solution: c_22^4 + + pol[0](0, 0) = 0.1393; // AMP Table 1, M solution: a_1^0 + pol[1](0, 0) = -0.02775; // AMP Table 1, M solution: a_1^1 + pol[2](0, 0) = 0.3952; // AMP Table 1, M solution: a_1^2 + pol[0](0, 1) = 3.241; // AMP Table 1, M solution: a_2^0 + pol[1](0, 1) = -3.432; // AMP Table 1, M solution: a_2^1 + pol[2](0, 1) = 1.141; // AMP Table 1, M solution: a_2^2 + + sP(0, 0) = -0.0074; // AMP Table 1, M solution: s_0 + sP(0, 1) = 0.9828; // AMP Table 1, M solution: s_1 + + if (kachaev){ + // change parameters according to Kachaev's prescription + c[4](0, 0) = 0; // was 0.1957; + c[4](1, 1) = 0; // was -0.3977; + a[0](0, 1) = 0; // was 0.0150 + a[0](1, 0) = 0; // was 0.0150 + // a[1] are the f's from the AMP paper + a[1](0, 0) = 0; + a[1](0, 1) = 0; + a[1](1, 0) = 0; + a[1](1, 1) = 0; + } + + const complex imag(0, 1); + double m = sqrt(s); - double g2 = 4.21 * g1; + if (fabs(s - sP(0, 1)) < 1e-6) { + m += 1e-6; + s = m * m; + } + + double mKm = (mK+mKs)/2.; + + const complex qPiPi = breakupMomentumComplex(s, mPi, mPi); + const complex qPi0Pi0 = breakupMomentumComplex(s, mPi0, mPi0); + const complex qKK = breakupMomentumComplex(s, mK, mK); + const complex qK0K0 = breakupMomentumComplex(s, mKs, mKs); + complex qKmKm = breakupMomentumComplex(s, mKm, mKm ); + + vector > rho(2); + if ( vesSheet ) { + if (qKmKm.imag() > 0) + qKmKm *= -1; + rho[0] = (2. * qPiPi) / m; + rho[1] = (2. * qKmKm) / m; + } else { + rho[0] = ((2. * qPiPi) / m + (2. * qPi0Pi0) / m) / 2.; + rho[1] = ((2. * qKK) / m + (2. * qK0K0) / m) / 2.; + } - double rho_pi = 2 * breakupMomentum(s, m1, m2) / m; - double rho_K = 0; - double rho_K_cont = 0; - if ( s > 4*mK*mK ) - rho_K = 2 * breakupMomentum(s, mK, mK) / m; + const double scale = (s / (4 * mKm * mKm)) - 1; + + vector > M(4); + for (unsigned int i = 0; i < 4; ++i) { + // translate 1x4 into 2x2 + int x = i/2; + int y = i%2; + + for (unsigned int j = 0; j < 2; ++j) { + const complex fa = 1. / (s - sP(0, j)); + M[i] += fa * a[j](x,y); + } + + for (unsigned int j = 0; j < 5; ++j) { + const complex sc = pow(scale, (int)j); + M[i] += sc * c[j](x,y); + } + } + + for (unsigned int i = 0; i < 4; ++i) { + int index = i/2; + if (i == 0 || i ==3) + M[i] -= imag*rho[index]; + } + + //cout << M[0] << " " << M[1] << " " << M[2] << " " << M[3] << endl; + + // modification: off-diagonal terms set to 0 + M[1] = 0; + M[2] = 0; + + //vector > T(4); + complex det = M[0]*M[3] - M[1]*M[2]; + //T[0] = M[3]/det; + //T[1] = -1*M[2]/det; + //T[2] = -1*M[1]/det; + //T[3] = M[0]/det; + complex alpha(alpha1, alpha2); + const complex amp = M[3]/det - alpha * M[2]/det; + + //if (_debug) + //printDebug << name() << "(m = " << maxPrecision(mass) << " GeV) = " + //<< maxPrecisionDouble(amp) << endl; + + return amp; +} + +complex +AMP_K1_my(double s, double m1, double m2) +{ + // parameter + //matrix > T(2, 2); + vector a(2, TMatrixD(2,2)); + vector c(5, TMatrixD(2,2)); + TMatrixD sP(1,2); + + sP(0, 0) = -0.0110; // AMP Table 1, K1 solution: s_0 + sP(0, 1) = 0.9247; // AMP Table 1, K1 solution: s_1 + + const double f[2] = {-0.2242, 0.5829}; // AMP Table 1, K1 solution: f_1^1 and f_2^1 + + c[0](0, 0) = 0.7347; // AMP Table 1, K1 solution: c_11^0 + c[1](0, 0) = -0.5266; // AMP Table 1, K1 solution: c_11^1 + c[2](0, 0) = 2.6151; // AMP Table 1, K1 solution: c_11^2 + c[3](0, 0) = -1.7747; // AMP Table 1, K1 solution: c_11^3 + c[4](0, 0) = 0.8031; // AMP Table 1, K1 solution: c_11^4 + c[0](0, 1) = c[0](1, 0) = -3.2762; // AMP Table 1, K1 solution: c_12^0 + c[1](0, 1) = c[1](1, 0) = -0.6662; // AMP Table 1, K1 solution: c_12^1 + c[2](0, 1) = c[2](1, 0) = 0.8778; // AMP Table 1, K1 solution: c_12^2 + c[3](0, 1) = c[3](1, 0) = -2.1190; // AMP Table 1, K1 solution: c_12^3 + c[4](0, 1) = c[4](1, 0) = 0.2319; // AMP Table 1, K1 solution: c_12^4 + c[0](1, 1) = -2.6785; // AMP Table 1, K1 solution: c_22^0 + c[1](1, 1) = 7.9951; // AMP Table 1, K1 solution: c_22^1 + c[2](1, 1) = 5.5763; // AMP Table 1, K1 solution: c_22^2 + c[3](1, 1) = -1.4956; // AMP Table 1, K1 solution: c_22^3 + c[4](1, 1) = 0 ; // AMP Table 1, K1 solution: c_22^4 + + const complex imag(0, 1); + + double m = sqrt(s); + + double mKm = (mK+mKs)/2.; + + const complex qPiPi = breakupMomentumComplex(s, mPi, mPi); + const complex qPi0Pi0 = breakupMomentumComplex(s, mPi0, mPi0); + const complex qKK = breakupMomentumComplex(s, mK, mK); + const complex qK0K0 = breakupMomentumComplex(s, mKs, mKs); + //complex qKmKm = breakupMomentumComplex(s, mKm, mKm ); + + vector > rho(2); + rho[0] = ((2. * qPiPi) / m + (2. * qPi0Pi0) / m) / 2.; + rho[1] = ((2. * qKK) / m + (2. * qK0K0) / m) / 2.; + + const double scale = (s / (4 * mKm * mKm)) - 1; + + vector > K(4); + //vector > Khat(4); + for (unsigned int i = 0; i < 4; ++i) { + // translate 1x4 into 2x2 + int x = i/2; + int y = i%2; + + K[i] = 1;//( s - sP(0,0) ) / (4*mK*mK); //wrong in paper, comes later + K[i] /= ( sP(0,1) - s ) * ( sP(0,1) - sP(0,0) ); + K[i] = f[x] * f[y]; + + //cout << K[0] << " " << K[1] << " " << K[2] << " " << K[3] << endl; + + for (unsigned int j = 0; j < 5; ++j) { + const complex sc = pow(scale, (int)j); + K[i] += sc * c[j](x,y); + } + + K[i] *= ( s - sP(0,0) ) / (4*mK*mK); + + //Khat[i] = K[i]/( s - sP(0,0) ); + } + + // invert + vector > M(4); + complex det_K = (K[0]*K[3] - K[1]*K[2]); + det_K = 1. / det_K; + M[0] += det_K * K[3]; + M[1] -= det_K * K[1]; + M[2] -= det_K * K[2]; + M[3] += det_K * K[0]; + + for (unsigned int i = 0; i < 4; ++i) { + int index = i/2; + if (i == 0 || i ==3) + M[i] -= imag*rho[index]; + } + + // invert + vector > T(4); + complex det_M = (M[0]*M[3] - M[1]*M[2]); + det_M = 1. / det_M; + T[0] += det_M * M[3]; + T[1] -= det_M * M[1]; + T[2] -= det_M * M[2]; + T[3] += det_M * M[0]; + + // modification: off-diagonal terms set to 0 + //T[1] = 0; + //T[2] = 0; + + const complex amp = T[0]; + + return amp; +} + +complex +AMP_Kred(double s, double m1, double m2, double alpha1, double alpha2) +{ + const complex zi(0, 1); + const double pi = 3.1415926; + //double pie = pi/180.; + //complex zie = zi*pie; + complex zieps = zi*0.00001; + //double pip = atan(1.)*4; + double empic2 = mPi*mPi*4.; + //double empi02 = mPi0*mPi0*4.; + //double emp = 0.137; // mPim??? + //double emrho2 = 0.77*0.77; // mRho??? + double emK02 = mKs*mKs*4.; + double emKc2 = mK*mK*4.; + //double emK2 = 0.4956*0.4956*4.; // mKm??? + //double ehadbot = 0.28; + //double ehadtop = 1.6; + + // parameter + vector c(3, TMatrixD(2,2)); + TMatrixD f(2,2); + + double s0 = 0.41*empic2/4.; // S0 + const double S[3] = {s0, 0.26261, 1.0811}; // S0, S1, S2 + //const double S[3] = {s0, 0.26261, alpha1}; // S0, S1, S2 + + f(0,0) = 0.38949; // F11 + f(0,1) = 0.33961; // F12 + f(1,0) = 0.24150; // F21 + f(1,1) = -0.78538; // F22 + + c[0](0, 0) = 0.14760; // C110 + c[1](0, 0) = 0.062181; // C111 + c[2](0, 0) = 0.029465; // C112 + c[0](0, 1) = c[0](1, 0) = 0.10914; // C120 + c[1](0, 1) = c[1](1, 0) = -0.17912; // C121 + c[2](0, 1) = c[2](1, 0) = 0.10758; // C122 + c[0](1, 1) = -0.27253; // C220 + c[1](1, 1) = 0.79442; // C221 + c[2](1, 1) = -0.49529; // C222 + + + complex zs = s + zieps; + complex zrho1 = sqrt(1. - empic2/zs); // does that work?? + //double rho1 = sqrt(1. - empic2/s); + //double rho2; + complex zrho2 = sqrt(1. - emK02/zs)*0.5 + sqrt(1. - emKc2/zs)*0.5; + complex zf1 = zrho1/pi * log((zrho1+1.)/(zrho1-1.)); + complex zf2 = zrho2/pi * log((zrho2+1.)/(zrho2-1.)); + /*if (s > emK02) + rho2 = abs(zrho2); else - rho_K_cont = sqrt(4*mK*mK/s - 1); + rho2 = 0;*/ + double q2 = s/emKc2 - 1; + double q4 = q2*q2; + double q[3] = {1,q2,q4}; + + vector AQL(4); + for (unsigned int i = 0; i < 4; ++i) { + // translate 1x4 into 2x2 + int x = i/2; + int y = i%2; + + AQL[i] = 0; + for (unsigned int j = 0; j < 3; ++j) + AQL[i] += q[j] * c[j](x,y); + } + + vector APL1(4); + APL1[0] = f(0,0)*f(0,0); + APL1[1] = f(0,0)*f(1,0); + APL1[2] = f(0,0)*f(1,0); + APL1[3] = f(1,0)*f(1,0); + for (unsigned int i = 0; i < 4; ++i) + APL1[i] /= (S[1] - s)*(S[1] - S[0]); + + vector APL2(4); + APL2[0] = f(0,1)*f(0,1); + APL2[1] = f(0,1)*f(1,1); + APL2[2] = f(0,1)*f(1,1); + APL2[3] = f(1,1)*f(1,1); + for (unsigned int i = 0; i < 4; ++i) + APL2[i] /= (S[2] - s)*(S[2] - S[0]); + + vector ALL(4); + for (unsigned int i = 0; i < 4; ++i){ + ALL[i] = APL1[i] + APL2[i] + AQL[i]; + ALL[i] *= (s - S[0])/emKc2; + } + double detLL = ALL[0]*ALL[3] - ALL[1]*ALL[2]; + + complex denominator = 1. + zf1 * ALL[0] + zf2 * ALL[3] + zf1 * zf2 * detLL; + + vector > T(4); + for (unsigned int i = 0; i < 4; ++i) + T[i] = ALL[i]; + T[0] += zf2 * detLL; + T[3] += zf1 * detLL; + + for (unsigned int i = 0; i < 4; ++i){ + T[i] /= denominator; + T[i] *= emKc2/(s - S[0]); + } + + complex alpha(alpha1, alpha2); + //complex alpha(0.09, -0.5); + const complex amp = T[0] + alpha * T[1]; - double numerator = 1.; + //if (sqrt(s) < 1.6) + return amp; + //else return zi; +} + +complex +Flatte(double s, double m1, double m2, double m0, double g1, double g2) +{ + double m = sqrt(s); + + const complex imag(0, 1); + + double rho_pi = 2. * breakupMomentum(s, m1, m2) / m; + complex rho_K = 2. * breakupMomentumComplex(s, mK, mK) / m; + + double numerator = m0 * g1*rho_pi; - complex denominator = complex(m0*m0 - s + g2*rho_K_cont, - g1*rho_pi - g2*rho_K); + complex denominator = m0*m0 - s - imag * m0 * ( g1*rho_pi + g2*rho_K); return numerator / denominator; } diff --git a/bw.h b/bw.h index 1ccdaaf..7f25748 100644 --- a/bw.h +++ b/bw.h @@ -17,6 +17,9 @@ const double hbarc = 0.1973269631; // GeV fm double blattWeisskopf(int L, double p); double breakupMomentum(double s, double m1, double m2); std::complex Flatte(double s, double m1, double m2, double m0, double g1); +std::complex AMP_M(double s, double m1, double m2, double alpha1, double alpha2); +std::complex AMP_K1(double s, double m1, double m2); +std::complex AMP_Kred(double s, double m1, double m2, double alpha1, double alpha2); std::complex BW(double s, double m1, double m2, int J, double m0, double Gamma0); std::complex BW_a2_pietap(double s); std::complex BW_a2_pietap_coupled(double s);