#include "HADRONS++/Current_Library/VA_0_PPP.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "ATOOLS/Org/MyStrStream.H"
#include "ATOOLS/Org/My_MPI.H"
using namespace HADRONS;
using namespace ATOOLS;
using namespace std;
VA_0_PPP::FF_Base::~FF_Base()
{
}
VA_0_PPP::VA_0_PPP(const ATOOLS::Flavour_Vector& flavs,
const std::vector<int>& indices, const std::string& name) :
Current_Base(flavs, indices, name)
{
int nPion_0(0), nPion_ch(0), nKaon_0(0), nKaon_ch(0), nKaon_S (0), nKaon_L(0);
// count number of pions, kaons and calc. mass^2
for( int i=0; i<3; i++ ) {
if( m_flavs[p_i[i]].Kfcode() == kf_pi_plus ) nPion_ch++;
if( m_flavs[p_i[i]].Kfcode() == kf_pi ) nPion_0++;
if( m_flavs[p_i[i]].Kfcode() == kf_K_plus ) nKaon_ch++;
if( m_flavs[p_i[i]].Kfcode() == kf_K ) nKaon_0++;
if( m_flavs[p_i[i]].Kfcode() == kf_K_L ) nKaon_L++;
if( m_flavs[p_i[i]].Kfcode() == kf_K_S ) nKaon_S++;
m_ms[i] = sqr( m_flavs[p_i[i]].HadMass() );
}
// sanity check
if (nPion_ch+nPion_0+nKaon_ch+nKaon_L+nKaon_S+nKaon_0 != 3) {
msg_Error()<<"ERROR in HADRONS::VA_0_PPP constructor\n"
<<" number of three outgoing pseudoscalars != 3 ?!.\n"
<<" Don't know, what to do. Will abort."<<endl;
Abort();
}
// define mode number
m_mode = nPion_ch*1000 + nPion_0*100 + nKaon_ch*10 + (nKaon_L+nKaon_S+nKaon_0);
m_kaon_mode = 100*nKaon_L + 10*nKaon_S + nKaon_0;
}
void VA_0_PPP::SetModelParameters( struct GeneralModel _md )
{
m_Vud = _md("Vud", Tools::Vud);
m_Vus = _md("Vus", Tools::Vus);
double fpi = _md("fpi", 0.1307 );
m_B123 = 0.;
double A123;
switch( m_mode ) {
case 3000 : /* pi- pi- pi+ mode */
A123 = m_Vud*SQRT_05;
m_B123 = 0.;
break;
case 1200 : /* pi0 pi0 pi- mode */
A123 = m_Vud*SQRT_05;
m_B123 = 0.;
break;
case 1020 : /* K- pi- K+ */
m_B123 = 2.;
A123 = -0.5*m_Vud;
break;
case 1002 : /* K0 pi- K0b */
m_B123 = -2.;
A123 = -0.5*m_Vud;
break;
case 111 : /* K- pi0 K0 */
A123 = 1.5*m_Vud*SQRT_05;
m_B123 = -2./3.;
break;
case 210 : /* pi0 pi0 K- mode */
A123 = m_Vus/4.;
m_B123 = 1.;
break;
case 2010 : /* K- pi- pi+ */
m_B123 = -1.;
A123 = -0.5*m_Vus;
break;
case 1101 : /* pi- K0b pi0 */
m_B123 = -2/3.;
A123 = 1.5*m_Vus*SQRT_05;
break;
case 30 : /* K- K- K+ */
A123 = m_Vus;
break;
case 12 : /* K- K0b K0 */
A123 = -0.5*m_Vus;
break;
default : msg_Error()<<"Warning in HADRONS::Tau_Decay_MEs.C in Tau_Three_Pseudo::GetA123() :\n"
<<" Obviously this three pseudoscalar channel (code "<<m_mode<<")\n"
<<" doesn't have a global A123. Maybe it is not implemented yet.\n"
<<" Take A123=1., will continue and hope for the best."<<endl;
A123 = 1.;
break;
}
m_global = 4./(3.*fpi)*A123;
switch( int(_md("FORM_FACTOR", 1)) ) {
case 2 : p_ff = new RChT(m_mode,m_kaon_mode,_md,m_ms);
break;
case 1 : p_ff = new KS95(m_mode,m_kaon_mode,_md,m_ms);
break;
case 3 : p_ff = new KS(m_mode,m_kaon_mode,_md,m_ms);
break;
}
}
void VA_0_PPP::Calc(const ATOOLS::Vec4D_Vector& moms, bool m_anti)
{
Vec4D p1( moms[p_i[0]] ),
p2( moms[p_i[1]] ),
p3( moms[p_i[2]] );
Vec4D Q( p1+p2+p3 );
double s = (p1+p3).Abs2(),
t = (p2+p3).Abs2();
double Q2 = Q.Abs2();
double dot1 = Q*(p1-p3),
dot2 = Q*(p2-p3);
double d1 = dot1/Q2,
d2 = dot2/Q2;
Complex F1 = FormFactor( 1, Q2, s, t );
Complex F2 = FormFactor( 2, Q2, s, t );
Complex FS = FormFactor( 3, Q2, s, t );
Complex FV = m_B123*FormFactor( 4, Q2, s, t );
// F1 = Complex(0.,0.);
// F2 = Complex(0.,0.);
// FS = Complex(0.,0.);
Insert( m_global*( FS+F1*(1.-d1)-F2*d2 )*Vec4C(p1)
+ m_global*( FS-F1*d1+F2*(1.-d2) )*Vec4C(p2)
+ m_global*( FS-F1*(1.+d1)-F2*(1.+d2) )*Vec4C(p3)
+ m_global*Complex(0.,1.)*FV*Vec4C(cross(p1,p2,p3)) , 0);
// for( int h=0; h<4; h++ ) { // helicity comb. (nutau,tau)
// _ampls_tensor->push_back(
// ( F.X( m_nutau, m_pseudo_1, 0, h, m_cR, m_cL ) * ( FS + F1*(1.-d1) - F2*d2 )
// + F.X( m_nutau, m_pseudo_2, 0, h, m_cR, m_cL ) * ( FS - F1*d1 + F2*(1.-d2) )
// + F.X( m_nutau, m_pseudo_3, 0, h, m_cR, m_cL ) * ( FS - F1*(1.+d1) - F2*(1.+d2) )
// + Complex(0.,1.)*F.X( m_nutau, cross(p1,p2,p3), 0, h, m_cR, m_cL ) * FV
// ) * m_global
// );
// }
}
Complex VA_0_PPP::FormFactor( int j, double Q2, double s, double t )
{
return( p_ff->FormFactor(j,Q2,s,t) );
}
VA_0_PPP::FF_Base::FF_Base(int mode, int kaon_mode, GeneralModel _md, double * _ms)
: m_mode (mode), m_kaon_mode (kaon_mode)
{
// set resonances and parameters
m_ms[0] = _ms[0];
m_ms[1] = _ms[1];
m_ms[2] = _ms[2];
m_deltas = 0;
kf_code resA, resAA, resV[2], resVV[2];
switch( m_mode ) {
case 3000 : /* pi- pi- pi+ mode */
resV[0] = kf_rho_770;
resV[1] = kf_rho_770;
break;
case 1200 : /* pi0 pi0 pi- mode */
resV[0] = kf_rho_770_plus;
resV[1] = kf_rho_770_plus;
break;
case 1020 : /* K- pi- K+ */
resV[0] = kf_rho_770;
resV[1] = kf_K_star_892;
break;
case 1002 : /* K0 pi- K0b */
resV[0] = kf_rho_770;
resV[1] = kf_K_star_892_plus;
break;
case 111 : /* K- pi0 K0 */
resV[0] = kf_rho_770_plus;
resV[1] = kf_K_star_892;
break;
case 210 : /* pi0 pi0 K- mode */
m_deltas= 1;
resV[0] = kf_K_star_892_plus;
resV[1] = kf_K_star_892_plus;
break;
case 2010 : /* K- pi- pi+ */
m_deltas= 1;
resV[0] = kf_K_star_892;
resV[1] = kf_rho_770;
break;
case 1101 : /* pi- K0b pi0 */
m_deltas= 1;
resV[0] = kf_rho_770_plus;
resV[1] = kf_K_star_892;
break;
case 30 : /* K- K- K+ */
m_deltas= 1;
resV[0] = kf_rho_770;
resV[1] = kf_rho_770;
break;
case 12 : /* K- K0b K0 */
m_deltas= 1;
resV[0] = kf_rho_770_plus;
resV[1] = kf_rho_770_plus;
break;
default:
THROW(fatal_error, "Three pseudoscalar mode not recognized: "+ToString(m_mode));
break;
}
// higher resonances
for( int i=0; i<2; ++i ) {
if( resV[i] == kf_rho_770 ) resVV[i] = kf_rho_1450;
if( resV[i] == kf_rho_770_plus ) resVV[i] = kf_rho_1450_plus;
if( resV[i] == kf_K_star_892 ) resVV[i] = kf_K_star_1410;
if( resV[i] == kf_K_star_892_plus ) resVV[i] = kf_K_star_1410_plus;
}
// axial resonance
int running = int( _md("RUNNING_WIDTH", 3 ) ); // running width
double mA, mAA, wA, wAA;
if( m_deltas ) {
running = 0;
resA = kf_K_1_1400_plus;
resAA = kf_K_1_1270_plus;
mA = _md("Mass_"+Flavour(resA).IDName(), 1.402 ); // mass of axial resonance
mAA = _md("Mass_"+Flavour(resAA).IDName(), 1.270 ); // mass of axial resonance'
wA = _md("Width_"+Flavour(resA).IDName(), 0.174 ); // width of axial resonance
wAA = _md("Width_"+Flavour(resAA).IDName(), 0.090 ); // width of axial resonance'
m_alpha = _md("alpha_"+Flavour(resAA).IDName(), 1.3 ); // weight factor for A'
}
else {
resA = kf_a_1_1260_plus;
resAA = kf_a_1_1260_plus;
mA = _md("Mass_"+Flavour(resA).IDName(), 1.251 ); // mass of axial resonance
mAA = _md("Mass_"+Flavour(resAA).IDName(), 1.251 ); // mass of axial resonance'
wA = _md("Width_"+Flavour(resA).IDName(), 0.599 ); // width of axial resonance
wAA = _md("Width_"+Flavour(resAA).IDName(), 0.599 ); // width of axial resonance'
m_alpha = _md("alpha_"+Flavour(resAA).IDName(), 0. ); // weight factor for A'
}
m_A = ResonanceFlavour( resA, mA, wA, running&1);
m_AA = ResonanceFlavour( resAA, mAA, wAA, running&1);
ResonanceFlavour help_rho = ResonanceFlavour(
kf_rho_770_plus,
_md("Mass_rho(770)+", _md("Mass_rho(770)", 0.773 )),
_md("Width_rho(770)+", _md("Width_rho(770)", 0.145 )),
running&2
);
ResonanceFlavour help_rhoP = ResonanceFlavour(
kf_rho_1450_plus,
_md("Mass_rho(1450)+", _md("Mass_rho(1450)", 1.370 )),
_md("Width_rho(1450)+", _md("Width_rho(1450)", 0.510 )),
running&2
);
double help_beta = _md("beta_rho(1450)+", _md("beta_rho(1450)", -0.145) );
m_A.InitialiseThreeBodyResonance(help_rho, help_rhoP, help_beta);
m_AA.InitialiseThreeBodyResonance(help_rho, help_rhoP, help_beta);
double mV[2], mVV[2], wV[2], wVV[2];
for( int i=0; i<2; ++i ) {
if( resV[i] == kf_rho_770 || resV[i] == kf_rho_770_plus ) {
mV[i] = _md("Mass_"+Flavour(resV[i]).IDName(), 0.773 ); // mass^2 of vector resonance ij
wV[i] = _md("Width_"+Flavour(resV[i]).IDName(), 0.145 ); // width^2 of vector resonance ij
mVV[i] = _md("Mass_"+Flavour(resVV[i]).IDName(), 1.370 ); // mass^2 of vector resonance' ij
wVV[i] = _md("Width_"+Flavour(resVV[i]).IDName(),0.510 ); // width^2 of vector resonance' ij
m_Beta[i] = _md("beta_"+Flavour(resVV[i]).IDName(), -0.145 ); // weight factor for Vij'
}
else if( resV[i] == kf_K_star_892 || resV[i] == kf_K_star_892_plus ) {
mV[i] = _md("Mass_"+Flavour(resV[i]).IDName(), 0.8921 ); // mass^2 of vector resonance ij
wV[i] = _md("Width_"+Flavour(resV[i]).IDName(), 0.0513 ); // width^2 of vector resonance ij
mVV[i] = _md("Mass_"+Flavour(resVV[i]).IDName(), 1.412 ); // mass^2 of vector resonance' ij
wVV[i] = _md("Width_"+Flavour(resVV[i]).IDName(),0.227 ); // width^2 of vector resonance' ij
m_Beta[i] = _md("beta_"+Flavour(resVV[i]).IDName(), -0.135 ); // weight factor for Vij'
}
else {
mV[i] = _md("Mass_"+Flavour(resV[i]).IDName(), Flavour(resV[i]).HadMass()); // mass^2 of vector resonance ij
wV[i] = _md("Width_"+Flavour(resV[i]).IDName(), Flavour(resV[i]).Width()); // width^2 of vector resonance ij
mVV[i] = _md("Mass_"+Flavour(resVV[i]).IDName(), Flavour(resVV[i]).HadMass()); // mass^2 of vector resonance' ij
wVV[i] = _md("Width_"+Flavour(resVV[i]).IDName(),Flavour(resVV[i]).Width()); // width^2 of vector resonance' ij
m_Beta[i] = _md("beta_"+Flavour(resVV[i]).IDName(), 0. ); // weight factor for Vij'
}
m_V[i] = ResonanceFlavour( resV[i], mV[i], wV[i], running&2 );
m_VV[i] = ResonanceFlavour( resVV[i], mVV[i], wVV[i], running&2 );
}
m_fpi2 = sqr( _md("fpi", 0.1307) ); // pion decay constant
}
VA_0_PPP::RChT::RChT(int mode, int kaon_mode, GeneralModel _md, double * _ms)
: FF_Base(mode,kaon_mode,_md,_ms)
{
// set resonances and parameters
if( m_mode != 1200 &&
m_mode != 3000 &&
m_mode != 1020 ) {
msg_Error()<<"Error: The mode "<<m_mode<<endl
<<" hasn't been implemented yet (RChT). Please use KS model."
<<" Don't know what to do. Will abort"<<endl;
Abort();
}
m_MO = _md("Mass_omega(782)", Flavour(kf_omega_782).HadMass());
m_MO2 = sqr(m_MO);
m_GO = _md("Width_omega(782)", Flavour(kf_omega_782).Width());
m_l0 = _md("lambda0", 1.); // fit parameter lambda0
m_gammaR = _md("gamma_rho(770)", 1.); // global factor for rho width
m_m = Flavour( kf_pi_plus ).HadMass(); // pion mass
m_m2 = sqr(m_m); // pion mass^2
m_mK2 = sqr( Flavour( kf_K_plus ).HadMass() ); // Kaon mass^2
m_exp_alpha = _md("exp_alpha", 2.45); // exponent in off-shell GA
m_l1 = _md("lambda1", 0.5); // fit parameter
m_l2 = _md("lambda2", 0.); // fit parameter
m_lsum = m_l1 + m_l2;
// constraints due to short-distance behaviour
m_FV2 = 2*m_fpi2; // vector coupling
m_FV = sqrt(m_FV2);
m_GV = m_FV/2.;
m_FA2 = m_FV2 - m_fpi2; // axial coupling
m_FA = sqrt(m_FA2);
}
double VA_0_PPP::RChT::MassWidthVector( int a, double s )
{
if(m_V[0].Running()) {
switch(m_mode) {
case 3000:
case 1200: return MassWidthVector( s );
case 1020: return (a==0)? MassWidthVector(s) : m_V[1].MassWidth();
}
// default:
msg_Error()<<"Warning: this form factor (RChT) for the three-pseudoe mode "<<m_mode<<"\n"
<<" hasn't been implemented yet. Please use KS model."<<endl;
}
return( m_V[a].MassWidth() );
}
double VA_0_PPP::RChT::MassWidthVector( double s )
{
double MVGV (0.);
if( s>4.*m_m2 ) MVGV += pow( 1.-4.*m_m2/s, 1.5 );
if( s>4.*m_mK2 ) MVGV += pow( 1.-4.*m_mK2/s, 1.5 ) / 2.;
MVGV *= m_gammaR*m_V[0].Mass2()*s/(96.*M_PI*m_fpi2);
return MVGV;
}
double VA_0_PPP::RChT::MassWidthAxial( double Q2 )
{
if( !m_deltas && m_A.Running() )
return( m_A.OffShellMassWidth(Q2)*pow(m_A.Mass2()/Q2,m_exp_alpha-2.) );
return m_A.MassWidth();
}
Complex VA_0_PPP::RChT::FormFactor( int j, double Q2, double s, double t )
{
switch( m_mode ) {
case 1200:
case 3000: { // 3pion mode
if (j==1 || j==2) { // axial contributions
double u = Q2-s-t+Mass2(1-1)+Mass2(2-1)+Mass2(3-1);
double x = (j==1)? s : t;
double y = (j==1)? t : s;
double MVGV_x = MassWidthVector(x);
double MVGV_y = MassWidthVector(y);
double MAGA = MassWidthAxial(Q2);
double F_Q2_x = x/2./Q2 - m_l0*m_m2/Q2;
double F_Q2_y = y/2./Q2 - m_l0*m_m2/Q2;
double MV2 = m_V[0].Mass2();
Complex alpha = 1. - 3./2.*x/Complex(x-MV2, MVGV_x);
Complex beta = -3./2.*x/Complex(x-MV2, MVGV_x)
+ F_Q2_x*(2.*Q2+x-u)/Complex(x-MV2, MVGV_x)
+ F_Q2_y*(u-x)/Complex(y-MV2, MVGV_y);
return alpha - Q2/Complex(Q2-m_A.Mass2(),MAGA)*beta;
}
else { // pseudoscalar, vector
return Complex(0.,0.);
}
}
case 1020: { // K- pi- K+ mode
double u = Q2-s-t+Mass2(1-1)+Mass2(2-1)+Mass2(3-1);
switch(j) {
case 1:
case 2: { // axial contributions
double x = (j==1)? s : t;
double y = (j==1)? t : s;
int a = j-1; // resonance assoc. with x
int b = (a==1)? 0 : 1; // resonance assoc. with y
Complex FAchi = Complex(1.,0.);
Complex FA1r = 0.5 * m_FV*m_GV/m_fpi2 * (
1./Complex(m_V[a].Mass2()-x,-1.*MassWidthVector(a,x))*
( 3.*x + Mass2(2)-Mass2(a) + (1.-2.*m_GV/m_FV)*(2.*Q2-2.*x-u+Mass2(a)-Mass2(b)) )
+1./Complex(m_V[b].Mass2()-y,-1.*MassWidthVector(b,y))*
( 2.*(Mass2(b)-Mass2(2)) + (1.-2.*m_GV/m_FV)*(u-x+Mass2(2)-Mass2(b)) )
);
Complex FA2r = -sqrt(2.) * m_FA*m_GV/m_fpi2 * Q2/Complex(m_A.Mass2()-Q2,-1.*MassWidthAxial(Q2)) * (
1./Complex(m_V[a].Mass2()-x,-1.*MassWidthVector(a,x))*
( m_lsum*(-3.*x+Mass2(a)-Mass2(2)) + FFunc(x,Mass2(b),Q2)*(2.*Q2+x-u+Mass2(2)-Mass2(b)) )
+1./Complex(m_V[b].Mass2()-y,-1.*MassWidthVector(b,y))*
( 2.*m_lsum*(Mass2(2)+Mass2(b)) + FFunc(y,Mass2(a),Q2)*(u-x+Mass2(b)-Mass2(2)) )
);
return FAchi + FA1r + FA2r;
}
case 3: { // pseudoscalar contribution
Complex FSchi = 3./2.* Mass2(1)/(Q2-Mass2(1)) *
( 1. + (Mass2(2)-u)/Q2 );
Complex FS1r = 3./2.*sqr(m_GV)/m_fpi2 * Mass2(1)/(Q2*(Q2-Mass2(1))) *
( s*(t-u)/Complex(m_V[0].Mass2()-s,-1.*MassWidthVector(0,s))
+ (t*(s-u)+(Q2-Mass2(0))*(Mass2(1)-Mass2(2)))/Complex(m_V[1].Mass2()-t,-1.*MassWidthVector(1,t))
);
return FSchi + FS1r;
}
case 4: { // vector contribution
double c1c2c5 = 0.;
// double c1c2c52c6 = -3./(32.*sqr(M_PI)) * m_V[0].Mass()/(sqrt(2.)*m_FV);
double c1c2c52c6 = -3./(96.*sqr(M_PI)) * m_V[0].Mass()/(sqrt(2.)*m_GV);
double c1c28c3c5 = 0.;
double c4 = 0.;
// double d3 = -3./(64.*sqr(M_PI)) * m_V[0].Mass2()/m_FV2 + m_fpi2/(8.*m_FV2);
double d3 = -3./(192.*sqr(M_PI)) * m_V[0].Mass2()/(m_FV*m_GV);
double d18d2d3 = m_fpi2/(8.*m_FV2);
double g12g2g3 = 0.;
double g2 = 3./(192.*sqr(M_PI)) * m_V[0].Mass()/(sqrt(2.)*m_FV);
Complex FVchi = 3./(2.*sqr(M_PI)*m_fpi2);
Complex FV1r1 = 6.*m_GV/(sqrt(2.)*m_fpi2*m_V[0].Mass()) * (
1./Complex(m_MO2-s,-1.*m_MO*m_GO) * (
c1c2c5*Q2 - c1c2c52c6*s + c1c28c3c5*Mass2(1) ) +
1./Complex(m_V[1].Mass2()-t,-1.*MassWidthVector(1,t)) * (
c1c2c5*Q2 - c1c2c52c6*t + c1c28c3c5*Mass2(0) + 8.*c4*(Mass2(0)-Mass2(1)) )
);
Complex FV1r2 = -12.*m_FV/(sqrt(2.)*m_V[0].Mass()*m_fpi2) *
1./Complex(m_V[0].Mass2()-Q2,-1*MassWidthVector(0,Q2)) * (
g12g2g3*(s+t)-2.*g2*Q2 );
Complex FV2r = -6.*m_FV*m_GV/m_fpi2 *
1./Complex(m_V[0].Mass2()-Q2,-1.*MassWidthVector(0,Q2)) * (
1./Complex(m_MO2-s,-1.*m_MO*m_GO) * (
d3*(Q2+s) + d18d2d3*Mass2(1) ) +
1./Complex(m_V[1].Mass2()-t,-1.*MassWidthVector(1,t)) * (
d3*(Q2+t) + d18d2d3*Mass2(0) )
);
return FVchi + FV1r1 + FV1r2 + FV2r;
}
}
}
}
return (j>3)? Complex(0.,0.) : Complex(1.,0.);
}
double VA_0_PPP::RChT::FFunc( double a, double b, double c)
{
return m_l2 + m_l1*a/c - m_l0*b/c;
}
// Parameterisation
// DECKER, FINKEMEIER, MIRKES hep-ph/9310270
VA_0_PPP::KS::KS(int mode, int kaon_mode, GeneralModel _md, double * _ms)
: FF_Base(mode,kaon_mode,_md,_ms)
{
bool anomaly (false);
// special constants for this particular parameterisation
switch( m_mode ) {
case 3000 : /* pi- pi- pi+ mode */
m_X123 = 2.*Mass2(0);
m_ms123 = Mass2(0);
m_G123 = 1;
break;
case 1200 : /* pi0 pi0 pi- mode */
m_X123 = Mass2(2);
m_ms123 = Mass2(2);
m_G123 = 1;
break;
case 1020 : /* K- pi- K+ */
m_X123 = Mass2(1) + Mass2(0);
m_ms123 = Mass2(1);
m_G123 = 1;
anomaly = true;
break;
case 1002 : /* K0 pi- K0b */
m_X123 = Mass2(1) + Mass2(0);
m_ms123 = Mass2(1);
m_G123 = 1;
anomaly = true;
break;
case 111 : /* K- pi0 K0 */
m_X123 = 0.;
m_ms123 = Mass2(1);
m_G123 = 0;
break;
case 210 : /* pi0 pi0 K- mode */
m_X123 = -2.*(Mass2(1)+Mass2(2));
m_ms123 = Mass2(2);
m_G123 = 1;
break;
case 2010 : /* K- pi- pi+ */
m_X123 = Mass2(1)+Mass2(0);
m_ms123 = Mass2(0);
m_G123 = 1;
anomaly = true;
break;
case 1101 : /* pi- K0b pi0 */
m_X123 = 0.;
m_ms123 = Mass2(1);
m_G123 = 0;
anomaly = true;
break;
case 30 : /* K- K- K+ */
m_X123 = 2.*Mass2(2);
m_ms123 = Mass2(2);
m_G123 = 1;
break;
case 12 : /* K- K0b K0 */
m_X123 = 2.*Mass2(0);
m_ms123 = Mass2(0);
m_G123 = 1;
break;
default:
msg_Error()<<METHOD<<" Warning: Three pseudoscalar mode not recognized. "
<<"mode="<<m_mode<<endl;
break;
}
kf_code resAnoV = (m_deltas)? kf_K_star_892_plus : kf_rho_770_plus;
kf_code resAnoVV = (m_deltas)? kf_K_star_1410_plus : kf_rho_1450_plus;
kf_code resAnoVVV = (m_deltas)? kf_K_star_1680_plus : kf_rho_1700_plus;
if( anomaly ) {
int running = int( _md("RUNNING_WIDTH", 3 ) ); // running width
double MV = _md("Mass_anomaly_"+Flavour(resAnoV).IDName(), Flavour(resAnoV).HadMass() ); // mass V
double MVV = _md("Mass_anomaly_"+Flavour(resAnoVV).IDName(), Flavour(resAnoVV).HadMass() ); // mass V'
double MVVV = _md("Mass_anomaly_"+Flavour(resAnoVVV).IDName(), Flavour(resAnoVVV).HadMass() ); // mass V'
double GV = _md("Width_anomaly_"+Flavour(resAnoV).IDName(), Flavour(resAnoV).Width() ); // width V
double GVV = _md("Width_anomaly_"+Flavour(resAnoVV).IDName(), Flavour(resAnoVV).Width() ); // width V'
double GVVV = _md("Width_anomaly_"+Flavour(resAnoVVV).IDName(), Flavour(resAnoVVV).Width() ); // width V'
m_AnoV = ResonanceFlavour( resAnoV, MV, GV, running&2 );
m_AnoVV = ResonanceFlavour( resAnoVV, MVV, GVV, running&2 );
m_AnoVVV = ResonanceFlavour( resAnoVVV, MVVV, GVVV, running&2 );
m_AlphaV = _md("alpha_anomaly_K*(892)", _md("alpha_anomaly_K*(892)+", 0. )); // alpha for K*
m_BetaV[0] = _md("beta_anomaly_"+Flavour(resAnoVV).IDName(), 0. ); // beta for V'
m_BetaV[1] = _md("gamma_anomaly_"+Flavour(resAnoVVV).IDName(), 0. ); // gamma for V''
}
else {
m_AnoV = ResonanceFlavour( kf_rho_770_plus, 0., 0., 0 );
m_AnoVV = ResonanceFlavour( kf_rho_770_plus, 0., 0., 0 );
m_AnoVVV = ResonanceFlavour( kf_rho_770_plus, 0., 0., 0 );
m_AlphaV = 0.;
m_BetaV[0] = 0.;
m_BetaV[1] = 0.;
}
}
// methods for axial and scalar FF
Complex VA_0_PPP::KS::BW_V( int a, double s )
{
if( !m_G123 && a==1 ) return Complex(0.,0.); // BW_V23 = 0 if G123=0
return m_V[a].BreitWigner(s);
}
Complex VA_0_PPP::KS::BW_VV( int a, double s )
{
if( !m_G123 && a==1 ) return Complex(0.,0.); // BW_V23' = 0 if G123=0
return m_VV[a].BreitWigner(s);
}
Complex VA_0_PPP::KS::BW_A( double s )
{
if (!IsZero(m_alpha)) {
return( ( m_A.BreitWigner(s)+m_alpha*m_AA.BreitWigner(s) ) / (1.+m_alpha) );
}
return m_A.BreitWigner(s);
}
Complex VA_0_PPP::KS::Tvector1( int a, int b, double x )
{
Complex ret =
BW_V(a,x) *( 1.-(Mass2(a)-Mass2(b))/(3.*m_V[a].Mass2()) )
+ m_Beta[a]*( BW_VV(a,x)*( 1.-(Mass2(a)-Mass2(b))/(3.*m_VV[a].Mass2()) ) );
return ret/(1.+m_Beta[a]);
}
Complex VA_0_PPP::KS::Tvector2( int a, int b, double x )
{
Complex ret = BW_V(a,x)/m_V[a].Mass2() + m_Beta[a]*BW_VV(a,x)/m_V[a].Mass2();
return ret*( 2.*(Mass2(a)-Mass2(b)) )/( 3.*(1.+m_Beta[a]) );
}
Complex VA_0_PPP::KS::TSvector( int a, int b, int c, double Q2, double s, double t )
{
Complex ret =
BW_V(a,s) * ( m_ms123*( Q2-2.*t-s+2.*Mass2(a)+Mass2(c) )
- (Mass2(a)-Mass2(b))/m_V[a].Mass2()*( m_ms123*(Q2+s-Mass2(c))-Q2*(s-m_V[a].Mass2()) ) )
+ m_Beta[a]*BW_VV(a,s) * ( m_ms123*( Q2-2.*t-s+2.*Mass2(a)+Mass2(c) )
- (Mass2(a)-Mass2(b))/m_VV[a].Mass2()*( m_ms123*(Q2+s-Mass2(c))-Q2*(s-m_VV[a].Mass2()) ) );
return ret/(1.+m_Beta[a]);
}
Complex VA_0_PPP::KS::Tgen( int a, int b, int c, double s, double t)
{
return( Tvector1(a-1,c-1,s) + Tvector2(b-1,c-1,t) );
}
// methods for vector FF
Complex VA_0_PPP::KS::T_V( double q2 )
{
Complex ret(0.,0.);
ret = m_AnoV.BreitWigner(q2)
+ m_BetaV[0] * m_AnoVV.BreitWigner(q2)
+ m_BetaV[1] * m_AnoVVV.BreitWigner(q2);
ret /= 1.+m_BetaV[0]+m_BetaV[1];
return ret;
}
Complex VA_0_PPP::KS::T_V13V23( double s, double t )
{
Complex ret(0.,0.);
int n_rho (1-1), n_k (2-1); // what is the rho and K* resonance
double x (s), y (t);
if( m_mode==2010 ) {
n_rho = 2-1; n_k = 1-1; // swap them in Kpipi channel
x=t; y=s;
}
ret = ( BW_V(n_rho,x) + m_Beta[n_rho]*BW_VV(n_rho,x) )/( 1.+m_Beta[n_rho] );
ret += m_AlphaV * BW_V(n_k,y);
ret /= 1.+m_AlphaV;
return ret;
}
// handling of form factors
Complex VA_0_PPP::KS::FormFactor( int j, double Q2, double s, double t )
{
Complex FF(0.,0.);
switch( j ) {
case 1 : { FF = BW_A(Q2)*Tgen(1,2,3,s,t); // axial
break; }
case 2 : { FF = BW_A(Q2)*Tgen(2,1,3,t,s); // axial
break; }
case 3 : { FF = m_X123 + BW_A(Q2)*(Q2-m_A.Mass2())/(m_A.Mass2()*Q2) // scalar
*( TSvector(1-1,3-1,2-1,Q2,s,t) + TSvector(2-1,3-1,1-1,Q2,t,s) );
FF /= 2.*(Q2-m_ms123);
break; }
case 4 : { FF = 3./(8.*sqr(M_PI)*m_fpi2)*T_V(Q2)*T_V13V23(s,t); // vector
break; }
}
return FF;
}
// Parameterisation
// FINKEMEIER, MIRKES hep-ph/9503474
VA_0_PPP::KS95::KS95(int mode, int kaon_mode, GeneralModel _md, double * _ms)
: FF_Base(mode,kaon_mode,_md,_ms)
{
kf_code resV, resVV, resVVV;
int running = int( _md("RUNNING_WIDTH", 3 ) ); // running width
double MV, MVV, MVVV, GV, GVV, GVVV;
if( m_deltas==0 ) {
resV = kf_rho_770_plus;
resVV = kf_rho_1450_plus;
resVVV = kf_rho_1700_plus;
MV = _md("q2_Mass_"+Flavour(resV).IDName(), 0.7769 ); // mass V
MVV = _md("q2_Mass_"+Flavour(resVV).IDName(), 1.363 ); // mass V'
MVVV = _md("q2_Mass_"+Flavour(resVVV).IDName(), 1.700 ); // mass V"
GV = _md("q2_Width_"+Flavour(resV).IDName(), 0.1490 ); // width V
GVV = _md("q2_Width_"+Flavour(resVV).IDName(), 0.310 ); // width V'
GVVV = _md("q2_Width_"+Flavour(resVVV).IDName(), 0.235 ); // width V"
m_beta_VectorFF = _md("q2_beta_"+Flavour(resVV).IDName(), -0.108 ); // beta for V'
m_gamma_VectorFF = _md("q2_gamma_"+Flavour(resVVV).IDName(), 0.020 ); // gamma for V"
}
else {
resV = kf_K_star_892_plus;
resVV = kf_K_star_1410_plus;
resVVV = kf_K_star_1680_plus;
MV = _md("q2_Mass_"+Flavour(resV).IDName(), 0.892 ); // mass V
MVV = _md("q2_Mass_"+Flavour(resVV).IDName(), 1.412 ); // mass V'
MVVV = _md("q2_Mass_"+Flavour(resVVV).IDName(), 1.714 ); // mass V"
GV = _md("q2_Width_"+Flavour(resV).IDName(), 0.051 ); // width V
GVV = _md("q2_Width_"+Flavour(resVV).IDName(), 0.227 ); // width V'
GVVV = _md("q2_Width_"+Flavour(resVVV).IDName(), 0.323 ); // width V"
m_beta_VectorFF = _md("q2_beta_"+Flavour(resVV).IDName(), -0.25 ); // beta for V'
m_gamma_VectorFF = _md("q2_gamma_"+Flavour(resVVV).IDName(), 0.038 ); // gamma for V"
}
m_V_VectorFF = ResonanceFlavour( resV, MV, GV, running&2 );
m_VV_VectorFF = ResonanceFlavour( resVV, MVV, GVV, running&2 );
m_VVV_VectorFF = ResonanceFlavour( resVVV, MVVV, GVVV, running&2 );
// omega - phi - resonances
kf_code resO = kf_omega_782;
kf_code resP = kf_phi_1020;
double MO = _md("G3_Mass_"+Flavour(resO).IDName(), 0.782 ); // mass omega
double MP = _md("G3_Mass_"+Flavour(resP).IDName(), 1.020 ); // mass phi
double GO = _md("G3_Width_"+Flavour(resO).IDName(), 0.00843 ); // mass omega
double GP = _md("G3_Width_"+Flavour(resP).IDName(), 0.00443 ); // mass phi
m_Omega_VectorFF = ResonanceFlavour( resO, MO, GO, 0 );
m_Phi_VectorFF = ResonanceFlavour( resP, MP, GP, 0 );
m_eps_VectorFF = _md("G3_epsilon_"+Flavour(resP).IDName(), 0.05 ); // epsilon for phi
}
Complex VA_0_PPP::KS95::Axial(double s)
{
if (!IsZero(m_alpha)) {
return( ( m_A.BreitWigner(s)+m_alpha*m_AA.BreitWigner(s) ) / (1.+m_alpha) );
}
return m_A.BreitWigner(s);
}
Complex VA_0_PPP::KS95::TVector(int i, double s)
{
return (m_V[i].BreitWigner(s)+m_Beta[i]*m_VV[i].BreitWigner(s))/(1.+m_Beta[i]);
}
Complex VA_0_PPP::KS95::Vector(int i, double s, double u)
{
switch(m_mode) {
case 1101 : { // pi- K0b pi0
if(m_kaon_mode==3) { // K0 K0
switch(i) {
case 0 : return TVector(0,s);
case 1 : return TVector(1,s);
}
}
else if(m_kaon_mode==20) { // KS KS
switch(i) {
case 0 : return 0.5*TVector(1,u);
case 1 : return -0.5*(TVector(1,s)+TVector(1,u));
}
}
else if(m_kaon_mode==200) { // KL KL
switch(i) {
case 0 : return 0.5*TVector(1,u);
case 1 : return -0.5*(TVector(1,s)+TVector(1,u));
}
switch(i) {
case 0 : return -0.5*TVector(1,u);
case 1 : return 0.5*(TVector(1,s)+TVector(1,u));
}
}
else if(m_kaon_mode==110) { // KS KL
switch(i) {
case 0 : return 0.5*(2.*TVector(0,s)+TVector(1,u));
case 1 : return 0.5*(TVector(1,s)-TVector(1,u));
}
}
}
case 111 : { // K- pi0 K0
switch(i) {
case 0 : return 2./3.*TVector(0,s) + 1./3.*TVector(1,u);
case 1 : return 1./3.*TVector(1,s) - 1./3.*TVector(1,u);
}
}
}
return TVector(i,s);
}
Complex VA_0_PPP::KS95::TOmega(double s)
{
return (m_Omega_VectorFF.BreitWigner(s) + m_eps_VectorFF*m_Phi_VectorFF.BreitWigner(s)) / (1.+m_eps_VectorFF);
}
Complex VA_0_PPP::KS95::G3(double q2, double s, double t, double u )
{
Complex vector1
= m_V_VectorFF.BreitWigner(q2)
+ m_beta_VectorFF*m_VV_VectorFF.BreitWigner(q2)
+ m_gamma_VectorFF*m_VVV_VectorFF.BreitWigner(q2);
vector1 /= 1.+m_beta_VectorFF+m_gamma_VectorFF;
Complex vector(0.,0.);
switch(m_mode) {
case 3000 :
case 1200 : return Complex(0.,0.);
case 210 : // pi0 pi0 K-
{
vector = TVector(2-1,t) - TVector(1-1,s);
}
case 2010 : // K- pi- pi+
{
vector = TVector(2-1,t) + TVector(1-1,s);
}
case 1101 : // pi- K0b pi0
{
vector = 2.*TVector(1-1,s) + TVector(2-1,t) + TVector(2-1,u);
vector /= 3.;
}
case 1020 : // K- pi- K+
{
vector = sqrt(2.)*TOmega(s) + TVector(2-1,t);
vector *= (sqrt(2.)-1.);
}
case 1002 : // K0 pi- K0b
{
if(m_kaon_mode==2) { // K0 K0
vector = sqrt(2.)*TOmega(s) + TVector(2-1,t);
vector *= (sqrt(2.)-1.);
}
else if(m_kaon_mode==20) { // KS KS
vector = TVector(2-1,t)-TVector(2-1,u);
vector *= -0.5*(sqrt(2.)-1.);
}
else if(m_kaon_mode==200) { // KL KL
vector = TVector(2-1,t)-TVector(2-1,u);
vector *= -0.5*(sqrt(2.)-1.);
}
else if(m_kaon_mode==110) { // KS KL
vector = 2.*sqrt(2.)*TOmega(s)+TVector(2-1,t)+TVector(2-1,u);
vector *= 0.5*(sqrt(2.)-1.);
}
}
case 111 : // K- pi0 K0
{
vector = TVector(1-1,s) - TVector(2-1,t);
vector *= (sqrt(2.)-1.);
}
}
return vector1*vector;
}
// handling of form factors
Complex VA_0_PPP::KS95::FormFactor( int j, double Q2, double s, double t )
{
double u = Q2 + Mass2(1-1) + Mass2(2-1) + Mass2(3-1) - s - t;
switch( j ) {
case 1 : // axial
return Axial(Q2)*Vector(1-1,s,u);
case 2 : // axial
return Axial(Q2)*Vector(2-1,t,u);
case 3 : // scalar
break;
case 4 : // vector
return 6./(8.*sqr(M_PI)*m_fpi2)*G3(Q2,s,t,u);
break;
}
Complex FF(0.,0.);
return FF;
}
DEFINE_CURRENT_GETTER(HADRONS::VA_0_PPP,"VA_0_PPP")
void ATOOLS::Getter<HADRONS::Current_Base,HADRONS::ME_Parameters,HADRONS::VA_0_PPP>::
PrintInfo(std::ostream &st,const size_t width) const {
st<<"Example: $ 0 \\rightarrow \\pi \\pi K $ \n\n"
<<"Order: special, see tau decays. \n\n"
<<"Available form factors: \n "
<<" \\begin{itemize} \n"
<<" \\item {\\tt FORM\\_FACTOR = 1 :} Kuehn-Santamaria \n"
<<" \\item {\\tt FORM\\_FACTOR = 2 :} Resonance Chiral Theory \n"
<<" \\end{itemize} \n"
<<"Reference: https://sherpa.hepforge.org/olddokuwiki/data/media/publications/theses/diplom\\__laubrich.pdf \n"
<<std::endl;
}