#include "HADRONS++/Current_Library/VA_0_PiPiPiPi3Charged.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/My_File.H" #include "ATOOLS/Org/My_MPI.H" #include "ATOOLS/Org/Shell_Tools.H" using namespace HADRONS; using namespace ATOOLS; using namespace std; void VA_0_PiPiPiPi3Charged::SetModelParameters( GeneralModel _md ) { double Vud = _md("Vud", Tools::Vud); m_global = Vud; p_lorenz = new Novo(_md); } void VA_0_PiPiPiPi3Charged::LorenzBase::SetPrivates( const ATOOLS::Vec4D * _p,int * i) { p_i = i; p_p = _p; m_q = p_p[p_i[0]]+p_p[p_i[1]]+p_p[p_i[2]]+p_p[p_i[3]]; // = q m_q2 = m_q.Abs2(); } // Novosibirsk parameterisation // see hep-ph/0201149 for details and // hep-ph/0312240 for errata VA_0_PiPiPiPi3Charged::Novo::Novo( GeneralModel _md ) : LorenzBase() { m_mpi2 = sqr( Flavour(kf_pi_plus).HadMass() ); m_rho = SimpleResonanceFlavour( string("kf_rho_770_plus"), _md("Mass_rho(770)+", _md("Mass_rho(770)", 0.7761 )), _md("Width_rho(770)+", _md("Width_rho(770)", 0.1445 )) ); m_sigma = SimpleResonanceFlavour( string("kf_f_0_600"), _md("Mass_sigma", 0.800), _md("Width_sigma", 0.800) ); m_omega = SimpleResonanceFlavour( string("kf_omega_782"), _md("Mass_omega(782)", 0.782), _md("Width_omega(782)", 0.00841) ); m_a1 = SimpleResonanceFlavour( string("a(1)(1260)+"), _md("Mass_a(1)(1260)+", 1.230 ), _md("Width_a(1)(1260)+", 0.450 ) ); m_Lambda2 = _md("Lambda2", 1.2); double z_abs = _md("z_abs", 1.3998721 ); double z_phase = _md("z_phase", 0.43585036 ); m_z = Complex( z_abs*cos(z_phase), z_abs*sin(z_phase) ); // G(q2) function (read as histogram) // if file does not exist if (!FileExists(rpa->gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/G_pi-pi+pi-pi0.dat")) { THROW(fatal_error, rpa->gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/G_pi-pi+pi-pi0.dat not found."); } else { // read table and create histogram msg_Tracking()<<"HADRONS::VA_0_PiPiPiPi3Charged::Novo::Novo(...) \n" <<" Read G_{pi-pi+pi-pi0}(q2)."<gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/G_pi-pi+pi-pi0.dat"); infile.Open(); std::string found_file_name = infile.FileContent(); p_G = new Histogram("PhaseSpaceFunctions/G_pi-pi+pi-pi0.dat",0,found_file_name); } // Gomega(q2) function (read as histogram) // if file does not exist if (!FileExists(rpa->gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/Gomega_pi-pi+pi-pi0.dat")) { THROW(fatal_error, rpa->gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/Gomega_pi-pi+pi-pi0.dat not found."); } else { // read table and create histogram msg_Tracking()<<"HADRONS::VA_0_PiPiPiPi3Charged::Novo::Novo(...) \n" <<" Read Gomega_{pi+pi0pi-pi-}(q2)."<gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/Gomega_pi-pi+pi-pi0.dat"); infile.Open(); std::string found_file_name = infile.FileContent(); p_Go = new Histogram("PhaseSpaceFunctions/Gomega_pi-pi+pi-pi0.dat",0,found_file_name); } // RunningWidth_a1(q2) function (read as histogram) // if file does not exist if (!FileExists(rpa->gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/RunningWidth_a1_4pi-channel.dat")) { THROW(fatal_error, rpa->gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/RunningWidth_a1_4pi-channel.dat not found"); } else { // read table and create histogram msg_Tracking()<<"HADRONS::VA_0_PiPiPiPi3Charged::Novo::Novo(...) \n" <<" Read a1's running width (q2)."<gen.Variable("SHERPA_SHARE_PATH")+"/PhaseSpaceFunctions/RunningWidth_a1_4pi-channel.dat"); infile.Open(); std::string found_file_name = infile.FileContent(); p_a1width = new Histogram("PhaseSpaceFunctions/RunningWidth_a1_4pi-channel.dat",0,found_file_name); } // scale histograms if MASS_SCALE is set in decay channel file // if not set -> scaled by tau_mass/tau_mass = 1 double tau_mass = Flavour(kf_tau).HadMass(); double mass_scale = double(_md("MASS_SCALE", tau_mass)) / tau_mass; p_G->ScaleHistogramWidth(mass_scale); p_Go->ScaleHistogramWidth(mass_scale); p_a1width->ScaleHistogramWidth(mass_scale * mass_scale); } // global form factor double VA_0_PiPiPiPi3Charged::Novo::G( double q2 ) { if(sqrt(q2)<0.6) return 0.; double val(0.); p_G->Extrapolate( sqrt(q2)+p_G->BinSize(), &val, 0 ); // shift +BinSize to get right value return val*76.565*sqrt(0.7179*sqrt(q2)-0.27505)/(sqr(m_rho.Mass2())*sqrt(q2)); // note: magic correction from TAUOLA (hep-ph/0312240) } double VA_0_PiPiPiPi3Charged::Novo::Go( double q2 ) { if(sqrt(q2)<0.6) return 0.; double val(0.); p_Go->Extrapolate( sqrt(q2)+p_Go->BinSize(), &val, 0 ); // shift +BinSize to get right value return val*886.84*sqrt(0.70983*sqrt(q2)-0.26689)/(sqr(m_rho.Mass2())*sqrt(q2)); // note: magic correction from TAUOLA (hep-ph/0312240) } // a1 form factor double VA_0_PiPiPiPi3Charged::Novo::F2_a1( double s ) { return sqr( (1.+m_a1.Mass2()/m_Lambda2)/(1.+s/m_Lambda2) ); } // a1 propagator Complex VA_0_PiPiPiPi3Charged::Novo::Da1( double s ) { double val(0); p_a1width->Extrapolate( s+p_a1width->BinSize(), &val, 0 ); // shift +BinSize to get right value return Complex( s/m_a1.Mass2()-1., sqrt(s)*val/m_a1.Mass2() ); } // rho h function double VA_0_PiPiPiPi3Charged::Novo::hrho( double s ) { double epsilon = 1.e-8; if( s>4*m_mpi2 ) { double root = sqrt(1.-4.*m_mpi2/s); return root*log((1.+root)/(1.-root))*(s-4.*m_mpi2)/M_PI; } else if( s>epsilon ) return 0.; else return -8.*m_mpi2/M_PI; } // rho propagator Complex VA_0_PiPiPiPi3Charged::Novo::Drho( double s ) { double real, imag; double g2 = pow( m_rho.Mass2()-4.*m_mpi2, 1.5 )/sqrt(m_rho.Mass2()); double root = sqrt(1.-4.*m_mpi2/m_rho.Mass2()); double dhrhods = root/M_PI * ( root + (1.+2.*m_mpi2/m_rho.Mass2())*log((1.+root)/(1.-root)) ); double dm = ( hrho(s) - hrho(m_rho.Mass2()) - (s-m_rho.Mass2())*dhrhods ) / g2; real = s-m_rho.Mass2() - m_rho.Mass()*m_rho.Width()*dm; if( s-4*m_mpi2 <= 0. ) imag = 0.; else { double g1 = pow( s-4.*m_mpi2, 1.5 )/sqrt(s); imag = m_rho.Mass()*m_rho.Width()*g1/g2; } double dm0 = ( hrho(0.) - hrho(m_rho.Mass2()) + m_rho.Mass2()*dhrhods ) / g2; double normalisation = m_rho.Mass2() + m_rho.Mass()*m_rho.Width()*dm0; return Complex(real/normalisation,imag/normalisation); } // sigma propagator Complex VA_0_PiPiPiPi3Charged::Novo::Dsigma( double s ) { double real, imag; real = s-m_sigma.Mass2(); if( s-4*m_mpi2 <= 0. ) imag = 0.; else { double g1 = sqrt( 1.-4.*m_mpi2/s ); double g2 = sqrt( 1.-4.*m_mpi2/m_sigma.Mass2() ); imag = m_sigma.Mass()*m_sigma.Width()*g1/g2; } double normalisation = m_sigma.Mass2(); return Complex(real/normalisation,imag/normalisation); } // omega form factor double VA_0_PiPiPiPi3Charged::Novo::F2_o( double s ) { return 1.; } // omega propagator Complex VA_0_PiPiPiPi3Charged::Novo::Do( double s ) { double real, imag; real = s-m_omega.Mass2(); double q (sqrt(s)), g; if( q<1.0 ) { q -= m_omega.Mass(); g = 1+q*(17.560+q*(141.110+q*(894.884+q*(4977.35+q*(7610.66-42524.4*q))))); if( g<0. ) g = 0.; } else { g = -1333.26+q*(4860.19+q*(-6000.81+2504.97*q)); } imag = m_omega.Mass()*m_omega.Width()*g; double normalisation = m_omega.Mass2(); return Complex(real/normalisation,imag/normalisation); } // t1 function Vec4C VA_0_PiPiPiPi3Charged::Novo::t1( int a, int b, int c, int d ) { // subtract 1 for C++ convention a--; b--; c--; d--; Vec4D qa1 = m_q - p_p[p_i[a]]; // virtuality of a1 double qa12 = qa1.Abs2(); Complex global = Complex(-1.*F2_a1(qa12), 0.) / Da1(qa12) / Drho((p_p[p_i[c]]+p_p[p_i[d]]).Abs2()); double dot0 = qa1*m_q; double dot1 = qa1*p_p[p_i[c]]; double dot2 = qa1*p_p[p_i[d]]; double dot3 = (m_q*p_p[p_i[d]]) * (p_p[p_i[a]]*p_p[p_i[c]]); double dot4 = (m_q*p_p[p_i[c]]) * (p_p[p_i[d]]*p_p[p_i[a]]); return global * ( dot0*(p_p[p_i[d]]*dot1 - p_p[p_i[c]]*dot2) + qa1*(dot3-dot4) ); } // t2 function Vec4C VA_0_PiPiPiPi3Charged::Novo::t2( int a, int b, int c, int d ) { // subtract 1 for C++ convention a--; b--; c--; d--; Vec4D qa1 = m_q - p_p[p_i[a]]; // virtuality of a1 double qa12 = qa1.Abs2(); Complex global = m_z * F2_a1(qa12) / Da1(qa12) / Dsigma((p_p[p_i[c]]+p_p[p_i[d]]).Abs2()); double dot0 = m_q*qa1; double dot1 = m_q*p_p[p_i[b]]; return global * ( p_p[p_i[b]]*dot0 - qa1*dot1 ) * qa12; } // t3 function Vec4C VA_0_PiPiPiPi3Charged::Novo::t3( int a, int b, int c, int d ) { // subtract 1 for C++ convention a--; b--; c--; d--; Vec4D qo = m_q - p_p[p_i[a]]; // virtuality of omega double qo2 = qo.Abs2(); Complex global = F2_o(qo2) / Do(qo2) / Drho((p_p[p_i[c]]+p_p[p_i[d]]).Abs2()); double dotq2 = m_q*p_p[p_i[b]]; double dotq3 = m_q*p_p[p_i[c]]; double dotq4 = m_q*p_p[p_i[d]]; double dot12 = p_p[p_i[a]]*p_p[p_i[b]]; double dot13 = p_p[p_i[a]]*p_p[p_i[c]]; double dot14 = p_p[p_i[a]]*p_p[p_i[d]]; return global * ( p_p[p_i[b]]*(dotq3*dot14-dotq4*dot13) - dotq2*(p_p[p_i[c]]*dot14-p_p[p_i[d]]*dot13) + dot12*(p_p[p_i[c]]*dotq4-p_p[p_i[d]]*dotq3) ); } Vec4C VA_0_PiPiPiPi3Charged::Novo::operator()() { // a1 -> rho pi current Vec4C J_a1rpi = G(m_q2) * ( t1(1,2,3,4) + t1(3,2,1,4) + t1(1,3,2,4) + t1(3,1,2,4) + t1(4,3,1,2) + t1(4,1,3,2) ); // a1 -> sigma pi current Vec4C J_a1spi = G(m_q2) * ( t2(4,3,1,2) + t2(4,1,3,2) - t2(1,4,3,2) - t2(3,4,1,2) ); // omega -> rho pi current Vec4C J_omrpi = Go(m_q2) * ( t3(1,2,3,4) + t3(3,2,1,4) - t3(1,3,2,4) - t3(3,1,2,4) - t3(1,4,3,2) - t3(3,4,1,2) ); // sum (factor 2! due to 2 identical particles) return (J_a1rpi+J_a1spi+J_omrpi)/sqrt(2.); } // CLEO parameterisation // see hep-ex/9908024 and CERN-TH.6793/93 for details VA_0_PiPiPiPi3Charged::KS::KS( GeneralModel _md ) : LorenzBase() { for (int i=1; i<=3; i++ ) { m_r[i] = m_q-p_p[p_i[i]]; m_s[i] = (p_p[p_i[0]]+p_p[p_i[i]]).Abs2(); } // redefinition of variables m_s[0] = (p_p[p_i[1]]+p_p[p_i[3]]).Abs2(); // = t4 // unused variables yet m_r[0] = ATOOLS::Vec4D(0.,0.,0.,0.); m_fpi2 = sqr(_md("fpi",0.0924))*2.; // redefine fpi m_grop = _md("grop", 12.924); m_Go3p = _md("Go3p", 1476.); m_mpi2 = sqr( Flavour(kf_pi_plus).HadMass() ); m_mpi02 = sqr( Flavour(kf_pi).HadMass() ); double MR = _md("Mass_rho(770)+", Flavour(kf_rho_770_plus).HadMass() ); double MRR = _md("Mass_rho(1450)+", Flavour(kf_rho_1450_plus).HadMass() ); double MRRR = _md("Mass_rho(1700)+", Flavour(kf_rho_1700_plus).HadMass() ); double MO = _md("Mass_omega(782)", Flavour(kf_omega_782).HadMass() ); double MF = _md("Mass_f(0)(980)", Flavour(kf_f_0_980).HadMass() ); double MS = _md("Mass_sigma", Flavour(kf_f_0_980).HadMass() ); double MA = _md("Mass_a(1)(1260)+", Flavour(kf_a_1_1260_plus).HadMass()); double GR = _md("Width_rho(770)+", Flavour(kf_rho_770_plus).Width() ); double GRR = _md("Width_rho(1450)+", Flavour(kf_rho_1450_plus).Width() ); double GRRR = _md("Width_rho(1700)+", Flavour(kf_rho_1700_plus).Width() ); double GO = _md("Width_omega(782)", Flavour(kf_omega_782).Width() ); double GF = _md("Width_f(0)(980)", Flavour(kf_f_0_980).Width() ); double GS = _md("Width_sigma", Flavour(kf_f_0_980).Width() ); double GA = _md("Width_a(1)(1260)+", Flavour(kf_a_1_1260_plus).Width()); m_Rho = ResonanceFlavour( kf_rho_770_plus, MR, GR, 1 ); m_RR = ResonanceFlavour( kf_rho_1450_plus, MRR, GRR, 1 ); m_RRR = ResonanceFlavour( kf_rho_1700_plus, MRRR, GRRR, 1 ); m_O = ResonanceFlavour( kf_omega_782, MO, GO, 0 ); m_F = ResonanceFlavour( kf_f_0_980, MF, GF, 1 ); m_S = ResonanceFlavour( kf_f_0_980, MS, GS, 1 ); m_A = ResonanceFlavour( kf_a_1_1260_plus, MA, GA, 1 ); m_beta = _md("beta", 0.); m_gamma = _md("gamma", 0.); m_sigma = _md("sigma", 0.); m_A.InitialiseThreeBodyResonance(m_Rho, m_RR, m_beta); m_Frho = _md("frho", 0.266)*m_Rho.Mass2(); // m_R[0]=m_R[1] = 0.; // m_R[2] = -2.; // m_R[3]=m_R[4] = 1.; m_R[0] = 0.; m_R[1] = -2.; m_R[2]=m_R[3] = 1.; Complex sum (0.,0.); char helps[20]; double absol, phase; for (int i=0; i<4; i++) { snprintf( helps,20,"alpha_%i",i ); absol = _md(helps+string("_abs"), 0. ); phase = _md(helps+string("_phase"), 0. ); m_Alpha[i] = Complex( absol*cos(phase), absol*sin(phase) ); sum += m_Alpha[i]; } m_SumAlpha = sum; for (int i=0; i<4; i++) { snprintf( helps, 20, "beta_omega_pi_%i", i ); absol = _md( helps+string("_abs"), 0. ); phase = _md( helps+string("_phase"), 0. ); m_Beta_opi[i] = Complex( absol*cos(phase), absol*sin(phase) ); snprintf( helps, 20, "beta_a1_pi_%i", i ); absol = _md( helps+string("_abs"), 0. ); phase = _md( helps+string("_phase"), 0. ); m_Beta_api[i] = Complex( absol*cos(phase), absol*sin(phase) ); snprintf( helps, 20, "beta_sigma_rho_%i", i ); absol = _md( helps+string("_abs"), 0. ); phase = _md( helps+string("_phase"), 0. ); m_Beta_srh[i] = Complex( absol*cos(phase), absol*sin(phase) ); snprintf( helps, 20, "beta_f0_rho_%i", i ); absol = _md( helps+string("_abs"), 0. ); phase = _md( helps+string("_phase"), 0. ); m_Beta_frh[i] = Complex( absol*cos(phase), absol*sin(phase) ); } } Complex VA_0_PiPiPiPi3Charged::KS::Fk( double x, Complex * _beta ) { Complex BW_R = m_Rho.BreitWigner(x); Complex BW_RR = m_RR.BreitWigner(x); Complex BW_RRR = m_RRR.BreitWigner(x); return( (_beta[0] + _beta[1]*BW_R + _beta[2]*BW_RR + _beta[3]*BW_RRR)/(_beta[0]+_beta[1]+_beta[2]+_beta[3]) ); } Complex VA_0_PiPiPiPi3Charged::KS::Trho( double x ) { Complex BW_R = m_Rho.BreitWigner(x); Complex BW_RR = m_RR.BreitWigner(x); Complex BW_RRR = m_RRR.BreitWigner(x); return( (BW_R + m_beta*BW_RR + m_gamma*BW_RRR)/(1.+m_beta+m_gamma) ); } Complex VA_0_PiPiPiPi3Charged::KS::TTrho( double x ) { Complex BW_R = m_Rho.BreitWignerAlt(x); Complex BW_RR = m_RR.BreitWignerAlt(x); return( BW_R + m_sigma*BW_RR ); } double VA_0_PiPiPiPi3Charged::KS::Dots( int k, int l ) // k=2,3; l=0,1,2,3 (l!=k) // special dot products used for anomalous part of J_omega { int pre = l-1; // predecessor of l if (pre==-1) pre = 3; if (pre==k) pre -= 1; int suc = l+1; // successor of l if (suc==k) suc += 1; if (suc==4) suc = 0; return( (m_r[k]*p_p[p_i[pre]])*(p_p[p_i[suc]]*p_p[p_i[k]]) - (m_r[k]*p_p[p_i[suc]])*(p_p[p_i[pre]]*p_p[p_i[k]]) ); } Vec4C VA_0_PiPiPiPi3Charged::KS::OmegaPi() { // chiral part Vec4C J_chi( 0.0,0.0,0.0,0.0 ); Vec4C sum1( 0.0,0.0,0.0,0.0 ); Vec4C sum2( 0.0,0.0,0.0,0.0 ); for (int k=1; k<=3; k++) { sum2 = Vec4C( 0.0,0.0,0.0,0.0 ); for (int l=1; l<=3; l++) if (l!=k) { sum2 += Vec4C(m_q-2.*p_p[p_i[l]]) * ((m_r[l]*(p_p[p_i[k]]-p_p[p_i[0]]))/m_r[l].Abs2()); } sum1 += m_R[k]*Trho(m_s[k]) * (Vec4C(p_p[p_i[k]])-Vec4C(p_p[p_i[0]])-sum2); } J_chi = 2.*sqrt(3.)/m_fpi2*Trho(m_q2) * sum1; // anomalous part sum1 = Vec4C( 0.0,0.0,0.0,0.0 ); for( int k=2; k<=3; k++ ) { sum2 = Vec4C( 0.0,0.0,0.0,0.0 ); for (int l=0; l<=3; l++) if (l!=k) { sum2 += Vec4C(p_p[p_i[l]]) * Dots(k,l); } sum1 += m_O.BreitWignerAlt(m_r[k].Abs2())*sum2; } Vec4C J_a = m_Go3p*m_Frho*m_grop*TTrho(m_q2) * sum1; // total return( (J_chi + J_a)*Fk(m_q2,m_Beta_opi) ); } Vec4C VA_0_PiPiPiPi3Charged::KS::AonePi() { Vec4C term1(0.,0.,0.,0.), term2(0.,0.,0.,0.); Complex A, B, C; Vec4D P, Q, R; // 1st term P = p_p[p_i[0]]-p_p[p_i[2]]; Q = p_p[p_i[0]]-p_p[p_i[3]]; R = m_r[1]; A = m_Rho.BreitWigner(m_s[2]); B = m_Rho.BreitWigner(m_s[3]); C = m_A.BreitWigner(R.Abs2()); term1 = A*Vec4C(p_p[p_i[0]]-p_p[p_i[2]]) + B*Vec4C(p_p[p_i[0]]-p_p[p_i[3]]) - Vec4C(m_q-p_p[p_i[1]])*(( A*(R*P) + B*(R*Q) )/R.Abs2()) - Vec4C(m_q)*(( A*(m_q*P) + B*(m_q*Q) - (m_q*Q)*(A*(R*P)+B*(R*Q))/R.Abs2() )/m_q2); term1 *= C; // 2nd and 3rd term int ind; Vec4C help; for (int k=2; k<=3; k++) { ind = (k==2)? 3 : 2; P = p_p[p_i[0]]-p_p[p_i[1]]; Q = p_p[p_i[ind]]-p_p[p_i[1]]; R = m_r[k]; A = m_Rho.BreitWigner(m_s[1]); double t = (ind==2)? (p_p[p_i[1]]+p_p[p_i[2]]).Abs2() : m_s[0]; // s[0] = t[3] B = m_Rho.BreitWigner(t); C = m_A.BreitWigner(R.Abs2()); help = A*Vec4C(p_p[p_i[0]]-p_p[p_i[1]]) + B*Vec4C(p_p[p_i[1]]-p_p[p_i[ind]]) - Vec4C(m_q-p_p[p_i[k]])*(( A*(R*P) + B*(R*Q) )/R.Abs2()) - Vec4C(m_q)*(( A*(m_q*P) + B*(m_q*Q) - (m_q*Q)*(A*(R*P)+B*(R*Q))/R.Abs2() )/m_q2); term2 += C*help; } // total return (term1-term2)*Fk(m_q2,m_Beta_api); } Vec4C VA_0_PiPiPiPi3Charged::KS::SigmaRho() { int ind; Vec4C term( 0.0,0.0,0.0,0.0 ); Complex BW_S, BW_R; for (int k=2; k<=3; k++) { ind = (k==2)? 3 : 2; BW_S = m_S.BreitWigner(m_s[k]); double t = (ind==2)? (p_p[p_i[1]]+p_p[p_i[2]]).Abs2() : m_s[0]; // s[0] = t[3] BW_R = m_Rho.BreitWigner(t); term += BW_S*BW_R * Vec4C( p_p[p_i[1]] - p_p[p_i[ind]] + (m_q*(p_p[p_i[ind]]-p_p[p_i[1]]))*m_q ); } return term*Fk(m_q2,m_Beta_srh); } Vec4C VA_0_PiPiPiPi3Charged::KS::FzeroRho() { int ind; Vec4C term( 0.0,0.0,0.0,0.0 ); Complex BW_S, BW_R; for (int k=2; k<=3; k++) { ind = (k==2)? 3 : 2; BW_S = m_F.BreitWigner(m_s[k]); double t = (ind==2)? (p_p[p_i[1]]+p_p[p_i[2]]).Abs2() : m_s[0]; // s[0] = t[3] BW_R = m_Rho.BreitWigner(t); term += BW_S*BW_R * Vec4C( p_p[p_i[1]] - p_p[p_i[ind]] + (m_q*(p_p[p_i[ind]]-p_p[p_i[1]]))*m_q ); } return term*Fk(m_q2,m_Beta_frh); } Vec4C VA_0_PiPiPiPi3Charged::KS::operator()() { Vec4C help(0.,0.,0.,0.); help += m_Alpha[0] * OmegaPi(); help += m_Alpha[1] * AonePi(); help += m_Alpha[2] * SigmaRho(); help += m_Alpha[3] * FzeroRho(); return help/m_SumAlpha; } // General framework void VA_0_PiPiPiPi3Charged::Calc(const ATOOLS::Vec4D_Vector& moms, bool m_anti) { Vec4C help(0.,0.,0.,0.); // Lorentz structure p_lorenz->SetPrivates(&moms.front(), &p_i.front()); help = (*p_lorenz)(); Insert( help*m_global , 0); } DEFINE_CURRENT_GETTER(HADRONS::VA_0_PiPiPiPi3Charged,"VA_0_PiPiPiPi3Charged") void ATOOLS::Getter:: PrintInfo(std::ostream &st,const size_t width) const { st<<"Example: $ 0 \\rightarrow \\pi^+ \\pi^+ \\pi^- \\pi^0 $ \n\n" <<"Order: 0 = $\\pi^\\mp$, 1 = $\\pi^0$, 2 = $\\pi^\\pm$, 3 = $\\pi^\\pm$ " <<"Available form factors: \n" <<" \\begin{itemize} \n" <<" \\item Kuehn-Santamaria (default) \n" <<" \\end{itemize} \n" <<"Reference: https://sherpa.hepforge.org/olddokuwiki/data/media/publications/theses/diplom\\__laubrich.pdf \n" <