#include "PHOTONS++/MEs/Scalar_To_Vector_Lepton_Neutrino.H" #include "ATOOLS/Math/Poincare.H" #include "ATOOLS/Phys/Flavour.H" #include "ATOOLS/Phys/Particle.H" #include "METOOLS/Main/XYZFuncs.H" #include "METOOLS/Main/Polarization_Tools.H" #include "METOOLS/Loops/Master_Integrals.H" #include "METOOLS/Loops/PV_Integrals.H" using namespace PHOTONS; using namespace ATOOLS; using namespace METOOLS; using namespace std; Scalar_To_Vector_Lepton_Neutrino::Scalar_To_Vector_Lepton_Neutrino (const Particle_Vector_Vector& pvv) : PHOTONS_ME_Base(pvv), Dipole_FF(pvv) { m_name = "Scalar_To_Vector_Lepton_Neutrino"; m_flavs[0] = pvv[1][0]->Flav(); // neutral IS scalar m_masses[0] = pvv[1][0]->FinalMass(); m_flavs[3] = pvv[3][0]->Flav(); // neutrino (calcs with m_nu = 0) m_masses[3] = 0.; // switch ordering if necessary (FS charged scalar at pos 1) m_switch = pvv[2][0]->Flav().IsLepton(); // m_switch == true if first charged final state particle is lepton if (m_switch == false) { m_flavs[1] = pvv[2][0]->Flav(); // charged FS vector m_masses[1] = pvv[2][0]->FinalMass(); m_flavs[2] = pvv[2][1]->Flav(); // lepton m_masses[2] = pvv[2][1]->FinalMass(); } else { m_flavs[1] = pvv[2][1]->Flav(); // charged FS vector m_masses[1] = pvv[2][1]->FinalMass(); m_flavs[2] = pvv[2][0]->Flav(); // lepton m_masses[2] = pvv[2][0]->FinalMass(); } for (unsigned int i=4; i<9; i++) { m_flavs[i] = Flavour(kf_photon); m_masses[i] = 0.; } m_cL = 1.; m_cR = 0.; // operate in Form-Factor-Model, point-like otherwise m_ffmodel = false; m_fpluszero = 1.; m_fplusprimezero = 1.; } Scalar_To_Vector_Lepton_Neutrino::~Scalar_To_Vector_Lepton_Neutrino() { } void Scalar_To_Vector_Lepton_Neutrino::BoostOriginalPVVToMultipoleCMS() { // m_pvv_one already in multipole CMS // m_pvv_zero in arbitrary frame -> boost m_olddipole into its CMS // and rotate m_olddipole.at(0) into -z direction Vec4D sum(0.,0.,0.,0.); for (unsigned int i=0; iMomentum(); } Vec4D p1 = m_olddipole[0]->Momentum(); p_boost = new Poincare(sum); p_boost->Boost(p1); p_rot = new Poincare(p1,Vec4D(0.,0.,0.,1.)); for (unsigned int i=0; iMomentum(); p_boost->Boost(vec); p_rot->Rotate(vec); m_olddipole[i]->SetMomentum(vec); } for (unsigned int i=0; iMomentum(); p_boost->Boost(vec); p_rot->Rotate(vec); m_oldspectator[i]->SetMomentum(vec); } } void Scalar_To_Vector_Lepton_Neutrino::FillMomentumArrays (const Particle_Vector_Vector& pvv_one) { // m_moms0 - no photon Poincare boost(m_pvv_zero[1][0]->Momentum()); Vec4D vec(0.,0.,0.,0.); // neutral IS scalar vec = m_pvv_zero[1][0]->Momentum(); boost.Boost(vec); m_moms0[0] = vec; // neutrino vec = m_pvv_zero[3][0]->Momentum(); boost.Boost(vec); m_moms0[3] = vec; if (m_switch == false) { // charged FS vector vec = m_pvv_zero[2][0]->Momentum(); boost.Boost(vec); m_moms0[1] = vec; // lepton vec = m_pvv_zero[2][1]->Momentum(); boost.Boost(vec); m_moms0[2] = vec; } else { // charged FS scalar vec = m_pvv_zero[2][1]->Momentum(); boost.Boost(vec); m_moms0[1] = vec; // lepton vec = m_pvv_zero[2][0]->Momentum(); boost.Boost(vec); m_moms0[2] = vec; } // m_moms1 - project multiphoton state onto one photon phase space // do reconstruction procedure again pretending only one photon was generated // not necessary if only one photon if (pvv_one[4].size() == 1) { m_moms1[0][0] = pvv_one[1][0]->Momentum(); m_moms1[0][3] = pvv_one[3][0]->Momentum(); if (m_switch == false) { m_moms1[0][1] = pvv_one[2][0]->Momentum(); m_moms1[0][2] = pvv_one[2][1]->Momentum(); } else { m_moms1[0][1] = pvv_one[2][1]->Momentum(); m_moms1[0][2] = pvv_one[2][0]->Momentum(); } m_moms1[0][4] = pvv_one[4][0]->Momentum(); } else { Dipole_FF::DefineDipole(); BoostOriginalPVVToMultipoleCMS(); for (unsigned int i=0; iMomentum(); m_moms1[i][2] = m_newdipole[1]->Momentum(); } else { m_moms1[i][1] = m_newdipole[1]->Momentum(); m_moms1[i][2] = m_newdipole[0]->Momentum(); } m_moms1[i][3] = m_newspectator[0]->Momentum(); m_moms1[i][4] = m_softphotons[0]->Momentum(); m_moms1[i][0] = m_moms1[i][1]+m_moms1[i][2]+m_moms1[i][3]+m_moms1[i][4]; m_softphotons.clear(); } } } double Scalar_To_Vector_Lepton_Neutrino::Smod(unsigned int kk) { m_moms = m_moms1[kk]; Vec4D k = m_moms[4]; Vec4D pi = m_moms[1]; Vec4D pj = m_moms[2]; double Zi = -1.; // both charged particles of opposite charge double Zj = +1.; int ti = +1; // always both final state int tj = +1; return m_alpha/(4.*M_PI*M_PI)*Zi*Zj*ti*tj*(pi/(pi*k)-pj/(pj*k)).Abs2(); } Complex Scalar_To_Vector_Lepton_Neutrino::InfraredSubtractedME_0_0() { m_moms = m_moms0; Vec4C epsV = conj(Polarization_Vector(m_moms[1]).at(m_spins[1])); double t = (m_moms[0]-m_moms[1]).Abs2(); Vec4C auxvec = Contraction(AuxiliaryTensor(m_moms[0],m_moms[1],t),2,epsV); XYZFunc XYZ(4,m_moms,m_flavs,false); if (m_ffmodel == true) { return 0; } else { return m_sqrt2*m_GF*XYZ.X(3,m_spins[3],auxvec,2,m_spins[2],m_cR,m_cL); } return 0; } Complex Scalar_To_Vector_Lepton_Neutrino::InfraredSubtractedME_0_1() { m_moms = m_moms0; return 0; } Complex Scalar_To_Vector_Lepton_Neutrino::InfraredSubtractedME_0_2() { return 0; } Complex Scalar_To_Vector_Lepton_Neutrino::InfraredSubtractedME_1_05(unsigned int i) { m_moms = m_moms1[i]; // set to set of momenta to be used Vec4C epsV = conj(Polarization_Vector(m_moms[1]).at(m_spins[1])); Vec4C epsP = conj(Polarization_Vector(m_moms[4]).at(m_spins[4])); double mV = m_masses[1]; // charged FS vector mass/propagator pole double ml = m_masses[2]; // lepton mass/propagator pole Vec4D qV = m_moms[1]+m_moms[4]; // vector propagator momentum Vec4D ql = m_moms[2]+m_moms[4]; // lepton propagator momentum m_moms[5] = m_moms[6] = ql; // enter those into m_moms m_flavs[5] = m_flavs[2]; // set to corresponding particle/antiparticle m_flavs[6] = m_flavs[2].Bar(); Vec4D k = m_moms[4]; XYZFunc XYZ(7,m_moms,m_flavs,false); m_flavs[5] = m_flavs[6] = Flavour(kf_none); // Form-Factor-Model Complex r1 = 0; Complex r2 = 0; Complex r3 = 0; Complex r4 = 0; if (m_ffmodel == true) { return 0; } // point-like else { // A double tA = (m_moms[0]-m_moms[1]).Abs2(); Vec4C auxvecA = Contraction(AuxiliaryTensor(m_moms[0],m_moms[1],tA),2,epsV); for (unsigned int s=0; s<=1; s++) { r1 = r1 + XYZ.X(3,m_spins[3],auxvecA,5,s,m_cR,m_cL) * XYZ.X(5,s,epsP,2,m_spins[2],1.,1.); r2 = r2 + XYZ.X(3,m_spins[3],auxvecA,6,s,m_cR,m_cL) * XYZ.X(6,s,epsP,2,m_spins[2],1.,1.); } r1 = -(0.5*m_sqrt2*m_e*m_GF)/(ql.Abs2()-ml*ml) * (1-ml/ql.Abs()) * r1; r2 = -(0.5*m_sqrt2*m_e*m_GF)/(ql.Abs2()-ml*ml) * (1+ml/ql.Abs()) * r2; // B double tB = (m_moms[0]-qV).Abs2(); Vec4C auxvecB = Contraction(AuxiliaryTensor(m_moms[0],qV,tB),2,epsP)*((m_moms[1]+2*k)*epsV) +Contraction(AuxiliaryTensor(m_moms[0],qV,tB),2,(m_moms[1]-k))*(epsV*epsP) +Contraction(AuxiliaryTensor(m_moms[0],qV,tB),2,epsV)*((-2*m_moms[1]-k)*epsP) -1./(mV*mV)*Contraction(AuxiliaryTensor(m_moms[0],qV,tB),2,qV) *( (qV*epsP)*((m_moms[1]+2*k)*epsV) +(epsV*epsP)*(qV*(m_moms[1]-k)) +(qV*epsV)*((-2*m_moms[1]-k)*epsP) ); r3 = (m_sqrt2*m_e*m_GF)/(qV.Abs2()-mV*mV) * XYZ.X(3,m_spins[3],auxvecB,2,m_spins[2],m_cR,m_cL); return r1+r2+r3; } return 0; } Complex Scalar_To_Vector_Lepton_Neutrino::InfraredSubtractedME_1_15(unsigned int i) { return 0; } Complex Scalar_To_Vector_Lepton_Neutrino::InfraredSubtractedME_2_1(unsigned int i, unsigned int j) { return 0; } double Scalar_To_Vector_Lepton_Neutrino::GetBeta_0_0() { double sum = 0; for (unsigned int i=0; i<=1; i++) { // spin neutrino for (unsigned int j=0; j<=1; j++) { // spin lepton for (unsigned int m=0; m<=2; m++) { // spin FS vector for (unsigned int k=0; k<=0; k++) { // spin IS scalar m_spins[0] = k; m_spins[1] = m; m_spins[2] = j; m_spins[3] = i; Complex M_0_0 = InfraredSubtractedME_0_0(); sum = sum + (M_0_0*conj(M_0_0)).real(); } } } } // spin avaraging gives factor 1 return sum; } double Scalar_To_Vector_Lepton_Neutrino::GetBeta_0_1() { return 0; // double sum = 0; // for (unsigned int i=0; i<=1; i++) { // spin neutrino // for (unsigned int j=0; j<=1; j++) { // spin lepton // for (unsigned int m=0; m<=2; m++) { // spin FS vector // for (unsigned int k=0; k<=0; k++) { // spin IS scalar // m_spins[0] = k; // m_spins[1] = m; // m_spins[2] = j; // m_spins[3] = i; // Complex M_0_0 = InfraredSubtractedME_0_0(); // Complex M_0_1 = InfraredSubtractedME_0_1(); // sum = sum + (M_0_0*conj(M_0_1)+M_0_1*conj(M_0_0)).real(); // } // } // } // } // // spin avaraging gives factor 1 // return sum; } double Scalar_To_Vector_Lepton_Neutrino::GetBeta_0_2() { return 0; } double Scalar_To_Vector_Lepton_Neutrino::GetBeta_1_1(unsigned int a) { double sum = 0; for (unsigned int i=0; i<=1; i++) { // spin neutrino for (unsigned int j=0; j<=1; j++) { // spin lepton for (unsigned int m=0; m<=2; m++) { // spin FS vector for (unsigned int k=0; k<=0; k++) { // spin IS scalar for (unsigned int l=0; l<=1; l++) { // spin gamma m_spins[0] = k; m_spins[1] = m; m_spins[2] = j; m_spins[3] = i; m_spins[4] = l; Complex M_1_05 = InfraredSubtractedME_1_05(a); sum = sum + (M_1_05*conj(M_1_05)).real(); } } } } } // spin avaraging gives factor 1 sum = 1./(16.*M_PI*M_PI*M_PI)*sum - Smod(a)*GetBeta_0_0(); return sum; } double Scalar_To_Vector_Lepton_Neutrino::GetBeta_1_2(unsigned int i) { return 0; } double Scalar_To_Vector_Lepton_Neutrino::GetBeta_2_2(unsigned int i, unsigned int j) { return 0; } double Scalar_To_Vector_Lepton_Neutrino::FormFactorAp(double t) { return -0.0851648 - 0.00291091*t - 0.0000786814*t*t; } double Scalar_To_Vector_Lepton_Neutrino::FormFactorAm(double t) { return 0.; } double Scalar_To_Vector_Lepton_Neutrino::FormFactorF(double t) { return 4.8048 + 0.0914977*t + 0.00174239*t*t; } double Scalar_To_Vector_Lepton_Neutrino::FormFactorG(double t) { return 0.104396 + 0.00400149*t + 0.000117909*t*t; } Lorentz_Ten2D Scalar_To_Vector_Lepton_Neutrino::AuxiliaryTensor(Vec4D p1, Vec4D p2, double t) { // convert t to dimensionless quantity tt=t/[MeV] double tt = t*(1000.*1000.); return FormFactorG(tt)*EpsilonTensorContraction(p1,p2) -FormFactorF(tt)*MetricTensor() -FormFactorAp(tt)*BuildTensor(p1+p2,p1-p2) -FormFactorAm(tt)*BuildTensor(p1-p2,p1+p2); } DECLARE_PHOTONS_ME_GETTER(Scalar_To_Vector_Lepton_Neutrino, "Scalar_To_Vector_Lepton_Neutrino") PHOTONS_ME_Base *ATOOLS::Getter:: operator()(const Particle_Vector_Vector &pvv) const { if ( (pvv.size() == 4) && (pvv[0].size() == 0) && (pvv[1].size() == 1) && pvv[1][0]->Flav().IsScalar() && pvv[1][0]->Flav().IsHadron() && (((pvv[2].size() == 2) && pvv[2][0]->Flav().IsVector() && pvv[2][0]->Flav().IsHadron() && pvv[2][1]->Flav().IsLepton()) || ((pvv[2].size() == 2) && pvv[2][1]->Flav().IsVector() && pvv[2][1]->Flav().IsHadron() && pvv[2][0]->Flav().IsLepton()) ) && (pvv[3].size() == 1) && pvv[3][0]->Flav().IsLepton()) return new Scalar_To_Vector_Lepton_Neutrino(pvv); return NULL; }